1 /* tsum -- test file for the list summation function 2 3 Copyright 2004-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 #ifndef MPFR_NCANCEL 26 #define MPFR_NCANCEL 10 27 #endif 28 29 #include <time.h> 30 31 /* return the cpu time in milliseconds */ 32 static int 33 cputime (void) 34 { 35 return clock () / (CLOCKS_PER_SEC / 1000); 36 } 37 38 static mpfr_prec_t 39 get_prec_max (mpfr_t *t, int n) 40 { 41 mpfr_exp_t e, min, max; 42 int i; 43 44 min = MPFR_EMAX_MAX; 45 max = MPFR_EMIN_MIN; 46 for (i = 0; i < n; i++) 47 if (MPFR_IS_PURE_FP (t[i])) 48 { 49 e = MPFR_GET_EXP (t[i]); 50 if (e > max) 51 max = e; 52 e -= MPFR_GET_PREC (t[i]); 53 if (e < min) 54 min = e; 55 } 56 57 return min > max ? MPFR_PREC_MIN /* no pure FP values */ 58 : max - min + __gmpfr_ceil_log2 (n); 59 } 60 61 static void 62 get_exact_sum (mpfr_t sum, mpfr_t *tab, int n) 63 { 64 int i; 65 66 mpfr_set_prec (sum, get_prec_max (tab, n)); 67 mpfr_set_ui (sum, 0, MPFR_RNDN); 68 for (i = 0; i < n; i++) 69 if (mpfr_add (sum, sum, tab[i], MPFR_RNDN)) 70 { 71 printf ("FIXME: get_exact_sum is buggy.\n"); 72 exit (1); 73 } 74 } 75 76 static void 77 generic_tests (void) 78 { 79 mpfr_t exact_sum, sum1, sum2; 80 mpfr_t *t; 81 mpfr_ptr *p; 82 mpfr_prec_t precmax = 444; 83 int i, m, nmax = 500; 84 int rnd_mode; 85 86 t = (mpfr_t *) tests_allocate (nmax * sizeof(mpfr_t)); 87 p = (mpfr_ptr *) tests_allocate (nmax * sizeof(mpfr_ptr)); 88 for (i = 0; i < nmax; i++) 89 { 90 mpfr_init2 (t[i], precmax); 91 p[i] = t[i]; 92 } 93 mpfr_inits2 (precmax, exact_sum, sum1, sum2, (mpfr_ptr) 0); 94 95 for (m = 0; m < 20000; m++) 96 { 97 int non_uniform, n; 98 mpfr_prec_t prec; 99 100 non_uniform = randlimb () % 10; 101 n = (randlimb () % nmax) + 1; 102 prec = MPFR_PREC_MIN + (randlimb () % (precmax - MPFR_PREC_MIN + 1)); 103 mpfr_set_prec (sum1, prec); 104 mpfr_set_prec (sum2, prec); 105 106 for (i = 0; i < n; i++) 107 { 108 mpfr_set_prec (t[i], MPFR_PREC_MIN + 109 (randlimb () % (precmax - MPFR_PREC_MIN + 1))); 110 mpfr_urandomb (t[i], RANDS); 111 if (m % 8 != 0 && (m % 8 == 1 || (randlimb () & 1))) 112 mpfr_neg (t[i], t[i], MPFR_RNDN); 113 if (non_uniform && MPFR_NOTZERO (t[i])) 114 mpfr_set_exp (t[i], randlimb () % 1000); 115 /* putchar ("-0+"[SIGN (mpfr_sgn (t[i])) + 1]); */ 116 } 117 /* putchar ('\n'); */ 118 get_exact_sum (exact_sum, t, n); 119 RND_LOOP (rnd_mode) 120 if (rnd_mode == MPFR_RNDF) 121 { 122 int inex; 123 124 inex = mpfr_set (sum1, exact_sum, MPFR_RNDD); 125 mpfr_sum (sum2, p, n, MPFR_RNDF); 126 if (! mpfr_equal_p (sum1, sum2) && 127 (inex == 0 || 128 (mpfr_nextabove (sum1), ! mpfr_equal_p (sum1, sum2)))) 129 { 130 printf ("generic_tests failed on m = %d, MPFR_RNDF\n", m); 131 printf ("Exact sum = "); 132 mpfr_dump (exact_sum); 133 printf ("Got "); 134 mpfr_dump (sum2); 135 exit (1); 136 } 137 } 138 else 139 { 140 int inex1, inex2; 141 142 inex1 = mpfr_set (sum1, exact_sum, (mpfr_rnd_t) rnd_mode); 143 inex2 = mpfr_sum (sum2, p, n, (mpfr_rnd_t) rnd_mode); 144 if (!(mpfr_equal_p (sum1, sum2) && SAME_SIGN (inex1, inex2))) 145 { 146 printf ("generic_tests failed on m = %d, %s\n", m, 147 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd_mode)); 148 printf ("Expected "); 149 mpfr_dump (sum1); 150 printf ("with inex = %d\n", inex1); 151 printf ("Got "); 152 mpfr_dump (sum2); 153 printf ("with inex = %d\n", inex2); 154 exit (1); 155 } 156 } 157 } 158 159 for (i = 0; i < nmax; i++) 160 mpfr_clear (t[i]); 161 mpfr_clears (exact_sum, sum1, sum2, (mpfr_ptr) 0); 162 tests_free (t, nmax * sizeof(mpfr_t)); 163 tests_free (p, nmax * sizeof(mpfr_ptr)); 164 } 165 166 /* glibc free() error or segmentation fault when configured 167 * with GMP 6.0.0 built with "--disable-alloca ABI=32". 168 * GCC's address sanitizer shows a heap-buffer-overflow. 169 * Fixed in r9369 (before the merge into the trunk). The problem was due 170 * to the fact that mpn functions do not accept a zero size argument, and 171 * since mpn_add_1 is here a macro in GMP, there's no assertion even when 172 * GMP was built with assertion checking (--enable-assert). 173 */ 174 static 175 void check_simple (void) 176 { 177 mpfr_t tab[3], r; 178 mpfr_ptr tabp[3]; 179 int i; 180 181 mpfr_init2 (r, 16); 182 for (i = 0; i < 3; i++) 183 { 184 mpfr_init2 (tab[i], 16); 185 mpfr_set_ui (tab[i], 1, MPFR_RNDN); 186 tabp[i] = tab[i]; 187 } 188 189 i = mpfr_sum (r, tabp, 3, MPFR_RNDN); 190 if (mpfr_cmp_ui (r, 3) || i != 0) 191 { 192 printf ("Error in check_simple\n"); 193 exit (1); 194 } 195 196 mpfr_clears (tab[0], tab[1], tab[2], r, (mpfr_ptr) 0); 197 } 198 199 static 200 void check_special (void) 201 { 202 mpfr_t tab[3], r; 203 mpfr_ptr tabp[3]; 204 int i; 205 206 mpfr_inits2 (53, tab[0], tab[1], tab[2], r, (mpfr_ptr) 0); 207 tabp[0] = tab[0]; 208 tabp[1] = tab[1]; 209 tabp[2] = tab[2]; 210 211 i = mpfr_sum (r, tabp, 0, MPFR_RNDN); 212 if (!MPFR_IS_ZERO (r) || !MPFR_IS_POS (r) || i != 0) 213 { 214 printf ("Special case n==0 failed!\n"); 215 exit (1); 216 } 217 218 mpfr_set_ui (tab[0], 42, MPFR_RNDN); 219 i = mpfr_sum (r, tabp, 1, MPFR_RNDN); 220 if (mpfr_cmp_ui (r, 42) || i != 0) 221 { 222 printf ("Special case n==1 failed!\n"); 223 exit (1); 224 } 225 226 mpfr_set_ui (tab[1], 17, MPFR_RNDN); 227 MPFR_SET_NAN (tab[2]); 228 i = mpfr_sum (r, tabp, 3, MPFR_RNDN); 229 if (!MPFR_IS_NAN (r) || i != 0) 230 { 231 printf ("Special case NAN failed!\n"); 232 exit (1); 233 } 234 235 MPFR_SET_INF (tab[2]); 236 MPFR_SET_POS (tab[2]); 237 i = mpfr_sum (r, tabp, 3, MPFR_RNDN); 238 if (!MPFR_IS_INF (r) || !MPFR_IS_POS (r) || i != 0) 239 { 240 printf ("Special case +INF failed!\n"); 241 exit (1); 242 } 243 244 MPFR_SET_INF (tab[2]); 245 MPFR_SET_NEG (tab[2]); 246 i = mpfr_sum (r, tabp, 3, MPFR_RNDN); 247 if (!MPFR_IS_INF (r) || !MPFR_IS_NEG (r) || i != 0) 248 { 249 printf ("Special case -INF failed!\n"); 250 exit (1); 251 } 252 253 MPFR_SET_ZERO (tab[1]); 254 i = mpfr_sum (r, tabp, 2, MPFR_RNDN); 255 if (mpfr_cmp_ui (r, 42) || i != 0) 256 { 257 printf ("Special case 42+0 failed!\n"); 258 exit (1); 259 } 260 261 MPFR_SET_NAN (tab[0]); 262 i = mpfr_sum (r, tabp, 3, MPFR_RNDN); 263 if (!MPFR_IS_NAN (r) || i != 0) 264 { 265 printf ("Special case NAN+0+-INF failed!\n"); 266 exit (1); 267 } 268 269 mpfr_set_inf (tab[0], 1); 270 mpfr_set_ui (tab[1], 59, MPFR_RNDN); 271 mpfr_set_inf (tab[2], -1); 272 i = mpfr_sum (r, tabp, 3, MPFR_RNDN); 273 if (!MPFR_IS_NAN (r) || i != 0) 274 { 275 printf ("Special case +INF + 59 +-INF failed!\n"); 276 exit (1); 277 } 278 279 mpfr_clears (tab[0], tab[1], tab[2], r, (mpfr_ptr) 0); 280 } 281 282 #define NC 7 283 #define NS 6 284 285 static void 286 check_more_special (void) 287 { 288 char *str[NC] = { "NaN", "+Inf", "-Inf", "+0", "-0", "+1", "-1" }; 289 int i, r, k[NS]; 290 mpfr_t c[NC], s[NS], sum; 291 mpfr_ptr p[NS]; 292 int inex; 293 294 for (i = 0; i < NC; i++) 295 { 296 int ret; 297 mpfr_init2 (c[i], 8); 298 ret = mpfr_set_str (c[i], str[i], 0, MPFR_RNDN); 299 MPFR_ASSERTN (ret == 0); 300 } 301 for (i = 0; i < NS; i++) 302 mpfr_init2 (s[i], 8); 303 mpfr_init2 (sum, 8); 304 305 RND_LOOP(r) 306 { 307 i = 0; 308 while (1) 309 { 310 while (i < NS) 311 { 312 p[i] = c[0]; 313 mpfr_set_nan (s[i]); 314 k[i++] = 0; 315 } 316 inex = mpfr_sum (sum, p, NS, (mpfr_rnd_t) r); 317 if (! SAME_VAL (sum, s[NS-1]) || inex != 0) 318 { 319 printf ("Error in check_more_special on %s", 320 mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 321 for (i = 0; i < NS; i++) 322 printf (" %d", k[i]); 323 printf (" with\n"); 324 for (i = 0; i < NS; i++) 325 { 326 printf (" p[%d] = %s = ", i, str[k[i]]); 327 mpfr_dump (p[i]); 328 } 329 printf ("Expected "); 330 mpfr_dump (s[NS-1]); 331 printf ("with inex = 0\n"); 332 printf ("Got "); 333 mpfr_dump (sum); 334 printf ("with inex = %d\n", inex); 335 exit (1); 336 } 337 while (k[--i] == NC-1) 338 if (i == 0) 339 goto next_rnd; 340 p[i] = c[++k[i]]; 341 if (i == 0) 342 mpfr_set (s[i], p[i], (mpfr_rnd_t) r); 343 else 344 mpfr_add (s[i], s[i-1], p[i], (mpfr_rnd_t) r); 345 i++; 346 } 347 next_rnd: ; 348 } 349 350 for (i = 0; i < NC; i++) 351 mpfr_clear (c[i]); 352 for (i = 0; i < NS; i++) 353 mpfr_clear (s[i]); 354 mpfr_clear (sum); 355 } 356 357 /* i * 2^(46+h) + j * 2^(45+h) + k * 2^(44+h) + f * 2^(-2), 358 with -1 <= i, j, k <= 1, i != 0, -3 <= f <= 3, and 359 * prec set up so that ulp(exact sum) = 2^0, then 360 * prec set up so that ulp(exact sum) = 2^(44+h) when possible, 361 i.e. when prec >= MPFR_PREC_MIN. 362 ------ 363 Some explanations: 364 ulp(exact sum) = 2^q means EXP(exact sum) - prec = q where prec is 365 the precision of the output. Thus ulp(exact sum) = 2^0 is achieved 366 by setting prec = EXP(s3), where s3 is the exact sum (computed with 367 mpfr_add's and sufficient precision). Then ulp(exact sum) = 2^(44+h) 368 is achieved by subtracting 44+h from prec. The loop on prec does 369 this. Since EXP(s3) <= 47+h, prec <= 3 at the second iteration, 370 thus there will be at most 2 iterations. Whether a second iteration 371 is done or not depends on EXP(s3), i.e. the values of the parameters, 372 and the value of MPFR_PREC_MIN. */ 373 static void 374 check1 (int h) 375 { 376 mpfr_t sum1, sum2, s1, s2, s3, t[4]; 377 mpfr_ptr p[4]; 378 int i, j, k, f, prec, r, inex1, inex2; 379 380 mpfr_init2 (sum1, 47 + h); 381 mpfr_init2 (sum2, 47 + h); 382 mpfr_init2 (s1, 3); 383 mpfr_init2 (s2, 3); 384 mpfr_init2 (s3, 49 + h); 385 for (i = 0; i < 4; i++) 386 { 387 mpfr_init2 (t[i], 2); 388 p[i] = t[i]; 389 } 390 391 for (i = -1; i <= 1; i += 2) 392 { 393 mpfr_set_si_2exp (t[0], i, 46 + h, MPFR_RNDN); 394 for (j = -1; j <= 1; j++) 395 { 396 mpfr_set_si_2exp (t[1], j, 45 + h, MPFR_RNDN); 397 inex1 = mpfr_add (s1, t[0], t[1], MPFR_RNDN); 398 MPFR_ASSERTN (inex1 == 0); 399 for (k = -1; k <= 1; k++) 400 { 401 mpfr_set_si_2exp (t[2], k, 44 + h, MPFR_RNDN); 402 inex1 = mpfr_add (s2, s1, t[2], MPFR_RNDN); 403 MPFR_ASSERTN (inex1 == 0); 404 for (f = -3; f <= 3; f++) 405 { 406 mpfr_set_si_2exp (t[3], f, -2, MPFR_RNDN); 407 inex1 = mpfr_add (s3, s2, t[3], MPFR_RNDN); 408 MPFR_ASSERTN (inex1 == 0); 409 for (prec = mpfr_get_exp (s3); 410 prec >= MPFR_PREC_MIN; 411 prec -= 44 + h) 412 { 413 mpfr_set_prec (sum1, prec); 414 mpfr_set_prec (sum2, prec); 415 RND_LOOP_NO_RNDF (r) 416 { 417 inex1 = mpfr_set (sum1, s3, (mpfr_rnd_t) r); 418 inex2 = mpfr_sum (sum2, p, 4, (mpfr_rnd_t) r); 419 MPFR_ASSERTN (mpfr_check (sum1)); 420 MPFR_ASSERTN (mpfr_check (sum2)); 421 if (!(mpfr_equal_p (sum1, sum2) && 422 SAME_SIGN (inex1, inex2))) 423 { 424 printf ("Error in check1 on %s, prec = %d, " 425 "i = %d, j = %d, k = %d, f = %d, " 426 "h = %d\n", 427 mpfr_print_rnd_mode ((mpfr_rnd_t) r), 428 prec, i, j, k, f, h); 429 printf ("Expected "); 430 mpfr_dump (sum1); 431 printf ("with inex = %d\n", inex1); 432 printf ("Got "); 433 mpfr_dump (sum2); 434 printf ("with inex = %d\n", inex2); 435 exit (1); 436 } 437 } 438 } 439 } 440 } 441 } 442 } 443 444 for (i = 0; i < 4; i++) 445 mpfr_clear (t[i]); 446 mpfr_clears (sum1, sum2, s1, s2, s3, (mpfr_ptr) 0); 447 } 448 449 /* With N = 2 * GMP_NUMB_BITS: 450 i * 2^N + j + k * 2^(-1) + f1 * 2^(-N) + f2 * 2^(-N), 451 with i = -1 or 1, j = 0 or i, -1 <= k <= 1, -1 <= f1 <= 1, -1 <= f2 <= 1 452 ulp(exact sum) = 2^0. */ 453 static void 454 check2 (void) 455 { 456 mpfr_t sum1, sum2, s1, s2, s3, s4, t[5]; 457 mpfr_ptr p[5]; 458 int i, j, k, f1, f2, prec, r, inex1, inex2; 459 460 #define N (2 * GMP_NUMB_BITS) 461 462 mpfr_init2 (sum1, N+1); 463 mpfr_init2 (sum2, N+1); 464 mpfr_init2 (s1, N+1); 465 mpfr_init2 (s2, N+2); 466 mpfr_init2 (s3, 2*N+1); 467 mpfr_init2 (s4, 2*N+1); 468 for (i = 0; i < 5; i++) 469 { 470 mpfr_init2 (t[i], 2); 471 p[i] = t[i]; 472 } 473 474 for (i = -1; i <= 1; i += 2) 475 { 476 mpfr_set_si_2exp (t[0], i, N, MPFR_RNDN); 477 for (j = 0; j != 2*i; j += i) 478 { 479 mpfr_set_si (t[1], j, MPFR_RNDN); 480 inex1 = mpfr_add (s1, t[0], t[1], MPFR_RNDN); 481 MPFR_ASSERTN (inex1 == 0); 482 for (k = -1; k <= 1; k++) 483 { 484 mpfr_set_si_2exp (t[2], k, -1, MPFR_RNDN); 485 inex1 = mpfr_add (s2, s1, t[2], MPFR_RNDN); 486 MPFR_ASSERTN (inex1 == 0); 487 for (f1 = -1; f1 <= 1; f1++) 488 { 489 mpfr_set_si_2exp (t[3], f1, -N, MPFR_RNDN); 490 inex1 = mpfr_add (s3, s2, t[3], MPFR_RNDN); 491 MPFR_ASSERTN (inex1 == 0); 492 for (f2 = -1; f2 <= 1; f2++) 493 { 494 mpfr_set_si_2exp (t[4], f2, -N, MPFR_RNDN); 495 inex1 = mpfr_add (s4, s3, t[4], MPFR_RNDN); 496 MPFR_ASSERTN (inex1 == 0); 497 prec = mpfr_get_exp (s4); 498 mpfr_set_prec (sum1, prec); 499 mpfr_set_prec (sum2, prec); 500 RND_LOOP_NO_RNDF (r) 501 { 502 inex1 = mpfr_set (sum1, s4, (mpfr_rnd_t) r); 503 inex2 = mpfr_sum (sum2, p, 5, (mpfr_rnd_t) r); 504 MPFR_ASSERTN (mpfr_check (sum1)); 505 MPFR_ASSERTN (mpfr_check (sum2)); 506 if (!(mpfr_equal_p (sum1, sum2) && 507 SAME_SIGN (inex1, inex2))) 508 { 509 printf ("Error in check2 on %s, prec = %d, " 510 "i = %d, j = %d, k = %d, f1 = %d, " 511 "f2 = %d\n", 512 mpfr_print_rnd_mode ((mpfr_rnd_t) r), 513 prec, i, j, k, f1, f2); 514 printf ("Expected "); 515 mpfr_dump (sum1); 516 printf ("with inex = %d\n", inex1); 517 printf ("Got "); 518 mpfr_dump (sum2); 519 printf ("with inex = %d\n", inex2); 520 exit (1); 521 } 522 } 523 } 524 } 525 } 526 } 527 } 528 529 for (i = 0; i < 5; i++) 530 mpfr_clear (t[i]); 531 mpfr_clears (sum1, sum2, s1, s2, s3, s4, (mpfr_ptr) 0); 532 } 533 534 /* t[i] = (2^17 - 1) * 2^(17*(i-8)) for 0 <= i <= 16. 535 * t[17] = 2^(17*9+1) * j for -4 <= j <= 4. 536 * t[18] = 2^(-1) * k for -1 <= k <= 1. 537 * t[19] = 2^(-17*8) * m for -3 <= m <= 3. 538 * prec = MPFR_PREC_MIN and 17*9+4 539 */ 540 static void 541 check3 (void) 542 { 543 mpfr_t sum1, sum2, s1, s2, s3, s4, t[20]; 544 mpfr_ptr p[20]; 545 mpfr_flags_t flags1, flags2; 546 int i, s, j, k, m, q, r, inex1, inex2; 547 int prec[2] = { MPFR_PREC_MIN, 17*9+4 }; 548 549 mpfr_init2 (s1, 17*17); 550 mpfr_init2 (s2, 17*17+4); 551 mpfr_init2 (s3, 17*17+4); 552 mpfr_init2 (s4, 17*17+5); 553 mpfr_set_ui (s1, 0, MPFR_RNDN); 554 for (i = 0; i < 20; i++) 555 { 556 mpfr_init2 (t[i], 20); 557 p[i] = t[i]; 558 if (i < 17) 559 { 560 mpfr_set_ui_2exp (t[i], 0x1ffff, 17*(i-8), MPFR_RNDN); 561 inex1 = mpfr_add (s1, s1, t[i], MPFR_RNDN); 562 MPFR_ASSERTN (inex1 == 0); 563 } 564 } 565 566 for (s = 1; s >= -1; s -= 2) 567 { 568 for (j = -4; j <= 4; j++) 569 { 570 mpfr_set_si_2exp (t[17], j, 17*9+1, MPFR_RNDN); 571 inex1 = mpfr_add (s2, s1, t[17], MPFR_RNDN); 572 MPFR_ASSERTN (inex1 == 0); 573 for (k = -1; k <= 1; k++) 574 { 575 mpfr_set_si_2exp (t[18], k, -1, MPFR_RNDN); 576 inex1 = mpfr_add (s3, s2, t[18], MPFR_RNDN); 577 MPFR_ASSERTN (inex1 == 0); 578 for (m = -3; m <= 3; m++) 579 { 580 mpfr_set_si_2exp (t[19], m, -17*8, MPFR_RNDN); 581 inex1 = mpfr_add (s4, s3, t[19], MPFR_RNDN); 582 MPFR_ASSERTN (inex1 == 0); 583 for (q = 0; q < 2; q++) 584 { 585 mpfr_inits2 (prec[q], sum1, sum2, (mpfr_ptr) 0); 586 RND_LOOP_NO_RNDF (r) 587 { 588 mpfr_clear_flags (); 589 inex1 = mpfr_set (sum1, s4, (mpfr_rnd_t) r); 590 flags1 = __gmpfr_flags; 591 mpfr_clear_flags (); 592 inex2 = mpfr_sum (sum2, p, 20, (mpfr_rnd_t) r); 593 flags2 = __gmpfr_flags; 594 MPFR_ASSERTN (mpfr_check (sum1)); 595 MPFR_ASSERTN (mpfr_check (sum2)); 596 if (!(mpfr_equal_p (sum1, sum2) && 597 SAME_SIGN (inex1, inex2) && 598 flags1 == flags2)) 599 { 600 printf ("Error in check3 on %s, " 601 "s = %d, j = %d, k = %d, m = %d\n", 602 mpfr_print_rnd_mode ((mpfr_rnd_t) r), 603 s, j, k, m); 604 printf ("Expected "); 605 mpfr_dump (sum1); 606 printf ("with inex = %d and flags =", inex1); 607 flags_out (flags1); 608 printf ("Got "); 609 mpfr_dump (sum2); 610 printf ("with inex = %d and flags =", inex2); 611 flags_out (flags2); 612 exit (1); 613 } 614 } 615 mpfr_clears (sum1, sum2, (mpfr_ptr) 0); 616 } /* q */ 617 } /* m */ 618 } /* k */ 619 } /* j */ 620 for (i = 0; i < 17; i++) 621 mpfr_neg (t[i], t[i], MPFR_RNDN); 622 mpfr_neg (s1, s1, MPFR_RNDN); 623 } /* s */ 624 625 for (i = 0; i < 20; i++) 626 mpfr_clear (t[i]); 627 mpfr_clears (s1, s2, s3, s4, (mpfr_ptr) 0); 628 } 629 630 /* Test of s * (q * 2^(n-1) - 2^k) + h + i * 2^(-2) + j * 2^(-2) 631 * with h = -1 or 1, -1 <= i odd <= j <= 3, 2 <= q <= 3, s = -1 or 1, 632 * prec n-k. 633 * On a 64-bit machine: 634 * MPFR_RNDN, tmd=2, rbit=0, sst=0, negative is checked with the inputs 635 * -3*2^58, 2^5, -1, 2^(-2), 3*2^(-2) 636 * MPFR_RNDN, tmd=2, rbit=0, sst=1, negative is checked with the inputs 637 * -3*2^58, 2^5, -1, 3*2^(-2), 3*2^(-2) 638 * 639 * Note: This test detects an error in a result when "sq + 3" is replaced 640 * by "sq + 2" (11th argument of the first sum_raw invocation) and the 641 * corresponding assertion d >= 3 is removed, confirming that one cannot 642 * decrease this proved error bound. 643 */ 644 static void 645 check4 (void) 646 { 647 mpfr_t sum1, sum2, s1, s2, s3, s4, t[5]; 648 mpfr_ptr p[5]; 649 int h, i, j, k, n, q, r, s, prec, inex1, inex2; 650 651 mpfr_inits2 (257, sum1, sum2, s1, s2, s3, s4, (mpfr_ptr) 0); 652 for (i = 0; i < 5; i++) 653 { 654 mpfr_init2 (t[i], 2); 655 p[i] = t[i]; 656 } 657 658 /* No GNU style for the many nested loops... */ 659 for (k = 1; k <= 64; k++) { 660 mpfr_set_si_2exp (t[0], -1, k, MPFR_RNDN); 661 for (n = k + MPFR_PREC_MIN; n <= k + 65; n++) { 662 prec = n - k; 663 mpfr_set_prec (sum1, prec); 664 mpfr_set_prec (sum2, prec); 665 for (q = 2; q <= 3; q++) { 666 mpfr_set_si_2exp (t[1], q, n - 1, MPFR_RNDN); 667 inex1 = mpfr_add (s1, t[0], t[1], MPFR_RNDN); 668 MPFR_ASSERTN (inex1 == 0); 669 for (s = -1; s <= 1; s += 2) { 670 mpfr_neg (t[0], t[0], MPFR_RNDN); 671 mpfr_neg (t[1], t[1], MPFR_RNDN); 672 mpfr_neg (s1, s1, MPFR_RNDN); 673 for (h = -1; h <= 1; h += 2) { 674 mpfr_set_si (t[2], h, MPFR_RNDN); 675 inex1 = mpfr_add (s2, s1, t[2], MPFR_RNDN); 676 MPFR_ASSERTN (inex1 == 0); 677 for (i = -1; i <= 3; i += 2) { 678 mpfr_set_si_2exp (t[3], i, -2, MPFR_RNDN); 679 inex1 = mpfr_add (s3, s2, t[3], MPFR_RNDN); 680 MPFR_ASSERTN (inex1 == 0); 681 for (j = i; j <= 3; j++) { 682 mpfr_set_si_2exp (t[4], j, -2, MPFR_RNDN); 683 inex1 = mpfr_add (s4, s3, t[4], MPFR_RNDN); 684 MPFR_ASSERTN (inex1 == 0); 685 RND_LOOP_NO_RNDF (r) { 686 inex1 = mpfr_set (sum1, s4, (mpfr_rnd_t) r); 687 inex2 = mpfr_sum (sum2, p, 5, (mpfr_rnd_t) r); 688 MPFR_ASSERTN (mpfr_check (sum1)); 689 MPFR_ASSERTN (mpfr_check (sum2)); 690 if (!(mpfr_equal_p (sum1, sum2) && 691 SAME_SIGN (inex1, inex2))) 692 { 693 printf ("Error in check4 on %s, " 694 "k = %d, n = %d (prec %d), " 695 "q = %d, s = %d, h = %d, i = %d, j = %d\n", 696 mpfr_print_rnd_mode ((mpfr_rnd_t) r), 697 k, n, prec, q, s, h, i, j); 698 printf ("Expected "); 699 mpfr_dump (sum1); 700 printf ("with inex = %d\n", inex1); 701 printf ("Got "); 702 mpfr_dump (sum2); 703 printf ("with inex = %d\n", inex2); 704 exit (1); 705 } 706 } 707 } 708 } 709 } 710 } 711 } 712 } 713 } 714 715 for (i = 0; i < 5; i++) 716 mpfr_clear (t[i]); 717 mpfr_clears (sum1, sum2, s1, s2, s3, s4, (mpfr_ptr) 0); 718 } 719 720 /* bug reported by Joseph S. Myers on 2013-10-27 721 https://sympa.inria.fr/sympa/arc/mpfr/2013-10/msg00015.html */ 722 static void 723 bug20131027 (void) 724 { 725 mpfr_t sum, t[4]; 726 mpfr_ptr p[4]; 727 char *s[4] = { 728 "0x1p1000", 729 "-0x0.fffffffffffff80000000000000001p1000", 730 "-0x1p947", 731 "0x1p880" 732 }; 733 int i, r; 734 735 mpfr_init2 (sum, 53); 736 737 for (i = 0; i < 4; i++) 738 { 739 mpfr_init2 (t[i], i == 0 ? 53 : 1000); 740 mpfr_set_str (t[i], s[i], 0, MPFR_RNDN); 741 p[i] = t[i]; 742 } 743 744 RND_LOOP(r) 745 { 746 int expected_sign = (mpfr_rnd_t) r == MPFR_RNDD ? -1 : 1; 747 int inex; 748 749 inex = mpfr_sum (sum, p, 4, (mpfr_rnd_t) r); 750 751 if (MPFR_NOTZERO (sum) || MPFR_SIGN (sum) != expected_sign || inex != 0) 752 { 753 printf ("mpfr_sum incorrect in bug20131027 for %s:\n" 754 "expected %c0 with inex = 0, got ", 755 mpfr_print_rnd_mode ((mpfr_rnd_t) r), 756 expected_sign > 0 ? '+' : '-'); 757 mpfr_dump (sum); 758 printf ("with inex = %d\n", inex); 759 exit (1); 760 } 761 } 762 763 for (i = 0; i < 4; i++) 764 mpfr_clear (t[i]); 765 mpfr_clear (sum); 766 } 767 768 /* Occurs in branches/new-sum/src/sum.c@9344 on a 64-bit machine. */ 769 static void 770 bug20150327 (void) 771 { 772 mpfr_t sum1, sum2, t[3]; 773 mpfr_ptr p[3]; 774 char *s[3] = { "0.10000111110101000010101011100001", "1E-100", "0.1E95" }; 775 int i, r; 776 777 mpfr_inits2 (58, sum1, sum2, (mpfr_ptr) 0); 778 779 for (i = 0; i < 3; i++) 780 { 781 mpfr_init2 (t[i], 64); 782 mpfr_set_str (t[i], s[i], 2, MPFR_RNDN); 783 p[i] = t[i]; 784 } 785 786 RND_LOOP_NO_RNDF (r) 787 { 788 int inex1, inex2; 789 790 mpfr_set (sum1, t[2], MPFR_RNDN); 791 inex1 = -1; 792 if (MPFR_IS_LIKE_RNDU ((mpfr_rnd_t) r, 1)) 793 { 794 mpfr_nextabove (sum1); 795 inex1 = 1; 796 } 797 798 inex2 = mpfr_sum (sum2, p, 3, (mpfr_rnd_t) r); 799 800 if (!(mpfr_equal_p (sum1, sum2) && SAME_SIGN (inex1, inex2))) 801 { 802 printf ("mpfr_sum incorrect in bug20150327 for %s:\n", 803 mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 804 printf ("Expected "); 805 mpfr_dump (sum1); 806 printf ("with inex = %d\n", inex1); 807 printf ("Got "); 808 mpfr_dump (sum2); 809 printf ("with inex = %d\n", inex2); 810 exit (1); 811 } 812 } 813 814 for (i = 0; i < 3; i++) 815 mpfr_clear (t[i]); 816 mpfr_clears (sum1, sum2, (mpfr_ptr) 0); 817 } 818 819 /* TODO: A test with more inputs (but can't be compared to mpfr_add). */ 820 static void 821 check_extreme (void) 822 { 823 mpfr_t u, v, w, x, y; 824 mpfr_ptr t[2]; 825 int i, inex1, inex2, r; 826 827 t[0] = u; 828 t[1] = v; 829 830 mpfr_inits2 (32, u, v, w, x, y, (mpfr_ptr) 0); 831 mpfr_setmin (u, mpfr_get_emax ()); 832 mpfr_setmax (v, mpfr_get_emin ()); 833 mpfr_setmin (w, mpfr_get_emax () - 40); 834 RND_LOOP_NO_RNDF (r) 835 for (i = 0; i < 2; i++) 836 { 837 mpfr_set_prec (x, 64); 838 inex1 = mpfr_add (x, u, w, MPFR_RNDN); 839 MPFR_ASSERTN (inex1 == 0); 840 inex1 = mpfr_prec_round (x, 32, (mpfr_rnd_t) r); 841 inex2 = mpfr_sum (y, t, 2, (mpfr_rnd_t) r); 842 if (!(mpfr_equal_p (x, y) && SAME_SIGN (inex1, inex2))) 843 { 844 printf ("Error in check_extreme (%s, i = %d)\n", 845 mpfr_print_rnd_mode ((mpfr_rnd_t) r), i); 846 printf ("Expected "); 847 mpfr_dump (x); 848 printf ("with inex = %d\n", inex1); 849 printf ("Got "); 850 mpfr_dump (y); 851 printf ("with inex = %d\n", inex2); 852 exit (1); 853 } 854 mpfr_neg (v, v, MPFR_RNDN); 855 mpfr_neg (w, w, MPFR_RNDN); 856 } 857 mpfr_clears (u, v, w, x, y, (mpfr_ptr) 0); 858 } 859 860 /* Generic random tests with cancellations. 861 * 862 * In summary, we do 4000 tests of the following form: 863 * 1. We set the first MPFR_NCANCEL members of an array to random values, 864 * with a random exponent taken in 4 ranges, depending on the value of 865 * i % 4 (see code below). 866 * 2. For each of the next MPFR_NCANCEL iterations: 867 * A. we randomly permute some terms of the array (to make sure that a 868 * particular order doesn't have an influence on the result); 869 * B. we compute the sum in a random rounding mode; 870 * C. if this sum is zero, we end the current test (there is no longer 871 * anything interesting to test); 872 * D. we check that this sum is below some bound (chosen as infinite 873 * for the first iteration of (2), i.e. this test will be useful 874 * only for the next iterations, after cancellations); 875 * E. we put the opposite of this sum in the array, the goal being to 876 * introduce a chain of cancellations; 877 * F. we compute the bound for the next iteration, derived from (E). 878 * 3. We do another iteration like (2), but with reusing a random element 879 * of the array. This last test allows one to check the support of 880 * reused arguments. Before this support (r10467), it triggers an 881 * assertion failure with (almost?) all seeds, and if assertions are 882 * not checked, tsum fails in most cases but not all. 883 */ 884 static void 885 cancel (void) 886 { 887 mpfr_t x[2 * MPFR_NCANCEL], bound; 888 mpfr_ptr px[2 * MPFR_NCANCEL]; 889 int i, j, k, n; 890 891 mpfr_init2 (bound, 2); 892 893 /* With 4000 tests, tsum will fail in most cases without support of 894 reused arguments (before r10467). */ 895 for (i = 0; i < 4000; i++) 896 { 897 mpfr_set_inf (bound, 1); 898 for (n = 0; n <= numberof (x); n++) 899 { 900 mpfr_prec_t p; 901 mpfr_rnd_t rnd; 902 903 if (n < numberof (x)) 904 { 905 px[n] = x[n]; 906 p = MPFR_PREC_MIN + (randlimb () % 256); 907 mpfr_init2 (x[n], p); 908 k = n; 909 } 910 else 911 { 912 /* Reuse of a random member of the array. */ 913 k = randlimb () % n; 914 } 915 916 if (n < MPFR_NCANCEL) 917 { 918 mpfr_exp_t e; 919 920 MPFR_ASSERTN (n < numberof (x)); 921 e = (i & 1) ? 0 : mpfr_get_emin (); 922 tests_default_random (x[n], 256, e, 923 ((i & 2) ? e + 2000 : mpfr_get_emax ()), 924 0); 925 } 926 else 927 { 928 /* random permutation with n random transpositions */ 929 for (j = 0; j < n; j++) 930 { 931 int k1, k2; 932 933 k1 = randlimb () % (n-1); 934 k2 = randlimb () % (n-1); 935 mpfr_swap (x[k1], x[k2]); 936 } 937 938 rnd = RND_RAND (); 939 940 #if MPFR_DEBUG 941 printf ("mpfr_sum cancellation test\n"); 942 for (j = 0; j < n; j++) 943 { 944 printf (" x%d[%3ld] = ", j, mpfr_get_prec(x[j])); 945 mpfr_out_str (stdout, 16, 0, x[j], MPFR_RNDN); 946 printf ("\n"); 947 } 948 printf (" rnd = %s, output prec = %ld\n", 949 mpfr_print_rnd_mode (rnd), mpfr_get_prec (x[n])); 950 #endif 951 952 mpfr_sum (x[k], px, n, rnd); 953 954 if (mpfr_zero_p (x[k])) 955 { 956 if (k == n) 957 n++; 958 break; 959 } 960 961 if (mpfr_cmpabs (x[k], bound) > 0) 962 { 963 printf ("Error in cancel on i = %d, n = %d\n", i, n); 964 printf ("Expected bound: "); 965 mpfr_dump (bound); 966 printf ("x[%d]: ", k); 967 mpfr_dump (x[k]); 968 exit (1); 969 } 970 971 if (k != n) 972 break; 973 974 /* For the bound, use MPFR_RNDU due to possible underflow. 975 It would be nice to add some specific underflow checks, 976 though there are already ones in check_underflow(). */ 977 mpfr_set_ui_2exp (bound, 1, 978 mpfr_get_exp (x[n]) - p - (rnd == MPFR_RNDN), 979 MPFR_RNDU); 980 /* The next sum will be <= bound in absolute value 981 (the equality can be obtained in all rounding modes 982 since the sum will be rounded). */ 983 984 mpfr_neg (x[n], x[n], MPFR_RNDN); 985 } 986 } 987 988 while (--n >= 0) 989 mpfr_clear (x[n]); 990 } 991 992 mpfr_clear (bound); 993 } 994 995 #ifndef NOVFL 996 # define NOVFL 30 997 #endif 998 999 static void 1000 check_overflow (void) 1001 { 1002 mpfr_t sum1, sum2, x, y; 1003 mpfr_ptr t[2 * NOVFL]; 1004 mpfr_exp_t emin, emax; 1005 int i, r; 1006 1007 emin = mpfr_get_emin (); 1008 emax = mpfr_get_emax (); 1009 set_emin (MPFR_EMIN_MIN); 1010 set_emax (MPFR_EMAX_MAX); 1011 1012 mpfr_inits2 (32, sum1, sum2, x, y, (mpfr_ptr) 0); 1013 mpfr_setmax (x, mpfr_get_emax ()); 1014 mpfr_neg (y, x, MPFR_RNDN); 1015 1016 for (i = 0; i < 2 * NOVFL; i++) 1017 t[i] = i < NOVFL ? x : y; 1018 1019 /* Two kinds of test: 1020 * i = 1: overflow. 1021 * i = 2: intermediate overflow (exact sum is 0). 1022 */ 1023 for (i = 1; i <= 2; i++) 1024 RND_LOOP(r) 1025 { 1026 int inex1, inex2; 1027 1028 inex1 = mpfr_add (sum1, x, i == 1 ? x : y, (mpfr_rnd_t) r); 1029 inex2 = mpfr_sum (sum2, t, i * NOVFL, (mpfr_rnd_t) r); 1030 MPFR_ASSERTN (mpfr_check (sum1)); 1031 MPFR_ASSERTN (mpfr_check (sum2)); 1032 if (!(mpfr_equal_p (sum1, sum2) && SAME_SIGN (inex1, inex2))) 1033 { 1034 printf ("Error in check_overflow on %s, i = %d\n", 1035 mpfr_print_rnd_mode ((mpfr_rnd_t) r), i); 1036 printf ("Expected "); 1037 mpfr_dump (sum1); 1038 printf ("with inex = %d\n", inex1); 1039 printf ("Got "); 1040 mpfr_dump (sum2); 1041 printf ("with inex = %d\n", inex2); 1042 exit (1); 1043 } 1044 } 1045 1046 mpfr_clears (sum1, sum2, x, y, (mpfr_ptr) 0); 1047 1048 set_emin (emin); 1049 set_emax (emax); 1050 } 1051 1052 #ifndef NUNFL 1053 # define NUNFL 9 1054 #endif 1055 1056 static void 1057 check_underflow (void) 1058 { 1059 mpfr_t sum1, sum2, t[NUNFL]; 1060 mpfr_ptr p[NUNFL]; 1061 mpfr_prec_t precmax = 444; 1062 mpfr_exp_t emin, emax; 1063 unsigned int ex_flags, flags; 1064 int c, i; 1065 1066 emin = mpfr_get_emin (); 1067 emax = mpfr_get_emax (); 1068 set_emin (MPFR_EMIN_MIN); 1069 set_emax (MPFR_EMAX_MAX); 1070 1071 ex_flags = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT; 1072 1073 mpfr_init2 (sum1, MPFR_PREC_MIN); 1074 mpfr_init2 (sum2, precmax); 1075 1076 for (i = 0; i < NUNFL; i++) 1077 { 1078 mpfr_init2 (t[i], precmax); 1079 p[i] = t[i]; 1080 } 1081 1082 for (c = 0; c < 8; c++) 1083 { 1084 mpfr_prec_t fprec; 1085 int n, neg, r; 1086 1087 fprec = MPFR_PREC_MIN + (randlimb () % (precmax - MPFR_PREC_MIN + 1)); 1088 n = 3 + (randlimb () % (NUNFL - 2)); 1089 MPFR_ASSERTN (n <= NUNFL); 1090 1091 mpfr_set_prec (sum2, (randlimb () & 1) ? MPFR_PREC_MIN : precmax); 1092 mpfr_set_prec (t[0], fprec + 64); 1093 mpfr_set_zero (t[0], 1); 1094 1095 for (i = 1; i < n; i++) 1096 { 1097 int inex; 1098 1099 mpfr_set_prec (t[i], MPFR_PREC_MIN + 1100 (randlimb () % (fprec - MPFR_PREC_MIN + 1))); 1101 do 1102 mpfr_urandomb (t[i], RANDS); 1103 while (MPFR_IS_ZERO (t[i])); 1104 mpfr_set_exp (t[i], MPFR_EMIN_MIN); 1105 inex = mpfr_sub (t[0], t[0], t[i], MPFR_RNDN); 1106 MPFR_ASSERTN (inex == 0); 1107 } 1108 1109 neg = randlimb () & 1; 1110 if (neg) 1111 mpfr_nextbelow (t[0]); 1112 else 1113 mpfr_nextabove (t[0]); 1114 1115 RND_LOOP(r) 1116 { 1117 int inex1, inex2; 1118 1119 mpfr_set_zero (sum1, 1); 1120 if (neg) 1121 mpfr_nextbelow (sum1); 1122 else 1123 mpfr_nextabove (sum1); 1124 inex1 = mpfr_div_2ui (sum1, sum1, 2, (mpfr_rnd_t) r); 1125 1126 mpfr_clear_flags (); 1127 inex2 = mpfr_sum (sum2, p, n, (mpfr_rnd_t) r); 1128 flags = __gmpfr_flags; 1129 1130 MPFR_ASSERTN (mpfr_check (sum1)); 1131 MPFR_ASSERTN (mpfr_check (sum2)); 1132 1133 if (flags != ex_flags) 1134 { 1135 printf ("Bad flags in check_underflow on %s, c = %d\n", 1136 mpfr_print_rnd_mode ((mpfr_rnd_t) r), c); 1137 printf ("Expected flags:"); 1138 flags_out (ex_flags); 1139 printf ("Got flags: "); 1140 flags_out (flags); 1141 printf ("sum = "); 1142 mpfr_dump (sum2); 1143 exit (1); 1144 } 1145 1146 if (!(mpfr_equal_p (sum1, sum2) && SAME_SIGN (inex1, inex2))) 1147 { 1148 printf ("Error in check_underflow on %s, c = %d\n", 1149 mpfr_print_rnd_mode ((mpfr_rnd_t) r), c); 1150 printf ("Expected "); 1151 mpfr_dump (sum1); 1152 printf ("with inex = %d\n", inex1); 1153 printf ("Got "); 1154 mpfr_dump (sum2); 1155 printf ("with inex = %d\n", inex2); 1156 exit (1); 1157 } 1158 } 1159 } 1160 1161 for (i = 0; i < NUNFL; i++) 1162 mpfr_clear (t[i]); 1163 mpfr_clears (sum1, sum2, (mpfr_ptr) 0); 1164 1165 set_emin (emin); 1166 set_emax (emax); 1167 } 1168 1169 static void 1170 check_coverage (void) 1171 { 1172 #ifdef MPFR_COV_CHECK 1173 int r, i, j, k, p, q; 1174 int err = 0; 1175 1176 RND_LOOP_NO_RNDF (r) 1177 for (i = 0; i < 1 + ((mpfr_rnd_t) r == MPFR_RNDN); i++) 1178 for (j = 0; j < 2; j++) 1179 for (k = 0; k < 3; k++) 1180 for (p = 0; p < 2; p++) 1181 for (q = 0; q < 2; q++) 1182 if (!__gmpfr_cov_sum_tmd[r][i][j][k][p][q]) 1183 { 1184 printf ("TMD not tested on %s, tmd=%d, rbit=%d, sst=%d," 1185 " %s, sq %s MPFR_PREC_MIN\n", 1186 mpfr_print_rnd_mode ((mpfr_rnd_t) r), i+1, j, k-1, 1187 p ? "pos" : "neg", q ? ">" : "=="); 1188 err = 1; 1189 } 1190 1191 if (err) 1192 exit (1); 1193 #endif 1194 } 1195 1196 static int 1197 mpfr_sum_naive (mpfr_t s, mpfr_t *x, int n, mpfr_rnd_t rnd) 1198 { 1199 int ret, i; 1200 switch (n) 1201 { 1202 case 0: 1203 return mpfr_set_ui (s, 0, rnd); 1204 case 1: 1205 return mpfr_set (s, x[0], rnd); 1206 default: 1207 ret = mpfr_add (s, x[0], x[1], rnd); 1208 for (i = 2; i < n; i++) 1209 ret = mpfr_add (s, s, x[i], rnd); 1210 return ret; 1211 } 1212 } 1213 1214 /* check adding n random numbers, iterated k times */ 1215 static void 1216 check_random (int n, int k, mpfr_prec_t prec, mpfr_rnd_t rnd) 1217 { 1218 mpfr_t *x, s, ref_s; 1219 mpfr_ptr *y; 1220 int i, st, ret = 0, ref_ret = 0; 1221 gmp_randstate_t state; 1222 1223 gmp_randinit_default (state); 1224 mpfr_init2 (s, prec); 1225 mpfr_init2 (ref_s, prec); 1226 x = (mpfr_t *) tests_allocate (n * sizeof (mpfr_t)); 1227 y = (mpfr_ptr *) tests_allocate (n * sizeof (mpfr_ptr)); 1228 for (i = 0; i < n; i++) 1229 { 1230 y[i] = x[i]; 1231 mpfr_init2 (x[i], prec); 1232 mpfr_urandom (x[i], state, rnd); 1233 } 1234 1235 st = cputime (); 1236 for (i = 0; i < k; i++) 1237 ref_ret = mpfr_sum_naive (ref_s, x, n, rnd); 1238 printf ("mpfr_sum_naive took %dms\n", cputime () - st); 1239 1240 st = cputime (); 1241 for (i = 0; i < k; i++) 1242 ret = mpfr_sum (s, y, n, rnd); 1243 printf ("mpfr_sum took %dms\n", cputime () - st); 1244 1245 if (n <= 2) 1246 { 1247 MPFR_ASSERTN (mpfr_cmp (ref_s, s) == 0); 1248 MPFR_ASSERTN (ref_ret == ret); 1249 } 1250 1251 for (i = 0; i < n; i++) 1252 mpfr_clear (x[i]); 1253 tests_free (x, n * sizeof (mpfr_t)); 1254 tests_free (y, n * sizeof (mpfr_ptr)); 1255 mpfr_clear (s); 1256 mpfr_clear (ref_s); 1257 gmp_randclear (state); 1258 } 1259 1260 /* This bug appears when porting sum.c for MPFR 3.1.4 (0.11E826 is returned), 1261 but does not appear in the trunk. It was due to buggy MPFR_IS_LIKE_RNDD 1262 and MPFR_IS_LIKE_RNDU macros. The fix was done in r9295 in the trunk and 1263 it has been merged in r10234 in the 3.1 branch. Note: the bug would have 1264 been caught by generic_tests anyway, but a simple testcase is easier for 1265 checking with log messages (MPFR_LOG_ALL=1). */ 1266 static void 1267 bug20160315 (void) 1268 { 1269 mpfr_t r, t[4]; 1270 mpfr_ptr p[4]; 1271 char *s[4] = { "0.10E20", "0", "0.11E382", "0.10E826" }; 1272 int i; 1273 1274 mpfr_init2 (r, 2); 1275 for (i = 0; i < 4; i++) 1276 { 1277 mpfr_init2 (t[i], 2); 1278 mpfr_set_str_binary (t[i], s[i]); 1279 p[i] = t[i]; 1280 } 1281 mpfr_sum (r, p, 4, MPFR_RNDN); 1282 if (! mpfr_equal_p (r, t[3])) 1283 { 1284 printf ("Error in bug20160315.\n"); 1285 printf ("Expected "); 1286 mpfr_dump (t[3]); 1287 printf ("Got "); 1288 mpfr_dump (r); 1289 exit (1); 1290 } 1291 for (i = 0; i < 4; i++) 1292 mpfr_clear (t[i]); 1293 mpfr_clear (r); 1294 } 1295 1296 int 1297 main (int argc, char *argv[]) 1298 { 1299 int h; 1300 1301 if (argc == 5) 1302 { 1303 check_random (atoi (argv[1]), atoi (argv[2]), atoi (argv[3]), 1304 (mpfr_rnd_t) atoi (argv[4])); 1305 return 0; 1306 } 1307 1308 tests_start_mpfr (); 1309 1310 if (argc != 1) 1311 { 1312 fprintf (stderr, "Usage: tsum\n tsum n k prec rnd\n"); 1313 exit (1); 1314 } 1315 1316 check_simple (); 1317 check_special (); 1318 check_more_special (); 1319 for (h = 0; h <= 64; h++) 1320 check1 (h); 1321 check2 (); 1322 check3 (); 1323 check4 (); 1324 bug20131027 (); 1325 bug20150327 (); 1326 bug20160315 (); 1327 generic_tests (); 1328 check_extreme (); 1329 cancel (); 1330 check_overflow (); 1331 check_underflow (); 1332 1333 check_coverage (); 1334 tests_end_mpfr (); 1335 return 0; 1336 } 1337