1/* Test file for mpfr_sub. 2 3Copyright 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_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 53static void 54check_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 318static void 319bug_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 */ 355static void 356check_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 399static void 400check_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 */ 489static void 490bug20101017 (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 */ 550static void 551check_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", 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", 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", 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 638int 639main (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