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