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
cputime(void)33 cputime (void)
34 {
35 return clock () / (CLOCKS_PER_SEC / 1000);
36 }
37
38 static mpfr_prec_t
get_prec_max(mpfr_t * t,int n)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
get_exact_sum(mpfr_ptr sum,mpfr_t * tab,int n)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
generic_tests(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
check_simple(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
check_special(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
check_more_special(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
check1(int h)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
check2(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
check3(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
check4(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
bug20131027(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
bug20150327(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
check_extreme(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
cancel(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
check_overflow(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
check_underflow(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
check_coverage(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
mpfr_sum_naive(mpfr_ptr s,mpfr_t * x,int n,mpfr_rnd_t rnd)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
check_random(int n,int k,mpfr_prec_t prec,mpfr_rnd_t rnd)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
bug20160315(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
main(int argc,char * argv[])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