1 /* Test file for mpfr_sub. 2 3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramel 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 <stdio.h> 24 #include <stdlib.h> 25 26 #include "mpfr-test.h" 27 28 #ifdef CHECK_EXTERNAL 29 static int 30 test_sub (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_sub (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_sub mpfr_sub 51 #endif 52 53 static void 54 check_diverse (void) 55 { 56 mpfr_t x, y, z; 57 int inexact; 58 59 mpfr_init (x); 60 mpfr_init (y); 61 mpfr_init (z); 62 63 /* check corner case cancel=0, but add_exp=1 */ 64 mpfr_set_prec (x, 2); 65 mpfr_set_prec (y, 4); 66 mpfr_set_prec (z, 2); 67 mpfr_setmax (y, __gmpfr_emax); 68 mpfr_set_str_binary (z, "0.1E-10"); /* tiny */ 69 test_sub (x, y, z, MPFR_RNDN); /* should round to 2^emax, i.e. overflow */ 70 if (!mpfr_inf_p (x) || mpfr_sgn (x) < 0) 71 { 72 printf ("Error in mpfr_sub(a,b,c,RNDN) for b=maxfloat, prec(a)<prec(b), c tiny\n"); 73 exit (1); 74 } 75 76 /* other coverage test */ 77 mpfr_set_prec (x, 2); 78 mpfr_set_prec (y, 2); 79 mpfr_set_prec (z, 2); 80 mpfr_set_ui (y, 1, MPFR_RNDN); 81 mpfr_set_si (z, -2, MPFR_RNDN); 82 test_sub (x, y, z, MPFR_RNDD); 83 if (mpfr_cmp_ui (x, 3)) 84 { 85 printf ("Error in mpfr_sub(1,-2,RNDD)\n"); 86 exit (1); 87 } 88 89 mpfr_set_prec (x, 288); 90 mpfr_set_prec (y, 288); 91 mpfr_set_prec (z, 288); 92 mpfr_set_str_binary (y, "0.111000110011000001000111101010111011110011101001101111111110000011100101000001001010110010101010011001010100000001110011110001010101101010001011101110100100001011110100110000101101100011010001001011011010101010000010001101001000110010010111111011110001111101001000101101001100101100101000E80"); 93 mpfr_set_str_binary (z, "0.100001111111101001011010001100110010100111001110000110011101001011010100001000000100111011010110110010000000000010101101011000010000110001110010100001100101011100100100001011000100011110000001010101000100011101001000010111100000111000111011001000100100011000100000010010111000000100100111E-258"); 94 inexact = test_sub (x, y, z, MPFR_RNDN); 95 if (inexact <= 0) 96 { 97 printf ("Wrong inexact flag for prec=288\n"); 98 exit (1); 99 } 100 101 mpfr_set_prec (x, 32); 102 mpfr_set_prec (y, 63); 103 mpfr_set_prec (z, 63); 104 mpfr_set_str_binary (x, "0.101101111011011100100100100111E31"); 105 mpfr_set_str_binary (y, "0.111110010010100100110101101010001001100101110001000101110111111E-1"); 106 test_sub (z, x, y, MPFR_RNDN); 107 mpfr_set_str_binary (y, "0.1011011110110111001001001001101100000110110101101100101001011E31"); 108 if (mpfr_cmp (z, y)) 109 { 110 printf ("Error in mpfr_sub (5)\n"); 111 printf ("expected "); mpfr_print_binary (y); puts (""); 112 printf ("got "); mpfr_print_binary (z); puts (""); 113 exit (1); 114 } 115 116 mpfr_set_prec (y, 63); 117 mpfr_set_prec (z, 63); 118 mpfr_set_str_binary (y, "0.1011011110110111001001001001101100000110110101101100101001011E31"); 119 mpfr_sub_ui (z, y, 1541116494, MPFR_RNDN); 120 mpfr_set_str_binary (y, "-0.11111001001010010011010110101E-1"); 121 if (mpfr_cmp (z, y)) 122 { 123 printf ("Error in mpfr_sub (7)\n"); 124 printf ("expected "); mpfr_print_binary (y); puts (""); 125 printf ("got "); mpfr_print_binary (z); puts (""); 126 exit (1); 127 } 128 129 mpfr_set_prec (y, 63); 130 mpfr_set_prec (z, 63); 131 mpfr_set_str_binary (y, "0.1011011110110111001001001001101100000110110101101100101001011E31"); 132 mpfr_sub_ui (z, y, 1541116494, MPFR_RNDN); 133 mpfr_set_str_binary (y, "-0.11111001001010010011010110101E-1"); 134 if (mpfr_cmp (z, y)) 135 { 136 printf ("Error in mpfr_sub (6)\n"); 137 printf ("expected "); mpfr_print_binary (y); puts (""); 138 printf ("got "); mpfr_print_binary (z); puts (""); 139 exit (1); 140 } 141 142 mpfr_set_prec (x, 32); 143 mpfr_set_prec (y, 32); 144 mpfr_set_str_binary (x, "0.10110111101001110100100101111000E0"); 145 mpfr_set_str_binary (y, "0.10001100100101000100110111000100E0"); 146 if ((inexact = test_sub (x, x, y, MPFR_RNDN))) 147 { 148 printf ("Wrong inexact flag (2): got %d instead of 0\n", inexact); 149 exit (1); 150 } 151 152 mpfr_set_prec (x, 32); 153 mpfr_set_prec (y, 32); 154 mpfr_set_str_binary (x, "0.11111000110111011000100111011010E0"); 155 mpfr_set_str_binary (y, "0.10011111101111000100001000000000E-3"); 156 if ((inexact = test_sub (x, x, y, MPFR_RNDN))) 157 { 158 printf ("Wrong inexact flag (1): got %d instead of 0\n", inexact); 159 exit (1); 160 } 161 162 mpfr_set_prec (x, 33); 163 mpfr_set_prec (y, 33); 164 mpfr_set_ui (x, 1, MPFR_RNDN); 165 mpfr_set_str_binary (y, "-0.1E-32"); 166 mpfr_add (x, x, y, MPFR_RNDN); 167 mpfr_set_str_binary (y, "0.111111111111111111111111111111111E0"); 168 if (mpfr_cmp (x, y)) 169 { 170 printf ("Error in mpfr_sub (1 - 1E-33) with prec=33\n"); 171 printf ("Expected "); mpfr_print_binary (y); puts (""); 172 printf ("got "); mpfr_print_binary (x); puts (""); 173 exit (1); 174 } 175 176 mpfr_set_prec (x, 32); 177 mpfr_set_prec (y, 33); 178 mpfr_set_ui (x, 1, MPFR_RNDN); 179 mpfr_set_str_binary (y, "-0.1E-32"); 180 mpfr_add (x, x, y, MPFR_RNDN); 181 if (mpfr_cmp_ui (x, 1)) 182 { 183 printf ("Error in mpfr_sub (1 - 1E-33) with prec=32\n"); 184 printf ("Expected 1.0, got "); mpfr_print_binary (x); puts (""); 185 exit (1); 186 } 187 188 mpfr_set_prec (x, 65); 189 mpfr_set_prec (y, 65); 190 mpfr_set_prec (z, 64); 191 mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101"); 192 mpfr_set_str_binary (y, "0.1110111011110001110111011111111111101000011001011100101100101100"); 193 test_sub (z, x, y, MPFR_RNDZ); 194 if (mpfr_cmp_ui (z, 1)) 195 { 196 printf ("Error in mpfr_sub (1)\n"); 197 exit (1); 198 } 199 test_sub (z, x, y, MPFR_RNDU); 200 mpfr_set_str_binary (x, "1.000000000000000000000000000000000000000000000000000000000000001"); 201 if (mpfr_cmp (z, x)) 202 { 203 printf ("Error in mpfr_sub (2)\n"); 204 printf ("Expected "); mpfr_print_binary (x); puts (""); 205 printf ("Got "); mpfr_print_binary (z); puts (""); 206 exit (1); 207 } 208 mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101"); 209 test_sub (z, x, y, MPFR_RNDN); 210 if (mpfr_cmp_ui (z, 1)) 211 { 212 printf ("Error in mpfr_sub (3)\n"); 213 exit (1); 214 } 215 inexact = test_sub (z, x, y, MPFR_RNDA); 216 mpfr_set_str_binary (x, "1.000000000000000000000000000000000000000000000000000000000000001"); 217 if (mpfr_cmp (z, x) || inexact <= 0) 218 { 219 printf ("Error in mpfr_sub (4)\n"); 220 exit (1); 221 } 222 mpfr_set_prec (x, 66); 223 mpfr_set_str_binary (x, "1.11101110111100011101110111111111111010000110010111001011001010111"); 224 test_sub (z, x, y, MPFR_RNDN); 225 if (mpfr_cmp_ui (z, 1)) 226 { 227 printf ("Error in mpfr_sub (5)\n"); 228 exit (1); 229 } 230 231 /* check in-place operations */ 232 mpfr_set_ui (x, 1, MPFR_RNDN); 233 test_sub (x, x, x, MPFR_RNDN); 234 if (mpfr_cmp_ui(x, 0)) 235 { 236 printf ("Error for mpfr_sub (x, x, x, MPFR_RNDN) with x=1.0\n"); 237 exit (1); 238 } 239 240 mpfr_set_prec (x, 53); 241 mpfr_set_prec (y, 53); 242 mpfr_set_prec (z, 53); 243 mpfr_set_str1 (x, "1.229318102e+09"); 244 mpfr_set_str1 (y, "2.32221184180698677665e+05"); 245 test_sub (z, x, y, MPFR_RNDN); 246 if (mpfr_cmp_str1 (z, "1229085880.815819263458251953125")) 247 { 248 printf ("Error in mpfr_sub (1.22e9 - 2.32e5)\n"); 249 printf ("expected 1229085880.815819263458251953125, got "); 250 mpfr_out_str(stdout, 10, 0, z, MPFR_RNDN); 251 putchar('\n'); 252 exit (1); 253 } 254 255 mpfr_set_prec (x, 112); 256 mpfr_set_prec (y, 98); 257 mpfr_set_prec (z, 54); 258 mpfr_set_str_binary (x, "0.11111100100000000011000011100000101101010001000111E-401"); 259 mpfr_set_str_binary (y, "0.10110000100100000101101100011111111011101000111000101E-464"); 260 test_sub (z, x, y, MPFR_RNDN); 261 if (mpfr_cmp (z, x)) { 262 printf ("mpfr_sub(z, x, y) failed for prec(x)=112, prec(y)=98\n"); 263 printf ("expected "); mpfr_print_binary (x); puts (""); 264 printf ("got "); mpfr_print_binary (z); puts (""); 265 exit (1); 266 } 267 268 mpfr_set_prec (x, 33); 269 mpfr_set_ui (x, 1, MPFR_RNDN); 270 mpfr_div_2exp (x, x, 32, MPFR_RNDN); 271 mpfr_sub_ui (x, x, 1, MPFR_RNDN); 272 273 mpfr_set_prec (x, 5); 274 mpfr_set_prec (y, 5); 275 mpfr_set_str_binary (x, "1e-12"); 276 mpfr_set_ui (y, 1, MPFR_RNDN); 277 test_sub (x, y, x, MPFR_RNDD); 278 mpfr_set_str_binary (y, "0.11111"); 279 if (mpfr_cmp (x, y)) 280 { 281 printf ("Error in mpfr_sub (x, y, x, MPFR_RNDD) for x=2^(-12), y=1\n"); 282 exit (1); 283 } 284 285 mpfr_set_prec (x, 24); 286 mpfr_set_prec (y, 24); 287 mpfr_set_str_binary (x, "-0.100010000000000000000000E19"); 288 mpfr_set_str_binary (y, "0.100000000000000000000100E15"); 289 mpfr_add (x, x, y, MPFR_RNDD); 290 mpfr_set_str_binary (y, "-0.1E19"); 291 if (mpfr_cmp (x, y)) 292 { 293 printf ("Error in mpfr_add (2)\n"); 294 exit (1); 295 } 296 297 mpfr_set_prec (x, 2); 298 mpfr_set_prec (y, 10); 299 mpfr_set_prec (z, 10); 300 mpfr_set_ui (y, 0, MPFR_RNDN); 301 mpfr_set_str_binary (z, "0.10001"); 302 if (test_sub (x, y, z, MPFR_RNDN) <= 0) 303 { 304 printf ("Wrong inexact flag in x=mpfr_sub(0,z) for prec(z)>prec(x)\n"); 305 exit (1); 306 } 307 if (test_sub (x, z, y, MPFR_RNDN) >= 0) 308 { 309 printf ("Wrong inexact flag in x=mpfr_sub(z,0) for prec(z)>prec(x)\n"); 310 exit (1); 311 } 312 313 mpfr_clear (x); 314 mpfr_clear (y); 315 mpfr_clear (z); 316 } 317 318 static void 319 bug_ddefour(void) 320 { 321 mpfr_t ex, ex1, ex2, ex3, tot, tot1; 322 323 mpfr_init2(ex, 53); 324 mpfr_init2(ex1, 53); 325 mpfr_init2(ex2, 53); 326 mpfr_init2(ex3, 53); 327 mpfr_init2(tot, 150); 328 mpfr_init2(tot1, 150); 329 330 mpfr_set_ui( ex, 1, MPFR_RNDN); 331 mpfr_mul_2exp( ex, ex, 906, MPFR_RNDN); 332 mpfr_log( tot, ex, MPFR_RNDN); 333 mpfr_set( ex1, tot, MPFR_RNDN); /* ex1 = high(tot) */ 334 test_sub( ex2, tot, ex1, MPFR_RNDN); /* ex2 = high(tot - ex1) */ 335 test_sub( tot1, tot, ex1, MPFR_RNDN); /* tot1 = tot - ex1 */ 336 mpfr_set( ex3, tot1, MPFR_RNDN); /* ex3 = high(tot - ex1) */ 337 338 if (mpfr_cmp(ex2, ex3)) 339 { 340 printf ("Error in ddefour test.\n"); 341 printf ("ex2="); mpfr_print_binary (ex2); puts (""); 342 printf ("ex3="); mpfr_print_binary (ex3); puts (""); 343 exit (1); 344 } 345 346 mpfr_clear (ex); 347 mpfr_clear (ex1); 348 mpfr_clear (ex2); 349 mpfr_clear (ex3); 350 mpfr_clear (tot); 351 mpfr_clear (tot1); 352 } 353 354 /* if u = o(x-y), v = o(u-x), w = o(v+y), then x-y = u-w */ 355 static void 356 check_two_sum (mpfr_prec_t p) 357 { 358 mpfr_t x, y, u, v, w; 359 mpfr_rnd_t rnd; 360 int inexact; 361 362 mpfr_init2 (x, p); 363 mpfr_init2 (y, p); 364 mpfr_init2 (u, p); 365 mpfr_init2 (v, p); 366 mpfr_init2 (w, p); 367 mpfr_urandomb (x, RANDS); 368 mpfr_urandomb (y, RANDS); 369 if (mpfr_cmpabs (x, y) < 0) 370 mpfr_swap (x, y); 371 rnd = MPFR_RNDN; 372 inexact = test_sub (u, x, y, rnd); 373 test_sub (v, u, x, rnd); 374 mpfr_add (w, v, y, rnd); 375 /* as u = (x-y) - w, we should have inexact and w of opposite signs */ 376 if (((inexact == 0) && mpfr_cmp_ui (w, 0)) || 377 ((inexact > 0) && (mpfr_cmp_ui (w, 0) <= 0)) || 378 ((inexact < 0) && (mpfr_cmp_ui (w, 0) >= 0))) 379 { 380 printf ("Wrong inexact flag for prec=%u, rnd=%s\n", (unsigned)p, 381 mpfr_print_rnd_mode (rnd)); 382 printf ("x="); mpfr_print_binary(x); puts (""); 383 printf ("y="); mpfr_print_binary(y); puts (""); 384 printf ("u="); mpfr_print_binary(u); puts (""); 385 printf ("v="); mpfr_print_binary(v); puts (""); 386 printf ("w="); mpfr_print_binary(w); puts (""); 387 printf ("inexact = %d\n", inexact); 388 exit (1); 389 } 390 mpfr_clear (x); 391 mpfr_clear (y); 392 mpfr_clear (u); 393 mpfr_clear (v); 394 mpfr_clear (w); 395 } 396 397 #define MAX_PREC 200 398 399 static void 400 check_inexact (void) 401 { 402 mpfr_t x, y, z, u; 403 mpfr_prec_t px, py, pu, pz; 404 int inexact, cmp; 405 mpfr_rnd_t rnd; 406 407 mpfr_init (x); 408 mpfr_init (y); 409 mpfr_init (z); 410 mpfr_init (u); 411 412 mpfr_set_prec (x, 2); 413 mpfr_set_ui (x, 6, MPFR_RNDN); 414 mpfr_div_2exp (x, x, 4, MPFR_RNDN); /* x = 6/16 */ 415 mpfr_set_prec (y, 2); 416 mpfr_set_si (y, -1, MPFR_RNDN); 417 mpfr_div_2exp (y, y, 4, MPFR_RNDN); /* y = -1/16 */ 418 inexact = test_sub (y, y, x, MPFR_RNDN); /* y = round(-7/16) = -1/2 */ 419 if (inexact >= 0) 420 { 421 printf ("Error: wrong inexact flag for -1/16 - (6/16)\n"); 422 exit (1); 423 } 424 425 for (px=2; px<MAX_PREC; px++) 426 { 427 mpfr_set_prec (x, px); 428 do 429 { 430 mpfr_urandomb (x, RANDS); 431 } 432 while (mpfr_cmp_ui (x, 0) == 0); 433 for (pu=2; pu<MAX_PREC; pu++) 434 { 435 mpfr_set_prec (u, pu); 436 do 437 { 438 mpfr_urandomb (u, RANDS); 439 } 440 while (mpfr_cmp_ui (u, 0) == 0); 441 { 442 py = 2 + (randlimb () % (MAX_PREC - 2)); 443 mpfr_set_prec (y, py); 444 /* warning: MPFR_EXP is undefined for 0 */ 445 pz = (mpfr_cmpabs (x, u) >= 0) ? MPFR_EXP(x) - MPFR_EXP(u) 446 : MPFR_EXP(u) - MPFR_EXP(x); 447 pz = pz + MAX(MPFR_PREC(x), MPFR_PREC(u)); 448 mpfr_set_prec (z, pz); 449 rnd = RND_RAND (); 450 if (test_sub (z, x, u, rnd)) 451 { 452 printf ("z <- x - u should be exact\n"); 453 exit (1); 454 } 455 { 456 rnd = RND_RAND (); 457 inexact = test_sub (y, x, u, rnd); 458 cmp = mpfr_cmp (y, z); 459 if (((inexact == 0) && (cmp != 0)) || 460 ((inexact > 0) && (cmp <= 0)) || 461 ((inexact < 0) && (cmp >= 0))) 462 { 463 printf ("Wrong inexact flag for rnd=%s\n", 464 mpfr_print_rnd_mode(rnd)); 465 printf ("expected %d, got %d\n", cmp, inexact); 466 printf ("x="); mpfr_print_binary (x); puts (""); 467 printf ("u="); mpfr_print_binary (u); puts (""); 468 printf ("y= "); mpfr_print_binary (y); puts (""); 469 printf ("x-u="); mpfr_print_binary (z); puts (""); 470 exit (1); 471 } 472 } 473 } 474 } 475 } 476 477 mpfr_clear (x); 478 mpfr_clear (y); 479 mpfr_clear (z); 480 mpfr_clear (u); 481 } 482 483 /* Bug found by Jakub Jelinek 484 * http://bugzilla.redhat.com/643657 485 * https://gforge.inria.fr/tracker/index.php?func=detail&aid=11301 486 * The consequence can be either an assertion failure (i = 2 in the 487 * testcase below, in debug mode) or an incorrectly rounded value. 488 */ 489 static void 490 bug20101017 (void) 491 { 492 mpfr_t a, b, c; 493 int inex; 494 int i; 495 496 mpfr_init2 (a, GMP_NUMB_BITS * 2); 497 mpfr_init2 (b, GMP_NUMB_BITS); 498 mpfr_init2 (c, GMP_NUMB_BITS); 499 500 /* a = 2^(2N) + k.2^(2N-1) + 2^N and b = 1 501 with N = GMP_NUMB_BITS and k = 0 or 1. 502 c = a - b should round to the same value as a. */ 503 504 for (i = 2; i <= 3; i++) 505 { 506 mpfr_set_ui_2exp (a, i, GMP_NUMB_BITS - 1, MPFR_RNDN); 507 mpfr_add_ui (a, a, 1, MPFR_RNDN); 508 mpfr_mul_2ui (a, a, GMP_NUMB_BITS, MPFR_RNDN); 509 mpfr_set_ui (b, 1, MPFR_RNDN); 510 inex = mpfr_sub (c, a, b, MPFR_RNDN); 511 mpfr_set (b, a, MPFR_RNDN); 512 if (! mpfr_equal_p (c, b)) 513 { 514 printf ("Error in bug20101017 for i = %d.\n", i); 515 printf ("Expected "); 516 mpfr_out_str (stdout, 16, 0, b, MPFR_RNDN); 517 putchar ('\n'); 518 printf ("Got "); 519 mpfr_out_str (stdout, 16, 0, c, MPFR_RNDN); 520 putchar ('\n'); 521 exit (1); 522 } 523 if (inex >= 0) 524 { 525 printf ("Error in bug20101017 for i = %d: bad inex value.\n", i); 526 printf ("Expected negative, got %d.\n", inex); 527 exit (1); 528 } 529 } 530 531 mpfr_set_prec (a, 64); 532 mpfr_set_prec (b, 129); 533 mpfr_set_prec (c, 2); 534 mpfr_set_str_binary (b, "0.100000000000000000000000000000000000000000000000000000000000000010000000000000000000000000000000000000000000000000000000000000001E65"); 535 mpfr_set_str_binary (c, "0.10E1"); 536 inex = mpfr_sub (a, b, c, MPFR_RNDN); 537 if (mpfr_cmp_ui_2exp (a, 1, 64) != 0 || inex >= 0) 538 { 539 printf ("Error in mpfr_sub for b-c for b=2^64+1+2^(-64), c=1\n"); 540 printf ("Expected result 2^64 with inex < 0\n"); 541 printf ("Got "); mpfr_print_binary (a); 542 printf (" with inex=%d\n", inex); 543 exit (1); 544 } 545 546 mpfr_clears (a, b, c, (mpfr_ptr) 0); 547 } 548 549 /* hard test of rounding */ 550 static void 551 check_rounding (void) 552 { 553 mpfr_t a, b, c, res; 554 mpfr_prec_t p; 555 long k, l; 556 int i; 557 558 #define MAXKL (2 * GMP_NUMB_BITS) 559 for (p = MPFR_PREC_MIN; p <= GMP_NUMB_BITS; p++) 560 { 561 mpfr_init2 (a, p); 562 mpfr_init2 (res, p); 563 mpfr_init2 (b, p + 1 + MAXKL); 564 mpfr_init2 (c, MPFR_PREC_MIN); 565 566 /* b = 2^p + 1 + 2^(-k), c = 2^(-l) */ 567 for (k = 0; k <= MAXKL; k++) 568 for (l = 0; l <= MAXKL; l++) 569 { 570 mpfr_set_ui_2exp (b, 1, p, MPFR_RNDN); 571 mpfr_add_ui (b, b, 1, MPFR_RNDN); 572 mpfr_mul_2ui (b, b, k, MPFR_RNDN); 573 mpfr_add_ui (b, b, 1, MPFR_RNDN); 574 mpfr_div_2ui (b, b, k, MPFR_RNDN); 575 mpfr_set_ui_2exp (c, 1, -l, MPFR_RNDN); 576 i = mpfr_sub (a, b, c, MPFR_RNDN); 577 /* b - c = 2^p + 1 + 2^(-k) - 2^(-l), should be rounded to 578 2^p for l <= k, and 2^p+2 for l < k */ 579 if (l <= k) 580 { 581 if (mpfr_cmp_ui_2exp (a, 1, p) != 0) 582 { 583 printf ("Wrong result in check_rounding\n"); 584 printf ("p=%lu k=%ld l=%ld\n", (unsigned long) p, k, l); 585 printf ("b="); mpfr_print_binary (b); puts (""); 586 printf ("c="); mpfr_print_binary (c); puts (""); 587 printf ("Expected 2^%lu\n", (unsigned long) p); 588 printf ("Got "); mpfr_print_binary (a); puts (""); 589 exit (1); 590 } 591 if (i >= 0) 592 { 593 printf ("Wrong ternary value in check_rounding\n"); 594 printf ("p=%lu k=%ld l=%ld\n", (unsigned long) p, k, l); 595 printf ("b="); mpfr_print_binary (b); puts (""); 596 printf ("c="); mpfr_print_binary (c); puts (""); 597 printf ("a="); mpfr_print_binary (a); puts (""); 598 printf ("Expected < 0, got %d\n", i); 599 exit (1); 600 } 601 } 602 else /* l < k */ 603 { 604 mpfr_set_ui_2exp (res, 1, p, MPFR_RNDN); 605 mpfr_add_ui (res, res, 2, MPFR_RNDN); 606 if (mpfr_cmp (a, res) != 0) 607 { 608 printf ("Wrong result in check_rounding\n"); 609 printf ("b="); mpfr_print_binary (b); puts (""); 610 printf ("c="); mpfr_print_binary (c); puts (""); 611 printf ("Expected "); mpfr_print_binary (res); puts (""); 612 printf ("Got "); mpfr_print_binary (a); puts (""); 613 exit (1); 614 } 615 if (i <= 0) 616 { 617 printf ("Wrong ternary value in check_rounding\n"); 618 printf ("b="); mpfr_print_binary (b); puts (""); 619 printf ("c="); mpfr_print_binary (c); puts (""); 620 printf ("Expected > 0, got %d\n", i); 621 exit (1); 622 } 623 } 624 } 625 626 mpfr_clear (a); 627 mpfr_clear (res); 628 mpfr_clear (b); 629 mpfr_clear (c); 630 } 631 } 632 633 #define TEST_FUNCTION test_sub 634 #define TWO_ARGS 635 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS) 636 #include "tgeneric.c" 637 638 int 639 main (void) 640 { 641 mpfr_prec_t p; 642 unsigned int i; 643 644 tests_start_mpfr (); 645 646 bug20101017 (); 647 check_rounding (); 648 check_diverse (); 649 check_inexact (); 650 bug_ddefour (); 651 for (p=2; p<200; p++) 652 for (i=0; i<50; i++) 653 check_two_sum (p); 654 test_generic (2, 800, 100); 655 656 tests_end_mpfr (); 657 return 0; 658 } 659