1 /* tsum -- test file for the list summation function 2 3 Copyright 2004-2023 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 https://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_ptr 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 || RAND_BOOL ())) 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 void 175 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 void 200 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 const 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 const 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 const char *s[3] = { 775 "0.10000111110101000010101011100001", 776 "1E-100", 777 "0.1E95" 778 }; 779 int i, r; 780 781 mpfr_inits2 (58, sum1, sum2, (mpfr_ptr) 0); 782 783 for (i = 0; i < 3; i++) 784 { 785 mpfr_init2 (t[i], 64); 786 mpfr_set_str (t[i], s[i], 2, MPFR_RNDN); 787 p[i] = t[i]; 788 } 789 790 RND_LOOP_NO_RNDF (r) 791 { 792 int inex1, inex2; 793 794 mpfr_set (sum1, t[2], MPFR_RNDN); 795 inex1 = -1; 796 if (MPFR_IS_LIKE_RNDU ((mpfr_rnd_t) r, 1)) 797 { 798 mpfr_nextabove (sum1); 799 inex1 = 1; 800 } 801 802 inex2 = mpfr_sum (sum2, p, 3, (mpfr_rnd_t) r); 803 804 if (!(mpfr_equal_p (sum1, sum2) && SAME_SIGN (inex1, inex2))) 805 { 806 printf ("mpfr_sum incorrect in bug20150327 for %s:\n", 807 mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 808 printf ("Expected "); 809 mpfr_dump (sum1); 810 printf ("with inex = %d\n", inex1); 811 printf ("Got "); 812 mpfr_dump (sum2); 813 printf ("with inex = %d\n", inex2); 814 exit (1); 815 } 816 } 817 818 for (i = 0; i < 3; i++) 819 mpfr_clear (t[i]); 820 mpfr_clears (sum1, sum2, (mpfr_ptr) 0); 821 } 822 823 /* TODO: A test with more inputs (but can't be compared to mpfr_add). */ 824 static void 825 check_extreme (void) 826 { 827 mpfr_t u, v, w, x, y; 828 mpfr_ptr t[2]; 829 int i, inex1, inex2, r; 830 831 t[0] = u; 832 t[1] = v; 833 834 mpfr_inits2 (32, u, v, w, x, y, (mpfr_ptr) 0); 835 mpfr_setmin (u, mpfr_get_emax ()); 836 mpfr_setmax (v, mpfr_get_emin ()); 837 mpfr_setmin (w, mpfr_get_emax () - 40); 838 RND_LOOP_NO_RNDF (r) 839 for (i = 0; i < 2; i++) 840 { 841 mpfr_set_prec (x, 64); 842 inex1 = mpfr_add (x, u, w, MPFR_RNDN); 843 MPFR_ASSERTN (inex1 == 0); 844 inex1 = mpfr_prec_round (x, 32, (mpfr_rnd_t) r); 845 inex2 = mpfr_sum (y, t, 2, (mpfr_rnd_t) r); 846 if (!(mpfr_equal_p (x, y) && SAME_SIGN (inex1, inex2))) 847 { 848 printf ("Error in check_extreme (%s, i = %d)\n", 849 mpfr_print_rnd_mode ((mpfr_rnd_t) r), i); 850 printf ("Expected "); 851 mpfr_dump (x); 852 printf ("with inex = %d\n", inex1); 853 printf ("Got "); 854 mpfr_dump (y); 855 printf ("with inex = %d\n", inex2); 856 exit (1); 857 } 858 mpfr_neg (v, v, MPFR_RNDN); 859 mpfr_neg (w, w, MPFR_RNDN); 860 } 861 mpfr_clears (u, v, w, x, y, (mpfr_ptr) 0); 862 } 863 864 /* Generic random tests with cancellations. 865 * 866 * In summary, we do 4000 tests of the following form: 867 * 1. We set the first MPFR_NCANCEL members of an array to random values, 868 * with a random exponent taken in 4 ranges, depending on the value of 869 * i % 4 (see code below). 870 * 2. For each of the next MPFR_NCANCEL iterations: 871 * A. we randomly permute some terms of the array (to make sure that a 872 * particular order doesn't have an influence on the result); 873 * B. we compute the sum in a random rounding mode; 874 * C. if this sum is zero, we end the current test (there is no longer 875 * anything interesting to test); 876 * D. we check that this sum is below some bound (chosen as infinite 877 * for the first iteration of (2), i.e. this test will be useful 878 * only for the next iterations, after cancellations); 879 * E. we put the opposite of this sum in the array, the goal being to 880 * introduce a chain of cancellations; 881 * F. we compute the bound for the next iteration, derived from (E). 882 * 3. We do another iteration like (2), but with reusing a random element 883 * of the array. This last test allows one to check the support of 884 * reused arguments. Before this support (r10467), it triggers an 885 * assertion failure with (almost?) all seeds, and if assertions are 886 * not checked, tsum fails in most cases but not all. 887 */ 888 static void 889 cancel (void) 890 { 891 mpfr_t x[2 * MPFR_NCANCEL], bound; 892 mpfr_ptr px[2 * MPFR_NCANCEL]; 893 int i, j, k, n; 894 895 mpfr_init2 (bound, 2); 896 897 /* With 4000 tests, tsum will fail in most cases without support of 898 reused arguments (before r10467). */ 899 for (i = 0; i < 4000; i++) 900 { 901 mpfr_set_inf (bound, 1); 902 for (n = 0; n <= numberof (x); n++) 903 { 904 mpfr_prec_t p; 905 mpfr_rnd_t rnd; 906 907 if (n < numberof (x)) 908 { 909 px[n] = x[n]; 910 p = MPFR_PREC_MIN + (randlimb () % 256); 911 mpfr_init2 (x[n], p); 912 k = n; 913 } 914 else 915 { 916 /* Reuse of a random member of the array. */ 917 k = randlimb () % n; 918 } 919 920 if (n < MPFR_NCANCEL) 921 { 922 mpfr_exp_t e; 923 924 MPFR_ASSERTN (n < numberof (x)); 925 e = (i & 1) ? 0 : mpfr_get_emin (); 926 tests_default_random (x[n], 256, e, 927 ((i & 2) ? e + 2000 : mpfr_get_emax ()), 928 0); 929 } 930 else 931 { 932 /* random permutation with n random transpositions */ 933 for (j = 0; j < n; j++) 934 { 935 int k1, k2; 936 937 k1 = randlimb () % (n-1); 938 k2 = randlimb () % (n-1); 939 mpfr_swap (x[k1], x[k2]); 940 } 941 942 rnd = RND_RAND (); 943 944 #if MPFR_DEBUG 945 printf ("mpfr_sum cancellation test\n"); 946 for (j = 0; j < n; j++) 947 { 948 printf (" x%d[%3ld] = ", j, mpfr_get_prec(x[j])); 949 mpfr_out_str (stdout, 16, 0, x[j], MPFR_RNDN); 950 printf ("\n"); 951 } 952 printf (" rnd = %s, output prec = %ld\n", 953 mpfr_print_rnd_mode (rnd), mpfr_get_prec (x[n])); 954 #endif 955 956 mpfr_sum (x[k], px, n, rnd); 957 958 if (mpfr_zero_p (x[k])) 959 { 960 if (k == n) 961 n++; 962 break; 963 } 964 965 if (mpfr_cmpabs (x[k], bound) > 0) 966 { 967 printf ("Error in cancel on i = %d, n = %d\n", i, n); 968 printf ("Expected bound: "); 969 mpfr_dump (bound); 970 printf ("x[%d]: ", k); 971 mpfr_dump (x[k]); 972 exit (1); 973 } 974 975 if (k != n) 976 break; 977 978 /* For the bound, use MPFR_RNDU due to possible underflow. 979 It would be nice to add some specific underflow checks, 980 though there are already ones in check_underflow(). */ 981 mpfr_set_ui_2exp (bound, 1, 982 mpfr_get_exp (x[n]) - p - (rnd == MPFR_RNDN), 983 MPFR_RNDU); 984 /* The next sum will be <= bound in absolute value 985 (the equality can be obtained in all rounding modes 986 since the sum will be rounded). */ 987 988 mpfr_neg (x[n], x[n], MPFR_RNDN); 989 } 990 } 991 992 while (--n >= 0) 993 mpfr_clear (x[n]); 994 } 995 996 mpfr_clear (bound); 997 } 998 999 #ifndef NOVFL 1000 # define NOVFL 30 1001 #endif 1002 1003 static void 1004 check_overflow (void) 1005 { 1006 mpfr_t sum1, sum2, x, y; 1007 mpfr_ptr t[2 * NOVFL]; 1008 mpfr_exp_t emin, emax; 1009 int i, r; 1010 1011 emin = mpfr_get_emin (); 1012 emax = mpfr_get_emax (); 1013 set_emin (MPFR_EMIN_MIN); 1014 set_emax (MPFR_EMAX_MAX); 1015 1016 mpfr_inits2 (32, sum1, sum2, x, y, (mpfr_ptr) 0); 1017 mpfr_setmax (x, mpfr_get_emax ()); 1018 mpfr_neg (y, x, MPFR_RNDN); 1019 1020 for (i = 0; i < 2 * NOVFL; i++) 1021 t[i] = i < NOVFL ? x : y; 1022 1023 /* Two kinds of test: 1024 * i = 1: overflow. 1025 * i = 2: intermediate overflow (exact sum is 0). 1026 */ 1027 for (i = 1; i <= 2; i++) 1028 RND_LOOP(r) 1029 { 1030 int inex1, inex2; 1031 1032 inex1 = mpfr_add (sum1, x, i == 1 ? x : y, (mpfr_rnd_t) r); 1033 inex2 = mpfr_sum (sum2, t, i * NOVFL, (mpfr_rnd_t) r); 1034 MPFR_ASSERTN (mpfr_check (sum1)); 1035 MPFR_ASSERTN (mpfr_check (sum2)); 1036 if (!(mpfr_equal_p (sum1, sum2) && SAME_SIGN (inex1, inex2))) 1037 { 1038 printf ("Error in check_overflow on %s, i = %d\n", 1039 mpfr_print_rnd_mode ((mpfr_rnd_t) r), i); 1040 printf ("Expected "); 1041 mpfr_dump (sum1); 1042 printf ("with inex = %d\n", inex1); 1043 printf ("Got "); 1044 mpfr_dump (sum2); 1045 printf ("with inex = %d\n", inex2); 1046 exit (1); 1047 } 1048 } 1049 1050 mpfr_clears (sum1, sum2, x, y, (mpfr_ptr) 0); 1051 1052 set_emin (emin); 1053 set_emax (emax); 1054 } 1055 1056 #ifndef NUNFL 1057 # define NUNFL 9 1058 #endif 1059 1060 static void 1061 check_underflow (void) 1062 { 1063 mpfr_t sum1, sum2, t[NUNFL]; 1064 mpfr_ptr p[NUNFL]; 1065 mpfr_prec_t precmax = 444; 1066 mpfr_exp_t emin, emax; 1067 unsigned int ex_flags, flags; 1068 int c, i; 1069 1070 emin = mpfr_get_emin (); 1071 emax = mpfr_get_emax (); 1072 set_emin (MPFR_EMIN_MIN); 1073 set_emax (MPFR_EMAX_MAX); 1074 1075 ex_flags = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT; 1076 1077 mpfr_init2 (sum1, MPFR_PREC_MIN); 1078 mpfr_init2 (sum2, precmax); 1079 1080 for (i = 0; i < NUNFL; i++) 1081 { 1082 mpfr_init2 (t[i], precmax); 1083 p[i] = t[i]; 1084 } 1085 1086 for (c = 0; c < 8; c++) 1087 { 1088 mpfr_prec_t fprec; 1089 int n, neg, r; 1090 1091 fprec = MPFR_PREC_MIN + (randlimb () % (precmax - MPFR_PREC_MIN + 1)); 1092 n = 3 + (randlimb () % (NUNFL - 2)); 1093 MPFR_ASSERTN (n <= NUNFL); 1094 1095 mpfr_set_prec (sum2, RAND_BOOL () ? MPFR_PREC_MIN : precmax); 1096 mpfr_set_prec (t[0], fprec + 64); 1097 mpfr_set_zero (t[0], 1); 1098 1099 for (i = 1; i < n; i++) 1100 { 1101 int inex; 1102 1103 mpfr_set_prec (t[i], MPFR_PREC_MIN + 1104 (randlimb () % (fprec - MPFR_PREC_MIN + 1))); 1105 do 1106 mpfr_urandomb (t[i], RANDS); 1107 while (MPFR_IS_ZERO (t[i])); 1108 mpfr_set_exp (t[i], MPFR_EMIN_MIN); 1109 inex = mpfr_sub (t[0], t[0], t[i], MPFR_RNDN); 1110 MPFR_ASSERTN (inex == 0); 1111 } 1112 1113 neg = RAND_BOOL (); 1114 if (neg) 1115 mpfr_nextbelow (t[0]); 1116 else 1117 mpfr_nextabove (t[0]); 1118 1119 RND_LOOP(r) 1120 { 1121 int inex1, inex2; 1122 1123 mpfr_set_zero (sum1, 1); 1124 if (neg) 1125 mpfr_nextbelow (sum1); 1126 else 1127 mpfr_nextabove (sum1); 1128 inex1 = mpfr_div_2ui (sum1, sum1, 2, (mpfr_rnd_t) r); 1129 1130 mpfr_clear_flags (); 1131 inex2 = mpfr_sum (sum2, p, n, (mpfr_rnd_t) r); 1132 flags = __gmpfr_flags; 1133 1134 MPFR_ASSERTN (mpfr_check (sum1)); 1135 MPFR_ASSERTN (mpfr_check (sum2)); 1136 1137 if (flags != ex_flags) 1138 { 1139 printf ("Bad flags in check_underflow on %s, c = %d\n", 1140 mpfr_print_rnd_mode ((mpfr_rnd_t) r), c); 1141 printf ("Expected flags:"); 1142 flags_out (ex_flags); 1143 printf ("Got flags: "); 1144 flags_out (flags); 1145 printf ("sum = "); 1146 mpfr_dump (sum2); 1147 exit (1); 1148 } 1149 1150 if (!(mpfr_equal_p (sum1, sum2) && SAME_SIGN (inex1, inex2))) 1151 { 1152 printf ("Error in check_underflow on %s, c = %d\n", 1153 mpfr_print_rnd_mode ((mpfr_rnd_t) r), c); 1154 printf ("Expected "); 1155 mpfr_dump (sum1); 1156 printf ("with inex = %d\n", inex1); 1157 printf ("Got "); 1158 mpfr_dump (sum2); 1159 printf ("with inex = %d\n", inex2); 1160 exit (1); 1161 } 1162 } 1163 } 1164 1165 for (i = 0; i < NUNFL; i++) 1166 mpfr_clear (t[i]); 1167 mpfr_clears (sum1, sum2, (mpfr_ptr) 0); 1168 1169 set_emin (emin); 1170 set_emax (emax); 1171 } 1172 1173 static void 1174 check_coverage (void) 1175 { 1176 #ifdef MPFR_COV_CHECK 1177 int r, i, j, k, p, q; 1178 int err = 0; 1179 1180 RND_LOOP_NO_RNDF (r) 1181 for (i = 0; i < 1 + ((mpfr_rnd_t) r == MPFR_RNDN); i++) 1182 for (j = 0; j < 2; j++) 1183 for (k = 0; k < 3; k++) 1184 for (p = 0; p < 2; p++) 1185 for (q = 0; q < 2; q++) 1186 if (!__gmpfr_cov_sum_tmd[r][i][j][k][p][q]) 1187 { 1188 printf ("TMD not tested on %s, tmd=%d, rbit=%d, sst=%d," 1189 " %s, sq %s MPFR_PREC_MIN\n", 1190 mpfr_print_rnd_mode ((mpfr_rnd_t) r), i+1, j, k-1, 1191 p ? "pos" : "neg", q ? ">" : "=="); 1192 err = 1; 1193 } 1194 1195 if (err) 1196 exit (1); 1197 #endif 1198 } 1199 1200 static int 1201 mpfr_sum_naive (mpfr_ptr s, mpfr_t *x, int n, mpfr_rnd_t rnd) 1202 { 1203 int ret, i; 1204 switch (n) 1205 { 1206 case 0: 1207 return mpfr_set_ui (s, 0, rnd); 1208 case 1: 1209 return mpfr_set (s, x[0], rnd); 1210 default: 1211 ret = mpfr_add (s, x[0], x[1], rnd); 1212 for (i = 2; i < n; i++) 1213 ret = mpfr_add (s, s, x[i], rnd); 1214 return ret; 1215 } 1216 } 1217 1218 /* check adding n random numbers, iterated k times */ 1219 static void 1220 check_random (int n, int k, mpfr_prec_t prec, mpfr_rnd_t rnd) 1221 { 1222 mpfr_t *x, s, ref_s; 1223 mpfr_ptr *y; 1224 int i, st, ret = 0, ref_ret = 0; 1225 gmp_randstate_t state; 1226 1227 gmp_randinit_default (state); 1228 mpfr_init2 (s, prec); 1229 mpfr_init2 (ref_s, prec); 1230 x = (mpfr_t *) tests_allocate (n * sizeof (mpfr_t)); 1231 y = (mpfr_ptr *) tests_allocate (n * sizeof (mpfr_ptr)); 1232 for (i = 0; i < n; i++) 1233 { 1234 y[i] = x[i]; 1235 mpfr_init2 (x[i], prec); 1236 mpfr_urandom (x[i], state, rnd); 1237 } 1238 1239 st = cputime (); 1240 for (i = 0; i < k; i++) 1241 ref_ret = mpfr_sum_naive (ref_s, x, n, rnd); 1242 printf ("mpfr_sum_naive took %dms\n", cputime () - st); 1243 1244 st = cputime (); 1245 for (i = 0; i < k; i++) 1246 ret = mpfr_sum (s, y, n, rnd); 1247 printf ("mpfr_sum took %dms\n", cputime () - st); 1248 1249 if (n <= 2) 1250 { 1251 MPFR_ASSERTN (mpfr_cmp (ref_s, s) == 0); 1252 MPFR_ASSERTN (ref_ret == ret); 1253 } 1254 1255 for (i = 0; i < n; i++) 1256 mpfr_clear (x[i]); 1257 tests_free (x, n * sizeof (mpfr_t)); 1258 tests_free (y, n * sizeof (mpfr_ptr)); 1259 mpfr_clear (s); 1260 mpfr_clear (ref_s); 1261 gmp_randclear (state); 1262 } 1263 1264 /* This bug appears when porting sum.c for MPFR 3.1.4 (0.11E826 is returned), 1265 but does not appear in the trunk. It was due to buggy MPFR_IS_LIKE_RNDD 1266 and MPFR_IS_LIKE_RNDU macros. The fix was done in r9295 in the trunk and 1267 it has been merged in r10234 in the 3.1 branch. Note: the bug would have 1268 been caught by generic_tests anyway, but a simple testcase is easier for 1269 checking with log messages (MPFR_LOG_ALL=1). */ 1270 static void 1271 bug20160315 (void) 1272 { 1273 mpfr_t r, t[4]; 1274 mpfr_ptr p[4]; 1275 const char *s[4] = { "0.10E20", "0", "0.11E382", "0.10E826" }; 1276 int i; 1277 1278 mpfr_init2 (r, 2); 1279 for (i = 0; i < 4; i++) 1280 { 1281 mpfr_init2 (t[i], 2); 1282 mpfr_set_str_binary (t[i], s[i]); 1283 p[i] = t[i]; 1284 } 1285 mpfr_sum (r, p, 4, MPFR_RNDN); 1286 if (! mpfr_equal_p (r, t[3])) 1287 { 1288 printf ("Error in bug20160315.\n"); 1289 printf ("Expected "); 1290 mpfr_dump (t[3]); 1291 printf ("Got "); 1292 mpfr_dump (r); 1293 exit (1); 1294 } 1295 for (i = 0; i < 4; i++) 1296 mpfr_clear (t[i]); 1297 mpfr_clear (r); 1298 } 1299 1300 int 1301 main (int argc, char *argv[]) 1302 { 1303 int h; 1304 1305 if (argc == 5) 1306 { 1307 check_random (atoi (argv[1]), atoi (argv[2]), atoi (argv[3]), 1308 (mpfr_rnd_t) atoi (argv[4])); 1309 return 0; 1310 } 1311 1312 tests_start_mpfr (); 1313 1314 if (argc != 1) 1315 { 1316 fprintf (stderr, "Usage: tsum\n tsum n k prec rnd\n"); 1317 exit (1); 1318 } 1319 1320 check_simple (); 1321 check_special (); 1322 check_more_special (); 1323 for (h = 0; h <= 64; h++) 1324 check1 (h); 1325 check2 (); 1326 check3 (); 1327 check4 (); 1328 bug20131027 (); 1329 bug20150327 (); 1330 bug20160315 (); 1331 generic_tests (); 1332 check_extreme (); 1333 cancel (); 1334 check_overflow (); 1335 check_underflow (); 1336 1337 check_coverage (); 1338 tests_end_mpfr (); 1339 return 0; 1340 } 1341