1/* tsum -- test file for the list summation function 2 3Copyright 2004-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba 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 20https://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 "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 */ 32static int 33cputime (void) 34{ 35 return clock () / (CLOCKS_PER_SEC / 1000); 36} 37 38static mpfr_prec_t 39get_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 61static void 62get_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 76static void 77generic_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 */ 174static void 175check_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 199static void 200check_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 285static void 286check_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. */ 373static void 374check1 (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. */ 453static void 454check2 (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 */ 540static void 541check3 (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 */ 644static void 645check4 (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 */ 722static void 723bug20131027 (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. */ 769static void 770bug20150327 (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). */ 824static void 825check_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 */ 888static void 889cancel (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 1003static void 1004check_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 1060static void 1061check_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 1173static void 1174check_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 1200static int 1201mpfr_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 */ 1219static void 1220check_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). */ 1270static void 1271bug20160315 (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 1300int 1301main (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