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