1/* Test file for mpfr_div. 2 3Copyright 1999, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#include <stdio.h> 24#include <stdlib.h> 25 26#include "mpfr-test.h" 27 28#ifdef CHECK_EXTERNAL 29static int 30test_div (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) 31{ 32 int res; 33 int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_number_p (c); 34 if (ok) 35 { 36 mpfr_print_raw (b); 37 printf (" "); 38 mpfr_print_raw (c); 39 } 40 res = mpfr_div (a, b, c, rnd_mode); 41 if (ok) 42 { 43 printf (" "); 44 mpfr_print_raw (a); 45 printf ("\n"); 46 } 47 return res; 48} 49#else 50#define test_div mpfr_div 51#endif 52 53#define check53(n, d, rnd, res) check4(n, d, rnd, 53, res) 54 55/* return 0 iff a and b are of the same sign */ 56static int 57inex_cmp (int a, int b) 58{ 59 if (a > 0) 60 return (b > 0) ? 0 : 1; 61 else if (a == 0) 62 return (b == 0) ? 0 : 1; 63 else 64 return (b < 0) ? 0 : 1; 65} 66 67static void 68check4 (const char *Ns, const char *Ds, mpfr_rnd_t rnd_mode, int p, 69 const char *Qs) 70{ 71 mpfr_t q, n, d; 72 73 mpfr_inits2 (p, q, n, d, (mpfr_ptr) 0); 74 mpfr_set_str1 (n, Ns); 75 mpfr_set_str1 (d, Ds); 76 test_div(q, n, d, rnd_mode); 77 if (mpfr_cmp_str (q, Qs, ((p==53) ? 10 : 2), MPFR_RNDN) ) 78 { 79 printf ("mpfr_div failed for n=%s, d=%s, p=%d, rnd_mode=%s\n", 80 Ns, Ds, p, mpfr_print_rnd_mode (rnd_mode)); 81 printf ("got ");mpfr_print_binary(q); 82 mpfr_set_str (q, Qs, ((p==53) ? 10 : 2), MPFR_RNDN); 83 printf("\nexpected "); mpfr_print_binary(q); 84 putchar('\n'); 85 exit (1); 86 } 87 mpfr_clears (q, n, d, (mpfr_ptr) 0); 88} 89 90static void 91check24 (const char *Ns, const char *Ds, mpfr_rnd_t rnd_mode, const char *Qs) 92{ 93 mpfr_t q, n, d; 94 95 mpfr_inits2 (24, q, n, d, (mpfr_ptr) 0); 96 97 mpfr_set_str1 (n, Ns); 98 mpfr_set_str1 (d, Ds); 99 test_div(q, n, d, rnd_mode); 100 if (mpfr_cmp_str1 (q, Qs) ) 101 { 102 printf ("mpfr_div failed for n=%s, d=%s, prec=24, rnd_mode=%s\n", 103 Ns, Ds, mpfr_print_rnd_mode(rnd_mode)); 104 printf ("expected quotient is %s, got ", Qs); 105 mpfr_out_str(stdout,10,0,q, MPFR_RNDN); putchar('\n'); 106 exit (1); 107 } 108 mpfr_clears (q, n, d, (mpfr_ptr) 0); 109} 110 111/* the following examples come from the paper "Number-theoretic Test 112 Generation for Directed Rounding" from Michael Parks, Table 2 */ 113static void 114check_float(void) 115{ 116 check24("70368760954880.0", "8388609.0", MPFR_RNDN, "8.388609e6"); 117 check24("140737479966720.0", "16777213.0", MPFR_RNDN, "8.388609e6"); 118 check24("70368777732096.0", "8388611.0", MPFR_RNDN, "8.388609e6"); 119 check24("105553133043712.0", "12582911.0", MPFR_RNDN, "8.38861e6"); 120 /* the exponent for the following example was forgotten in 121 the Arith'14 version of Parks' paper */ 122 check24 ("12582913.0", "12582910.0", MPFR_RNDN, "1.000000238"); 123 check24 ("105553124655104.0", "12582910.0", MPFR_RNDN, "8388610.0"); 124 check24("140737479966720.0", "8388609.0", MPFR_RNDN, "1.6777213e7"); 125 check24("70368777732096.0", "8388609.0", MPFR_RNDN, "8.388611e6"); 126 check24("105553133043712.0", "8388610.0", MPFR_RNDN, "1.2582911e7"); 127 check24("105553124655104.0", "8388610.0", MPFR_RNDN, "1.258291e7"); 128 129 check24("70368760954880.0", "8388609.0", MPFR_RNDZ, "8.388608e6"); 130 check24("140737479966720.0", "16777213.0", MPFR_RNDZ, "8.388609e6"); 131 check24("70368777732096.0", "8388611.0", MPFR_RNDZ, "8.388608e6"); 132 check24("105553133043712.0", "12582911.0", MPFR_RNDZ, "8.38861e6"); 133 check24("12582913.0", "12582910.0", MPFR_RNDZ, "1.000000238"); 134 check24 ("105553124655104.0", "12582910.0", MPFR_RNDZ, "8388610.0"); 135 check24("140737479966720.0", "8388609.0", MPFR_RNDZ, "1.6777213e7"); 136 check24("70368777732096.0", "8388609.0", MPFR_RNDZ, "8.38861e6"); 137 check24("105553133043712.0", "8388610.0", MPFR_RNDZ, "1.2582911e7"); 138 check24("105553124655104.0", "8388610.0", MPFR_RNDZ, "1.258291e7"); 139 140 check24("70368760954880.0", "8388609.0", MPFR_RNDU, "8.388609e6"); 141 check24("140737479966720.0", "16777213.0", MPFR_RNDU, "8.38861e6"); 142 check24("70368777732096.0", "8388611.0", MPFR_RNDU, "8.388609e6"); 143 check24("105553133043712.0", "12582911.0", MPFR_RNDU, "8.388611e6"); 144 check24("12582913.0", "12582910.0", MPFR_RNDU, "1.000000357"); 145 check24 ("105553124655104.0", "12582910.0", MPFR_RNDU, "8388611.0"); 146 check24("140737479966720.0", "8388609.0", MPFR_RNDU, "1.6777214e7"); 147 check24("70368777732096.0", "8388609.0", MPFR_RNDU, "8.388611e6"); 148 check24("105553133043712.0", "8388610.0", MPFR_RNDU, "1.2582912e7"); 149 check24("105553124655104.0", "8388610.0", MPFR_RNDU, "1.2582911e7"); 150 151 check24("70368760954880.0", "8388609.0", MPFR_RNDD, "8.388608e6"); 152 check24("140737479966720.0", "16777213.0", MPFR_RNDD, "8.388609e6"); 153 check24("70368777732096.0", "8388611.0", MPFR_RNDD, "8.388608e6"); 154 check24("105553133043712.0", "12582911.0", MPFR_RNDD, "8.38861e6"); 155 check24("12582913.0", "12582910.0", MPFR_RNDD, "1.000000238"); 156 check24 ("105553124655104.0", "12582910.0", MPFR_RNDD, "8388610.0"); 157 check24("140737479966720.0", "8388609.0", MPFR_RNDD, "1.6777213e7"); 158 check24("70368777732096.0", "8388609.0", MPFR_RNDD, "8.38861e6"); 159 check24("105553133043712.0", "8388610.0", MPFR_RNDD, "1.2582911e7"); 160 check24("105553124655104.0", "8388610.0", MPFR_RNDD, "1.258291e7"); 161 162 check24("70368760954880.0", "8388609.0", MPFR_RNDA, "8.388609e6"); 163} 164 165static void 166check_double(void) 167{ 168 check53("0.0", "1.0", MPFR_RNDZ, "0.0"); 169 check53("-7.4988969224688591e63", "4.8816866450288732e306", MPFR_RNDD, 170 "-1.5361282826510687291e-243"); 171 check53("-1.33225773037748601769e+199", "3.63449540676937123913e+79", 172 MPFR_RNDZ, "-3.6655920045905428978e119"); 173 check53("9.89438396044940256501e-134", "5.93472984109987421717e-67",MPFR_RNDU, 174 "1.6672003992376663654e-67"); 175 check53("9.89438396044940256501e-134", "5.93472984109987421717e-67",MPFR_RNDA, 176 "1.6672003992376663654e-67"); 177 check53("9.89438396044940256501e-134", "-5.93472984109987421717e-67", 178 MPFR_RNDU, "-1.6672003992376663654e-67"); 179 check53("-4.53063926135729747564e-308", "7.02293374921793516813e-84", 180 MPFR_RNDD, "-6.4512060388748850857e-225"); 181 check53("6.25089225176473806123e-01","-2.35527154824420243364e-230", 182 MPFR_RNDD, "-2.6540006635008291192e229"); 183 check53("6.25089225176473806123e-01","-2.35527154824420243364e-230", 184 MPFR_RNDA, "-2.6540006635008291192e229"); 185 check53("6.52308934689126e15", "-1.62063546601505417497e273", MPFR_RNDN, 186 "-4.0250194961676020848e-258"); 187 check53("1.04636807108079349236e-189", "3.72295730823253012954e-292", 188 MPFR_RNDZ, "2.810583051186143125e102"); 189 /* problems found by Kevin under HP-PA */ 190 check53 ("2.861044553323177e-136", "-1.1120354257068143e+45", MPFR_RNDZ, 191 "-2.5727998292003016e-181"); 192 check53 ("-4.0559157245809205e-127", "-1.1237723844524865e+77", MPFR_RNDN, 193 "3.6091968273068081e-204"); 194 check53 ("-1.8177943561493235e-93", "-8.51233984260364e-104", MPFR_RNDU, 195 "2.1354814184595821e+10"); 196} 197 198static void 199check_64(void) 200{ 201 mpfr_t x,y,z; 202 203 mpfr_inits2 (64, x, y, z, (mpfr_ptr) 0); 204 205 mpfr_set_str_binary(x, "1.00100100110110101001010010101111000001011100100101010000000000E54"); 206 mpfr_set_str_binary(y, "1.00000000000000000000000000000000000000000000000000000000000000E584"); 207 test_div(z, x, y, MPFR_RNDU); 208 if (mpfr_cmp_str (z, "0.1001001001101101010010100101011110000010111001001010100000000000E-529", 2, MPFR_RNDN)) 209 { 210 printf("Error for tdiv for MPFR_RNDU and p=64\nx="); 211 mpfr_print_binary(x); 212 printf("\ny="); 213 mpfr_print_binary(y); 214 printf("\ngot "); 215 mpfr_print_binary(z); 216 printf("\nexpected 0.1001001001101101010010100101011110000010111001001010100000000000E-529\n"); 217 exit(1); 218 } 219 220 mpfr_clears (x, y, z, (mpfr_ptr) 0); 221} 222 223static void 224check_convergence (void) 225{ 226 mpfr_t x, y; int i, j; 227 228 mpfr_init2(x, 130); 229 mpfr_set_str_binary(x, "0.1011111101011010101000001010011111101000011100011101010011111011000011001010000000111100100111110011001010110100100001001000111001E6944"); 230 mpfr_init2(y, 130); 231 mpfr_set_ui(y, 5, MPFR_RNDN); 232 test_div(x, x, y, MPFR_RNDD); /* exact division */ 233 234 mpfr_set_prec(x, 64); 235 mpfr_set_prec(y, 64); 236 mpfr_set_str_binary(x, "0.10010010011011010100101001010111100000101110010010101E55"); 237 mpfr_set_str_binary(y, "0.1E585"); 238 test_div(x, x, y, MPFR_RNDN); 239 mpfr_set_str_binary(y, "0.10010010011011010100101001010111100000101110010010101E-529"); 240 if (mpfr_cmp (x, y)) 241 { 242 printf ("Error in mpfr_div for prec=64, rnd=MPFR_RNDN\n"); 243 printf ("got "); mpfr_print_binary(x); puts (""); 244 printf ("instead of "); mpfr_print_binary(y); puts (""); 245 exit(1); 246 } 247 248 for (i=32; i<=64; i+=32) 249 { 250 mpfr_set_prec(x, i); 251 mpfr_set_prec(y, i); 252 mpfr_set_ui(x, 1, MPFR_RNDN); 253 RND_LOOP(j) 254 { 255 mpfr_set_ui (y, 1, MPFR_RNDN); 256 test_div (y, x, y, (mpfr_rnd_t) j); 257 if (mpfr_cmp_ui (y, 1)) 258 { 259 printf ("mpfr_div failed for x=1.0, y=1.0, prec=%d rnd=%s\n", 260 i, mpfr_print_rnd_mode ((mpfr_rnd_t) j)); 261 printf ("got "); mpfr_print_binary(y); puts (""); 262 exit (1); 263 } 264 } 265 } 266 267 mpfr_clear (x); 268 mpfr_clear (y); 269} 270 271#define KMAX 10000 272 273/* given y = o(x/u), x, u, find the inexact flag by 274 multiplying y by u */ 275static int 276get_inexact (mpfr_t y, mpfr_t x, mpfr_t u) 277{ 278 mpfr_t xx; 279 int inex; 280 mpfr_init2 (xx, mpfr_get_prec (y) + mpfr_get_prec (u)); 281 mpfr_mul (xx, y, u, MPFR_RNDN); /* exact */ 282 inex = mpfr_cmp (xx, x); 283 mpfr_clear (xx); 284 return inex; 285} 286 287static void 288check_hard (void) 289{ 290 mpfr_t u, v, q, q2; 291 mpfr_prec_t precu, precv, precq; 292 int rnd; 293 int inex, inex2, i, j; 294 295 mpfr_init (q); 296 mpfr_init (q2); 297 mpfr_init (u); 298 mpfr_init (v); 299 300 for (precq = MPFR_PREC_MIN; precq <= 64; precq ++) 301 { 302 mpfr_set_prec (q, precq); 303 mpfr_set_prec (q2, precq + 1); 304 for (j = 0; j < 2; j++) 305 { 306 if (j == 0) 307 { 308 do 309 { 310 mpfr_urandomb (q2, RANDS); 311 } 312 while (mpfr_cmp_ui (q2, 0) == 0); 313 } 314 else /* use q2=1 */ 315 mpfr_set_ui (q2, 1, MPFR_RNDN); 316 for (precv = precq; precv <= 10 * precq; precv += precq) 317 { 318 mpfr_set_prec (v, precv); 319 do 320 { 321 mpfr_urandomb (v, RANDS); 322 } 323 while (mpfr_cmp_ui (v, 0) == 0); 324 for (precu = precq; precu <= 10 * precq; precu += precq) 325 { 326 mpfr_set_prec (u, precu); 327 mpfr_mul (u, v, q2, MPFR_RNDN); 328 mpfr_nextbelow (u); 329 for (i = 0; i <= 2; i++) 330 { 331 RND_LOOP(rnd) 332 { 333 inex = test_div (q, u, v, (mpfr_rnd_t) rnd); 334 inex2 = get_inexact (q, u, v); 335 if (inex_cmp (inex, inex2)) 336 { 337 printf ("Wrong inexact flag for rnd=%s: expected %d, got %d\n", 338 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), inex2, inex); 339 printf ("u= "); mpfr_dump (u); 340 printf ("v= "); mpfr_dump (v); 341 printf ("q= "); mpfr_dump (q); 342 mpfr_set_prec (q2, precq + precv); 343 mpfr_mul (q2, q, v, MPFR_RNDN); 344 printf ("q*v="); mpfr_dump (q2); 345 exit (1); 346 } 347 } 348 mpfr_nextabove (u); 349 } 350 } 351 } 352 } 353 } 354 355 mpfr_clear (q); 356 mpfr_clear (q2); 357 mpfr_clear (u); 358 mpfr_clear (v); 359} 360 361static void 362check_lowr (void) 363{ 364 mpfr_t x, y, z, z2, z3, tmp; 365 int k, c, c2; 366 367 368 mpfr_init2 (x, 1000); 369 mpfr_init2 (y, 100); 370 mpfr_init2 (tmp, 850); 371 mpfr_init2 (z, 10); 372 mpfr_init2 (z2, 10); 373 mpfr_init2 (z3, 50); 374 375 for (k = 1; k < KMAX; k++) 376 { 377 do 378 { 379 mpfr_urandomb (z, RANDS); 380 } 381 while (mpfr_cmp_ui (z, 0) == 0); 382 do 383 { 384 mpfr_urandomb (tmp, RANDS); 385 } 386 while (mpfr_cmp_ui (tmp, 0) == 0); 387 mpfr_mul (x, z, tmp, MPFR_RNDN); /* exact */ 388 c = test_div (z2, x, tmp, MPFR_RNDN); 389 390 if (c || mpfr_cmp (z2, z)) 391 { 392 printf ("Error in mpfr_div rnd=MPFR_RNDN\n"); 393 printf ("got "); mpfr_print_binary(z2); puts (""); 394 printf ("instead of "); mpfr_print_binary(z); puts (""); 395 printf ("inex flag = %d, expected 0\n", c); 396 exit (1); 397 } 398 } 399 400 /* x has still precision 1000, z precision 10, and tmp prec 850 */ 401 mpfr_set_prec (z2, 9); 402 for (k = 1; k < KMAX; k++) 403 { 404 mpfr_urandomb (z, RANDS); 405 do 406 { 407 mpfr_urandomb (tmp, RANDS); 408 } 409 while (mpfr_cmp_ui (tmp, 0) == 0); 410 mpfr_mul (x, z, tmp, MPFR_RNDN); /* exact */ 411 c = test_div (z2, x, tmp, MPFR_RNDN); 412 /* since z2 has one less bit that z, either the division is exact 413 if z is representable on 9 bits, or we have an even round case */ 414 415 c2 = get_inexact (z2, x, tmp); 416 if ((mpfr_cmp (z2, z) == 0 && c) || inex_cmp (c, c2)) 417 { 418 printf ("Error in mpfr_div rnd=MPFR_RNDN\n"); 419 printf ("got "); mpfr_print_binary(z2); puts (""); 420 printf ("instead of "); mpfr_print_binary(z); puts (""); 421 printf ("inex flag = %d, expected %d\n", c, c2); 422 exit (1); 423 } 424 else if (c == 2) 425 { 426 mpfr_nexttoinf (z); 427 if (mpfr_cmp(z2, z)) 428 { 429 printf ("Error in mpfr_div [even rnd?] rnd=MPFR_RNDN\n"); 430 printf ("Dividing "); 431 printf ("got "); mpfr_print_binary(z2); puts (""); 432 printf ("instead of "); mpfr_print_binary(z); puts (""); 433 printf ("inex flag = %d\n", 1); 434 exit (1); 435 } 436 } 437 else if (c == -2) 438 { 439 mpfr_nexttozero (z); 440 if (mpfr_cmp(z2, z)) 441 { 442 printf ("Error in mpfr_div [even rnd?] rnd=MPFR_RNDN\n"); 443 printf ("Dividing "); 444 printf ("got "); mpfr_print_binary(z2); puts (""); 445 printf ("instead of "); mpfr_print_binary(z); puts (""); 446 printf ("inex flag = %d\n", 1); 447 exit (1); 448 } 449 } 450 } 451 452 mpfr_set_prec(x, 1000); 453 mpfr_set_prec(y, 100); 454 mpfr_set_prec(tmp, 850); 455 mpfr_set_prec(z, 10); 456 mpfr_set_prec(z2, 10); 457 458 /* almost exact divisions */ 459 for (k = 1; k < KMAX; k++) 460 { 461 do 462 { 463 mpfr_urandomb (z, RANDS); 464 } 465 while (mpfr_cmp_ui (z, 0) == 0); 466 do 467 { 468 mpfr_urandomb (tmp, RANDS); 469 } 470 while (mpfr_cmp_ui (tmp, 0) == 0); 471 mpfr_mul(x, z, tmp, MPFR_RNDN); 472 mpfr_set(y, tmp, MPFR_RNDD); 473 mpfr_nexttoinf (x); 474 475 c = test_div(z2, x, y, MPFR_RNDD); 476 test_div(z3, x, y, MPFR_RNDD); 477 mpfr_set(z, z3, MPFR_RNDD); 478 479 if (c != -1 || mpfr_cmp(z2, z)) 480 { 481 printf ("Error in mpfr_div rnd=MPFR_RNDD\n"); 482 printf ("got "); mpfr_print_binary(z2); puts (""); 483 printf ("instead of "); mpfr_print_binary(z); puts (""); 484 printf ("inex flag = %d\n", c); 485 exit (1); 486 } 487 488 mpfr_set (y, tmp, MPFR_RNDU); 489 test_div (z3, x, y, MPFR_RNDU); 490 mpfr_set (z, z3, MPFR_RNDU); 491 c = test_div (z2, x, y, MPFR_RNDU); 492 if (c != 1 || mpfr_cmp (z2, z)) 493 { 494 printf ("Error in mpfr_div rnd=MPFR_RNDU\n"); 495 printf ("u="); mpfr_dump (x); 496 printf ("v="); mpfr_dump (y); 497 printf ("got "); mpfr_print_binary (z2); puts (""); 498 printf ("instead of "); mpfr_print_binary (z); puts (""); 499 printf ("inex flag = %d\n", c); 500 exit (1); 501 } 502 } 503 504 mpfr_clear (x); 505 mpfr_clear (y); 506 mpfr_clear (z); 507 mpfr_clear (z2); 508 mpfr_clear (z3); 509 mpfr_clear (tmp); 510} 511 512#define MAX_PREC 128 513 514static void 515check_inexact (void) 516{ 517 mpfr_t x, y, z, u; 518 mpfr_prec_t px, py, pu; 519 int inexact, cmp; 520 mpfr_rnd_t rnd; 521 522 mpfr_init (x); 523 mpfr_init (y); 524 mpfr_init (z); 525 mpfr_init (u); 526 527 mpfr_set_prec (x, 28); 528 mpfr_set_prec (y, 28); 529 mpfr_set_prec (z, 1023); 530 mpfr_set_str_binary (x, "0.1000001001101101111100010011E0"); 531 mpfr_set_str (z, "48284762641021308813686974720835219181653367326353400027913400579340343320519877153813133510034402932651132854764198688352364361009429039801248971901380781746767119334993621199563870113045276395603170432175354501451429471578325545278975153148347684600400321033502982713296919861760382863826626093689036010394", 10, MPFR_RNDN); 532 mpfr_div (x, x, z, MPFR_RNDN); 533 mpfr_set_str_binary (y, "0.1111001011001101001001111100E-1023"); 534 if (mpfr_cmp (x, y)) 535 { 536 printf ("Error in mpfr_div for prec=28, RNDN\n"); 537 printf ("Expected "); mpfr_dump (y); 538 printf ("Got "); mpfr_dump (x); 539 exit (1); 540 } 541 542 mpfr_set_prec (x, 53); 543 mpfr_set_str_binary (x, "0.11101100110010100011011000000100001111011111110010101E0"); 544 mpfr_set_prec (u, 127); 545 mpfr_set_str_binary (u, "0.1000001100110110110101110110101101111000110000001111111110000000011111001010110100110010111111111101000001011011101011101101000E-2"); 546 mpfr_set_prec (y, 95); 547 inexact = test_div (y, x, u, MPFR_RNDN); 548 if (inexact != (cmp = get_inexact (y, x, u))) 549 { 550 printf ("Wrong inexact flag (0): expected %d, got %d\n", cmp, inexact); 551 printf ("x="); mpfr_out_str (stdout, 10, 99, x, MPFR_RNDN); printf ("\n"); 552 printf ("u="); mpfr_out_str (stdout, 10, 99, u, MPFR_RNDN); printf ("\n"); 553 printf ("y="); mpfr_out_str (stdout, 10, 99, y, MPFR_RNDN); printf ("\n"); 554 exit (1); 555 } 556 557 mpfr_set_prec (x, 33); 558 mpfr_set_str_binary (x, "0.101111100011011101010011101100001E0"); 559 mpfr_set_prec (u, 2); 560 mpfr_set_str_binary (u, "0.1E0"); 561 mpfr_set_prec (y, 28); 562 if ((inexact = test_div (y, x, u, MPFR_RNDN) >= 0)) 563 { 564 printf ("Wrong inexact flag (1): expected -1, got %d\n", 565 inexact); 566 exit (1); 567 } 568 569 mpfr_set_prec (x, 129); 570 mpfr_set_str_binary (x, "0.111110101111001100000101011100101100110011011101010001000110110101100101000010000001110110100001101010001010100010001111001101010E-2"); 571 mpfr_set_prec (u, 15); 572 mpfr_set_str_binary (u, "0.101101000001100E-1"); 573 mpfr_set_prec (y, 92); 574 if ((inexact = test_div (y, x, u, MPFR_RNDN)) <= 0) 575 { 576 printf ("Wrong inexact flag for rnd=MPFR_RNDN(1): expected 1, got %d\n", 577 inexact); 578 mpfr_dump (x); 579 mpfr_dump (u); 580 mpfr_dump (y); 581 exit (1); 582 } 583 584 for (px=2; px<MAX_PREC; px++) 585 { 586 mpfr_set_prec (x, px); 587 mpfr_urandomb (x, RANDS); 588 for (pu=2; pu<=MAX_PREC; pu++) 589 { 590 mpfr_set_prec (u, pu); 591 do { mpfr_urandomb (u, RANDS); } while (mpfr_cmp_ui (u, 0) == 0); 592 { 593 py = MPFR_PREC_MIN + (randlimb () % (MAX_PREC - MPFR_PREC_MIN)); 594 mpfr_set_prec (y, py); 595 mpfr_set_prec (z, py + pu); 596 { 597 rnd = RND_RAND (); 598 inexact = test_div (y, x, u, rnd); 599 if (mpfr_mul (z, y, u, rnd)) 600 { 601 printf ("z <- y * u should be exact\n"); 602 exit (1); 603 } 604 cmp = mpfr_cmp (z, x); 605 if (((inexact == 0) && (cmp != 0)) || 606 ((inexact > 0) && (cmp <= 0)) || 607 ((inexact < 0) && (cmp >= 0))) 608 { 609 printf ("Wrong inexact flag for rnd=%s\n", 610 mpfr_print_rnd_mode(rnd)); 611 printf ("expected %d, got %d\n", cmp, inexact); 612 printf ("x="); mpfr_print_binary (x); puts (""); 613 printf ("u="); mpfr_print_binary (u); puts (""); 614 printf ("y="); mpfr_print_binary (y); puts (""); 615 printf ("y*u="); mpfr_print_binary (z); puts (""); 616 exit (1); 617 } 618 } 619 } 620 } 621 } 622 623 mpfr_clear (x); 624 mpfr_clear (y); 625 mpfr_clear (z); 626 mpfr_clear (u); 627} 628 629static void 630check_nan (void) 631{ 632 mpfr_t a, d, q; 633 mpfr_exp_t emax, emin; 634 int i; 635 636 mpfr_init2 (a, 100L); 637 mpfr_init2 (d, 100L); 638 mpfr_init2 (q, 100L); 639 640 /* 1/nan == nan */ 641 mpfr_set_ui (a, 1L, MPFR_RNDN); 642 MPFR_SET_NAN (d); 643 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 644 MPFR_ASSERTN (mpfr_nan_p (q)); 645 646 /* nan/1 == nan */ 647 MPFR_SET_NAN (a); 648 mpfr_set_ui (d, 1L, MPFR_RNDN); 649 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 650 MPFR_ASSERTN (mpfr_nan_p (q)); 651 652 /* +inf/1 == +inf */ 653 MPFR_SET_INF (a); 654 MPFR_SET_POS (a); 655 mpfr_set_ui (d, 1L, MPFR_RNDN); 656 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 657 MPFR_ASSERTN (mpfr_inf_p (q)); 658 MPFR_ASSERTN (mpfr_sgn (q) > 0); 659 660 /* 1/+inf == 0 */ 661 mpfr_set_ui (a, 1L, MPFR_RNDN); 662 MPFR_SET_INF (d); 663 MPFR_SET_POS (d); 664 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 665 MPFR_ASSERTN (mpfr_number_p (q)); 666 MPFR_ASSERTN (mpfr_sgn (q) == 0); 667 668 /* 0/0 == nan */ 669 mpfr_set_ui (a, 0L, MPFR_RNDN); 670 mpfr_set_ui (d, 0L, MPFR_RNDN); 671 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 672 MPFR_ASSERTN (mpfr_nan_p (q)); 673 674 /* +inf/+inf == nan */ 675 MPFR_SET_INF (a); 676 MPFR_SET_POS (a); 677 MPFR_SET_INF (d); 678 MPFR_SET_POS (d); 679 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 680 MPFR_ASSERTN (mpfr_nan_p (q)); 681 682 /* 1/+0 = +Inf */ 683 mpfr_set_ui (a, 1, MPFR_RNDZ); 684 mpfr_set_ui (d, 0, MPFR_RNDZ); 685 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 686 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0); 687 688 /* 1/-0 = -Inf */ 689 mpfr_set_ui (a, 1, MPFR_RNDZ); 690 mpfr_set_ui (d, 0, MPFR_RNDZ); 691 mpfr_neg (d, d, MPFR_RNDZ); 692 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 693 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0); 694 695 /* -1/+0 = -Inf */ 696 mpfr_set_si (a, -1, MPFR_RNDZ); 697 mpfr_set_ui (d, 0, MPFR_RNDZ); 698 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 699 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0); 700 701 /* -1/-0 = +Inf */ 702 mpfr_set_si (a, -1, MPFR_RNDZ); 703 mpfr_set_ui (d, 0, MPFR_RNDZ); 704 mpfr_neg (d, d, MPFR_RNDZ); 705 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */ 706 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0); 707 708 /* check overflow */ 709 emax = mpfr_get_emax (); 710 set_emax (1); 711 mpfr_set_ui (a, 1, MPFR_RNDZ); 712 mpfr_set_ui (d, 1, MPFR_RNDZ); 713 mpfr_div_2exp (d, d, 1, MPFR_RNDZ); 714 test_div (q, a, d, MPFR_RNDU); /* 1 / 0.5 = 2 -> overflow */ 715 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0); 716 set_emax (emax); 717 718 /* check underflow */ 719 emin = mpfr_get_emin (); 720 set_emin (-1); 721 mpfr_set_ui (a, 1, MPFR_RNDZ); 722 mpfr_div_2exp (a, a, 2, MPFR_RNDZ); 723 mpfr_set_prec (d, mpfr_get_prec (q) + 8); 724 for (i = -1; i <= 1; i++) 725 { 726 int sign; 727 728 /* Test 2^(-2) / (+/- (2 + eps)), with eps < 0, eps = 0, eps > 0. 729 -> underflow. 730 With div.c r5513, this test fails for eps > 0 in MPFR_RNDN. */ 731 mpfr_set_ui (d, 2, MPFR_RNDZ); 732 if (i < 0) 733 mpfr_nextbelow (d); 734 if (i > 0) 735 mpfr_nextabove (d); 736 for (sign = 0; sign <= 1; sign++) 737 { 738 mpfr_clear_flags (); 739 test_div (q, a, d, MPFR_RNDZ); /* result = 0 */ 740 MPFR_ASSERTN (mpfr_underflow_p ()); 741 MPFR_ASSERTN (sign ? MPFR_IS_NEG (q) : MPFR_IS_POS (q)); 742 MPFR_ASSERTN (MPFR_IS_ZERO (q)); 743 mpfr_clear_flags (); 744 test_div (q, a, d, MPFR_RNDN); /* result = 0 iff eps >= 0 */ 745 MPFR_ASSERTN (mpfr_underflow_p ()); 746 MPFR_ASSERTN (sign ? MPFR_IS_NEG (q) : MPFR_IS_POS (q)); 747 if (i < 0) 748 mpfr_nexttozero (q); 749 MPFR_ASSERTN (MPFR_IS_ZERO (q)); 750 mpfr_neg (d, d, MPFR_RNDN); 751 } 752 } 753 set_emin (emin); 754 755 mpfr_clear (a); 756 mpfr_clear (d); 757 mpfr_clear (q); 758} 759 760static void 761consistency (void) 762{ 763 mpfr_t x, y, z1, z2; 764 int i; 765 766 mpfr_inits (x, y, z1, z2, (mpfr_ptr) 0); 767 768 for (i = 0; i < 10000; i++) 769 { 770 mpfr_rnd_t rnd; 771 mpfr_prec_t px, py, pz, p; 772 int inex1, inex2; 773 774 rnd = RND_RAND (); 775 px = (randlimb () % 256) + 2; 776 py = (randlimb () % 128) + 2; 777 pz = (randlimb () % 256) + 2; 778 mpfr_set_prec (x, px); 779 mpfr_set_prec (y, py); 780 mpfr_set_prec (z1, pz); 781 mpfr_set_prec (z2, pz); 782 mpfr_urandomb (x, RANDS); 783 do 784 mpfr_urandomb (y, RANDS); 785 while (mpfr_zero_p (y)); 786 inex1 = mpfr_div (z1, x, y, rnd); 787 MPFR_ASSERTN (!MPFR_IS_NAN (z1)); 788 p = MAX (MAX (px, py), pz); 789 if (mpfr_prec_round (x, p, MPFR_RNDN) != 0 || 790 mpfr_prec_round (y, p, MPFR_RNDN) != 0) 791 { 792 printf ("mpfr_prec_round error for i = %d\n", i); 793 exit (1); 794 } 795 inex2 = mpfr_div (z2, x, y, rnd); 796 MPFR_ASSERTN (!MPFR_IS_NAN (z2)); 797 if (inex1 != inex2 || mpfr_cmp (z1, z2) != 0) 798 { 799 printf ("Consistency error for i = %d\n", i); 800 exit (1); 801 } 802 } 803 804 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0); 805} 806 807/* Reported by Carl Witty on 2007-06-03 */ 808static void 809test_20070603 (void) 810{ 811 mpfr_t n, d, q, c; 812 813 mpfr_init2 (n, 128); 814 mpfr_init2 (d, 128); 815 mpfr_init2 (q, 31); 816 mpfr_init2 (c, 31); 817 818 mpfr_set_str (n, "10384593717069655257060992206846485", 10, MPFR_RNDN); 819 mpfr_set_str (d, "10384593717069655257060992206847132", 10, MPFR_RNDN); 820 mpfr_div (q, n, d, MPFR_RNDU); 821 822 mpfr_set_ui (c, 1, MPFR_RNDN); 823 if (mpfr_cmp (q, c) != 0) 824 { 825 printf ("Error in test_20070603\nGot "); 826 mpfr_dump (q); 827 printf ("instead of "); 828 mpfr_dump (c); 829 exit (1); 830 } 831 832 /* same for 64-bit machines */ 833 mpfr_set_prec (n, 256); 834 mpfr_set_prec (d, 256); 835 mpfr_set_prec (q, 63); 836 mpfr_set_str (n, "822752278660603021077484591278675252491367930877209729029898240", 10, MPFR_RNDN); 837 mpfr_set_str (d, "822752278660603021077484591278675252491367930877212507873738752", 10, MPFR_RNDN); 838 mpfr_div (q, n, d, MPFR_RNDU); 839 if (mpfr_cmp (q, c) != 0) 840 { 841 printf ("Error in test_20070603\nGot "); 842 mpfr_dump (q); 843 printf ("instead of "); 844 mpfr_dump (c); 845 exit (1); 846 } 847 848 mpfr_clear (n); 849 mpfr_clear (d); 850 mpfr_clear (q); 851 mpfr_clear (c); 852} 853 854/* Bug found while adding tests for mpfr_cot */ 855static void 856test_20070628 (void) 857{ 858 mpfr_exp_t old_emax; 859 mpfr_t x, y; 860 int inex, err = 0; 861 862 old_emax = mpfr_get_emax (); 863 864 if (mpfr_set_emax (256)) 865 { 866 printf ("Can't change exponent range\n"); 867 exit (1); 868 } 869 870 mpfr_inits2 (53, x, y, (mpfr_ptr) 0); 871 mpfr_set_si (x, -1, MPFR_RNDN); 872 mpfr_set_si_2exp (y, 1, -256, MPFR_RNDN); 873 mpfr_clear_flags (); 874 inex = mpfr_div (x, x, y, MPFR_RNDD); 875 if (MPFR_SIGN (x) >= 0 || ! mpfr_inf_p (x)) 876 { 877 printf ("Error in test_20070628: expected -Inf, got\n"); 878 mpfr_dump (x); 879 err++; 880 } 881 if (inex >= 0) 882 { 883 printf ("Error in test_20070628: expected inex < 0, got %d\n", inex); 884 err++; 885 } 886 if (! mpfr_overflow_p ()) 887 { 888 printf ("Error in test_20070628: overflow flag is not set\n"); 889 err++; 890 } 891 mpfr_clears (x, y, (mpfr_ptr) 0); 892 mpfr_set_emax (old_emax); 893} 894 895#define TEST_FUNCTION test_div 896#define TWO_ARGS 897#define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS) 898#include "tgeneric.c" 899 900int 901main (int argc, char *argv[]) 902{ 903 tests_start_mpfr (); 904 905 check_inexact (); 906 check_hard (); 907 check_nan (); 908 check_lowr (); 909 check_float (); /* checks single precision */ 910 check_double (); 911 check_convergence (); 912 check_64 (); 913 914 check4("4.0","4.503599627370496e15", MPFR_RNDZ, 62, 915 "0.10000000000000000000000000000000000000000000000000000000000000E-49"); 916 check4("1.0","2.10263340267725788209e+187", MPFR_RNDU, 65, 917 "0.11010011111001101011111001100111110100000001101001111100111000000E-622"); 918 check4("2.44394909079968374564e-150", "2.10263340267725788209e+187",MPFR_RNDU, 919 65, 920 "0.11010011111001101011111001100111110100000001101001111100111000000E-1119"); 921 922 consistency (); 923 test_20070603 (); 924 test_20070628 (); 925 test_generic (2, 800, 50); 926 927 tests_end_mpfr (); 928 return 0; 929} 930