1/* Test file for mpfr_div (and some mpfr_div_ui, etc. tests). 2 3Copyright 1999, 2001-2023 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 20https://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 25static void 26check_equal (mpfr_srcptr a, mpfr_srcptr a2, const char *s, 27 mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t r) 28{ 29 if (SAME_VAL (a, a2)) 30 return; 31 if (r == MPFR_RNDF) /* RNDF might return different values */ 32 return; 33 printf ("Error in %s\n", mpfr_print_rnd_mode (r)); 34 printf ("b = "); 35 mpfr_dump (b); 36 printf ("c = "); 37 mpfr_dump (c); 38 printf ("mpfr_div result: "); 39 mpfr_dump (a); 40 printf ("%s result: ", s); 41 mpfr_dump (a2); 42 exit (1); 43} 44 45static int 46mpfr_all_div (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t r) 47{ 48 mpfr_t a2; 49 unsigned int oldflags, newflags; 50 int inex, inex2; 51 52 oldflags = __gmpfr_flags; 53 inex = mpfr_div (a, b, c, r); 54 55 /* this test makes no sense for RNDF, since it compares the ternary value 56 and the flags */ 57 if (a == b || a == c || r == MPFR_RNDF) 58 return inex; 59 60 newflags = __gmpfr_flags; 61 62 mpfr_init2 (a2, MPFR_PREC (a)); 63 64 if (mpfr_integer_p (b) && ! (MPFR_IS_ZERO (b) && MPFR_IS_NEG (b))) 65 { 66 /* b is an integer, but not -0 (-0 is rejected as 67 it becomes +0 when converted to an integer). */ 68 if (mpfr_fits_ulong_p (b, MPFR_RNDA)) 69 { 70 __gmpfr_flags = oldflags; 71 inex2 = mpfr_ui_div (a2, mpfr_get_ui (b, MPFR_RNDN), c, r); 72 if (!SAME_SIGN (inex2, inex)) 73 { 74 printf ("Error for ternary value (rnd=%s), mpfr_div %d, mpfr_ui_div %d\n", 75 mpfr_print_rnd_mode (r), inex, inex2); 76 exit (1); 77 } 78 if (__gmpfr_flags != newflags) 79 { 80 printf ("Error for flags, mpfr_div %d, mpfr_ui_div %d\n", 81 newflags, __gmpfr_flags); 82 exit (1); 83 } 84 check_equal (a, a2, "mpfr_ui_div", b, c, r); 85 } 86 if (mpfr_fits_slong_p (b, MPFR_RNDA)) 87 { 88 __gmpfr_flags = oldflags; 89 inex2 = mpfr_si_div (a2, mpfr_get_si (b, MPFR_RNDN), c, r); 90 MPFR_ASSERTN (SAME_SIGN (inex2, inex)); 91 MPFR_ASSERTN (__gmpfr_flags == newflags); 92 check_equal (a, a2, "mpfr_si_div", b, c, r); 93 } 94 } 95 96 if (mpfr_integer_p (c) && ! (MPFR_IS_ZERO (c) && MPFR_IS_NEG (c))) 97 { 98 /* c is an integer, but not -0 (-0 is rejected as 99 it becomes +0 when converted to an integer). */ 100 if (mpfr_fits_ulong_p (c, MPFR_RNDA)) 101 { 102 __gmpfr_flags = oldflags; 103 inex2 = mpfr_div_ui (a2, b, mpfr_get_ui (c, MPFR_RNDN), r); 104 MPFR_ASSERTN (SAME_SIGN (inex2, inex)); 105 MPFR_ASSERTN (__gmpfr_flags == newflags); 106 check_equal (a, a2, "mpfr_div_ui", b, c, r); 107 } 108 if (mpfr_fits_slong_p (c, MPFR_RNDA)) 109 { 110 __gmpfr_flags = oldflags; 111 inex2 = mpfr_div_si (a2, b, mpfr_get_si (c, MPFR_RNDN), r); 112 MPFR_ASSERTN (SAME_SIGN (inex2, inex)); 113 MPFR_ASSERTN (__gmpfr_flags == newflags); 114 check_equal (a, a2, "mpfr_div_si", b, c, r); 115 } 116 } 117 118 mpfr_clear (a2); 119 120 return inex; 121} 122 123#ifdef CHECK_EXTERNAL 124static int 125test_div (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) 126{ 127 int res; 128 int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_number_p (c); 129 if (ok) 130 { 131 mpfr_print_raw (b); 132 printf (" "); 133 mpfr_print_raw (c); 134 } 135 res = mpfr_all_div (a, b, c, rnd_mode); 136 if (ok) 137 { 138 printf (" "); 139 mpfr_print_raw (a); 140 printf ("\n"); 141 } 142 return res; 143} 144#else 145#define test_div mpfr_all_div 146#endif 147 148#define check53(n, d, rnd, res) check4(n, d, rnd, 53, res) 149 150/* return 0 iff a and b are of the same sign */ 151static int 152inex_cmp (int a, int b) 153{ 154 if (a > 0) 155 return (b > 0) ? 0 : 1; 156 else if (a == 0) 157 return (b == 0) ? 0 : 1; 158 else 159 return (b < 0) ? 0 : 1; 160} 161 162static void 163check4 (const char *Ns, const char *Ds, mpfr_rnd_t rnd_mode, int p, 164 const char *Qs) 165{ 166 mpfr_t q, n, d; 167 168 mpfr_inits2 (p, q, n, d, (mpfr_ptr) 0); 169 mpfr_set_str1 (n, Ns); 170 mpfr_set_str1 (d, Ds); 171 test_div(q, n, d, rnd_mode); 172 if (mpfr_cmp_str (q, Qs, ((p==53) ? 10 : 2), MPFR_RNDN) ) 173 { 174 printf ("mpfr_div failed for n=%s, d=%s, p=%d, rnd_mode=%s\n", 175 Ns, Ds, p, mpfr_print_rnd_mode (rnd_mode)); 176 printf ("got "); 177 mpfr_dump (q); 178 mpfr_set_str (q, Qs, ((p==53) ? 10 : 2), MPFR_RNDN); 179 printf ("expected "); 180 mpfr_dump (q); 181 exit (1); 182 } 183 mpfr_clears (q, n, d, (mpfr_ptr) 0); 184} 185 186static void 187check24 (const char *Ns, const char *Ds, mpfr_rnd_t rnd_mode, const char *Qs) 188{ 189 mpfr_t q, n, d; 190 191 mpfr_inits2 (24, q, n, d, (mpfr_ptr) 0); 192 193 mpfr_set_str1 (n, Ns); 194 mpfr_set_str1 (d, Ds); 195 test_div(q, n, d, rnd_mode); 196 if (mpfr_cmp_str1 (q, Qs) ) 197 { 198 printf ("mpfr_div failed for n=%s, d=%s, prec=24, rnd_mode=%s\n", 199 Ns, Ds, mpfr_print_rnd_mode(rnd_mode)); 200 printf ("expected quotient is %s, got ", Qs); 201 mpfr_out_str(stdout,10,0,q, MPFR_RNDN); putchar('\n'); 202 exit (1); 203 } 204 mpfr_clears (q, n, d, (mpfr_ptr) 0); 205} 206 207/* the following examples come from the paper "Number-theoretic Test 208 Generation for Directed Rounding" from Michael Parks, Table 2 */ 209static void 210check_float(void) 211{ 212 check24("70368760954880.0", "8388609.0", MPFR_RNDN, "8.388609e6"); 213 check24("140737479966720.0", "16777213.0", MPFR_RNDN, "8.388609e6"); 214 check24("70368777732096.0", "8388611.0", MPFR_RNDN, "8.388609e6"); 215 check24("105553133043712.0", "12582911.0", MPFR_RNDN, "8.38861e6"); 216 /* the exponent for the following example was forgotten in 217 the Arith'14 version of Parks' paper */ 218 check24 ("12582913.0", "12582910.0", MPFR_RNDN, "1.000000238"); 219 check24 ("105553124655104.0", "12582910.0", MPFR_RNDN, "8388610.0"); 220 check24("140737479966720.0", "8388609.0", MPFR_RNDN, "1.6777213e7"); 221 check24("70368777732096.0", "8388609.0", MPFR_RNDN, "8.388611e6"); 222 check24("105553133043712.0", "8388610.0", MPFR_RNDN, "1.2582911e7"); 223 check24("105553124655104.0", "8388610.0", MPFR_RNDN, "1.258291e7"); 224 225 check24("70368760954880.0", "8388609.0", MPFR_RNDZ, "8.388608e6"); 226 check24("140737479966720.0", "16777213.0", MPFR_RNDZ, "8.388609e6"); 227 check24("70368777732096.0", "8388611.0", MPFR_RNDZ, "8.388608e6"); 228 check24("105553133043712.0", "12582911.0", MPFR_RNDZ, "8.38861e6"); 229 check24("12582913.0", "12582910.0", MPFR_RNDZ, "1.000000238"); 230 check24 ("105553124655104.0", "12582910.0", MPFR_RNDZ, "8388610.0"); 231 check24("140737479966720.0", "8388609.0", MPFR_RNDZ, "1.6777213e7"); 232 check24("70368777732096.0", "8388609.0", MPFR_RNDZ, "8.38861e6"); 233 check24("105553133043712.0", "8388610.0", MPFR_RNDZ, "1.2582911e7"); 234 check24("105553124655104.0", "8388610.0", MPFR_RNDZ, "1.258291e7"); 235 236 check24("70368760954880.0", "8388609.0", MPFR_RNDU, "8.388609e6"); 237 check24("140737479966720.0", "16777213.0", MPFR_RNDU, "8.38861e6"); 238 check24("70368777732096.0", "8388611.0", MPFR_RNDU, "8.388609e6"); 239 check24("105553133043712.0", "12582911.0", MPFR_RNDU, "8.388611e6"); 240 check24("12582913.0", "12582910.0", MPFR_RNDU, "1.000000357"); 241 check24 ("105553124655104.0", "12582910.0", MPFR_RNDU, "8388611.0"); 242 check24("140737479966720.0", "8388609.0", MPFR_RNDU, "1.6777214e7"); 243 check24("70368777732096.0", "8388609.0", MPFR_RNDU, "8.388611e6"); 244 check24("105553133043712.0", "8388610.0", MPFR_RNDU, "1.2582912e7"); 245 check24("105553124655104.0", "8388610.0", MPFR_RNDU, "1.2582911e7"); 246 247 check24("70368760954880.0", "8388609.0", MPFR_RNDD, "8.388608e6"); 248 check24("140737479966720.0", "16777213.0", MPFR_RNDD, "8.388609e6"); 249 check24("70368777732096.0", "8388611.0", MPFR_RNDD, "8.388608e6"); 250 check24("105553133043712.0", "12582911.0", MPFR_RNDD, "8.38861e6"); 251 check24("12582913.0", "12582910.0", MPFR_RNDD, "1.000000238"); 252 check24 ("105553124655104.0", "12582910.0", MPFR_RNDD, "8388610.0"); 253 check24("140737479966720.0", "8388609.0", MPFR_RNDD, "1.6777213e7"); 254 check24("70368777732096.0", "8388609.0", MPFR_RNDD, "8.38861e6"); 255 check24("105553133043712.0", "8388610.0", MPFR_RNDD, "1.2582911e7"); 256 check24("105553124655104.0", "8388610.0", MPFR_RNDD, "1.258291e7"); 257 258 check24("70368760954880.0", "8388609.0", MPFR_RNDA, "8.388609e6"); 259} 260 261static void 262check_double(void) 263{ 264 check53("0.0", "1.0", MPFR_RNDZ, "0.0"); 265 check53("-7.4988969224688591e63", "4.8816866450288732e306", MPFR_RNDD, 266 "-1.5361282826510687291e-243"); 267 check53("-1.33225773037748601769e+199", "3.63449540676937123913e+79", 268 MPFR_RNDZ, "-3.6655920045905428978e119"); 269 check53("9.89438396044940256501e-134", "5.93472984109987421717e-67",MPFR_RNDU, 270 "1.6672003992376663654e-67"); 271 check53("9.89438396044940256501e-134", "5.93472984109987421717e-67",MPFR_RNDA, 272 "1.6672003992376663654e-67"); 273 check53("9.89438396044940256501e-134", "-5.93472984109987421717e-67", 274 MPFR_RNDU, "-1.6672003992376663654e-67"); 275 check53("-4.53063926135729747564e-308", "7.02293374921793516813e-84", 276 MPFR_RNDD, "-6.4512060388748850857e-225"); 277 check53("6.25089225176473806123e-01","-2.35527154824420243364e-230", 278 MPFR_RNDD, "-2.6540006635008291192e229"); 279 check53("6.25089225176473806123e-01","-2.35527154824420243364e-230", 280 MPFR_RNDA, "-2.6540006635008291192e229"); 281 check53("6.52308934689126e15", "-1.62063546601505417497e273", MPFR_RNDN, 282 "-4.0250194961676020848e-258"); 283 check53("1.04636807108079349236e-189", "3.72295730823253012954e-292", 284 MPFR_RNDZ, "2.810583051186143125e102"); 285 /* problems found by Kevin under HP-PA */ 286 check53 ("2.861044553323177e-136", "-1.1120354257068143e+45", MPFR_RNDZ, 287 "-2.5727998292003016e-181"); 288 check53 ("-4.0559157245809205e-127", "-1.1237723844524865e+77", MPFR_RNDN, 289 "3.6091968273068081e-204"); 290 check53 ("-1.8177943561493235e-93", "-8.51233984260364e-104", MPFR_RNDU, 291 "2.1354814184595821e+10"); 292} 293 294static void 295check_64(void) 296{ 297 mpfr_t x,y,z; 298 299 mpfr_inits2 (64, x, y, z, (mpfr_ptr) 0); 300 301 mpfr_set_str_binary(x, "1.00100100110110101001010010101111000001011100100101010000000000E54"); 302 mpfr_set_str_binary(y, "1.00000000000000000000000000000000000000000000000000000000000000E584"); 303 test_div(z, x, y, MPFR_RNDU); 304 if (mpfr_cmp_str (z, "0.1001001001101101010010100101011110000010111001001010100000000000E-529", 2, MPFR_RNDN)) 305 { 306 printf ("Error for tdiv for MPFR_RNDU and p=64\nx="); 307 mpfr_dump (x); 308 printf ("y="); 309 mpfr_dump (y); 310 printf ("got "); 311 mpfr_dump (z); 312 printf ("expected 0.1001001001101101010010100101011110000010111001001010100000000000E-529\n"); 313 exit (1); 314 } 315 316 mpfr_clears (x, y, z, (mpfr_ptr) 0); 317} 318 319static void 320check_convergence (void) 321{ 322 mpfr_t x, y; int i, j; 323 324 mpfr_init2(x, 130); 325 mpfr_set_str_binary(x, "0.1011111101011010101000001010011111101000011100011101010011111011000011001010000000111100100111110011001010110100100001001000111001E6944"); 326 mpfr_init2(y, 130); 327 mpfr_set_ui(y, 5, MPFR_RNDN); 328 test_div(x, x, y, MPFR_RNDD); /* exact division */ 329 330 mpfr_set_prec(x, 64); 331 mpfr_set_prec(y, 64); 332 mpfr_set_str_binary(x, "0.10010010011011010100101001010111100000101110010010101E55"); 333 mpfr_set_str_binary(y, "0.1E585"); 334 test_div(x, x, y, MPFR_RNDN); 335 mpfr_set_str_binary(y, "0.10010010011011010100101001010111100000101110010010101E-529"); 336 if (mpfr_cmp (x, y)) 337 { 338 printf ("Error in mpfr_div for prec=64, rnd=MPFR_RNDN\n"); 339 printf ("got "); mpfr_dump (x); 340 printf ("instead of "); mpfr_dump (y); 341 exit(1); 342 } 343 344 for (i=32; i<=64; i+=32) 345 { 346 mpfr_set_prec(x, i); 347 mpfr_set_prec(y, i); 348 mpfr_set_ui(x, 1, MPFR_RNDN); 349 RND_LOOP(j) 350 { 351 mpfr_set_ui (y, 1, MPFR_RNDN); 352 test_div (y, x, y, (mpfr_rnd_t) j); 353 if (mpfr_cmp_ui (y, 1)) 354 { 355 printf ("mpfr_div failed for x=1.0, y=1.0, prec=%d rnd=%s\n", 356 i, mpfr_print_rnd_mode ((mpfr_rnd_t) j)); 357 printf ("got "); mpfr_dump (y); 358 exit (1); 359 } 360 } 361 } 362 363 mpfr_clear (x); 364 mpfr_clear (y); 365} 366 367#define KMAX 10000 368 369/* given y = o(x/u), x, u, find the inexact flag by 370 multiplying y by u */ 371static int 372get_inexact (mpfr_ptr y, mpfr_ptr x, mpfr_ptr u) 373{ 374 mpfr_t xx; 375 int inex; 376 mpfr_init2 (xx, mpfr_get_prec (y) + mpfr_get_prec (u)); 377 mpfr_mul (xx, y, u, MPFR_RNDN); /* exact */ 378 inex = mpfr_cmp (xx, x); 379 mpfr_clear (xx); 380 return inex; 381} 382 383static void 384check_hard (void) 385{ 386 mpfr_t u, v, q, q2; 387 mpfr_prec_t precu, precv, precq; 388 int rnd; 389 int inex, inex2, i, j; 390 391 mpfr_init (q); 392 mpfr_init (q2); 393 mpfr_init (u); 394 mpfr_init (v); 395 396 for (precq = MPFR_PREC_MIN; precq <= 64; precq ++) 397 { 398 mpfr_set_prec (q, precq); 399 mpfr_set_prec (q2, precq + 1); 400 for (j = 0; j < 2; j++) 401 { 402 if (j == 0) 403 { 404 do 405 { 406 mpfr_urandomb (q2, RANDS); 407 } 408 while (mpfr_cmp_ui (q2, 0) == 0); 409 } 410 else /* use q2=1 */ 411 mpfr_set_ui (q2, 1, MPFR_RNDN); 412 for (precv = precq; precv <= 10 * precq; precv += precq) 413 { 414 mpfr_set_prec (v, precv); 415 do 416 { 417 mpfr_urandomb (v, RANDS); 418 } 419 while (mpfr_cmp_ui (v, 0) == 0); 420 for (precu = precq; precu <= 10 * precq; precu += precq) 421 { 422 mpfr_set_prec (u, precu); 423 mpfr_mul (u, v, q2, MPFR_RNDN); 424 mpfr_nextbelow (u); 425 for (i = 0; i <= 2; i++) 426 { 427 RND_LOOP_NO_RNDF (rnd) 428 { 429 inex = test_div (q, u, v, (mpfr_rnd_t) rnd); 430 inex2 = get_inexact (q, u, v); 431 if (inex_cmp (inex, inex2)) 432 { 433 printf ("Wrong inexact flag for rnd=%s: expected %d, got %d\n", 434 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), inex2, inex); 435 printf ("u= "); mpfr_dump (u); 436 printf ("v= "); mpfr_dump (v); 437 printf ("q= "); mpfr_dump (q); 438 mpfr_set_prec (q2, precq + precv); 439 mpfr_mul (q2, q, v, MPFR_RNDN); 440 printf ("q*v="); mpfr_dump (q2); 441 exit (1); 442 } 443 } 444 mpfr_nextabove (u); 445 } 446 } 447 } 448 } 449 } 450 451 mpfr_clear (q); 452 mpfr_clear (q2); 453 mpfr_clear (u); 454 mpfr_clear (v); 455} 456 457static void 458check_lowr (void) 459{ 460 mpfr_t x, y, z, z2, z3, tmp; 461 int k, c, c2; 462 463 464 mpfr_init2 (x, 1000); 465 mpfr_init2 (y, 100); 466 mpfr_init2 (tmp, 850); 467 mpfr_init2 (z, 10); 468 mpfr_init2 (z2, 10); 469 mpfr_init2 (z3, 50); 470 471 for (k = 1; k < KMAX; k++) 472 { 473 do 474 { 475 mpfr_urandomb (z, RANDS); 476 } 477 while (mpfr_cmp_ui (z, 0) == 0); 478 do 479 { 480 mpfr_urandomb (tmp, RANDS); 481 } 482 while (mpfr_cmp_ui (tmp, 0) == 0); 483 mpfr_mul (x, z, tmp, MPFR_RNDN); /* exact */ 484 c = test_div (z2, x, tmp, MPFR_RNDN); 485 486 if (c || mpfr_cmp (z2, z)) 487 { 488 printf ("Error in mpfr_div rnd=MPFR_RNDN\n"); 489 printf ("got "); mpfr_dump (z2); 490 printf ("instead of "); mpfr_dump (z); 491 printf ("inex flag = %d, expected 0\n", c); 492 exit (1); 493 } 494 } 495 496 /* x has still precision 1000, z precision 10, and tmp prec 850 */ 497 mpfr_set_prec (z2, 9); 498 for (k = 1; k < KMAX; k++) 499 { 500 mpfr_urandomb (z, RANDS); 501 do 502 { 503 mpfr_urandomb (tmp, RANDS); 504 } 505 while (mpfr_cmp_ui (tmp, 0) == 0); 506 mpfr_mul (x, z, tmp, MPFR_RNDN); /* exact */ 507 c = test_div (z2, x, tmp, MPFR_RNDN); 508 /* since z2 has one less bit that z, either the division is exact 509 if z is representable on 9 bits, or we have an even round case */ 510 511 c2 = get_inexact (z2, x, tmp); 512 if ((mpfr_cmp (z2, z) == 0 && c) || inex_cmp (c, c2)) 513 { 514 printf ("Error in mpfr_div rnd=MPFR_RNDN\n"); 515 printf ("got "); mpfr_dump (z2); 516 printf ("instead of "); mpfr_dump (z); 517 printf ("inex flag = %d, expected %d\n", c, c2); 518 exit (1); 519 } 520 else if (c == 2) 521 { 522 mpfr_nexttoinf (z); 523 if (mpfr_cmp(z2, z)) 524 { 525 printf ("Error in mpfr_div [even rnd?] rnd=MPFR_RNDN\n"); 526 printf ("Dividing "); 527 printf ("got "); mpfr_dump (z2); 528 printf ("instead of "); mpfr_dump (z); 529 printf ("inex flag = %d\n", 1); 530 exit (1); 531 } 532 } 533 else if (c == -2) 534 { 535 mpfr_nexttozero (z); 536 if (mpfr_cmp(z2, z)) 537 { 538 printf ("Error in mpfr_div [even rnd?] rnd=MPFR_RNDN\n"); 539 printf ("Dividing "); 540 printf ("got "); mpfr_dump (z2); 541 printf ("instead of "); mpfr_dump (z); 542 printf ("inex flag = %d\n", 1); 543 exit (1); 544 } 545 } 546 } 547 548 mpfr_set_prec(x, 1000); 549 mpfr_set_prec(y, 100); 550 mpfr_set_prec(tmp, 850); 551 mpfr_set_prec(z, 10); 552 mpfr_set_prec(z2, 10); 553 554 /* almost exact divisions */ 555 for (k = 1; k < KMAX; k++) 556 { 557 do 558 { 559 mpfr_urandomb (z, RANDS); 560 } 561 while (mpfr_cmp_ui (z, 0) == 0); 562 do 563 { 564 mpfr_urandomb (tmp, RANDS); 565 } 566 while (mpfr_cmp_ui (tmp, 0) == 0); 567 mpfr_mul(x, z, tmp, MPFR_RNDN); 568 mpfr_set(y, tmp, MPFR_RNDD); 569 mpfr_nexttoinf (x); 570 571 c = test_div(z2, x, y, MPFR_RNDD); 572 test_div(z3, x, y, MPFR_RNDD); 573 mpfr_set(z, z3, MPFR_RNDD); 574 575 if (c != -1 || mpfr_cmp(z2, z)) 576 { 577 printf ("Error in mpfr_div rnd=MPFR_RNDD\n"); 578 printf ("got "); mpfr_dump (z2); 579 printf ("instead of "); mpfr_dump (z); 580 printf ("inex flag = %d\n", c); 581 exit (1); 582 } 583 584 mpfr_set (y, tmp, MPFR_RNDU); 585 test_div (z3, x, y, MPFR_RNDU); 586 mpfr_set (z, z3, MPFR_RNDU); 587 c = test_div (z2, x, y, MPFR_RNDU); 588 if (c != 1 || mpfr_cmp (z2, z)) 589 { 590 printf ("Error in mpfr_div rnd=MPFR_RNDU\n"); 591 printf ("u="); mpfr_dump (x); 592 printf ("v="); mpfr_dump (y); 593 printf ("got "); mpfr_dump (z2); 594 printf ("instead of "); mpfr_dump (z); 595 printf ("inex flag = %d\n", c); 596 exit (1); 597 } 598 } 599 600 mpfr_clear (x); 601 mpfr_clear (y); 602 mpfr_clear (z); 603 mpfr_clear (z2); 604 mpfr_clear (z3); 605 mpfr_clear (tmp); 606} 607 608#define MAX_PREC 128 609 610static void 611check_inexact (void) 612{ 613 mpfr_t x, y, z, u; 614 mpfr_prec_t px, py, pu; 615 int inexact, cmp; 616 mpfr_rnd_t rnd; 617 618 mpfr_init (x); 619 mpfr_init (y); 620 mpfr_init (z); 621 mpfr_init (u); 622 623 mpfr_set_prec (x, 28); 624 mpfr_set_prec (y, 28); 625 mpfr_set_prec (z, 1023); 626 mpfr_set_str_binary (x, "0.1000001001101101111100010011E0"); 627 mpfr_set_str (z, "48284762641021308813686974720835219181653367326353400027913400579340343320519877153813133510034402932651132854764198688352364361009429039801248971901380781746767119334993621199563870113045276395603170432175354501451429471578325545278975153148347684600400321033502982713296919861760382863826626093689036010394", 10, MPFR_RNDN); 628 mpfr_div (x, x, z, MPFR_RNDN); 629 mpfr_set_str_binary (y, "0.1111001011001101001001111100E-1023"); 630 if (mpfr_cmp (x, y)) 631 { 632 printf ("Error in mpfr_div for prec=28, RNDN\n"); 633 printf ("Expected "); mpfr_dump (y); 634 printf ("Got "); mpfr_dump (x); 635 exit (1); 636 } 637 638 mpfr_set_prec (x, 53); 639 mpfr_set_str_binary (x, "0.11101100110010100011011000000100001111011111110010101E0"); 640 mpfr_set_prec (u, 127); 641 mpfr_set_str_binary (u, "0.1000001100110110110101110110101101111000110000001111111110000000011111001010110100110010111111111101000001011011101011101101000E-2"); 642 mpfr_set_prec (y, 95); 643 inexact = test_div (y, x, u, MPFR_RNDN); 644 if (inexact != (cmp = get_inexact (y, x, u))) 645 { 646 printf ("Wrong inexact flag (0): expected %d, got %d\n", cmp, inexact); 647 printf ("x="); mpfr_out_str (stdout, 10, 99, x, MPFR_RNDN); printf ("\n"); 648 printf ("u="); mpfr_out_str (stdout, 10, 99, u, MPFR_RNDN); printf ("\n"); 649 printf ("y="); mpfr_out_str (stdout, 10, 99, y, MPFR_RNDN); printf ("\n"); 650 exit (1); 651 } 652 653 mpfr_set_prec (x, 33); 654 mpfr_set_str_binary (x, "0.101111100011011101010011101100001E0"); 655 mpfr_set_prec (u, 2); 656 mpfr_set_str_binary (u, "0.1E0"); 657 mpfr_set_prec (y, 28); 658 inexact = test_div (y, x, u, MPFR_RNDN); 659 if (inexact >= 0) 660 { 661 printf ("Wrong inexact flag (1): expected -1, got %d\n", 662 inexact); 663 exit (1); 664 } 665 666 mpfr_set_prec (x, 129); 667 mpfr_set_str_binary (x, "0.111110101111001100000101011100101100110011011101010001000110110101100101000010000001110110100001101010001010100010001111001101010E-2"); 668 mpfr_set_prec (u, 15); 669 mpfr_set_str_binary (u, "0.101101000001100E-1"); 670 mpfr_set_prec (y, 92); 671 inexact = test_div (y, x, u, MPFR_RNDN); 672 if (inexact <= 0) 673 { 674 printf ("Wrong inexact flag for rnd=MPFR_RNDN(1): expected 1, got %d\n", 675 inexact); 676 mpfr_dump (x); 677 mpfr_dump (u); 678 mpfr_dump (y); 679 exit (1); 680 } 681 682 for (px=2; px<MAX_PREC; px++) 683 { 684 mpfr_set_prec (x, px); 685 mpfr_urandomb (x, RANDS); 686 for (pu=2; pu<=MAX_PREC; pu++) 687 { 688 mpfr_set_prec (u, pu); 689 do { mpfr_urandomb (u, RANDS); } while (mpfr_cmp_ui (u, 0) == 0); 690 { 691 py = MPFR_PREC_MIN + (randlimb () % (MAX_PREC - MPFR_PREC_MIN)); 692 mpfr_set_prec (y, py); 693 mpfr_set_prec (z, py + pu); 694 { 695 /* inexact is undefined for RNDF */ 696 rnd = RND_RAND_NO_RNDF (); 697 inexact = test_div (y, x, u, rnd); 698 if (mpfr_mul (z, y, u, rnd)) 699 { 700 printf ("z <- y * u should be exact\n"); 701 exit (1); 702 } 703 cmp = mpfr_cmp (z, x); 704 if (((inexact == 0) && (cmp != 0)) || 705 ((inexact > 0) && (cmp <= 0)) || 706 ((inexact < 0) && (cmp >= 0))) 707 { 708 printf ("Wrong inexact flag for rnd=%s\n", 709 mpfr_print_rnd_mode(rnd)); 710 printf ("expected %d, got %d\n", cmp, inexact); 711 printf ("x="); mpfr_dump (x); 712 printf ("u="); mpfr_dump (u); 713 printf ("y="); mpfr_dump (y); 714 printf ("y*u="); mpfr_dump (z); 715 exit (1); 716 } 717 } 718 } 719 } 720 } 721 722 mpfr_clear (x); 723 mpfr_clear (y); 724 mpfr_clear (z); 725 mpfr_clear (u); 726} 727 728static void 729check_special (void) 730{ 731 mpfr_t a, d, q; 732 mpfr_exp_t emax, emin; 733 int i; 734 735 mpfr_init2 (a, 100L); 736 mpfr_init2 (d, 100L); 737 mpfr_init2 (q, 100L); 738 739 /* 1/nan == nan */ 740 mpfr_set_ui (a, 1L, MPFR_RNDN); 741 MPFR_SET_NAN (d); 742 mpfr_clear_flags (); 743 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 744 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_NAN); 745 746 /* nan/1 == nan */ 747 MPFR_SET_NAN (a); 748 mpfr_set_ui (d, 1L, MPFR_RNDN); 749 mpfr_clear_flags (); 750 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 751 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_NAN); 752 753 /* +inf/1 == +inf */ 754 MPFR_SET_INF (a); 755 MPFR_SET_POS (a); 756 mpfr_set_ui (d, 1L, MPFR_RNDN); 757 mpfr_clear_flags (); 758 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 759 MPFR_ASSERTN (mpfr_inf_p (q)); 760 MPFR_ASSERTN (mpfr_sgn (q) > 0); 761 MPFR_ASSERTN (__gmpfr_flags == 0); 762 763 /* +inf/-1 == -inf */ 764 MPFR_SET_INF (a); 765 MPFR_SET_POS (a); 766 mpfr_set_si (d, -1, MPFR_RNDN); 767 mpfr_clear_flags (); 768 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 769 MPFR_ASSERTN (mpfr_inf_p (q)); 770 MPFR_ASSERTN (mpfr_sgn (q) < 0); 771 MPFR_ASSERTN (__gmpfr_flags == 0); 772 773 /* -inf/1 == -inf */ 774 MPFR_SET_INF (a); 775 MPFR_SET_NEG (a); 776 mpfr_set_ui (d, 1L, MPFR_RNDN); 777 mpfr_clear_flags (); 778 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 779 MPFR_ASSERTN (mpfr_inf_p (q)); 780 MPFR_ASSERTN (mpfr_sgn (q) < 0); 781 MPFR_ASSERTN (__gmpfr_flags == 0); 782 783 /* -inf/-1 == +inf */ 784 MPFR_SET_INF (a); 785 MPFR_SET_NEG (a); 786 mpfr_set_si (d, -1, MPFR_RNDN); 787 mpfr_clear_flags (); 788 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 789 MPFR_ASSERTN (mpfr_inf_p (q)); 790 MPFR_ASSERTN (mpfr_sgn (q) > 0); 791 MPFR_ASSERTN (__gmpfr_flags == 0); 792 793 /* 1/+inf == +0 */ 794 mpfr_set_ui (a, 1L, MPFR_RNDN); 795 MPFR_SET_INF (d); 796 MPFR_SET_POS (d); 797 mpfr_clear_flags (); 798 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 799 MPFR_ASSERTN (MPFR_IS_ZERO (q)); 800 MPFR_ASSERTN (MPFR_IS_POS (q)); 801 MPFR_ASSERTN (__gmpfr_flags == 0); 802 803 /* 1/-inf == -0 */ 804 mpfr_set_ui (a, 1L, MPFR_RNDN); 805 MPFR_SET_INF (d); 806 MPFR_SET_NEG (d); 807 mpfr_clear_flags (); 808 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 809 MPFR_ASSERTN (MPFR_IS_ZERO (q)); 810 MPFR_ASSERTN (MPFR_IS_NEG (q)); 811 MPFR_ASSERTN (__gmpfr_flags == 0); 812 813 /* -1/+inf == -0 */ 814 mpfr_set_si (a, -1, MPFR_RNDN); 815 MPFR_SET_INF (d); 816 MPFR_SET_POS (d); 817 mpfr_clear_flags (); 818 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 819 MPFR_ASSERTN (MPFR_IS_ZERO (q)); 820 MPFR_ASSERTN (MPFR_IS_NEG (q)); 821 MPFR_ASSERTN (__gmpfr_flags == 0); 822 823 /* -1/-inf == +0 */ 824 mpfr_set_si (a, -1, MPFR_RNDN); 825 MPFR_SET_INF (d); 826 MPFR_SET_NEG (d); 827 mpfr_clear_flags (); 828 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 829 MPFR_ASSERTN (MPFR_IS_ZERO (q)); 830 MPFR_ASSERTN (MPFR_IS_POS (q)); 831 MPFR_ASSERTN (__gmpfr_flags == 0); 832 833 /* 0/0 == nan */ 834 mpfr_set_ui (a, 0L, MPFR_RNDN); 835 mpfr_set_ui (d, 0L, MPFR_RNDN); 836 mpfr_clear_flags (); 837 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 838 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_NAN); 839 840 /* +inf/+inf == nan */ 841 MPFR_SET_INF (a); 842 MPFR_SET_POS (a); 843 MPFR_SET_INF (d); 844 MPFR_SET_POS (d); 845 mpfr_clear_flags (); 846 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 847 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_NAN); 848 849 /* 1/+0 = +inf */ 850 mpfr_set_ui (a, 1, MPFR_RNDZ); 851 mpfr_set_ui (d, 0, MPFR_RNDZ); 852 mpfr_clear_flags (); 853 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 854 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0); 855 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_DIVBY0); 856 857 /* 1/-0 = -inf */ 858 mpfr_set_ui (a, 1, MPFR_RNDZ); 859 mpfr_set_ui (d, 0, MPFR_RNDZ); 860 mpfr_neg (d, d, MPFR_RNDZ); 861 mpfr_clear_flags (); 862 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 863 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0); 864 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_DIVBY0); 865 866 /* -1/+0 = -inf */ 867 mpfr_set_si (a, -1, MPFR_RNDZ); 868 mpfr_set_ui (d, 0, MPFR_RNDZ); 869 mpfr_clear_flags (); 870 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 871 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0); 872 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_DIVBY0); 873 874 /* -1/-0 = +inf */ 875 mpfr_set_si (a, -1, MPFR_RNDZ); 876 mpfr_set_ui (d, 0, MPFR_RNDZ); 877 mpfr_neg (d, d, MPFR_RNDZ); 878 mpfr_clear_flags (); 879 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 880 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0); 881 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_DIVBY0); 882 883 /* +inf/+0 = +inf */ 884 MPFR_SET_INF (a); 885 MPFR_SET_POS (a); 886 mpfr_set_ui (d, 0, MPFR_RNDZ); 887 mpfr_clear_flags (); 888 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 889 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0); 890 MPFR_ASSERTN (__gmpfr_flags == 0); 891 892 /* +inf/-0 = -inf */ 893 MPFR_SET_INF (a); 894 MPFR_SET_POS (a); 895 mpfr_set_ui (d, 0, MPFR_RNDZ); 896 mpfr_neg (d, d, MPFR_RNDZ); 897 mpfr_clear_flags (); 898 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 899 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0); 900 MPFR_ASSERTN (__gmpfr_flags == 0); 901 902 /* -inf/+0 = -inf */ 903 MPFR_SET_INF (a); 904 MPFR_SET_NEG (a); 905 mpfr_set_ui (d, 0, MPFR_RNDZ); 906 mpfr_clear_flags (); 907 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 908 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0); 909 MPFR_ASSERTN (__gmpfr_flags == 0); 910 911 /* -inf/-0 = +inf */ 912 MPFR_SET_INF (a); 913 MPFR_SET_NEG (a); 914 mpfr_set_ui (d, 0, MPFR_RNDZ); 915 mpfr_neg (d, d, MPFR_RNDZ); 916 mpfr_clear_flags (); 917 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 918 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0); 919 MPFR_ASSERTN (__gmpfr_flags == 0); 920 921 /* check overflow */ 922 emax = mpfr_get_emax (); 923 set_emax (1); 924 mpfr_set_ui (a, 1, MPFR_RNDZ); 925 mpfr_set_ui (d, 1, MPFR_RNDZ); 926 mpfr_div_2ui (d, d, 1, MPFR_RNDZ); 927 mpfr_clear_flags (); 928 test_div (q, a, d, MPFR_RNDU); /* 1 / 0.5 = 2 -> overflow */ 929 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0); 930 MPFR_ASSERTN (__gmpfr_flags == (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT)); 931 set_emax (emax); 932 933 /* check underflow */ 934 emin = mpfr_get_emin (); 935 set_emin (-1); 936 mpfr_set_ui (a, 1, MPFR_RNDZ); 937 mpfr_div_2ui (a, a, 2, MPFR_RNDZ); 938 mpfr_set_prec (d, mpfr_get_prec (q) + 8); 939 for (i = -1; i <= 1; i++) 940 { 941 int sign; 942 943 /* Test 2^(-2) / (+/- (2 + eps)), with eps < 0, eps = 0, eps > 0. 944 -> underflow. 945 With div.c r5513, this test fails for eps > 0 in MPFR_RNDN. */ 946 mpfr_set_ui (d, 2, MPFR_RNDZ); 947 if (i < 0) 948 mpfr_nextbelow (d); 949 if (i > 0) 950 mpfr_nextabove (d); 951 for (sign = 0; sign <= 1; sign++) 952 { 953 mpfr_clear_flags (); 954 test_div (q, a, d, MPFR_RNDZ); /* result = 0 */ 955 MPFR_ASSERTN (__gmpfr_flags == 956 (MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT)); 957 MPFR_ASSERTN (sign ? MPFR_IS_NEG (q) : MPFR_IS_POS (q)); 958 MPFR_ASSERTN (MPFR_IS_ZERO (q)); 959 mpfr_clear_flags (); 960 test_div (q, a, d, MPFR_RNDN); /* result = 0 iff eps >= 0 */ 961 MPFR_ASSERTN (__gmpfr_flags == 962 (MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT)); 963 MPFR_ASSERTN (sign ? MPFR_IS_NEG (q) : MPFR_IS_POS (q)); 964 if (i < 0) 965 mpfr_nexttozero (q); 966 MPFR_ASSERTN (MPFR_IS_ZERO (q)); 967 mpfr_neg (d, d, MPFR_RNDN); 968 } 969 } 970 set_emin (emin); 971 972 mpfr_clear (a); 973 mpfr_clear (d); 974 mpfr_clear (q); 975} 976 977static void 978consistency (void) 979{ 980 mpfr_t x, y, z1, z2; 981 int i; 982 983 mpfr_inits (x, y, z1, z2, (mpfr_ptr) 0); 984 985 for (i = 0; i < 10000; i++) 986 { 987 mpfr_rnd_t rnd; 988 mpfr_prec_t px, py, pz, p; 989 int inex1, inex2; 990 991 /* inex is undefined for RNDF */ 992 rnd = RND_RAND_NO_RNDF (); 993 px = (randlimb () % 256) + 2; 994 py = (randlimb () % 128) + 2; 995 pz = (randlimb () % 256) + 2; 996 mpfr_set_prec (x, px); 997 mpfr_set_prec (y, py); 998 mpfr_set_prec (z1, pz); 999 mpfr_set_prec (z2, pz); 1000 mpfr_urandomb (x, RANDS); 1001 do 1002 mpfr_urandomb (y, RANDS); 1003 while (mpfr_zero_p (y)); 1004 inex1 = mpfr_div (z1, x, y, rnd); 1005 MPFR_ASSERTN (!MPFR_IS_NAN (z1)); 1006 p = MAX (MAX (px, py), pz); 1007 if (mpfr_prec_round (x, p, MPFR_RNDN) != 0 || 1008 mpfr_prec_round (y, p, MPFR_RNDN) != 0) 1009 { 1010 printf ("mpfr_prec_round error for i = %d\n", i); 1011 exit (1); 1012 } 1013 inex2 = mpfr_div (z2, x, y, rnd); 1014 MPFR_ASSERTN (!MPFR_IS_NAN (z2)); 1015 if (inex1 != inex2 || mpfr_cmp (z1, z2) != 0) 1016 { 1017 printf ("Consistency error for i = %d, rnd = %s\n", i, 1018 mpfr_print_rnd_mode (rnd)); 1019 printf ("inex1=%d inex2=%d\n", inex1, inex2); 1020 printf ("z1="); mpfr_dump (z1); 1021 printf ("z2="); mpfr_dump (z2); 1022 exit (1); 1023 } 1024 } 1025 1026 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0); 1027} 1028 1029/* Reported by Carl Witty on 2007-06-03 */ 1030static void 1031test_20070603 (void) 1032{ 1033 mpfr_t n, d, q, c; 1034 1035 mpfr_init2 (n, 128); 1036 mpfr_init2 (d, 128); 1037 mpfr_init2 (q, 31); 1038 mpfr_init2 (c, 31); 1039 1040 mpfr_set_str (n, "10384593717069655257060992206846485", 10, MPFR_RNDN); 1041 mpfr_set_str (d, "10384593717069655257060992206847132", 10, MPFR_RNDN); 1042 mpfr_div (q, n, d, MPFR_RNDU); 1043 1044 mpfr_set_ui (c, 1, MPFR_RNDN); 1045 if (mpfr_cmp (q, c) != 0) 1046 { 1047 printf ("Error in test_20070603\nGot "); 1048 mpfr_dump (q); 1049 printf ("instead of "); 1050 mpfr_dump (c); 1051 exit (1); 1052 } 1053 1054 /* same for 64-bit machines */ 1055 mpfr_set_prec (n, 256); 1056 mpfr_set_prec (d, 256); 1057 mpfr_set_prec (q, 63); 1058 mpfr_set_str (n, "822752278660603021077484591278675252491367930877209729029898240", 10, MPFR_RNDN); 1059 mpfr_set_str (d, "822752278660603021077484591278675252491367930877212507873738752", 10, MPFR_RNDN); 1060 mpfr_div (q, n, d, MPFR_RNDU); 1061 if (mpfr_cmp (q, c) != 0) 1062 { 1063 printf ("Error in test_20070603\nGot "); 1064 mpfr_dump (q); 1065 printf ("instead of "); 1066 mpfr_dump (c); 1067 exit (1); 1068 } 1069 1070 mpfr_clear (n); 1071 mpfr_clear (d); 1072 mpfr_clear (q); 1073 mpfr_clear (c); 1074} 1075 1076/* Bug found while adding tests for mpfr_cot */ 1077static void 1078test_20070628 (void) 1079{ 1080 mpfr_exp_t old_emax; 1081 mpfr_t x, y; 1082 int inex, err = 0; 1083 1084 old_emax = mpfr_get_emax (); 1085 set_emax (256); 1086 1087 mpfr_inits2 (53, x, y, (mpfr_ptr) 0); 1088 mpfr_set_si (x, -1, MPFR_RNDN); 1089 mpfr_set_si_2exp (y, 1, -256, MPFR_RNDN); 1090 mpfr_clear_flags (); 1091 inex = mpfr_div (x, x, y, MPFR_RNDD); 1092 if (MPFR_IS_POS (x) || ! mpfr_inf_p (x)) 1093 { 1094 printf ("Error in test_20070628: expected -Inf, got\n"); 1095 mpfr_dump (x); 1096 err++; 1097 } 1098 if (inex >= 0) 1099 { 1100 printf ("Error in test_20070628: expected inex < 0, got %d\n", inex); 1101 err++; 1102 } 1103 if (! mpfr_overflow_p ()) 1104 { 1105 printf ("Error in test_20070628: overflow flag is not set\n"); 1106 err++; 1107 } 1108 mpfr_clears (x, y, (mpfr_ptr) 0); 1109 set_emax (old_emax); 1110} 1111 1112/* Bug in mpfr_divhigh_n_basecase when all limbs of q (except the most 1113 significant one) are B-1 where B=2^GMP_NUMB_BITS. Since we truncate 1114 the divisor at each step, it might happen at some point that 1115 (np[n-1],np[n-2]) > (d1,d0), and not only the equality. 1116 Reported by Ricky Farr 1117 <https://sympa.inria.fr/sympa/arc/mpfr/2015-10/msg00023.html> 1118 To get a failure, a MPFR_DIVHIGH_TAB entry below the MPFR_DIV_THRESHOLD 1119 limit must have a value 0. With most mparam.h files, this cannot occur. To 1120 make the bug appear, one can configure MPFR with -DMPFR_TUNE_COVERAGE. */ 1121static void 1122test_20151023 (void) 1123{ 1124 mpfr_prec_t p; 1125 mpfr_t n, d, q, q0; 1126 int inex, i; 1127 1128 for (p = GMP_NUMB_BITS; p <= 2000; p++) 1129 { 1130 mpfr_init2 (n, 2*p); 1131 mpfr_init2 (d, p); 1132 mpfr_init2 (q, p); 1133 mpfr_init2 (q0, GMP_NUMB_BITS); 1134 1135 /* generate a random divisor of p bits */ 1136 do 1137 mpfr_urandomb (d, RANDS); 1138 while (mpfr_zero_p (d)); 1139 /* generate a random non-zero quotient of GMP_NUMB_BITS bits */ 1140 do 1141 mpfr_urandomb (q0, RANDS); 1142 while (mpfr_zero_p (q0)); 1143 /* zero-pad the quotient to p bits */ 1144 inex = mpfr_prec_round (q0, p, MPFR_RNDN); 1145 MPFR_ASSERTN(inex == 0); 1146 1147 for (i = 0; i < 3; i++) 1148 { 1149 /* i=0: try with the original quotient xxx000...000 1150 i=1: try with the original quotient minus one ulp 1151 i=2: try with the original quotient plus one ulp */ 1152 if (i == 1) 1153 mpfr_nextbelow (q0); 1154 else if (i == 2) 1155 { 1156 mpfr_nextabove (q0); 1157 mpfr_nextabove (q0); 1158 } 1159 1160 inex = mpfr_mul (n, d, q0, MPFR_RNDN); 1161 MPFR_ASSERTN(inex == 0); 1162 mpfr_nextabove (n); 1163 mpfr_div (q, n, d, MPFR_RNDN); 1164 if (! mpfr_equal_p (q, q0)) 1165 { 1166 printf ("Error in test_20151023 for p=%ld, rnd=RNDN, i=%d\n", 1167 (long) p, i); 1168 printf ("n="); mpfr_dump (n); 1169 printf ("d="); mpfr_dump (d); 1170 printf ("expected q0="); mpfr_dump (q0); 1171 printf ("got q="); mpfr_dump (q); 1172 exit (1); 1173 } 1174 1175 inex = mpfr_mul (n, d, q0, MPFR_RNDN); 1176 MPFR_ASSERTN(inex == 0); 1177 mpfr_nextbelow (n); 1178 mpfr_div (q, n, d, MPFR_RNDN); 1179 MPFR_ASSERTN(mpfr_cmp (q, q0) == 0); 1180 } 1181 1182 mpfr_clear (n); 1183 mpfr_clear (d); 1184 mpfr_clear (q); 1185 mpfr_clear (q0); 1186 } 1187} 1188 1189/* test a random division of p+extra bits divided by p+extra bits, 1190 with quotient of p bits only, where the p+extra bit approximation 1191 of the quotient is very near a rounding frontier. */ 1192static void 1193test_bad_aux (mpfr_prec_t p, mpfr_prec_t extra) 1194{ 1195 mpfr_t u, v, w, q0, q; 1196 1197 mpfr_init2 (u, p + extra); 1198 mpfr_init2 (v, p + extra); 1199 mpfr_init2 (w, p + extra); 1200 mpfr_init2 (q0, p); 1201 mpfr_init2 (q, p); 1202 do mpfr_urandomb (q0, RANDS); while (mpfr_zero_p (q0)); 1203 do mpfr_urandomb (v, RANDS); while (mpfr_zero_p (v)); 1204 1205 mpfr_set (w, q0, MPFR_RNDN); /* exact */ 1206 mpfr_nextabove (w); /* now w > q0 */ 1207 mpfr_mul (u, v, w, MPFR_RNDU); /* thus u > v*q0 */ 1208 mpfr_div (q, u, v, MPFR_RNDU); /* should have q > q0 */ 1209 MPFR_ASSERTN (mpfr_cmp (q, q0) > 0); 1210 mpfr_div (q, u, v, MPFR_RNDZ); /* should have q = q0 */ 1211 MPFR_ASSERTN (mpfr_cmp (q, q0) == 0); 1212 1213 mpfr_set (w, q0, MPFR_RNDN); /* exact */ 1214 mpfr_nextbelow (w); /* now w < q0 */ 1215 mpfr_mul (u, v, w, MPFR_RNDZ); /* thus u < v*q0 */ 1216 mpfr_div (q, u, v, MPFR_RNDZ); /* should have q < q0 */ 1217 MPFR_ASSERTN (mpfr_cmp (q, q0) < 0); 1218 mpfr_div (q, u, v, MPFR_RNDU); /* should have q = q0 */ 1219 MPFR_ASSERTN (mpfr_cmp (q, q0) == 0); 1220 1221 mpfr_clear (u); 1222 mpfr_clear (v); 1223 mpfr_clear (w); 1224 mpfr_clear (q0); 1225 mpfr_clear (q); 1226} 1227 1228static void 1229test_bad (void) 1230{ 1231 mpfr_prec_t p, extra; 1232 1233 for (p = MPFR_PREC_MIN; p <= 1024; p += 17) 1234 for (extra = 2; extra <= 64; extra++) 1235 test_bad_aux (p, extra); 1236} 1237 1238#define TEST_FUNCTION test_div 1239#define TWO_ARGS 1240#define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS) 1241#include "tgeneric.c" 1242 1243static void 1244test_extreme (void) 1245{ 1246 mpfr_t x, y, z; 1247 mpfr_exp_t emin, emax; 1248 mpfr_prec_t p[4] = { 8, 32, 64, 256 }; 1249 int xi, yi, zi, j, r; 1250 unsigned int flags, ex_flags; 1251 1252 emin = mpfr_get_emin (); 1253 emax = mpfr_get_emax (); 1254 1255 set_emin (MPFR_EMIN_MIN); 1256 set_emax (MPFR_EMAX_MAX); 1257 1258 for (xi = 0; xi < 4; xi++) 1259 { 1260 mpfr_init2 (x, p[xi]); 1261 mpfr_setmax (x, MPFR_EMAX_MAX); 1262 MPFR_ASSERTN (mpfr_check (x)); 1263 for (yi = 0; yi < 4; yi++) 1264 { 1265 mpfr_init2 (y, p[yi]); 1266 mpfr_setmin (y, MPFR_EMIN_MIN); 1267 for (j = 0; j < 2; j++) 1268 { 1269 MPFR_ASSERTN (mpfr_check (y)); 1270 for (zi = 0; zi < 4; zi++) 1271 { 1272 mpfr_init2 (z, p[zi]); 1273 RND_LOOP (r) 1274 { 1275 mpfr_clear_flags (); 1276 mpfr_div (z, x, y, (mpfr_rnd_t) r); 1277 flags = __gmpfr_flags; 1278 MPFR_ASSERTN (mpfr_check (z)); 1279 ex_flags = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT; 1280 if (flags != ex_flags) 1281 { 1282 printf ("Bad flags in test_extreme on z = a/b" 1283 " with %s and\n", 1284 mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 1285 printf ("a = "); 1286 mpfr_dump (x); 1287 printf ("b = "); 1288 mpfr_dump (y); 1289 printf ("Expected flags:"); 1290 flags_out (ex_flags); 1291 printf ("Got flags: "); 1292 flags_out (flags); 1293 printf ("z = "); 1294 mpfr_dump (z); 1295 exit (1); 1296 } 1297 mpfr_clear_flags (); 1298 mpfr_div (z, y, x, (mpfr_rnd_t) r); 1299 flags = __gmpfr_flags; 1300 MPFR_ASSERTN (mpfr_check (z)); 1301 ex_flags = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT; 1302 if (flags != ex_flags) 1303 { 1304 printf ("Bad flags in test_extreme on z = a/b" 1305 " with %s and\n", 1306 mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 1307 printf ("a = "); 1308 mpfr_dump (y); 1309 printf ("b = "); 1310 mpfr_dump (x); 1311 printf ("Expected flags:"); 1312 flags_out (ex_flags); 1313 printf ("Got flags: "); 1314 flags_out (flags); 1315 printf ("z = "); 1316 mpfr_dump (z); 1317 exit (1); 1318 } 1319 } 1320 mpfr_clear (z); 1321 } /* zi */ 1322 mpfr_nextabove (y); 1323 } /* j */ 1324 mpfr_clear (y); 1325 } /* yi */ 1326 mpfr_clear (x); 1327 } /* xi */ 1328 1329 set_emin (emin); 1330 set_emax (emax); 1331} 1332 1333static void 1334testall_rndf (mpfr_prec_t pmax) 1335{ 1336 mpfr_t a, b, c, d; 1337 mpfr_prec_t pa, pb, pc; 1338 1339 for (pa = MPFR_PREC_MIN; pa <= pmax; pa++) 1340 { 1341 mpfr_init2 (a, pa); 1342 mpfr_init2 (d, pa); 1343 for (pb = MPFR_PREC_MIN; pb <= pmax; pb++) 1344 { 1345 mpfr_init2 (b, pb); 1346 mpfr_set_ui (b, 1, MPFR_RNDN); 1347 while (mpfr_cmp_ui (b, 2) < 0) 1348 { 1349 for (pc = MPFR_PREC_MIN; pc <= pmax; pc++) 1350 { 1351 mpfr_init2 (c, pc); 1352 mpfr_set_ui (c, 1, MPFR_RNDN); 1353 while (mpfr_cmp_ui (c, 2) < 0) 1354 { 1355 mpfr_div (a, b, c, MPFR_RNDF); 1356 mpfr_div (d, b, c, MPFR_RNDD); 1357 if (!mpfr_equal_p (a, d)) 1358 { 1359 mpfr_div (d, b, c, MPFR_RNDU); 1360 if (!mpfr_equal_p (a, d)) 1361 { 1362 printf ("Error: mpfr_div(a,b,c,RNDF) does not " 1363 "match RNDD/RNDU\n"); 1364 printf ("b="); mpfr_dump (b); 1365 printf ("c="); mpfr_dump (c); 1366 printf ("a="); mpfr_dump (a); 1367 exit (1); 1368 } 1369 } 1370 mpfr_nextabove (c); 1371 } 1372 mpfr_clear (c); 1373 } 1374 mpfr_nextabove (b); 1375 } 1376 mpfr_clear (b); 1377 } 1378 mpfr_clear (a); 1379 mpfr_clear (d); 1380 } 1381} 1382 1383static void 1384test_mpfr_divsp2 (void) 1385{ 1386 mpfr_t u, v, q; 1387 1388 /* test to exercise r2 = v1 in mpfr_divsp2 */ 1389 mpfr_init2 (u, 128); 1390 mpfr_init2 (v, 128); 1391 mpfr_init2 (q, 83); 1392 1393 mpfr_set_str (u, "286677858044426991425771529092412636160", 10, MPFR_RNDN); 1394 mpfr_set_str (v, "241810647971575979588130185988987264768", 10, MPFR_RNDN); 1395 mpfr_div (q, u, v, MPFR_RNDN); 1396 mpfr_set_str (u, "5732952910203749289426944", 10, MPFR_RNDN); 1397 mpfr_div_2ui (u, u, 82, MPFR_RNDN); 1398 MPFR_ASSERTN(mpfr_equal_p (q, u)); 1399 1400 mpfr_clear (u); 1401 mpfr_clear (v); 1402 mpfr_clear (q); 1403} 1404 1405/* Assertion failure in r10769 with --enable-assert --enable-gmp-internals 1406 (same failure in tatan on a similar example). */ 1407static void 1408test_20160831 (void) 1409{ 1410 mpfr_t u, v, q; 1411 1412 mpfr_inits2 (124, u, v, q, (mpfr_ptr) 0); 1413 1414 mpfr_set_ui (u, 1, MPFR_RNDN); 1415 mpfr_set_str (v, "0x40000000000000005", 16, MPFR_RNDN); 1416 mpfr_div (q, u, v, MPFR_RNDN); 1417 mpfr_set_str (u, "0xfffffffffffffffecp-134", 16, MPFR_RNDN); 1418 MPFR_ASSERTN (mpfr_equal_p (q, u)); 1419 1420 mpfr_set_prec (u, 128); 1421 mpfr_set_prec (v, 128); 1422 mpfr_set_str (u, "186127091671619245460026015088243485690", 10, MPFR_RNDN); 1423 mpfr_set_str (v, "205987256581218236405412302590110119580", 10, MPFR_RNDN); 1424 mpfr_div (q, u, v, MPFR_RNDN); 1425 mpfr_set_str (u, "19217137613667309953639458782352244736", 10, MPFR_RNDN); 1426 mpfr_div_2ui (u, u, 124, MPFR_RNDN); 1427 MPFR_ASSERTN (mpfr_equal_p (q, u)); 1428 1429 mpfr_clears (u, v, q, (mpfr_ptr) 0); 1430} 1431 1432/* With r11138, on a 64-bit machine: 1433 div.c:128: MPFR assertion failed: qx >= __gmpfr_emin 1434*/ 1435static void 1436test_20170104 (void) 1437{ 1438 mpfr_t u, v, q; 1439 mpfr_exp_t emin; 1440 1441 emin = mpfr_get_emin (); 1442 set_emin (-42); 1443 1444 mpfr_init2 (u, 12); 1445 mpfr_init2 (v, 12); 1446 mpfr_init2 (q, 11); 1447 mpfr_set_str_binary (u, "0.111111111110E-29"); 1448 mpfr_set_str_binary (v, "0.111111111111E14"); 1449 mpfr_div (q, u, v, MPFR_RNDN); 1450 mpfr_clears (u, v, q, (mpfr_ptr) 0); 1451 1452 set_emin (emin); 1453} 1454 1455/* With r11140, on a 64-bit machine with GMP_CHECK_RANDOMIZE=1484406128: 1456 Consistency error for i = 2577 1457*/ 1458static void 1459test_20170105 (void) 1460{ 1461 mpfr_t x, y, z, t; 1462 1463 mpfr_init2 (x, 138); 1464 mpfr_init2 (y, 6); 1465 mpfr_init2 (z, 128); 1466 mpfr_init2 (t, 128); 1467 mpfr_set_str_binary (x, "0.100110111001001000101111010010011101111110111111110001110100000001110111010100111010100011101010110000010100000011100100110101101011000000E-6"); 1468 mpfr_set_str_binary (y, "0.100100E-2"); 1469 /* up to exponents, x/y is exactly 367625553447399614694201910705139062483, 1470 which has 129 bits, thus we are in the round-to-nearest-even case, and since 1471 the penultimate bit of x/y is 1, we should round upwards */ 1472 mpfr_set_str_binary (t, "0.10001010010010010000110110010110111111111100011011101010000000000110101000010001011110011011010000111010000000001100101101101010E-3"); 1473 mpfr_div (z, x, y, MPFR_RNDN); 1474 MPFR_ASSERTN(mpfr_equal_p (z, t)); 1475 mpfr_clears (x, y, z, t, (mpfr_ptr) 0); 1476} 1477 1478/* The real cause of the mpfr_ttanh failure from the non-regression test 1479 added in tests/ttanh.c@11993 was actually due to a bug in mpfr_div, as 1480 this can be seen by comparing the logs of the 3.1 branch and the trunk 1481 r11992 with MPFR_LOG_ALL=1 MPFR_LOG_PREC=50 on this particular test 1482 (this was noticed because adding this test to the 3.1 branch did not 1483 yield a failure like in the trunk, though the mpfr_ttanh code did not 1484 change until r11993). This was the bug actually fixed in r12002. 1485*/ 1486static void 1487test_20171219 (void) 1488{ 1489 mpfr_t x, y, z, t; 1490 1491 mpfr_inits2 (126, x, y, z, t, (mpfr_ptr) 0); 1492 mpfr_set_str_binary (x, "0.111000000000000111100000000000011110000000000001111000000000000111100000000000011110000000000001111000000000000111100000000000E1"); 1493 /* x = 36347266450315671364380109803814927 / 2^114 */ 1494 mpfr_add_ui (y, x, 2, MPFR_RNDN); 1495 /* y = 77885641318594292392624080437575695 / 2^114 */ 1496 mpfr_div (z, x, y, MPFR_RNDN); 1497 mpfr_set_ui_2exp (t, 3823, -13, MPFR_RNDN); 1498 MPFR_ASSERTN (mpfr_equal_p (z, t)); 1499 mpfr_clears (x, y, z, t, (mpfr_ptr) 0); 1500} 1501 1502#if !defined(MPFR_GENERIC_ABI) && GMP_NUMB_BITS == 64 1503/* exercise mpfr_div2_approx */ 1504static void 1505test_mpfr_div2_approx (unsigned long n) 1506{ 1507 mpfr_t x, y, z; 1508 1509 mpfr_init2 (x, 113); 1510 mpfr_init2 (y, 113); 1511 mpfr_init2 (z, 113); 1512 while (n--) 1513 { 1514 mpfr_urandomb (x, RANDS); 1515 mpfr_urandomb (y, RANDS); 1516 mpfr_div (z, x, y, MPFR_RNDN); 1517 } 1518 mpfr_clear (x); 1519 mpfr_clear (y); 1520 mpfr_clear (z); 1521} 1522#endif 1523 1524/* bug found in ttan with GMP_CHECK_RANDOMIZE=1514257254 */ 1525static void 1526bug20171218 (void) 1527{ 1528 mpfr_t s, c; 1529 mpfr_init2 (s, 124); 1530 mpfr_init2 (c, 124); 1531 mpfr_set_str_binary (s, "-0.1110000111100001111000011110000111100001111000011110000111100001111000011110000111100001111000011110000111100001111000011110E0"); 1532 mpfr_set_str_binary (c, "0.1111000011110000111100001111000011110000111100001111000011110000111100001111000011110000111100001111000011110000111100001111E-1"); 1533 mpfr_div (c, s, c, MPFR_RNDN); 1534 mpfr_set_str_binary (s, "-1.111000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"); 1535 MPFR_ASSERTN(mpfr_equal_p (c, s)); 1536 mpfr_clear (s); 1537 mpfr_clear (c); 1538} 1539 1540/* Extended test based on a bug found with flint-arb test suite with a 1541 32-bit ABI: https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=888459 1542 Division of the form: (1 - 2^(-pa)) / (1 - 2^(-pb)). 1543 The result is compared to the one obtained by increasing the precision of 1544 the divisor (without changing its value, so that both results should be 1545 equal). For all of these tests, a failure may occur in r12126 only with 1546 pb=GMP_NUMB_BITS and MPFR_RNDN (and some particular values of pa and pc). 1547 This bug was introduced by r9086, where mpfr_div uses mpfr_div_ui when 1548 the divisor has only one limb. 1549*/ 1550static void 1551bug20180126 (void) 1552{ 1553 mpfr_t a, b1, b2, c1, c2; 1554 int pa, i, j, pc, sa, sb, r, inex1, inex2; 1555 1556 for (pa = 100; pa < 800; pa += 11) 1557 for (i = 1; i <= 4; i++) 1558 for (j = -2; j <= 2; j++) 1559 { 1560 int pb = GMP_NUMB_BITS * i + j; 1561 1562 if (pb > pa) 1563 continue; 1564 1565 mpfr_inits2 (pa, a, b1, (mpfr_ptr) 0); 1566 mpfr_inits2 (pb, b2, (mpfr_ptr) 0); 1567 1568 mpfr_set_ui (a, 1, MPFR_RNDN); 1569 mpfr_nextbelow (a); /* 1 - 2^(-pa) */ 1570 mpfr_set_ui (b2, 1, MPFR_RNDN); 1571 mpfr_nextbelow (b2); /* 1 - 2^(-pb) */ 1572 inex1 = mpfr_set (b1, b2, MPFR_RNDN); 1573 MPFR_ASSERTN (inex1 == 0); 1574 1575 for (pc = 32; pc <= 320; pc += 32) 1576 { 1577 mpfr_inits2 (pc, c1, c2, (mpfr_ptr) 0); 1578 1579 for (sa = 0; sa < 2; sa++) 1580 { 1581 for (sb = 0; sb < 2; sb++) 1582 { 1583 RND_LOOP_NO_RNDF (r) 1584 { 1585 MPFR_ASSERTN (mpfr_equal_p (b1, b2)); 1586 inex1 = mpfr_div (c1, a, b1, (mpfr_rnd_t) r); 1587 inex2 = mpfr_div (c2, a, b2, (mpfr_rnd_t) r); 1588 1589 if (! mpfr_equal_p (c1, c2) || 1590 ! SAME_SIGN (inex1, inex2)) 1591 { 1592 printf ("Error in bug20180126 for " 1593 "pa=%d pb=%d pc=%d sa=%d sb=%d %s\n", 1594 pa, pb, pc, sa, sb, 1595 mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 1596 printf ("inex1 = %d, c1 = ", inex1); 1597 mpfr_dump (c1); 1598 printf ("inex2 = %d, c2 = ", inex2); 1599 mpfr_dump (c2); 1600 exit (1); 1601 } 1602 } 1603 1604 mpfr_neg (b1, b1, MPFR_RNDN); 1605 mpfr_neg (b2, b2, MPFR_RNDN); 1606 } /* sb */ 1607 1608 mpfr_neg (a, a, MPFR_RNDN); 1609 } /* sa */ 1610 1611 mpfr_clears (c1, c2, (mpfr_ptr) 0); 1612 } /* pc */ 1613 1614 mpfr_clears (a, b1, b2, (mpfr_ptr) 0); 1615 } /* j */ 1616} 1617 1618static void 1619coverage (mpfr_prec_t pmax) 1620{ 1621 mpfr_prec_t p; 1622 1623 for (p = MPFR_PREC_MIN; p <= pmax; p++) 1624 { 1625 int inex; 1626 mpfr_t q, u, v; 1627 1628 mpfr_init2 (q, p); 1629 mpfr_init2 (u, p); 1630 mpfr_init2 (v, p); 1631 1632 /* exercise case qx < emin */ 1633 mpfr_set_ui_2exp (u, 1, mpfr_get_emin (), MPFR_RNDN); 1634 mpfr_set_ui (v, 4, MPFR_RNDN); 1635 1636 mpfr_clear_flags (); 1637 /* u/v = 2^(emin-2), should be rounded to +0 for RNDN */ 1638 inex = mpfr_div (q, u, v, MPFR_RNDN); 1639 MPFR_ASSERTN(inex < 0); 1640 MPFR_ASSERTN(mpfr_zero_p (q) && mpfr_signbit (q) == 0); 1641 MPFR_ASSERTN(mpfr_underflow_p ()); 1642 1643 mpfr_clear_flags (); 1644 /* u/v = 2^(emin-2), should be rounded to 2^(emin-1) for RNDU */ 1645 inex = mpfr_div (q, u, v, MPFR_RNDU); 1646 MPFR_ASSERTN(inex > 0); 1647 MPFR_ASSERTN(mpfr_cmp_ui_2exp (q, 1, mpfr_get_emin () - 1) == 0); 1648 MPFR_ASSERTN(mpfr_underflow_p ()); 1649 1650 mpfr_clear_flags (); 1651 /* u/v = 2^(emin-2), should be rounded to +0 for RNDZ */ 1652 inex = mpfr_div (q, u, v, MPFR_RNDZ); 1653 MPFR_ASSERTN(inex < 0); 1654 MPFR_ASSERTN(mpfr_zero_p (q) && mpfr_signbit (q) == 0); 1655 MPFR_ASSERTN(mpfr_underflow_p ()); 1656 1657 if (p == 1) 1658 goto end_of_loop; 1659 1660 mpfr_set_ui_2exp (u, 1, mpfr_get_emin (), MPFR_RNDN); 1661 mpfr_nextbelow (u); /* u = (1-2^(-p))*2^emin */ 1662 mpfr_set_ui (v, 2, MPFR_RNDN); 1663 1664 mpfr_clear_flags (); 1665 /* u/v = (1-2^(-p))*2^(emin-1), will round to 2^(emin-1) for RNDN */ 1666 inex = mpfr_div (q, u, v, MPFR_RNDN); 1667 MPFR_ASSERTN(inex > 0); 1668 MPFR_ASSERTN(mpfr_cmp_ui_2exp (q, 1, mpfr_get_emin () - 1) == 0); 1669 MPFR_ASSERTN(mpfr_underflow_p ()); 1670 1671 mpfr_clear_flags (); 1672 /* u/v should round to 2^(emin-1) for RNDU */ 1673 inex = mpfr_div (q, u, v, MPFR_RNDU); 1674 MPFR_ASSERTN(inex > 0); 1675 MPFR_ASSERTN(mpfr_cmp_ui_2exp (q, 1, mpfr_get_emin () - 1) == 0); 1676 MPFR_ASSERTN(mpfr_underflow_p ()); 1677 1678 mpfr_clear_flags (); 1679 /* u/v should round to +0 for RNDZ */ 1680 inex = mpfr_div (q, u, v, MPFR_RNDZ); 1681 MPFR_ASSERTN(inex < 0); 1682 MPFR_ASSERTN(mpfr_zero_p (q) && mpfr_signbit (q) == 0); 1683 MPFR_ASSERTN(mpfr_underflow_p ()); 1684 1685 end_of_loop: 1686 mpfr_clear (q); 1687 mpfr_clear (u); 1688 mpfr_clear (v); 1689 } 1690} 1691 1692/* coverage for case usize >= n + n in Mulders' algorithm */ 1693static void 1694coverage2 (void) 1695{ 1696 mpfr_prec_t p; 1697 mpfr_t q, u, v, t, w; 1698 int inex, inex2; 1699 1700 p = MPFR_DIV_THRESHOLD * GMP_NUMB_BITS; 1701 mpfr_init2 (q, p); 1702 mpfr_init2 (u, 2 * p + 3 * GMP_NUMB_BITS); 1703 mpfr_init2 (v, p); 1704 do mpfr_urandomb (u, RANDS); while (mpfr_zero_p (u)); 1705 do mpfr_urandomb (v, RANDS); while (mpfr_zero_p (v)); 1706 inex = mpfr_div (q, u, v, MPFR_RNDN); 1707 mpfr_init2 (t, mpfr_get_prec (u)); 1708 mpfr_init2 (w, mpfr_get_prec (u)); 1709 inex2 = mpfr_mul (t, q, v, MPFR_RNDN); 1710 MPFR_ASSERTN(inex2 == 0); 1711 if (inex == 0) /* check q*v = u */ 1712 MPFR_ASSERTN(mpfr_equal_p (u, t)); 1713 else 1714 { 1715 if (inex > 0) 1716 mpfr_nextbelow (q); 1717 else 1718 mpfr_nextabove (q); 1719 inex2 = mpfr_mul (w, q, v, MPFR_RNDN); 1720 MPFR_ASSERTN(inex2 == 0); 1721 inex2 = mpfr_sub (t, t, u, MPFR_RNDN); 1722 MPFR_ASSERTN(inex2 == 0); 1723 inex2 = mpfr_sub (w, w, u, MPFR_RNDN); 1724 MPFR_ASSERTN(inex2 == 0); 1725 MPFR_ASSERTN(mpfr_cmpabs (t, w) <= 0); 1726 if (mpfr_cmpabs (t, w) == 0) /* even rule: significand of q should now 1727 be odd */ 1728 MPFR_ASSERTN(mpfr_min_prec (q) == mpfr_get_prec (q)); 1729 } 1730 1731 mpfr_clear (q); 1732 mpfr_clear (u); 1733 mpfr_clear (v); 1734 mpfr_clear (t); 1735 mpfr_clear (w); 1736} 1737 1738int 1739main (int argc, char *argv[]) 1740{ 1741 tests_start_mpfr (); 1742 1743 coverage (1024); 1744 coverage2 (); 1745 bug20180126 (); 1746 bug20171218 (); 1747 testall_rndf (9); 1748 test_20170105 (); 1749 check_inexact (); 1750 check_hard (); 1751 check_special (); 1752 check_lowr (); 1753 check_float (); /* checks single precision */ 1754 check_double (); 1755 check_convergence (); 1756 check_64 (); 1757 1758 check4("4.0","4.503599627370496e15", MPFR_RNDZ, 62, 1759 "0.10000000000000000000000000000000000000000000000000000000000000E-49"); 1760 check4("1.0","2.10263340267725788209e+187", MPFR_RNDU, 65, 1761 "0.11010011111001101011111001100111110100000001101001111100111000000E-622"); 1762 check4("2.44394909079968374564e-150", "2.10263340267725788209e+187",MPFR_RNDU, 1763 65, 1764 "0.11010011111001101011111001100111110100000001101001111100111000000E-1119"); 1765 1766 consistency (); 1767 test_20070603 (); 1768 test_20070628 (); 1769 test_20151023 (); 1770 test_20160831 (); 1771 test_20170104 (); 1772 test_20171219 (); 1773 test_generic (MPFR_PREC_MIN, 800, 50); 1774 test_bad (); 1775 test_extreme (); 1776 test_mpfr_divsp2 (); 1777#if !defined(MPFR_GENERIC_ABI) && GMP_NUMB_BITS == 64 1778 test_mpfr_div2_approx (1000000); 1779#endif 1780 1781 tests_end_mpfr (); 1782 return 0; 1783} 1784