tcmp2.c revision 1.1.1.4
1/* Test file for mpfr_cmp2. 2 3Copyright 1999-2003, 2006-2018 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#include "mpfr-test.h" 24 25/* set bit n of x to b, where bit 0 is the most significant one */ 26static void 27set_bit (mpfr_t x, unsigned int n, int b) 28{ 29 unsigned l; 30 mp_size_t xn; 31 32 xn = (MPFR_PREC(x) - 1) / mp_bits_per_limb; 33 l = n / mp_bits_per_limb; 34 n %= mp_bits_per_limb; 35 n = mp_bits_per_limb - 1 - n; 36 if (b) 37 MPFR_MANT(x)[xn - l] |= MPFR_LIMB_ONE << n; 38 else 39 MPFR_MANT(x)[xn - l] &= ~(MPFR_LIMB_ONE << n); 40} 41 42/* check that for x = 1.u 1 v 0^k low(x) 43 y = 1.u 0 v 1^k low(y) 44 mpfr_cmp2 (x, y) returns 1 + |u| + |v| + k for low(x) >= low(y), 45 and 1 + |u| + |v| + k + 1 otherwise */ 46static void 47worst_cases (void) 48{ 49 mpfr_t x, y; 50 unsigned int i, j, k, b, expected; 51 mpfr_prec_t l; 52 53 mpfr_init2 (x, 200); 54 mpfr_init2 (y, 200); 55 56 mpfr_set_ui (y, 1, MPFR_RNDN); 57 for (i = 1; i < MPFR_PREC(x); i++) 58 { 59 mpfr_set_ui (x, 1, MPFR_RNDN); 60 mpfr_div_2exp (y, y, 1, MPFR_RNDN); /* y = 1/2^i */ 61 62 l = 0; 63 if (mpfr_cmp2 (x, y, &l) <= 0 || l != 1) 64 { 65 printf ("Error in mpfr_cmp2:\nx="); 66 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN); 67 printf ("\ny="); 68 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN); 69 printf ("\ngot %lu instead of 1\n", (unsigned long) l); 70 exit (1); 71 } 72 73 mpfr_add (x, x, y, MPFR_RNDN); /* x = 1 + 1/2^i */ 74 l = 0; 75 if (mpfr_cmp2 (x, y, &l) <= 0 || l != 0) 76 { 77 printf ("Error in mpfr_cmp2:\nx="); 78 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN); 79 printf ("\ny="); 80 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN); 81 printf ("\ngot %lu instead of 0\n", (unsigned long) l); 82 exit (1); 83 } 84 } 85 86 for (i = 0; i < 64; i++) /* |u| = i */ 87 { 88 mpfr_urandomb (x, RANDS); 89 mpfr_set (y, x, MPFR_RNDN); 90 set_bit (x, i + 1, 1); 91 set_bit (y, i + 1, 0); 92 for (j = 0; j < 64; j++) /* |v| = j */ 93 { 94 b = randlimb () % 2; 95 set_bit (x, i + j + 2, b); 96 set_bit (y, i + j + 2, b); 97 for (k = 0; k < 64; k++) 98 { 99 if (k) 100 set_bit (x, i + j + k + 1, 0); 101 if (k) 102 set_bit (y, i + j + k + 1, 1); 103 set_bit (x, i + j + k + 2, 1); 104 set_bit (y, i + j + k + 2, 0); 105 l = 0; mpfr_cmp2 (x, y, &l); 106 expected = i + j + k + 1; 107 if (l != expected) 108 { 109 printf ("Error in mpfr_cmp2:\nx="); 110 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN); 111 printf ("\ny="); 112 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN); 113 printf ("\ngot %lu instead of %u\n", 114 (unsigned long) l, expected); 115 exit (1); 116 } 117 set_bit (x, i + j + k + 2, 0); 118 set_bit (x, i + j + k + 3, 0); 119 set_bit (y, i + j + k + 3, 1); 120 l = 0; mpfr_cmp2 (x, y, &l); 121 expected = i + j + k + 2; 122 if (l != expected) 123 { 124 printf ("Error in mpfr_cmp2:\nx="); 125 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN); 126 printf ("\ny="); 127 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN); 128 printf ("\ngot %lu instead of %u\n", 129 (unsigned long) l, expected); 130 exit (1); 131 } 132 } 133 } 134 } 135 136 mpfr_clear (x); 137 mpfr_clear (y); 138} 139 140static void 141tcmp2 (double x, double y, int i) 142{ 143 mpfr_t xx, yy; 144 mpfr_prec_t j; 145 146 if (i == -1) 147 { 148 if (x == y) 149 i = 53; 150 else 151 i = (int) (__gmpfr_floor_log2 (x) - __gmpfr_floor_log2 (x - y)); 152 } 153 mpfr_init2(xx, 53); mpfr_init2(yy, 53); 154 mpfr_set_d (xx, x, MPFR_RNDN); 155 mpfr_set_d (yy, y, MPFR_RNDN); 156 j = 0; 157 if (mpfr_cmp2 (xx, yy, &j) == 0) 158 { 159 if (x != y) 160 { 161 printf ("Error in mpfr_cmp2 for\nx="); 162 mpfr_out_str (stdout, 2, 0, xx, MPFR_RNDN); 163 printf ("\ny="); 164 mpfr_out_str (stdout, 2, 0, yy, MPFR_RNDN); 165 printf ("\ngot sign 0 for x != y\n"); 166 exit (1); 167 } 168 } 169 else if (j != (unsigned) i) 170 { 171 printf ("Error in mpfr_cmp2 for\nx="); 172 mpfr_out_str (stdout, 2, 0, xx, MPFR_RNDN); 173 printf ("\ny="); 174 mpfr_out_str (stdout, 2, 0, yy, MPFR_RNDN); 175 printf ("\ngot %lu instead of %d\n", (unsigned long) j, i); 176 exit (1); 177 } 178 mpfr_clear(xx); mpfr_clear(yy); 179} 180 181static void 182special (void) 183{ 184 mpfr_t x, y; 185 mpfr_prec_t j; 186 187 mpfr_init (x); mpfr_init (y); 188 189 /* bug found by Nathalie Revol, 21 March 2001 */ 190 mpfr_set_prec (x, 65); 191 mpfr_set_prec (y, 65); 192 mpfr_set_str_binary (x, "0.10000000000000000000000000000000000001110010010110100110011110000E1"); 193 mpfr_set_str_binary (y, "0.11100100101101001100111011111111110001101001000011101001001010010E-35"); 194 j = 0; 195 if (mpfr_cmp2 (x, y, &j) <= 0 || j != 1) 196 { 197 printf ("Error in mpfr_cmp2:\n"); 198 printf ("x="); 199 mpfr_dump (x); 200 printf ("y="); 201 mpfr_dump (y); 202 printf ("got %lu, expected 1\n", (unsigned long) j); 203 exit (1); 204 } 205 206 mpfr_set_prec(x, 127); mpfr_set_prec(y, 127); 207 mpfr_set_str_binary(x, "0.1011010000110111111000000101011110110001000101101011011110010010011110010000101101000010011001100110010000000010110000101000101E6"); 208 mpfr_set_str_binary(y, "0.1011010000110111111000000101011011111100011101000011001111000010100010100110110100110010011001100110010000110010010110000010110E6"); 209 j = 0; 210 if (mpfr_cmp2(x, y, &j) <= 0 || j != 32) 211 { 212 printf ("Error in mpfr_cmp2:\n"); 213 printf ("x="); mpfr_dump (x); 214 printf ("y="); mpfr_dump (y); 215 printf ("got %lu, expected 32\n", (unsigned long) j); 216 exit (1); 217 } 218 219 mpfr_set_prec (x, 128); mpfr_set_prec (y, 239); 220 mpfr_set_str_binary (x, "0.10001000110110000111011000101011111100110010010011001101000011111010010110001000000010100110100111111011011010101100100000000000E167"); 221 mpfr_set_str_binary (y, "0.10001000110110000111011000101011111100110010010011001101000011111010010110001000000010100110100111111011011010101100011111111111111111111111111111111111111111111111011111100101011100011001101000100111000010000000000101100110000111111000101E167"); 222 j = 0; 223 if (mpfr_cmp2(x, y, &j) <= 0 || j != 164) 224 { 225 printf ("Error in mpfr_cmp2:\n"); 226 printf ("x="); mpfr_dump (x); 227 printf ("y="); mpfr_dump (y); 228 printf ("got %lu, expected 164\n", (unsigned long) j); 229 exit (1); 230 } 231 232 /* bug found by Nathalie Revol, 29 March 2001 */ 233 mpfr_set_prec (x, 130); mpfr_set_prec (y, 130); 234 mpfr_set_str_binary (x, "0.1100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E2"); 235 mpfr_set_str_binary (y, "0.1011111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111100E2"); 236 j = 0; 237 if (mpfr_cmp2(x, y, &j) <= 0 || j != 127) 238 { 239 printf ("Error in mpfr_cmp2:\n"); 240 printf ("x="); mpfr_dump (x); 241 printf ("y="); mpfr_dump (y); 242 printf ("got %lu, expected 127\n", (unsigned long) j); 243 exit (1); 244 } 245 246 /* bug found by Nathalie Revol, 2 Apr 2001 */ 247 mpfr_set_prec (x, 65); mpfr_set_prec (y, 65); 248 mpfr_set_ui (x, 5, MPFR_RNDN); 249 mpfr_set_str_binary (y, "0.10011111111111111111111111111111111111111111111111111111111111101E3"); 250 j = 0; 251 if (mpfr_cmp2(x, y, &j) <= 0 || j != 63) 252 { 253 printf ("Error in mpfr_cmp2:\n"); 254 printf ("x="); mpfr_dump (x); 255 printf ("y="); mpfr_dump (y); 256 printf ("got %lu, expected 63\n", (unsigned long) j); 257 exit (1); 258 } 259 260 /* bug found by Nathalie Revol, 2 Apr 2001 */ 261 mpfr_set_prec (x, 65); mpfr_set_prec (y, 65); 262 mpfr_set_str_binary (x, "0.10011011111000101001110000000000000000000000000000000000000000000E-69"); 263 mpfr_set_str_binary (y, "0.10011011111000101001101111111111111111111111111111111111111111101E-69"); 264 j = 0; 265 if (mpfr_cmp2(x, y, &j) <= 0 || j != 63) 266 { 267 printf ("Error in mpfr_cmp2:\n"); 268 printf ("x="); mpfr_dump (x); 269 printf ("y="); mpfr_dump (y); 270 printf ("got %lu, expected 63\n", (unsigned long) j); 271 exit (1); 272 } 273 274 mpfr_set_prec (x, 65); 275 mpfr_set_prec (y, 32); 276 mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101"); 277 mpfr_set_str_binary (y, "0.11101110111100011101110111111111"); 278 if (mpfr_cmp2 (x, y, &j) <= 0 || j != 0) 279 { 280 printf ("Error in mpfr_cmp2 (1)\n"); 281 exit (1); 282 } 283 284 mpfr_set_prec (x, 2 * GMP_NUMB_BITS); 285 mpfr_set_prec (y, GMP_NUMB_BITS); 286 mpfr_set_ui (x, 1, MPFR_RNDN); /* x = 1 */ 287 mpfr_set_ui (y, 1, MPFR_RNDN); 288 mpfr_nextbelow (y); /* y = 1 - 2^(-GMP_NUMB_BITS) */ 289 mpfr_cmp2 (x, y, &j); 290 if (mpfr_cmp2 (x, y, &j) <= 0 || j != GMP_NUMB_BITS) 291 { 292 printf ("Error in mpfr_cmp2 (2)\n"); 293 exit (1); 294 } 295 296 mpfr_clear (x); 297 mpfr_clear (y); 298} 299 300/* Compare (m,kx) and (m,ky), where (m,k) means m fixed limbs followed by 301 k zero limbs. */ 302static void 303test_equal (void) 304{ 305 mpfr_t w, x, y; 306 int m, kx, ky, inex; 307 mpfr_prec_t j; 308 309 for (m = 1; m <= 4; m++) 310 { 311 mpfr_init2 (w, m * GMP_NUMB_BITS); 312 for (kx = 0; kx <= 4; kx++) 313 for (ky = 0; ky <= 4; ky++) 314 { 315 do mpfr_urandomb (w, RANDS); while (mpfr_zero_p (w)); 316 mpfr_init2 (x, (m + kx) * GMP_NUMB_BITS 317 - (kx == 0 ? 0 : randlimb () % GMP_NUMB_BITS)); 318 mpfr_init2 (y, (m + ky) * GMP_NUMB_BITS 319 - (ky == 0 ? 0 : randlimb () % GMP_NUMB_BITS)); 320 inex = mpfr_set (x, w, MPFR_RNDN); 321 MPFR_ASSERTN (inex == 0); 322 inex = mpfr_set (y, w, MPFR_RNDN); 323 MPFR_ASSERTN (inex == 0); 324 MPFR_ASSERTN (mpfr_equal_p (x, y)); 325 if (randlimb () & 1) 326 mpfr_neg (x, x, MPFR_RNDN); 327 if (randlimb () & 1) 328 mpfr_neg (y, y, MPFR_RNDN); 329 if (mpfr_cmp2 (x, y, &j) != 0) 330 { 331 printf ("Error in test_equal for m = %d, kx = %d, ky = %d\n", 332 m, kx, ky); 333 printf (" x = "); 334 mpfr_dump (x); 335 printf (" y = "); 336 mpfr_dump (y); 337 exit (1); 338 } 339 mpfr_clears (x, y, (mpfr_ptr) 0); 340 } 341 mpfr_clear (w); 342 } 343} 344 345int 346main (void) 347{ 348 int i; 349 long j; 350 double x, y, z; 351 352 tests_start_mpfr (); 353 mpfr_test_init (); 354 355 worst_cases (); 356 special (); 357 tcmp2 (5.43885304644369510000e+185, -1.87427265794105340000e-57, 1); 358 tcmp2 (1.06022698059744327881e+71, 1.05824655795525779205e+71, -1); 359 tcmp2 (1.0, 1.0, 53); 360 /* warning: cmp2 does not allow 0 as input */ 361 for (x = 0.5, i = 1; i < 54; i++) 362 { 363 tcmp2 (1.0, 1.0-x, i); 364 x /= 2.0; 365 } 366 for (x = 0.5, i = 1; i < 100; i++) 367 { 368 tcmp2 (1.0, x, 1); 369 x /= 2.0; 370 } 371 for (j = 0; j < 100000; j++) 372 { 373 x = DBL_RAND (); 374 y = DBL_RAND (); 375 if (x < y) 376 { 377 z = x; 378 x = y; 379 y = z; 380 } 381 if (y != 0.0) 382 tcmp2 (x, y, -1); 383 } 384 385 test_equal (); 386 387 tests_end_mpfr (); 388 389 return 0; 390} 391