1 /* Test file for mpfr_mul.
2
3 Copyright 1999, 2001-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 #ifdef CHECK_EXTERNAL
26 static int
test_mul(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode)27 test_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
28 {
29 int res;
30 int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_number_p (c);
31 if (ok)
32 {
33 mpfr_print_raw (b);
34 printf (" ");
35 mpfr_print_raw (c);
36 }
37 res = mpfr_mul (a, b, c, rnd_mode);
38 if (ok)
39 {
40 printf (" ");
41 mpfr_print_raw (a);
42 printf ("\n");
43 }
44 return res;
45 }
46 #else
47 #define test_mul mpfr_mul
48 #endif
49
50 /* checks that xs * ys gives the expected result res */
51 static void
check(const char * xs,const char * ys,mpfr_rnd_t rnd_mode,unsigned int px,unsigned int py,unsigned int pz,const char * res)52 check (const char *xs, const char *ys, mpfr_rnd_t rnd_mode,
53 unsigned int px, unsigned int py, unsigned int pz, const char *res)
54 {
55 mpfr_t xx, yy, zz;
56
57 mpfr_init2 (xx, px);
58 mpfr_init2 (yy, py);
59 mpfr_init2 (zz, pz);
60 mpfr_set_str1 (xx, xs);
61 mpfr_set_str1 (yy, ys);
62 test_mul(zz, xx, yy, rnd_mode);
63 if (mpfr_cmp_str1 (zz, res) )
64 {
65 printf ("(1) mpfr_mul failed for x=%s y=%s with rnd=%s\n",
66 xs, ys, mpfr_print_rnd_mode (rnd_mode));
67 printf ("correct is %s, mpfr_mul gives ", res);
68 mpfr_out_str (stdout, 10, 0, zz, MPFR_RNDN);
69 putchar ('\n');
70 exit (1);
71 }
72 mpfr_clear (xx);
73 mpfr_clear (yy);
74 mpfr_clear (zz);
75 }
76
77 static void
check53(const char * xs,const char * ys,mpfr_rnd_t rnd_mode,const char * zs)78 check53 (const char *xs, const char *ys, mpfr_rnd_t rnd_mode, const char *zs)
79 {
80 mpfr_t xx, yy, zz;
81
82 mpfr_inits2 (53, xx, yy, zz, (mpfr_ptr) 0);
83 mpfr_set_str1 (xx, xs);
84 mpfr_set_str1 (yy, ys);
85 test_mul (zz, xx, yy, rnd_mode);
86 if (mpfr_cmp_str1 (zz, zs) )
87 {
88 printf ("(2) mpfr_mul failed for x=%s y=%s with rnd=%s\n",
89 xs, ys, mpfr_print_rnd_mode(rnd_mode));
90 printf ("correct result is %s,\n mpfr_mul gives ", zs);
91 mpfr_out_str (stdout, 10, 0, zz, MPFR_RNDN);
92 putchar ('\n');
93 exit (1);
94 }
95 mpfr_clears (xx, yy, zz, (mpfr_ptr) 0);
96 }
97
98 /* checks that x*y gives the right result with 24 bits of precision */
99 static void
check24(const char * xs,const char * ys,mpfr_rnd_t rnd_mode,const char * zs)100 check24 (const char *xs, const char *ys, mpfr_rnd_t rnd_mode, const char *zs)
101 {
102 mpfr_t xx, yy, zz;
103
104 mpfr_inits2 (24, xx, yy, zz, (mpfr_ptr) 0);
105 mpfr_set_str1 (xx, xs);
106 mpfr_set_str1 (yy, ys);
107 test_mul (zz, xx, yy, rnd_mode);
108 if (mpfr_cmp_str1 (zz, zs) )
109 {
110 printf ("(3) mpfr_mul failed for x=%s y=%s with "
111 "rnd=%s\n", xs, ys, mpfr_print_rnd_mode(rnd_mode));
112 printf ("correct result is gives %s, mpfr_mul gives ", zs);
113 mpfr_out_str (stdout, 10, 0, zz, MPFR_RNDN);
114 putchar ('\n');
115 exit (1);
116 }
117 mpfr_clears (xx, yy, zz, (mpfr_ptr) 0);
118 }
119
120 /* the following examples come from the paper "Number-theoretic Test
121 Generation for Directed Rounding" from Michael Parks, Table 1 */
122 static void
check_float(void)123 check_float (void)
124 {
125 check24("8388609.0", "8388609.0", MPFR_RNDN, "70368760954880.0");
126 check24("16777213.0", "8388609.0", MPFR_RNDN, "140737479966720.0");
127 check24("8388611.0", "8388609.0", MPFR_RNDN, "70368777732096.0");
128 check24("12582911.0", "8388610.0", MPFR_RNDN, "105553133043712.0");
129 check24("12582914.0", "8388610.0", MPFR_RNDN, "105553158209536.0");
130 check24("13981013.0", "8388611.0", MPFR_RNDN, "117281279442944.0");
131 check24("11184811.0", "8388611.0", MPFR_RNDN, "93825028587520.0");
132 check24("11184810.0", "8388611.0", MPFR_RNDN, "93825020198912.0");
133 check24("13981014.0", "8388611.0", MPFR_RNDN, "117281287831552.0");
134
135 check24("8388609.0", "8388609.0", MPFR_RNDZ, "70368760954880.0");
136 check24("16777213.0", "8388609.0", MPFR_RNDZ, "140737471578112.0");
137 check24("8388611.0", "8388609.0", MPFR_RNDZ, "70368777732096.0");
138 check24("12582911.0", "8388610.0", MPFR_RNDZ, "105553124655104.0");
139 check24("12582914.0", "8388610.0", MPFR_RNDZ, "105553158209536.0");
140 check24("13981013.0", "8388611.0", MPFR_RNDZ, "117281271054336.0");
141 check24("11184811.0", "8388611.0", MPFR_RNDZ, "93825028587520.0");
142 check24("11184810.0", "8388611.0", MPFR_RNDZ, "93825011810304.0");
143 check24("13981014.0", "8388611.0", MPFR_RNDZ, "117281287831552.0");
144
145 check24("8388609.0", "8388609.0", MPFR_RNDU, "70368769343488.0");
146 check24("16777213.0", "8388609.0", MPFR_RNDU, "140737479966720.0");
147 check24("8388611.0", "8388609.0", MPFR_RNDU, "70368786120704.0");
148 check24("12582911.0", "8388610.0", MPFR_RNDU, "105553133043712.0");
149 check24("12582914.0", "8388610.0", MPFR_RNDU, "105553166598144.0");
150 check24("13981013.0", "8388611.0", MPFR_RNDU, "117281279442944.0");
151 check24("11184811.0", "8388611.0", MPFR_RNDU, "93825036976128.0");
152 check24("11184810.0", "8388611.0", MPFR_RNDU, "93825020198912.0");
153 check24("13981014.0", "8388611.0", MPFR_RNDU, "117281296220160.0");
154
155 check24("8388609.0", "8388609.0", MPFR_RNDD, "70368760954880.0");
156 check24("16777213.0", "8388609.0", MPFR_RNDD, "140737471578112.0");
157 check24("8388611.0", "8388609.0", MPFR_RNDD, "70368777732096.0");
158 check24("12582911.0", "8388610.0", MPFR_RNDD, "105553124655104.0");
159 check24("12582914.0", "8388610.0", MPFR_RNDD, "105553158209536.0");
160 check24("13981013.0", "8388611.0", MPFR_RNDD, "117281271054336.0");
161 check24("11184811.0", "8388611.0", MPFR_RNDD, "93825028587520.0");
162 check24("11184810.0", "8388611.0", MPFR_RNDD, "93825011810304.0");
163 check24("13981014.0", "8388611.0", MPFR_RNDD, "117281287831552.0");
164 }
165
166 /* check sign of result */
167 static void
check_sign(void)168 check_sign (void)
169 {
170 mpfr_t a, b;
171
172 mpfr_init2 (a, 53);
173 mpfr_init2 (b, 53);
174 mpfr_set_si (a, -1, MPFR_RNDN);
175 mpfr_set_ui (b, 2, MPFR_RNDN);
176 test_mul(a, b, b, MPFR_RNDN);
177 if (mpfr_cmp_ui (a, 4) )
178 {
179 printf ("2.0*2.0 gives \n");
180 mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
181 putchar ('\n');
182 exit (1);
183 }
184 mpfr_clear(a); mpfr_clear(b);
185 }
186
187 /* checks that the inexact return value is correct */
188 static void
check_exact(void)189 check_exact (void)
190 {
191 mpfr_t a, b, c, d;
192 mpfr_prec_t prec;
193 int i, inexact;
194 mpfr_rnd_t rnd;
195
196 mpfr_init (a);
197 mpfr_init (b);
198 mpfr_init (c);
199 mpfr_init (d);
200
201 mpfr_set_prec (a, 17);
202 mpfr_set_prec (b, 17);
203 mpfr_set_prec (c, 32);
204 mpfr_set_str_binary (a, "1.1000111011000100e-1");
205 mpfr_set_str_binary (b, "1.0010001111100111e-1");
206 if (test_mul (c, a, b, MPFR_RNDZ))
207 {
208 printf ("wrong return value (1)\n");
209 exit (1);
210 }
211
212 for (prec = MPFR_PREC_MIN; prec < 100; prec++)
213 {
214 mpfr_set_prec (a, prec);
215 mpfr_set_prec (b, prec);
216 /* for prec=1, ensure PREC(c) >= 1 */
217 mpfr_set_prec (c, 2 * prec - 2 + (prec == 1));
218 mpfr_set_prec (d, 2 * prec);
219 for (i = 0; i < 1000; i++)
220 {
221 mpfr_urandomb (a, RANDS);
222 mpfr_urandomb (b, RANDS);
223 rnd = RND_RAND ();
224 inexact = test_mul (c, a, b, rnd);
225 if (test_mul (d, a, b, rnd)) /* should be always exact */
226 {
227 printf ("unexpected inexact return value\n");
228 exit (1);
229 }
230 if ((inexact == 0) && mpfr_cmp (c, d) && rnd != MPFR_RNDF)
231 {
232 printf ("rnd=%s: inexact=0 but results differ\n",
233 mpfr_print_rnd_mode (rnd));
234 printf ("a=");
235 mpfr_out_str (stdout, 2, 0, a, rnd);
236 printf ("\nb=");
237 mpfr_out_str (stdout, 2, 0, b, rnd);
238 printf ("\nc=");
239 mpfr_out_str (stdout, 2, 0, c, rnd);
240 printf ("\nd=");
241 mpfr_out_str (stdout, 2, 0, d, rnd);
242 printf ("\n");
243 exit (1);
244 }
245 else if (inexact && (mpfr_cmp (c, d) == 0) && rnd != MPFR_RNDF)
246 {
247 printf ("inexact!=0 but results agree\n");
248 printf ("prec=%u rnd=%s a=", (unsigned int) prec,
249 mpfr_print_rnd_mode (rnd));
250 mpfr_out_str (stdout, 2, 0, a, rnd);
251 printf ("\nb=");
252 mpfr_out_str (stdout, 2, 0, b, rnd);
253 printf ("\nc=");
254 mpfr_out_str (stdout, 2, 0, c, rnd);
255 printf ("\nd=");
256 mpfr_out_str (stdout, 2, 0, d, rnd);
257 printf ("\n");
258 exit (1);
259 }
260 }
261 }
262
263 mpfr_clear (a);
264 mpfr_clear (b);
265 mpfr_clear (c);
266 mpfr_clear (d);
267 }
268
269 static void
check_max(void)270 check_max (void)
271 {
272 mpfr_t xx, yy, zz;
273 mpfr_exp_t emin;
274 int inex;
275
276 mpfr_init2(xx, 4);
277 mpfr_init2(yy, 4);
278 mpfr_init2(zz, 4);
279 mpfr_set_str1 (xx, "0.68750");
280 mpfr_mul_2si(xx, xx, MPFR_EMAX_DEFAULT/2, MPFR_RNDN);
281 mpfr_set_str1 (yy, "0.68750");
282 mpfr_mul_2si(yy, yy, MPFR_EMAX_DEFAULT - MPFR_EMAX_DEFAULT/2 + 1, MPFR_RNDN);
283 mpfr_clear_flags();
284 test_mul(zz, xx, yy, MPFR_RNDU);
285 if (!(mpfr_overflow_p() && MPFR_IS_INF(zz)))
286 {
287 printf("check_max failed (should be an overflow)\n");
288 exit(1);
289 }
290
291 mpfr_clear_flags();
292 test_mul(zz, xx, yy, MPFR_RNDD);
293 if (mpfr_overflow_p() || MPFR_IS_INF(zz))
294 {
295 printf("check_max failed (should NOT be an overflow)\n");
296 exit(1);
297 }
298 mpfr_set_str1 (xx, "0.93750");
299 mpfr_mul_2si(xx, xx, MPFR_EMAX_DEFAULT, MPFR_RNDN);
300 if (!(MPFR_IS_FP(xx) && MPFR_IS_FP(zz)))
301 {
302 printf("check_max failed (internal error)\n");
303 exit(1);
304 }
305 if (mpfr_cmp(xx, zz) != 0)
306 {
307 printf("check_max failed: got ");
308 mpfr_out_str(stdout, 2, 0, zz, MPFR_RNDZ);
309 printf(" instead of ");
310 mpfr_out_str(stdout, 2, 0, xx, MPFR_RNDZ);
311 printf("\n");
312 exit(1);
313 }
314
315 /* check underflow */
316 emin = mpfr_get_emin ();
317 set_emin (0);
318 mpfr_set_str_binary (xx, "0.1E0");
319 mpfr_set_str_binary (yy, "0.1E0");
320 test_mul (zz, xx, yy, MPFR_RNDN);
321 /* exact result is 0.1E-1, which should round to 0 */
322 MPFR_ASSERTN(mpfr_cmp_ui (zz, 0) == 0 && MPFR_IS_POS(zz));
323 set_emin (emin);
324
325 /* coverage test for mpfr_powerof2_raw */
326 emin = mpfr_get_emin ();
327 set_emin (0);
328 mpfr_set_prec (xx, mp_bits_per_limb + 1);
329 mpfr_set_str_binary (xx, "0.1E0");
330 mpfr_nextabove (xx);
331 mpfr_set_str_binary (yy, "0.1E0");
332 test_mul (zz, xx, yy, MPFR_RNDN);
333 /* exact result is just above 0.1E-1, which should round to minfloat */
334 MPFR_ASSERTN(mpfr_cmp (zz, yy) == 0);
335 set_emin (emin);
336
337 /* coverage test for mulders.c, case n > MUL_FFT_THRESHOLD */
338 mpfr_set_prec (xx, (MUL_FFT_THRESHOLD + 1) * GMP_NUMB_BITS);
339 mpfr_set_prec (yy, (MUL_FFT_THRESHOLD + 1) * GMP_NUMB_BITS);
340 mpfr_set_prec (zz, (MUL_FFT_THRESHOLD + 1) * GMP_NUMB_BITS);
341 mpfr_set_ui (xx, 1, MPFR_RNDN);
342 mpfr_nextbelow (xx);
343 mpfr_set_ui (yy, 1, MPFR_RNDN);
344 mpfr_nextabove (yy);
345 /* xx = 1 - 2^(-p), yy = 1 + 2^(1-p), where p = PREC(x),
346 thus xx * yy should be rounded to 1 */
347 inex = mpfr_mul (zz, xx, yy, MPFR_RNDN);
348 MPFR_ASSERTN(inex < 0);
349 MPFR_ASSERTN(mpfr_cmp_ui (zz, 1) == 0);
350
351 mpfr_clear(xx);
352 mpfr_clear(yy);
353 mpfr_clear(zz);
354 }
355
356 static void
check_min(void)357 check_min(void)
358 {
359 mpfr_t xx, yy, zz;
360
361 mpfr_init2(xx, 4);
362 mpfr_init2(yy, 4);
363 mpfr_init2(zz, 3);
364 mpfr_set_str1(xx, "0.9375");
365 mpfr_mul_2si(xx, xx, MPFR_EMIN_DEFAULT/2, MPFR_RNDN);
366 mpfr_set_str1(yy, "0.9375");
367 mpfr_mul_2si(yy, yy, MPFR_EMIN_DEFAULT - MPFR_EMIN_DEFAULT/2 - 1, MPFR_RNDN);
368 test_mul(zz, xx, yy, MPFR_RNDD);
369 if (MPFR_NOTZERO (zz))
370 {
371 printf("check_min failed: got ");
372 mpfr_out_str(stdout, 2, 0, zz, MPFR_RNDZ);
373 printf(" instead of 0\n");
374 exit(1);
375 }
376
377 test_mul(zz, xx, yy, MPFR_RNDU);
378 mpfr_set_str1 (xx, "0.5");
379 mpfr_mul_2si(xx, xx, MPFR_EMIN_DEFAULT, MPFR_RNDN);
380 if (mpfr_sgn(xx) <= 0)
381 {
382 printf("check_min failed (internal error)\n");
383 exit(1);
384 }
385 if (mpfr_cmp(xx, zz) != 0)
386 {
387 printf("check_min failed: got ");
388 mpfr_out_str(stdout, 2, 0, zz, MPFR_RNDZ);
389 printf(" instead of ");
390 mpfr_out_str(stdout, 2, 0, xx, MPFR_RNDZ);
391 printf("\n");
392 exit(1);
393 }
394
395 mpfr_clear(xx);
396 mpfr_clear(yy);
397 mpfr_clear(zz);
398 }
399
400 static void
check_nans(void)401 check_nans (void)
402 {
403 mpfr_t p, x, y;
404
405 mpfr_init2 (x, 123L);
406 mpfr_init2 (y, 123L);
407 mpfr_init2 (p, 123L);
408
409 /* nan * 0 == nan */
410 mpfr_set_nan (x);
411 mpfr_set_ui (y, 0L, MPFR_RNDN);
412 test_mul (p, x, y, MPFR_RNDN);
413 MPFR_ASSERTN (mpfr_nan_p (p));
414
415 /* 1 * nan == nan */
416 mpfr_set_ui (x, 1L, MPFR_RNDN);
417 mpfr_set_nan (y);
418 test_mul (p, x, y, MPFR_RNDN);
419 MPFR_ASSERTN (mpfr_nan_p (p));
420
421 /* 0 * +inf == nan */
422 mpfr_set_ui (x, 0L, MPFR_RNDN);
423 mpfr_set_nan (y);
424 test_mul (p, x, y, MPFR_RNDN);
425 MPFR_ASSERTN (mpfr_nan_p (p));
426
427 /* +1 * +inf == +inf */
428 mpfr_set_ui (x, 1L, MPFR_RNDN);
429 mpfr_set_inf (y, 1);
430 test_mul (p, x, y, MPFR_RNDN);
431 MPFR_ASSERTN (mpfr_inf_p (p));
432 MPFR_ASSERTN (mpfr_sgn (p) > 0);
433
434 /* -1 * +inf == -inf */
435 mpfr_set_si (x, -1L, MPFR_RNDN);
436 mpfr_set_inf (y, 1);
437 test_mul (p, x, y, MPFR_RNDN);
438 MPFR_ASSERTN (mpfr_inf_p (p));
439 MPFR_ASSERTN (mpfr_sgn (p) < 0);
440
441 mpfr_clear (x);
442 mpfr_clear (y);
443 mpfr_clear (p);
444 }
445
446 #define BUFSIZE 1552
447
448 static void
get_string(char * s,FILE * fp)449 get_string (char *s, FILE *fp)
450 {
451 int c, n = BUFSIZE;
452
453 while ((c = getc (fp)) != '\n')
454 {
455 if (c == EOF)
456 {
457 printf ("Error in get_string: end of file\n");
458 exit (1);
459 }
460 *(unsigned char *)s++ = c;
461 if (--n == 0)
462 {
463 printf ("Error in get_string: buffer is too small\n");
464 exit (1);
465 }
466 }
467 *s = '\0';
468 }
469
470 static void
check_regression(void)471 check_regression (void)
472 {
473 mpfr_t x, y, z;
474 int i;
475 FILE *fp;
476 char s[BUFSIZE];
477
478 mpfr_inits2 (6177, x, y, z, (mpfr_ptr) 0);
479 /* we read long strings from a file since ISO C90 does not support strings of
480 length > 509 */
481 fp = src_fopen ("tmul.dat", "r");
482 if (fp == NULL)
483 {
484 fprintf (stderr, "Error, cannot open tmul.dat in srcdir\n");
485 exit (1);
486 }
487 get_string (s, fp);
488 mpfr_set_str (y, s, 16, MPFR_RNDN);
489 get_string (s, fp);
490 mpfr_set_str (z, s, 16, MPFR_RNDN);
491 i = mpfr_mul (x, y, z, MPFR_RNDN);
492 get_string (s, fp);
493 if (mpfr_cmp_str (x, s, 16, MPFR_RNDN) != 0 || i != -1)
494 {
495 printf ("Regression test 1 failed (i=%d, expected -1)\nx=", i);
496 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN); putchar ('\n');
497 exit (1);
498 }
499 fclose (fp);
500
501 mpfr_set_prec (x, 606);
502 mpfr_set_prec (y, 606);
503 mpfr_set_prec (z, 606);
504
505 mpfr_set_str (y, "-f.ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff92daefc3f8052ca9f58736564d9e93e62d324@-1", 16, MPFR_RNDN);
506 mpfr_set_str (z, "-f.ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff92daefc3f8052ca9f58736564d9e93e62d324@-1", 16, MPFR_RNDN);
507 i = mpfr_mul (x, y, z, MPFR_RNDU);
508 mpfr_set_str (y, "f.ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff25b5df87f00a5953eb0e6cac9b3d27cc5a64c@-1", 16, MPFR_RNDN);
509 if (mpfr_cmp (x, y) || i <= 0)
510 {
511 printf ("Regression test (2) failed! (i=%d - Expected 1)\n", i);
512 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN); putchar ('\n');
513 exit (1);
514 }
515
516 mpfr_set_prec (x, 184);
517 mpfr_set_prec (y, 92);
518 mpfr_set_prec (z, 1023);
519
520 mpfr_set_str (y, "6.9b8c8498882770d8038c3b0@-1", 16, MPFR_RNDN);
521 mpfr_set_str (z, "7.44e24b986e7fb296f1e936ce749fec3504cbf0d5ba769466b1c9f1578115efd5d29b4c79271191a920a99280c714d3a657ad6e3afbab77ffce9d697e9bb9110e26d676069afcea8b69f1d1541f2365042d80a97c21dcccd8ace4f1bb58b49922003e738e6f37bb82ef653cb2e87f763974e6ae50ae54e7724c38b80653e3289@255", 16, MPFR_RNDN);
522 i = mpfr_mul (x, y, z, MPFR_RNDU);
523 mpfr_set_prec (y, 184);
524 mpfr_set_str (y, "3.0080038f2ac5054e3e71ccbb95f76aaab2221715025a28@255",
525 16, MPFR_RNDN);
526 if (mpfr_cmp (x, y) || i <= 0)
527 {
528 printf ("Regression test (4) failed! (i=%d - expected 1)\n", i);
529 printf ("Ref: 3.0080038f2ac5054e3e71ccbb95f76aaab2221715025a28@255\n"
530 "Got: ");
531 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN);
532 printf ("\n");
533 exit (1);
534 }
535
536 mpfr_set_prec (x, 908);
537 mpfr_set_prec (y, 908);
538 mpfr_set_prec (z, 908);
539 mpfr_set_str (y, "-f.fffffffffffffffffffffffffffffffffffffffffffffffffffffff"
540 "fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff"
541 "ffffffffffffffffffffffffffffffffffffffffffffffffffffff99be91f83ec6f0ed28a3d42"
542 "e6e9a327230345ea6@-1", 16, MPFR_RNDN);
543 mpfr_set_str (z, "-f.fffffffffffffffffffffffffffffffffffffffffffffffffffffff"
544 "fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff"
545 "ffffffffffffffffffffffffffffffffffffffffffffffffffffff99be91f83ec6f0ed28a3d42"
546 "e6e9a327230345ea6@-1", 16, MPFR_RNDN);
547 i = mpfr_mul (x, y, z, MPFR_RNDU);
548 mpfr_set_str (y, "f.ffffffffffffffffffffffffffffffffffffffffffffffffffffffff"
549 "fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff"
550 "fffffffffffffffffffffffffffffffffffffffffffffffffffff337d23f07d8de1da5147a85c"
551 "dd3464e46068bd4d@-1", 16, MPFR_RNDN);
552 if (mpfr_cmp (x, y) || i <= 0)
553 {
554 printf ("Regression test (5) failed! (i=%d - expected 1)\n", i);
555 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN);
556 printf ("\n");
557 exit (1);
558 }
559
560
561 mpfr_set_prec (x, 50);
562 mpfr_set_prec (y, 40);
563 mpfr_set_prec (z, 53);
564 mpfr_set_str (y, "4.1ffffffff8", 16, MPFR_RNDN);
565 mpfr_set_str (z, "4.2000000ffe0000@-4", 16, MPFR_RNDN);
566 i = mpfr_mul (x, y, z, MPFR_RNDN);
567 if (mpfr_cmp_str (x, "1.104000041d6c0@-3", 16, MPFR_RNDN) != 0
568 || i <= 0)
569 {
570 printf ("Regression test (6) failed! (i=%d - expected 1)\nx=", i);
571 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN);
572 printf ("\nMore prec=");
573 mpfr_set_prec (x, 93);
574 mpfr_mul (x, y, z, MPFR_RNDN);
575 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN);
576 printf ("\n");
577 exit (1);
578 }
579
580 mpfr_set_prec (x, 439);
581 mpfr_set_prec (y, 393);
582 mpfr_set_str (y, "-1.921fb54442d18469898cc51701b839a252049c1114cf98e804177d"
583 "4c76273644a29410f31c6809bbdf2a33679a748636600",
584 16, MPFR_RNDN);
585 /* the following call to mpfr_mul with identical arguments is intentional */
586 i = mpfr_mul (x, y, y, MPFR_RNDU);
587 if (mpfr_cmp_str (x, "2.77a79937c8bbcb495b89b36602306b1c2159a8ff834288a19a08"
588 "84094f1cda3dc426da61174c4544a173de83c2500f8bfea2e0569e3698",
589 16, MPFR_RNDN) != 0
590 || i <= 0)
591 {
592 printf ("Regression test (7) failed! (i=%d - expected 1)\nx=", i);
593 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN);
594 printf ("\n");
595 exit (1);
596 }
597
598 mpfr_set_prec (x, 1023);
599 mpfr_set_prec (y, 1023);
600 mpfr_set_prec (z, 511);
601 mpfr_set_ui (x, 17, MPFR_RNDN);
602 mpfr_set_ui (y, 42, MPFR_RNDN);
603 i = mpfr_mul (z, x, y, MPFR_RNDN);
604 if (mpfr_cmp_ui (z, 17*42) != 0 || i != 0)
605 {
606 printf ("Regression test (8) failed! (i=%d - expected 0)\nz=", i);
607 mpfr_out_str (stdout, 16, 0, z, MPFR_RNDN);
608 printf ("\n");
609 exit (1);
610 }
611
612 mpfr_clears (x, y, z, (mpfr_ptr) 0);
613 }
614
615 #define TEST_FUNCTION test_mul
616 #define TWO_ARGS
617 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS)
618 #include "tgeneric.c"
619
620 /* multiplies x by 53-bit approximation of Pi */
621 static int
mpfr_mulpi(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t r)622 mpfr_mulpi (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r)
623 {
624 mpfr_t z;
625 int inex;
626
627 mpfr_init2 (z, 53);
628 mpfr_set_str_binary (z, "11.001001000011111101101010100010001000010110100011");
629 inex = mpfr_mul (y, x, z, r);
630 mpfr_clear (z);
631 return inex;
632 }
633
634 static void
valgrind20110503(void)635 valgrind20110503 (void)
636 {
637 mpfr_t a, b, c;
638
639 mpfr_init2 (a, 2);
640 mpfr_init2 (b, 2005);
641 mpfr_init2 (c, 2);
642
643 mpfr_set_ui (b, 5, MPFR_RNDN);
644 mpfr_nextabove (b);
645 mpfr_set_ui (c, 1, MPFR_RNDN);
646 mpfr_mul (a, b, c, MPFR_RNDZ);
647 /* After the call to mpfr_mulhigh_n, valgrind complains:
648 Conditional jump or move depends on uninitialised value(s) */
649
650 mpfr_clears (a, b, c, (mpfr_ptr) 0);
651 }
652
653 static void
testall_rndf(mpfr_prec_t pmax)654 testall_rndf (mpfr_prec_t pmax)
655 {
656 mpfr_t a, b, c, d;
657 mpfr_prec_t pa, pb, pc;
658
659 for (pa = MPFR_PREC_MIN; pa <= pmax; pa++)
660 {
661 mpfr_init2 (a, pa);
662 mpfr_init2 (d, pa);
663 for (pb = MPFR_PREC_MIN; pb <= pmax; pb++)
664 {
665 mpfr_init2 (b, pb);
666 mpfr_set_ui (b, 1, MPFR_RNDN);
667 while (mpfr_cmp_ui (b, 2) < 0)
668 {
669 for (pc = MPFR_PREC_MIN; pc <= pmax; pc++)
670 {
671 mpfr_init2 (c, pc);
672 mpfr_set_ui (c, 1, MPFR_RNDN);
673 while (mpfr_cmp_ui (c, 2) < 0)
674 {
675 mpfr_mul (a, b, c, MPFR_RNDF);
676 mpfr_mul (d, b, c, MPFR_RNDD);
677 if (!mpfr_equal_p (a, d))
678 {
679 mpfr_mul (d, b, c, MPFR_RNDU);
680 if (!mpfr_equal_p (a, d))
681 {
682 printf ("Error: mpfr_mul(a,b,c,RNDF) does not "
683 "match RNDD/RNDU\n");
684 printf ("b="); mpfr_dump (b);
685 printf ("c="); mpfr_dump (c);
686 printf ("a="); mpfr_dump (a);
687 exit (1);
688 }
689 }
690 mpfr_nextabove (c);
691 }
692 mpfr_clear (c);
693 }
694 mpfr_nextabove (b);
695 }
696 mpfr_clear (b);
697 }
698 mpfr_clear (a);
699 mpfr_clear (d);
700 }
701 }
702
703 /* Check underflow flag corresponds to *after* rounding.
704 *
705 * More precisely, we want to test mpfr_mul on inputs b and c such that
706 * EXP(b*c) < emin but EXP(round(b*c, p, rnd)) = emin. Thus an underflow
707 * must not be generated.
708 */
709 static void
test_underflow(mpfr_prec_t pmax)710 test_underflow (mpfr_prec_t pmax)
711 {
712 mpfr_exp_t emin;
713 mpfr_prec_t p;
714 mpfr_t a0, a, b, c;
715 int inex;
716
717 mpfr_init2 (a0, MPFR_PREC_MIN);
718 emin = mpfr_get_emin ();
719 mpfr_setmin (a0, emin); /* 0.5 * 2^emin */
720
721 /* for RNDN, we want b*c < 0.5 * 2^emin but RNDN(b*c, p) = 0.5 * 2^emin,
722 thus b*c >= (0.5 - 1/4 * ulp_p(0.5)) * 2^emin */
723 for (p = MPFR_PREC_MIN; p <= pmax; p++)
724 {
725 mpfr_init2 (a, p + 1);
726 mpfr_init2 (b, p + 10);
727 mpfr_init2 (c, p + 10);
728 do mpfr_urandomb (b, RANDS); while (MPFR_IS_ZERO (b));
729 inex = mpfr_set_ui_2exp (a, 1, -1, MPFR_RNDZ); /* a = 0.5 */
730 MPFR_ASSERTN (inex == 0);
731 mpfr_nextbelow (a); /* 0.5 - 1/2*ulp_{p+1}(0.5) = 0.5 - 1/4*ulp_p(0.5) */
732 inex = mpfr_div (c, a, b, MPFR_RNDU);
733 /* 0.5 - 1/4 * ulp_p(0.5) = a <= b*c < 0.5 */
734 mpfr_mul_2si (b, b, emin / 2, MPFR_RNDZ);
735 mpfr_mul_2si (c, c, (emin - 1) / 2, MPFR_RNDZ);
736 /* now (0.5 - 1/4 * ulp_p(0.5)) * 2^emin <= b*c < 0.5 * 2^emin,
737 thus b*c should be rounded to 0.5 * 2^emin */
738 mpfr_set_prec (a, p);
739 mpfr_clear_underflow ();
740 mpfr_mul (a, b, c, MPFR_RNDN);
741 if (mpfr_underflow_p () || ! mpfr_equal_p (a, a0))
742 {
743 printf ("Error, b*c incorrect or underflow flag incorrectly set"
744 " for emin=%" MPFR_EXP_FSPEC "d, rnd=%s\n",
745 (mpfr_eexp_t) emin, mpfr_print_rnd_mode (MPFR_RNDN));
746 printf ("b="); mpfr_dump (b);
747 printf ("c="); mpfr_dump (c);
748 printf ("a="); mpfr_dump (a);
749 mpfr_set_prec (a, mpfr_get_prec (b) + mpfr_get_prec (c));
750 mpfr_mul_2ui (b, b, 1, MPFR_RNDN);
751 inex = mpfr_mul (a, b, c, MPFR_RNDN);
752 MPFR_ASSERTN (inex == 0);
753 printf ("Exact 2*a="); mpfr_dump (a);
754 exit (1);
755 }
756 mpfr_clear (a);
757 mpfr_clear (b);
758 mpfr_clear (c);
759 }
760
761 /* for RNDU, we want b*c < 0.5*2^emin but RNDU(b*c, p) = 0.5*2^emin thus
762 b*c > (0.5 - 1/2 * ulp_p(0.5)) * 2^emin */
763 for (p = MPFR_PREC_MIN; p <= pmax; p++)
764 {
765 mpfr_init2 (a, p);
766 mpfr_init2 (b, p + 10);
767 mpfr_init2 (c, p + 10);
768 do mpfr_urandomb (b, RANDS); while (MPFR_IS_ZERO (b));
769 inex = mpfr_set_ui_2exp (a, 1, -1, MPFR_RNDZ); /* a = 0.5 */
770 MPFR_ASSERTN (inex == 0);
771 mpfr_nextbelow (a); /* 0.5 - 1/2 * ulp_p(0.5) */
772 inex = mpfr_div (c, a, b, MPFR_RNDU);
773 /* 0.5 - 1/2 * ulp_p(0.5) <= b*c < 0.5 */
774 mpfr_mul_2si (b, b, emin / 2, MPFR_RNDZ);
775 mpfr_mul_2si (c, c, (emin - 1) / 2, MPFR_RNDZ);
776 if (inex == 0)
777 mpfr_nextabove (c); /* ensures b*c > (0.5 - 1/2 * ulp_p(0.5)) * 2^emin.
778 Warning: for p=1, 0.5 - 1/2 * ulp_p(0.5)
779 = 0.25, thus b*c > 2^(emin-2), which should
780 also be rounded up with p=1 to 0.5 * 2^emin
781 with an unbounded exponent range. */
782 mpfr_clear_underflow ();
783 mpfr_mul (a, b, c, MPFR_RNDU);
784 if (mpfr_underflow_p () || ! mpfr_equal_p (a, a0))
785 {
786 printf ("Error, b*c incorrect or underflow flag incorrectly set"
787 " for emin=%" MPFR_EXP_FSPEC "d, rnd=%s\n",
788 (mpfr_eexp_t) emin, mpfr_print_rnd_mode (MPFR_RNDU));
789 printf ("b="); mpfr_dump (b);
790 printf ("c="); mpfr_dump (c);
791 printf ("a="); mpfr_dump (a);
792 mpfr_set_prec (a, mpfr_get_prec (b) + mpfr_get_prec (c));
793 mpfr_mul_2ui (b, b, 1, MPFR_RNDN);
794 inex = mpfr_mul (a, b, c, MPFR_RNDN);
795 MPFR_ASSERTN (inex == 0);
796 printf ("Exact 2*a="); mpfr_dump (a);
797 exit (1);
798 }
799 mpfr_clear (a);
800 mpfr_clear (b);
801 mpfr_clear (c);
802 }
803
804 mpfr_clear (a0);
805 }
806
807 /* checks special case where no underflow should occur */
808 static void
bug20161209(void)809 bug20161209 (void)
810 {
811 mpfr_exp_t emin;
812 mpfr_t x, y, z;
813
814 emin = mpfr_get_emin ();
815 set_emin (-1);
816
817 /* test for mpfr_mul_1 for 64-bit limb, mpfr_mul_2 for 32-bit limb */
818 mpfr_init2 (x, 53);
819 mpfr_init2 (y, 53);
820 mpfr_init2 (z, 53);
821 mpfr_set_str_binary (x, "0.101000001E-1"); /* x = 321/2^10 */
822 mpfr_set_str_binary (y, "0.110011000010100101111000011011000111011000001E-1");
823 /* y = 28059810762433/2^46 */
824 /* x * y = (2^53+1)/2^56 = 0.001...000[1]000..., and should round to 0.25 */
825 mpfr_mul (z, x, y, MPFR_RNDN);
826 MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
827
828 /* test for mpfr_mul_2 for 64-bit limb */
829 mpfr_set_prec (x, 65);
830 mpfr_set_prec (y, 65);
831 mpfr_set_prec (z, 65);
832 mpfr_set_str_binary (x, "0.101101000010010110100001E-1"); /* 11806113/2^25 */
833 mpfr_set_str_binary (y, "0.101101011110010101011010001111111001100001E-1");
834 /* y = 3124947910241/2^43 */
835 /* x * y = (2^65+1)/2^68 = 0.001...000[1]000..., and should round to 0.25 */
836 mpfr_mul (z, x, y, MPFR_RNDN);
837 MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
838
839 /* test for the generic code */
840 mpfr_set_prec (x, 54);
841 mpfr_set_prec (y, 55);
842 mpfr_set_prec (z, 53);
843 mpfr_set_str_binary (x, "0.101000001E-1");
844 mpfr_set_str_binary (y, "0.110011000010100101111000011011000111011000001E-1");
845 mpfr_mul (z, x, y, MPFR_RNDN);
846 MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
847
848 /* another test for the generic code */
849 mpfr_set_prec (x, 66);
850 mpfr_set_prec (y, 67);
851 mpfr_set_prec (z, 65);
852 mpfr_set_str_binary (x, "0.101101000010010110100001E-1");
853 mpfr_set_str_binary (y, "0.101101011110010101011010001111111001100001E-1");
854 mpfr_mul (z, x, y, MPFR_RNDN);
855 MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
856
857 mpfr_clear (x);
858 mpfr_clear (y);
859 mpfr_clear (z);
860 set_emin (emin);
861 }
862
863 /* test for case in mpfr_mul_1() where:
864 ax = __gmpfr_emin - 1
865 ap[0] == ~mask
866 rnd_mode = MPFR_RNDZ.
867 Whatever the values of rb and sb, we should round to zero (underflow). */
868 static void
bug20161209a(void)869 bug20161209a (void)
870 {
871 mpfr_exp_t emin;
872 mpfr_t x, y, z;
873
874 emin = mpfr_get_emin ();
875 set_emin (-1);
876
877 mpfr_init2 (x, 53);
878 mpfr_init2 (y, 53);
879 mpfr_init2 (z, 53);
880
881 /* case rb = sb = 0 */
882 mpfr_set_str_binary (x, "0.11010010100110000110110011111E-1");
883 mpfr_set_str_binary (y, "0.1001101110011000110100001");
884 /* x = 441650591/2^30, y = 20394401/2^25 */
885 /* x * y = (2^53-1)/2^55 = 0.00111...111[0]000..., and should round to 0 */
886 mpfr_mul (z, x, y, MPFR_RNDZ);
887 MPFR_ASSERTN(mpfr_zero_p (z));
888
889 /* case rb = 1, sb = 0 */
890 mpfr_set_str_binary (x, "0.111111111000000000000000000111111111E-1");
891 mpfr_set_str_binary (y, "0.1000000001000000001");
892 /* x = 68585259519/2^37, y = 262657/2^19 */
893 /* x * y = (2^54-1)/2^56 = 0.00111...111[1]000..., and should round to 0 */
894 mpfr_mul (z, x, y, MPFR_RNDZ);
895 MPFR_ASSERTN(mpfr_zero_p (z));
896
897 /* case rb = 0, sb = 1 */
898 mpfr_set_str_binary (x, "0.110010011001011110001100100001000001E-1");
899 mpfr_set_str_binary (y, "0.10100010100010111101");
900 /* x = 541144371852^37, y = 665789/2^20 */
901 /* x * y = (2^55-3)/2^57 = 0.00111...111[0]100..., and should round to 0 */
902 mpfr_mul (z, x, y, MPFR_RNDZ);
903 MPFR_ASSERTN(mpfr_zero_p (z));
904
905 /* case rb = sb = 1 */
906 mpfr_set_str_binary (x, "0.10100110001001001010001111110010100111E-1");
907 mpfr_set_str_binary (y, "0.110001010011101001");
908 /* x = 178394823847/2^39, y = 201961/2^18 */
909 /* x * y = (2^55-1)/2^57 = 0.00111...111[1]100..., and should round to 0 */
910 mpfr_mul (z, x, y, MPFR_RNDZ);
911 MPFR_ASSERTN(mpfr_zero_p (z));
912
913 /* similar test for mpfr_mul_2 (we only check rb = sb = 1 here) */
914 mpfr_set_prec (x, 65);
915 mpfr_set_prec (y, 65);
916 mpfr_set_prec (z, 65);
917 /* 2^67-1 = 193707721 * 761838257287 */
918 mpfr_set_str_binary (x, "0.1011100010111011111011001001E-1");
919 mpfr_set_str_binary (y, "0.1011000101100001000110010100010010000111");
920 mpfr_mul (z, x, y, MPFR_RNDZ);
921 MPFR_ASSERTN(mpfr_zero_p (z));
922
923 mpfr_clear (x);
924 mpfr_clear (y);
925 mpfr_clear (z);
926 set_emin (emin);
927 }
928
929 /* bug for RNDF */
930 static void
bug20170602(void)931 bug20170602 (void)
932 {
933 mpfr_t x, u, y, yd, yu;
934
935 mpfr_init2 (x, 493);
936 mpfr_init2 (u, 493);
937 mpfr_init2 (y, 503);
938 mpfr_init2 (yd, 503);
939 mpfr_init2 (yu, 503);
940 mpfr_set_str_binary (x, "0.1111100000000000001111111110000000001111111111111000000000000000000011111111111111111111111000000000000000000001111111111111111111111111111111111111111000000000011111111111111111111000000000011111111111111000000000000001110000000000000000000000000000000000000000011111111111110011111111111100000000000000011111111111111111110000000011111111111111111110011111111111110000000000001111111111111111000000000000000000000000000000000000111111111111111111111111111111011111111111111111111111111111100E44");
941 mpfr_set_str_binary (u, "0.1110000000000000001111111111111111111111111111111111111000000000111111111111111111111111111111000000000000000000001111111000000000000000011111111111111111111111111111111111111111111111111111111000000000000000011111111111111000000011111111111111111110000000000000001111111111111111111111111111111111111110000000000001111111111111111111111111111111111111000000000000000000000000000000000001111111111111000000000000000000001111111111100000000000000011111111111111111111111111111111111111111111111E35");
942 mpfr_mul (y, x, u, MPFR_RNDF);
943 mpfr_mul (yd, x, u, MPFR_RNDD);
944 mpfr_mul (yu, x, u, MPFR_RNDU);
945 if (mpfr_cmp (y, yd) != 0 && mpfr_cmp (y, yu) != 0)
946 {
947 printf ("RNDF is neither RNDD nor RNDU\n");
948 printf ("x"); mpfr_dump (x);
949 printf ("u"); mpfr_dump (u);
950 printf ("y(RNDF)="); mpfr_dump (y);
951 printf ("y(RNDD)="); mpfr_dump (yd);
952 printf ("y(RNDU)="); mpfr_dump (yu);
953 exit (1);
954 }
955 mpfr_clear (x);
956 mpfr_clear (u);
957 mpfr_clear (y);
958 mpfr_clear (yd);
959 mpfr_clear (yu);
960 }
961
962 /* Test for 1 to 3 limbs. */
963 static void
small_prec(void)964 small_prec (void)
965 {
966 mpfr_exp_t emin, emax;
967 mpfr_t x, y, z1, z2, zz;
968 int xq, yq, zq;
969 mpfr_rnd_t rnd;
970 mpfr_flags_t flags1, flags2;
971 int inex1, inex2;
972 int i, j, r, s, ediff;
973
974 emin = mpfr_get_emin ();
975 emax = mpfr_get_emax ();
976
977 /* The mpfr_mul implementation doesn't extend the exponent range,
978 so that it is OK to reduce it here for the test to make sure that
979 mpfr_mul_2si can be used. */
980 set_emin (-1000);
981 set_emax (1000);
982
983 mpfr_inits2 (3 * GMP_NUMB_BITS, x, y, z1, z2, (mpfr_ptr) 0);
984 mpfr_init2 (zz, 6 * GMP_NUMB_BITS);
985 for (i = 0; i < 3; i++)
986 for (j = 0; j < 10000; j++)
987 {
988 xq = i * GMP_NUMB_BITS + 1 + randlimb () % GMP_NUMB_BITS;
989 mpfr_set_prec (x, xq);
990 yq = i * GMP_NUMB_BITS + 1 + randlimb () % GMP_NUMB_BITS;
991 mpfr_set_prec (y, yq);
992 zq = i * GMP_NUMB_BITS + 1 + randlimb () % (GMP_NUMB_BITS-1);
993 mpfr_set_prec (z1, zq);
994 mpfr_set_prec (z2, zq);
995 s = r = randlimb () & 0x7f;
996 do mpfr_urandomb (x, RANDS); while (mpfr_zero_p (x));
997 if (s & 1)
998 mpfr_neg (x, x, MPFR_RNDN);
999 s >>= 1;
1000 if (s & 1)
1001 {
1002 do mpfr_urandomb (y, RANDS); while (mpfr_zero_p (y));
1003 }
1004 else
1005 {
1006 mpfr_ui_div (y, 1, x, MPFR_RNDN);
1007 mpfr_set_exp (y, 0);
1008 }
1009 s >>= 1;
1010 if (s & 1)
1011 mpfr_neg (y, y, MPFR_RNDN);
1012 s >>= 1;
1013 rnd = RND_RAND_NO_RNDF ();
1014 inex1 = mpfr_mul (zz, x, y, MPFR_RNDN);
1015 MPFR_ASSERTN (inex1 == 0);
1016 if (s == 0)
1017 {
1018 ediff = __gmpfr_emin - MPFR_EXP (x);
1019 mpfr_set_exp (x, __gmpfr_emin);
1020 }
1021 else if (s == 1)
1022 {
1023 ediff = __gmpfr_emax - MPFR_EXP (x) + 1;
1024 mpfr_set_exp (x, __gmpfr_emax);
1025 mpfr_mul_2ui (y, y, 1, MPFR_RNDN);
1026 }
1027 else
1028 ediff = 0;
1029 mpfr_clear_flags ();
1030 inex1 = mpfr_mul_2si (z1, zz, ediff, rnd);
1031 flags1 = __gmpfr_flags;
1032 mpfr_clear_flags ();
1033 inex2 = mpfr_mul (z2, x, y, rnd);
1034 flags2 = __gmpfr_flags;
1035 if (!(mpfr_equal_p (z1, z2) &&
1036 SAME_SIGN (inex1, inex2)
1037 && flags1 == flags2
1038 ))
1039 {
1040 printf ("Error in small_prec on i = %d, j = %d\n", i, j);
1041 printf ("r = 0x%x, xq = %d, yq = %d, zq = %d, rnd = %s\n",
1042 r, xq, yq, zq, mpfr_print_rnd_mode (rnd));
1043 printf ("x = ");
1044 mpfr_dump (x);
1045 printf ("y = ");
1046 mpfr_dump (y);
1047 printf ("ediff = %d\n", ediff);
1048 printf ("zz = ");
1049 mpfr_dump (zz);
1050 printf ("Expected ");
1051 mpfr_dump (z1);
1052 printf ("with inex = %d and flags =", inex1);
1053 flags_out (flags1);
1054 printf ("Got ");
1055 mpfr_dump (z2);
1056 printf ("with inex = %d and flags =", inex2);
1057 flags_out (flags2);
1058 exit (1);
1059 }
1060 }
1061 mpfr_clears (x, y, z1, z2, zz, (mpfr_ptr) 0);
1062
1063 set_emin (emin);
1064 set_emax (emax);
1065 }
1066
1067 /* check ax < __gmpfr_emin with rnd_mode == MPFR_RNDN, rb = 0 and sb <> 0 */
1068 static void
test_underflow2(void)1069 test_underflow2 (void)
1070 {
1071 mpfr_t x, y, z;
1072 mpfr_exp_t emin;
1073
1074 emin = mpfr_get_emin ();
1075 set_emin (0);
1076
1077 mpfr_init2 (x, 24);
1078 mpfr_init2 (y, 24);
1079 mpfr_init2 (z, 24);
1080
1081 mpfr_set_ui_2exp (x, 12913, -14, MPFR_RNDN);
1082 mpfr_set_ui_2exp (y, 5197, -13, MPFR_RNDN);
1083 mpfr_clear_underflow ();
1084 /* x*y = 0.0111111111111111111111111[01] thus underflow */
1085 mpfr_mul (z, y, x, MPFR_RNDN);
1086 MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -1) == 0);
1087 MPFR_ASSERTN(mpfr_underflow_p ());
1088
1089 mpfr_set_prec (x, 65);
1090 mpfr_set_prec (y, 65);
1091 mpfr_set_prec (z, 65);
1092 mpfr_set_str_binary (x, "0.10011101000110111011111100101001111");
1093 mpfr_set_str_binary (y, "0.110100001001000111000011011110011");
1094 mpfr_clear_underflow ();
1095 /* x*y = 0.011111111111111111111111111111111111111111111111111111111111111111[01] thus underflow */
1096 mpfr_mul (z, y, x, MPFR_RNDN);
1097 MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -1) == 0);
1098 MPFR_ASSERTN(mpfr_underflow_p ());
1099
1100 mpfr_clear (y);
1101 mpfr_clear (x);
1102 mpfr_clear (z);
1103
1104 set_emin (emin);
1105 }
1106
1107 static void
coverage(mpfr_prec_t pmax)1108 coverage (mpfr_prec_t pmax)
1109 {
1110 mpfr_t a, b, c;
1111 mpfr_prec_t p;
1112 int inex;
1113
1114 for (p = MPFR_PREC_MIN; p <= pmax; p++)
1115 {
1116 mpfr_init2 (a, p);
1117 mpfr_init2 (b, p);
1118 mpfr_init2 (c, p);
1119
1120 /* exercise case b*c = 2^(emin-2), which is just in the middle
1121 between 0 and the smallest positive number 0.5*2^emin */
1122 mpfr_set_ui_2exp (b, 1, mpfr_get_emin (), MPFR_RNDN);
1123 mpfr_set_ui_2exp (c, 1, -2, MPFR_RNDN);
1124 mpfr_clear_flags ();
1125 inex = mpfr_mul (a, b, c, MPFR_RNDN);
1126 MPFR_ASSERTN(inex < 0);
1127 MPFR_ASSERTN(mpfr_zero_p (a) && mpfr_signbit (a) == 0);
1128 MPFR_ASSERTN(mpfr_underflow_p ());
1129
1130 if (p == 1)
1131 goto end_of_loop;
1132
1133 /* case b*c > 2^(emin-2): b = (1-2^(-p))*2^emin,
1134 c = 0.25*(1+2^(1-p)), thus b*c = (1+2^(-p)-2^(1-2p))*2^(emin-2)
1135 should be rounded to 2^(emin-1) for RNDN */
1136 mpfr_nextbelow (b);
1137 mpfr_nextabove (c);
1138 mpfr_clear_flags ();
1139 inex = mpfr_mul (a, b, c, MPFR_RNDN);
1140 MPFR_ASSERTN(inex > 0);
1141 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, mpfr_get_emin () - 1) == 0);
1142 MPFR_ASSERTN(mpfr_underflow_p ());
1143
1144 /* b = (1-2^(1-p))*2^emin, c = 0.25*(1+2^(1-p)),
1145 thus b*c = (1-2^(2-2p))*2^(emin-2) should be rounded to 0 */
1146 mpfr_nextbelow (b);
1147 mpfr_clear_flags ();
1148 inex = mpfr_mul (a, b, c, MPFR_RNDN);
1149 MPFR_ASSERTN(inex < 0);
1150 MPFR_ASSERTN(mpfr_zero_p (a) && mpfr_signbit (a) == 0);
1151 MPFR_ASSERTN(mpfr_underflow_p ());
1152
1153 /* special case where b*c is in [nextbelow(0.5*2^emin),0.5*2^emin[ */
1154 if ((p % 2) == 0)
1155 {
1156 /* the middle of the interval [nextbelow(0.5*2^emin),0.5*2^emin[
1157 is (1-2^(-p-1))*2^(emin-1)
1158 = (1-2^(-p/2))*(1+2^(-p/2))*2^(emin-1) */
1159 mpfr_set_si_2exp (b, -1, -p/2, MPFR_RNDN);
1160 mpfr_add_ui (b, b, 1, MPFR_RNDN);
1161 mpfr_set_si_2exp (c, 1, -p/2, MPFR_RNDN);
1162 mpfr_add_ui (c, c, 1, MPFR_RNDN);
1163 MPFR_ASSERTN(mpfr_get_emin () < 0);
1164 mpfr_mul_2si (b, b, (mpfr_get_emin () - 1) / 2, MPFR_RNDN);
1165 mpfr_mul_2si (c, c, (mpfr_get_emin () - 2) / 2, MPFR_RNDN);
1166 mpfr_clear_flags ();
1167 inex = mpfr_mul (a, b, c, MPFR_RNDN);
1168 MPFR_ASSERTN(inex > 0);
1169 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, mpfr_get_emin () - 1) == 0);
1170 MPFR_ASSERTN(mpfr_underflow_p ());
1171 mpfr_clear_flags ();
1172 inex = mpfr_mul (a, b, c, MPFR_RNDU);
1173 MPFR_ASSERTN(inex > 0);
1174 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, mpfr_get_emin () - 1) == 0);
1175 MPFR_ASSERTN(mpfr_underflow_p ());
1176 mpfr_clear_flags ();
1177 inex = mpfr_mul (a, b, c, MPFR_RNDD);
1178 MPFR_ASSERTN(inex < 0);
1179 MPFR_ASSERTN(mpfr_zero_p (a) && mpfr_signbit (a) == 0);
1180 MPFR_ASSERTN(mpfr_underflow_p ());
1181 }
1182 else /* p is odd:
1183 b = (1-2^(-(p+1)/2))*2^...
1184 c = (1+2^(-(p+1)/2))*2^... */
1185 {
1186 mpfr_set_si_2exp (b, -1, -(p+1)/2, MPFR_RNDN);
1187 mpfr_add_ui (b, b, 1, MPFR_RNDN);
1188 mpfr_set_si_2exp (c, 1, -(p+1)/2, MPFR_RNDN);
1189 mpfr_add_ui (c, c, 1, MPFR_RNDN);
1190 MPFR_ASSERTN(mpfr_get_emin () < 0);
1191 mpfr_mul_2si (b, b, (mpfr_get_emin () - 1) / 2, MPFR_RNDN);
1192 mpfr_mul_2si (c, c, (mpfr_get_emin () - 2) / 2, MPFR_RNDN);
1193 mpfr_clear_flags ();
1194 inex = mpfr_mul (a, b, c, MPFR_RNDN);
1195 MPFR_ASSERTN(inex > 0);
1196 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, mpfr_get_emin () - 1) == 0);
1197 MPFR_ASSERTN(!mpfr_underflow_p ());
1198 mpfr_clear_flags ();
1199 inex = mpfr_mul (a, b, c, MPFR_RNDU);
1200 MPFR_ASSERTN(inex > 0);
1201 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, mpfr_get_emin () - 1) == 0);
1202 MPFR_ASSERTN(!mpfr_underflow_p ());
1203 mpfr_clear_flags ();
1204 inex = mpfr_mul (a, b, c, MPFR_RNDD);
1205 MPFR_ASSERTN(inex < 0);
1206 MPFR_ASSERTN(mpfr_zero_p (a) && mpfr_signbit (a) == 0);
1207 MPFR_ASSERTN(mpfr_underflow_p ());
1208 }
1209
1210 if (p <= 2) /* for p=2, 1+2^(-ceil((p+1)/2)) = 1 + 2^(-2) is not
1211 exactly representable */
1212 goto end_of_loop;
1213
1214 /* b = 1-2^(-ceil((p+1)/2))
1215 c = 1+2^(-ceil((p+1)/2))
1216 For p odd, b*c = 1-2^(p+1) should round to 1;
1217 for p even, b*c = 1-2^(p+2) should round to 1 too. */
1218 mpfr_set_si_2exp (b, -1, -(p+2)/2, MPFR_RNDN);
1219 mpfr_add_ui (b, b, 1, MPFR_RNDN);
1220 mpfr_set_si_2exp (c, 1, -(p+2)/2, MPFR_RNDN);
1221 mpfr_add_ui (c, c, 1, MPFR_RNDN);
1222 inex = mpfr_mul (a, b, c, MPFR_RNDN);
1223 MPFR_ASSERTN(inex > 0);
1224 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0);
1225 /* For RNDU, b*c should round to 1 */
1226 inex = mpfr_mul (a, b, c, MPFR_RNDU);
1227 MPFR_ASSERTN(inex > 0);
1228 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0);
1229 /* For RNDD, b*c should round to 1-2^(-p) */
1230 inex = mpfr_mul (a, b, c, MPFR_RNDD);
1231 MPFR_ASSERTN(inex < 0);
1232 mpfr_nextabove (a);
1233 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0);
1234
1235 /* same as above, but near emax, to exercise the case where a carry
1236 produces an overflow */
1237 mpfr_set_si_2exp (b, -1, -(p+2)/2, MPFR_RNDN);
1238 mpfr_add_ui (b, b, 1, MPFR_RNDN);
1239 mpfr_mul_2si (b, b, mpfr_get_emax (), MPFR_RNDN);
1240 mpfr_set_si_2exp (c, 1, -(p+2)/2, MPFR_RNDN);
1241 mpfr_add_ui (c, c, 1, MPFR_RNDN);
1242 /* b*c should round to 2^emax */
1243 mpfr_clear_flags ();
1244 inex = mpfr_mul (a, b, c, MPFR_RNDN);
1245 MPFR_ASSERTN(inex > 0);
1246 MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0);
1247 MPFR_ASSERTN(mpfr_overflow_p ());
1248 /* idem for RNDU */
1249 mpfr_clear_flags ();
1250 inex = mpfr_mul (a, b, c, MPFR_RNDU);
1251 MPFR_ASSERTN(inex > 0);
1252 MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0);
1253 MPFR_ASSERTN(mpfr_overflow_p ());
1254 /* For RNDD, b*c should round to (1-2^(-p))*2^emax */
1255 mpfr_clear_flags ();
1256 inex = mpfr_mul (a, b, c, MPFR_RNDD);
1257 MPFR_ASSERTN(inex < 0);
1258 MPFR_ASSERTN(!mpfr_inf_p (a));
1259 MPFR_ASSERTN(!mpfr_overflow_p ());
1260 mpfr_nextabove (a);
1261 MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0);
1262
1263 end_of_loop:
1264 mpfr_clear (a);
1265 mpfr_clear (b);
1266 mpfr_clear (c);
1267 }
1268 }
1269
1270 /* check special underflow case for precision = 64 */
1271 static void
coverage2(void)1272 coverage2 (void)
1273 {
1274 mpfr_t a, b, c;
1275 int inex;
1276 mpfr_exp_t emin;
1277
1278 emin = mpfr_get_emin (); /* save emin */
1279 set_emin (0);
1280
1281 mpfr_init2 (a, 64);
1282 mpfr_init2 (b, 64);
1283 mpfr_init2 (c, 64);
1284
1285 mpfr_set_str_binary (b, "1111110110100001011100100000100000110110001100100010011010011001E-64"); /* 18276014142440744601/2^64 */
1286 mpfr_set_str_binary (c, "1000000100110010000111000100010010010001000100101010111101010100E-64"); /* 9309534460545511252/2^64 */
1287 /* since 1/2-2^-66 < b0*c0 < 1/2, b0*c0 should be rounded to 1/2 */
1288 inex = mpfr_mul (a, b, c, MPFR_RNDN);
1289 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, -1) == 0);
1290 MPFR_ASSERTN(inex > 0);
1291
1292 mpfr_clear (a);
1293 mpfr_clear (b);
1294 mpfr_clear (c);
1295
1296 set_emin (emin); /* restore emin */
1297 }
1298
1299 int
main(int argc,char * argv[])1300 main (int argc, char *argv[])
1301 {
1302 tests_start_mpfr ();
1303
1304 coverage (1024);
1305 coverage2 ();
1306 testall_rndf (9);
1307 check_nans ();
1308 check_exact ();
1309 check_float ();
1310
1311 check53("6.9314718055994530941514e-1", "0.0", MPFR_RNDZ, "0.0");
1312 check53("0.0", "6.9314718055994530941514e-1", MPFR_RNDZ, "0.0");
1313 check_sign();
1314 check53("-4.165000000e4", "-0.00004801920768307322868063274915", MPFR_RNDN,
1315 "2.0");
1316 check53("2.71331408349172961467e-08", "-6.72658901114033715233e-165",
1317 MPFR_RNDZ, "-1.8251348697787782844e-172");
1318 check53("2.71331408349172961467e-08", "-6.72658901114033715233e-165",
1319 MPFR_RNDA, "-1.8251348697787786e-172");
1320 check53("0.31869277231188065", "0.88642843322303122", MPFR_RNDZ,
1321 "2.8249833483992453642e-1");
1322 check("8.47622108205396074254e-01", "3.24039313247872939883e-01", MPFR_RNDU,
1323 28, 45, 2, "0.375");
1324 check("8.47622108205396074254e-01", "3.24039313247872939883e-01", MPFR_RNDA,
1325 28, 45, 2, "0.375");
1326 check("2.63978122803639081440e-01", "6.8378615379333496093e-1", MPFR_RNDN,
1327 34, 23, 31, "0.180504585267044603");
1328 check("1.0", "0.11835170935876249132", MPFR_RNDU, 6, 41, 36,
1329 "0.1183517093595583");
1330 check53("67108865.0", "134217729.0", MPFR_RNDN, "9.007199456067584e15");
1331 check("1.37399642157394197284e-01", "2.28877275604219221350e-01", MPFR_RNDN,
1332 49, 15, 32, "0.0314472340833162888");
1333 check("4.03160720978664954828e-01", "5.854828e-1"
1334 /*"5.85483042917246621073e-01"*/, MPFR_RNDZ,
1335 51, 22, 32, "0.2360436821472831");
1336 check("3.90798504668055102229e-14", "9.85394674650308388664e-04", MPFR_RNDN,
1337 46, 22, 12, "0.385027296503914762e-16");
1338 check("4.58687081072827851358e-01", "2.20543551472118792844e-01", MPFR_RNDN,
1339 49, 3, 2, "0.09375");
1340 check_max();
1341 check_min();
1342 small_prec ();
1343
1344 check_regression ();
1345 test_generic (MPFR_PREC_MIN, 500, 100);
1346
1347 data_check ("data/mulpi", mpfr_mulpi, "mpfr_mulpi");
1348
1349 valgrind20110503 ();
1350 test_underflow (128);
1351 bug20161209 ();
1352 bug20161209a ();
1353 bug20170602 ();
1354 test_underflow2 ();
1355
1356 tests_end_mpfr ();
1357 return 0;
1358 }
1359