tmul.c revision 1.1.1.3.4.1 1 1.1 mrg /* Test file for mpfr_mul.
2 1.1 mrg
3 1.1.1.3.4.1 christos 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.4.1 christos 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.4.1 christos mpfr_out_str (stdout, 10, 0, zz, MPFR_RNDN);
69 1.1.1.3.4.1 christos putchar ('\n');
70 1.1 mrg exit (1);
71 1.1 mrg }
72 1.1.1.3.4.1 christos mpfr_clear (xx);
73 1.1.1.3.4.1 christos mpfr_clear (yy);
74 1.1.1.3.4.1 christos 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.4.1 christos mpfr_out_str (stdout, 10, 0, zz, MPFR_RNDN);
92 1.1.1.3.4.1 christos 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.4.1 christos mpfr_out_str (stdout, 10, 0, zz, MPFR_RNDN);
114 1.1.1.3.4.1 christos 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.4.1 christos mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
181 1.1.1.3.4.1 christos 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.4.1 christos 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.4.1 christos /* for prec=1, ensure PREC(c) >= 1 */
217 1.1.1.3.4.1 christos 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.4.1 christos if ((inexact == 0) && mpfr_cmp (c, d) && rnd != MPFR_RNDF)
231 1.1 mrg {
232 1.1.1.3.4.1 christos printf ("rnd=%s: inexact=0 but results differ\n",
233 1.1.1.3.4.1 christos mpfr_print_rnd_mode (rnd));
234 1.1.1.3.4.1 christos printf ("a=");
235 1.1.1.3.4.1 christos mpfr_out_str (stdout, 2, 0, a, rnd);
236 1.1.1.3.4.1 christos printf ("\nb=");
237 1.1.1.3.4.1 christos mpfr_out_str (stdout, 2, 0, b, rnd);
238 1.1.1.3.4.1 christos printf ("\nc=");
239 1.1.1.3.4.1 christos mpfr_out_str (stdout, 2, 0, c, rnd);
240 1.1.1.3.4.1 christos printf ("\nd=");
241 1.1.1.3.4.1 christos mpfr_out_str (stdout, 2, 0, d, rnd);
242 1.1.1.3.4.1 christos printf ("\n");
243 1.1 mrg exit (1);
244 1.1 mrg }
245 1.1.1.3.4.1 christos 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.4.1 christos 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.4.1 christos static void
638 1.1.1.3.4.1 christos testall_rndf (mpfr_prec_t pmax)
639 1.1.1.3.4.1 christos {
640 1.1.1.3.4.1 christos mpfr_t a, b, c, d;
641 1.1.1.3.4.1 christos mpfr_prec_t pa, pb, pc;
642 1.1.1.3.4.1 christos
643 1.1.1.3.4.1 christos for (pa = MPFR_PREC_MIN; pa <= pmax; pa++)
644 1.1.1.3.4.1 christos {
645 1.1.1.3.4.1 christos mpfr_init2 (a, pa);
646 1.1.1.3.4.1 christos mpfr_init2 (d, pa);
647 1.1.1.3.4.1 christos for (pb = MPFR_PREC_MIN; pb <= pmax; pb++)
648 1.1.1.3.4.1 christos {
649 1.1.1.3.4.1 christos mpfr_init2 (b, pb);
650 1.1.1.3.4.1 christos mpfr_set_ui (b, 1, MPFR_RNDN);
651 1.1.1.3.4.1 christos while (mpfr_cmp_ui (b, 2) < 0)
652 1.1.1.3.4.1 christos {
653 1.1.1.3.4.1 christos for (pc = MPFR_PREC_MIN; pc <= pmax; pc++)
654 1.1.1.3.4.1 christos {
655 1.1.1.3.4.1 christos mpfr_init2 (c, pc);
656 1.1.1.3.4.1 christos mpfr_set_ui (c, 1, MPFR_RNDN);
657 1.1.1.3.4.1 christos while (mpfr_cmp_ui (c, 2) < 0)
658 1.1.1.3.4.1 christos {
659 1.1.1.3.4.1 christos mpfr_mul (a, b, c, MPFR_RNDF);
660 1.1.1.3.4.1 christos mpfr_mul (d, b, c, MPFR_RNDD);
661 1.1.1.3.4.1 christos if (!mpfr_equal_p (a, d))
662 1.1.1.3.4.1 christos {
663 1.1.1.3.4.1 christos mpfr_mul (d, b, c, MPFR_RNDU);
664 1.1.1.3.4.1 christos if (!mpfr_equal_p (a, d))
665 1.1.1.3.4.1 christos {
666 1.1.1.3.4.1 christos printf ("Error: mpfr_mul(a,b,c,RNDF) does not "
667 1.1.1.3.4.1 christos "match RNDD/RNDU\n");
668 1.1.1.3.4.1 christos printf ("b="); mpfr_dump (b);
669 1.1.1.3.4.1 christos printf ("c="); mpfr_dump (c);
670 1.1.1.3.4.1 christos printf ("a="); mpfr_dump (a);
671 1.1.1.3.4.1 christos exit (1);
672 1.1.1.3.4.1 christos }
673 1.1.1.3.4.1 christos }
674 1.1.1.3.4.1 christos mpfr_nextabove (c);
675 1.1.1.3.4.1 christos }
676 1.1.1.3.4.1 christos mpfr_clear (c);
677 1.1.1.3.4.1 christos }
678 1.1.1.3.4.1 christos mpfr_nextabove (b);
679 1.1.1.3.4.1 christos }
680 1.1.1.3.4.1 christos mpfr_clear (b);
681 1.1.1.3.4.1 christos }
682 1.1.1.3.4.1 christos mpfr_clear (a);
683 1.1.1.3.4.1 christos mpfr_clear (d);
684 1.1.1.3.4.1 christos }
685 1.1.1.3.4.1 christos }
686 1.1.1.3.4.1 christos
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.4.1 christos /* checks special case where no underflow should occur */
792 1.1.1.3.4.1 christos static void
793 1.1.1.3.4.1 christos bug20161209 (void)
794 1.1.1.3.4.1 christos {
795 1.1.1.3.4.1 christos mpfr_exp_t emin;
796 1.1.1.3.4.1 christos mpfr_t x, y, z;
797 1.1.1.3.4.1 christos
798 1.1.1.3.4.1 christos emin = mpfr_get_emin ();
799 1.1.1.3.4.1 christos set_emin (-1);
800 1.1.1.3.4.1 christos
801 1.1.1.3.4.1 christos /* test for mpfr_mul_1 for 64-bit limb, mpfr_mul_2 for 32-bit limb */
802 1.1.1.3.4.1 christos mpfr_init2 (x, 53);
803 1.1.1.3.4.1 christos mpfr_init2 (y, 53);
804 1.1.1.3.4.1 christos mpfr_init2 (z, 53);
805 1.1.1.3.4.1 christos mpfr_set_str_binary (x, "0.101000001E-1"); /* x = 321/2^10 */
806 1.1.1.3.4.1 christos mpfr_set_str_binary (y, "0.110011000010100101111000011011000111011000001E-1");
807 1.1.1.3.4.1 christos /* y = 28059810762433/2^46 */
808 1.1.1.3.4.1 christos /* x * y = (2^53+1)/2^56 = 0.001...000[1]000..., and should round to 0.25 */
809 1.1.1.3.4.1 christos mpfr_mul (z, x, y, MPFR_RNDN);
810 1.1.1.3.4.1 christos MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
811 1.1.1.3.4.1 christos
812 1.1.1.3.4.1 christos /* test for mpfr_mul_2 for 64-bit limb */
813 1.1.1.3.4.1 christos mpfr_set_prec (x, 65);
814 1.1.1.3.4.1 christos mpfr_set_prec (y, 65);
815 1.1.1.3.4.1 christos mpfr_set_prec (z, 65);
816 1.1.1.3.4.1 christos mpfr_set_str_binary (x, "0.101101000010010110100001E-1"); /* 11806113/2^25 */
817 1.1.1.3.4.1 christos mpfr_set_str_binary (y, "0.101101011110010101011010001111111001100001E-1");
818 1.1.1.3.4.1 christos /* y = 3124947910241/2^43 */
819 1.1.1.3.4.1 christos /* x * y = (2^65+1)/2^68 = 0.001...000[1]000..., and should round to 0.25 */
820 1.1.1.3.4.1 christos mpfr_mul (z, x, y, MPFR_RNDN);
821 1.1.1.3.4.1 christos MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
822 1.1.1.3.4.1 christos
823 1.1.1.3.4.1 christos /* test for the generic code */
824 1.1.1.3.4.1 christos mpfr_set_prec (x, 54);
825 1.1.1.3.4.1 christos mpfr_set_prec (y, 55);
826 1.1.1.3.4.1 christos mpfr_set_prec (z, 53);
827 1.1.1.3.4.1 christos mpfr_set_str_binary (x, "0.101000001E-1");
828 1.1.1.3.4.1 christos mpfr_set_str_binary (y, "0.110011000010100101111000011011000111011000001E-1");
829 1.1.1.3.4.1 christos mpfr_mul (z, x, y, MPFR_RNDN);
830 1.1.1.3.4.1 christos MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
831 1.1.1.3.4.1 christos
832 1.1.1.3.4.1 christos /* another test for the generic code */
833 1.1.1.3.4.1 christos mpfr_set_prec (x, 66);
834 1.1.1.3.4.1 christos mpfr_set_prec (y, 67);
835 1.1.1.3.4.1 christos mpfr_set_prec (z, 65);
836 1.1.1.3.4.1 christos mpfr_set_str_binary (x, "0.101101000010010110100001E-1");
837 1.1.1.3.4.1 christos mpfr_set_str_binary (y, "0.101101011110010101011010001111111001100001E-1");
838 1.1.1.3.4.1 christos mpfr_mul (z, x, y, MPFR_RNDN);
839 1.1.1.3.4.1 christos MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
840 1.1.1.3.4.1 christos
841 1.1.1.3.4.1 christos mpfr_clear (x);
842 1.1.1.3.4.1 christos mpfr_clear (y);
843 1.1.1.3.4.1 christos mpfr_clear (z);
844 1.1.1.3.4.1 christos set_emin (emin);
845 1.1.1.3.4.1 christos }
846 1.1.1.3.4.1 christos
847 1.1.1.3.4.1 christos /* test for case in mpfr_mul_1() where:
848 1.1.1.3.4.1 christos ax = __gmpfr_emin - 1
849 1.1.1.3.4.1 christos ap[0] == ~mask
850 1.1.1.3.4.1 christos rnd_mode = MPFR_RNDZ.
851 1.1.1.3.4.1 christos Whatever the values of rb and sb, we should round to zero (underflow). */
852 1.1.1.3.4.1 christos static void
853 1.1.1.3.4.1 christos bug20161209a (void)
854 1.1.1.3.4.1 christos {
855 1.1.1.3.4.1 christos mpfr_exp_t emin;
856 1.1.1.3.4.1 christos mpfr_t x, y, z;
857 1.1.1.3.4.1 christos
858 1.1.1.3.4.1 christos emin = mpfr_get_emin ();
859 1.1.1.3.4.1 christos set_emin (-1);
860 1.1.1.3.4.1 christos
861 1.1.1.3.4.1 christos mpfr_init2 (x, 53);
862 1.1.1.3.4.1 christos mpfr_init2 (y, 53);
863 1.1.1.3.4.1 christos mpfr_init2 (z, 53);
864 1.1.1.3.4.1 christos
865 1.1.1.3.4.1 christos /* case rb = sb = 0 */
866 1.1.1.3.4.1 christos mpfr_set_str_binary (x, "0.11010010100110000110110011111E-1");
867 1.1.1.3.4.1 christos mpfr_set_str_binary (y, "0.1001101110011000110100001");
868 1.1.1.3.4.1 christos /* x = 441650591/2^30, y = 20394401/2^25 */
869 1.1.1.3.4.1 christos /* x * y = (2^53-1)/2^55 = 0.00111...111[0]000..., and should round to 0 */
870 1.1.1.3.4.1 christos mpfr_mul (z, x, y, MPFR_RNDZ);
871 1.1.1.3.4.1 christos MPFR_ASSERTN(mpfr_zero_p (z));
872 1.1.1.3.4.1 christos
873 1.1.1.3.4.1 christos /* case rb = 1, sb = 0 */
874 1.1.1.3.4.1 christos mpfr_set_str_binary (x, "0.111111111000000000000000000111111111E-1");
875 1.1.1.3.4.1 christos mpfr_set_str_binary (y, "0.1000000001000000001");
876 1.1.1.3.4.1 christos /* x = 68585259519/2^37, y = 262657/2^19 */
877 1.1.1.3.4.1 christos /* x * y = (2^54-1)/2^56 = 0.00111...111[1]000..., and should round to 0 */
878 1.1.1.3.4.1 christos mpfr_mul (z, x, y, MPFR_RNDZ);
879 1.1.1.3.4.1 christos MPFR_ASSERTN(mpfr_zero_p (z));
880 1.1.1.3.4.1 christos
881 1.1.1.3.4.1 christos /* case rb = 0, sb = 1 */
882 1.1.1.3.4.1 christos mpfr_set_str_binary (x, "0.110010011001011110001100100001000001E-1");
883 1.1.1.3.4.1 christos mpfr_set_str_binary (y, "0.10100010100010111101");
884 1.1.1.3.4.1 christos /* x = 541144371852^37, y = 665789/2^20 */
885 1.1.1.3.4.1 christos /* x * y = (2^55-3)/2^57 = 0.00111...111[0]100..., and should round to 0 */
886 1.1.1.3.4.1 christos mpfr_mul (z, x, y, MPFR_RNDZ);
887 1.1.1.3.4.1 christos MPFR_ASSERTN(mpfr_zero_p (z));
888 1.1.1.3.4.1 christos
889 1.1.1.3.4.1 christos /* case rb = sb = 1 */
890 1.1.1.3.4.1 christos mpfr_set_str_binary (x, "0.10100110001001001010001111110010100111E-1");
891 1.1.1.3.4.1 christos mpfr_set_str_binary (y, "0.110001010011101001");
892 1.1.1.3.4.1 christos /* x = 178394823847/2^39, y = 201961/2^18 */
893 1.1.1.3.4.1 christos /* x * y = (2^55-1)/2^57 = 0.00111...111[1]100..., and should round to 0 */
894 1.1.1.3.4.1 christos mpfr_mul (z, x, y, MPFR_RNDZ);
895 1.1.1.3.4.1 christos MPFR_ASSERTN(mpfr_zero_p (z));
896 1.1.1.3.4.1 christos
897 1.1.1.3.4.1 christos /* similar test for mpfr_mul_2 (we only check rb = sb = 1 here) */
898 1.1.1.3.4.1 christos mpfr_set_prec (x, 65);
899 1.1.1.3.4.1 christos mpfr_set_prec (y, 65);
900 1.1.1.3.4.1 christos mpfr_set_prec (z, 65);
901 1.1.1.3.4.1 christos /* 2^67-1 = 193707721 * 761838257287 */
902 1.1.1.3.4.1 christos mpfr_set_str_binary (x, "0.1011100010111011111011001001E-1");
903 1.1.1.3.4.1 christos mpfr_set_str_binary (y, "0.1011000101100001000110010100010010000111");
904 1.1.1.3.4.1 christos mpfr_mul (z, x, y, MPFR_RNDZ);
905 1.1.1.3.4.1 christos MPFR_ASSERTN(mpfr_zero_p (z));
906 1.1.1.3.4.1 christos
907 1.1.1.3.4.1 christos mpfr_clear (x);
908 1.1.1.3.4.1 christos mpfr_clear (y);
909 1.1.1.3.4.1 christos mpfr_clear (z);
910 1.1.1.3.4.1 christos set_emin (emin);
911 1.1.1.3.4.1 christos }
912 1.1.1.3.4.1 christos
913 1.1.1.3.4.1 christos /* bug for RNDF */
914 1.1.1.3.4.1 christos static void
915 1.1.1.3.4.1 christos bug20170602 (void)
916 1.1.1.3.4.1 christos {
917 1.1.1.3.4.1 christos mpfr_t x, u, y, yd, yu;
918 1.1.1.3.4.1 christos
919 1.1.1.3.4.1 christos mpfr_init2 (x, 493);
920 1.1.1.3.4.1 christos mpfr_init2 (u, 493);
921 1.1.1.3.4.1 christos mpfr_init2 (y, 503);
922 1.1.1.3.4.1 christos mpfr_init2 (yd, 503);
923 1.1.1.3.4.1 christos mpfr_init2 (yu, 503);
924 1.1.1.3.4.1 christos mpfr_set_str_binary (x, "0.1111100000000000001111111110000000001111111111111000000000000000000011111111111111111111111000000000000000000001111111111111111111111111111111111111111000000000011111111111111111111000000000011111111111111000000000000001110000000000000000000000000000000000000000011111111111110011111111111100000000000000011111111111111111110000000011111111111111111110011111111111110000000000001111111111111111000000000000000000000000000000000000111111111111111111111111111111011111111111111111111111111111100E44");
925 1.1.1.3.4.1 christos mpfr_set_str_binary (u, "0.1110000000000000001111111111111111111111111111111111111000000000111111111111111111111111111111000000000000000000001111111000000000000000011111111111111111111111111111111111111111111111111111111000000000000000011111111111111000000011111111111111111110000000000000001111111111111111111111111111111111111110000000000001111111111111111111111111111111111111000000000000000000000000000000000001111111111111000000000000000000001111111111100000000000000011111111111111111111111111111111111111111111111E35");
926 1.1.1.3.4.1 christos mpfr_mul (y, x, u, MPFR_RNDF);
927 1.1.1.3.4.1 christos mpfr_mul (yd, x, u, MPFR_RNDD);
928 1.1.1.3.4.1 christos mpfr_mul (yu, x, u, MPFR_RNDU);
929 1.1.1.3.4.1 christos if (mpfr_cmp (y, yd) != 0 && mpfr_cmp (y, yu) != 0)
930 1.1.1.3.4.1 christos {
931 1.1.1.3.4.1 christos printf ("RNDF is neither RNDD nor RNDU\n");
932 1.1.1.3.4.1 christos printf ("x"); mpfr_dump (x);
933 1.1.1.3.4.1 christos printf ("u"); mpfr_dump (u);
934 1.1.1.3.4.1 christos printf ("y(RNDF)="); mpfr_dump (y);
935 1.1.1.3.4.1 christos printf ("y(RNDD)="); mpfr_dump (yd);
936 1.1.1.3.4.1 christos printf ("y(RNDU)="); mpfr_dump (yu);
937 1.1.1.3.4.1 christos exit (1);
938 1.1.1.3.4.1 christos }
939 1.1.1.3.4.1 christos mpfr_clear (x);
940 1.1.1.3.4.1 christos mpfr_clear (u);
941 1.1.1.3.4.1 christos mpfr_clear (y);
942 1.1.1.3.4.1 christos mpfr_clear (yd);
943 1.1.1.3.4.1 christos mpfr_clear (yu);
944 1.1.1.3.4.1 christos }
945 1.1.1.3.4.1 christos
946 1.1.1.3.4.1 christos /* Test for 1 to 3 limbs. */
947 1.1.1.3.4.1 christos static void
948 1.1.1.3.4.1 christos small_prec (void)
949 1.1.1.3.4.1 christos {
950 1.1.1.3.4.1 christos mpfr_exp_t emin, emax;
951 1.1.1.3.4.1 christos mpfr_t x, y, z1, z2, zz;
952 1.1.1.3.4.1 christos int xq, yq, zq;
953 1.1.1.3.4.1 christos mpfr_rnd_t rnd;
954 1.1.1.3.4.1 christos mpfr_flags_t flags1, flags2;
955 1.1.1.3.4.1 christos int inex1, inex2;
956 1.1.1.3.4.1 christos int i, j, r, s, ediff;
957 1.1.1.3.4.1 christos
958 1.1.1.3.4.1 christos emin = mpfr_get_emin ();
959 1.1.1.3.4.1 christos emax = mpfr_get_emax ();
960 1.1.1.3.4.1 christos
961 1.1.1.3.4.1 christos /* The mpfr_mul implementation doesn't extend the exponent range,
962 1.1.1.3.4.1 christos so that it is OK to reduce it here for the test to make sure that
963 1.1.1.3.4.1 christos mpfr_mul_2si can be used. */
964 1.1.1.3.4.1 christos set_emin (-1000);
965 1.1.1.3.4.1 christos set_emax (1000);
966 1.1.1.3.4.1 christos
967 1.1.1.3.4.1 christos mpfr_inits2 (3 * GMP_NUMB_BITS, x, y, z1, z2, (mpfr_ptr) 0);
968 1.1.1.3.4.1 christos mpfr_init2 (zz, 6 * GMP_NUMB_BITS);
969 1.1.1.3.4.1 christos for (i = 0; i < 3; i++)
970 1.1.1.3.4.1 christos for (j = 0; j < 10000; j++)
971 1.1.1.3.4.1 christos {
972 1.1.1.3.4.1 christos xq = i * GMP_NUMB_BITS + 1 + randlimb () % GMP_NUMB_BITS;
973 1.1.1.3.4.1 christos mpfr_set_prec (x, xq);
974 1.1.1.3.4.1 christos yq = i * GMP_NUMB_BITS + 1 + randlimb () % GMP_NUMB_BITS;
975 1.1.1.3.4.1 christos mpfr_set_prec (y, yq);
976 1.1.1.3.4.1 christos zq = i * GMP_NUMB_BITS + 1 + randlimb () % (GMP_NUMB_BITS-1);
977 1.1.1.3.4.1 christos mpfr_set_prec (z1, zq);
978 1.1.1.3.4.1 christos mpfr_set_prec (z2, zq);
979 1.1.1.3.4.1 christos s = r = randlimb () & 0x7f;
980 1.1.1.3.4.1 christos do mpfr_urandomb (x, RANDS); while (mpfr_zero_p (x));
981 1.1.1.3.4.1 christos if (s & 1)
982 1.1.1.3.4.1 christos mpfr_neg (x, x, MPFR_RNDN);
983 1.1.1.3.4.1 christos s >>= 1;
984 1.1.1.3.4.1 christos if (s & 1)
985 1.1.1.3.4.1 christos {
986 1.1.1.3.4.1 christos do mpfr_urandomb (y, RANDS); while (mpfr_zero_p (y));
987 1.1.1.3.4.1 christos }
988 1.1.1.3.4.1 christos else
989 1.1.1.3.4.1 christos {
990 1.1.1.3.4.1 christos mpfr_ui_div (y, 1, x, MPFR_RNDN);
991 1.1.1.3.4.1 christos mpfr_set_exp (y, 0);
992 1.1.1.3.4.1 christos }
993 1.1.1.3.4.1 christos s >>= 1;
994 1.1.1.3.4.1 christos if (s & 1)
995 1.1.1.3.4.1 christos mpfr_neg (y, y, MPFR_RNDN);
996 1.1.1.3.4.1 christos s >>= 1;
997 1.1.1.3.4.1 christos rnd = RND_RAND_NO_RNDF ();
998 1.1.1.3.4.1 christos inex1 = mpfr_mul (zz, x, y, MPFR_RNDN);
999 1.1.1.3.4.1 christos MPFR_ASSERTN (inex1 == 0);
1000 1.1.1.3.4.1 christos if (s == 0)
1001 1.1.1.3.4.1 christos {
1002 1.1.1.3.4.1 christos ediff = __gmpfr_emin - MPFR_EXP (x);
1003 1.1.1.3.4.1 christos mpfr_set_exp (x, __gmpfr_emin);
1004 1.1.1.3.4.1 christos }
1005 1.1.1.3.4.1 christos else if (s == 1)
1006 1.1.1.3.4.1 christos {
1007 1.1.1.3.4.1 christos ediff = __gmpfr_emax - MPFR_EXP (x) + 1;
1008 1.1.1.3.4.1 christos mpfr_set_exp (x, __gmpfr_emax);
1009 1.1.1.3.4.1 christos mpfr_mul_2ui (y, y, 1, MPFR_RNDN);
1010 1.1.1.3.4.1 christos }
1011 1.1.1.3.4.1 christos else
1012 1.1.1.3.4.1 christos ediff = 0;
1013 1.1.1.3.4.1 christos mpfr_clear_flags ();
1014 1.1.1.3.4.1 christos inex1 = mpfr_mul_2si (z1, zz, ediff, rnd);
1015 1.1.1.3.4.1 christos flags1 = __gmpfr_flags;
1016 1.1.1.3.4.1 christos mpfr_clear_flags ();
1017 1.1.1.3.4.1 christos inex2 = mpfr_mul (z2, x, y, rnd);
1018 1.1.1.3.4.1 christos flags2 = __gmpfr_flags;
1019 1.1.1.3.4.1 christos if (!(mpfr_equal_p (z1, z2) &&
1020 1.1.1.3.4.1 christos SAME_SIGN (inex1, inex2) &&
1021 1.1.1.3.4.1 christos flags1 == flags2))
1022 1.1.1.3.4.1 christos {
1023 1.1.1.3.4.1 christos printf ("Error in small_prec on i = %d, j = %d\n", i, j);
1024 1.1.1.3.4.1 christos printf ("r = 0x%x, xq = %d, yq = %d, zq = %d, rnd = %s\n",
1025 1.1.1.3.4.1 christos r, xq, yq, zq, mpfr_print_rnd_mode (rnd));
1026 1.1.1.3.4.1 christos printf ("x = ");
1027 1.1.1.3.4.1 christos mpfr_dump (x);
1028 1.1.1.3.4.1 christos printf ("y = ");
1029 1.1.1.3.4.1 christos mpfr_dump (y);
1030 1.1.1.3.4.1 christos printf ("ediff = %d\n", ediff);
1031 1.1.1.3.4.1 christos printf ("zz = ");
1032 1.1.1.3.4.1 christos mpfr_dump (zz);
1033 1.1.1.3.4.1 christos printf ("Expected ");
1034 1.1.1.3.4.1 christos mpfr_dump (z1);
1035 1.1.1.3.4.1 christos printf ("with inex = %d and flags =", inex1);
1036 1.1.1.3.4.1 christos flags_out (flags1);
1037 1.1.1.3.4.1 christos printf ("Got ");
1038 1.1.1.3.4.1 christos mpfr_dump (z2);
1039 1.1.1.3.4.1 christos printf ("with inex = %d and flags =", inex2);
1040 1.1.1.3.4.1 christos flags_out (flags2);
1041 1.1.1.3.4.1 christos exit (1);
1042 1.1.1.3.4.1 christos }
1043 1.1.1.3.4.1 christos }
1044 1.1.1.3.4.1 christos mpfr_clears (x, y, z1, z2, zz, (mpfr_ptr) 0);
1045 1.1.1.3.4.1 christos
1046 1.1.1.3.4.1 christos set_emin (emin);
1047 1.1.1.3.4.1 christos set_emax (emax);
1048 1.1.1.3.4.1 christos }
1049 1.1.1.3.4.1 christos
1050 1.1.1.3.4.1 christos /* check ax < __gmpfr_emin with rnd_mode == MPFR_RNDN, rb = 0 and sb <> 0 */
1051 1.1.1.3.4.1 christos static void
1052 1.1.1.3.4.1 christos test_underflow2 (void)
1053 1.1.1.3.4.1 christos {
1054 1.1.1.3.4.1 christos mpfr_t x, y, z;
1055 1.1.1.3.4.1 christos mpfr_exp_t emin;
1056 1.1.1.3.4.1 christos
1057 1.1.1.3.4.1 christos emin = mpfr_get_emin ();
1058 1.1.1.3.4.1 christos mpfr_set_emin (0);
1059 1.1.1.3.4.1 christos
1060 1.1.1.3.4.1 christos mpfr_init2 (x, 24);
1061 1.1.1.3.4.1 christos mpfr_init2 (y, 24);
1062 1.1.1.3.4.1 christos mpfr_init2 (z, 24);
1063 1.1.1.3.4.1 christos
1064 1.1.1.3.4.1 christos mpfr_set_ui_2exp (x, 12913, -14, MPFR_RNDN);
1065 1.1.1.3.4.1 christos mpfr_set_ui_2exp (y, 5197, -13, MPFR_RNDN);
1066 1.1.1.3.4.1 christos mpfr_clear_underflow ();
1067 1.1.1.3.4.1 christos /* x*y = 0.0111111111111111111111111[01] thus underflow */
1068 1.1.1.3.4.1 christos mpfr_mul (z, y, x, MPFR_RNDN);
1069 1.1.1.3.4.1 christos MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -1) == 0);
1070 1.1.1.3.4.1 christos MPFR_ASSERTN(mpfr_underflow_p ());
1071 1.1.1.3.4.1 christos
1072 1.1.1.3.4.1 christos mpfr_set_prec (x, 65);
1073 1.1.1.3.4.1 christos mpfr_set_prec (y, 65);
1074 1.1.1.3.4.1 christos mpfr_set_prec (z, 65);
1075 1.1.1.3.4.1 christos mpfr_set_str_binary (x, "0.10011101000110111011111100101001111");
1076 1.1.1.3.4.1 christos mpfr_set_str_binary (y, "0.110100001001000111000011011110011");
1077 1.1.1.3.4.1 christos mpfr_clear_underflow ();
1078 1.1.1.3.4.1 christos /* x*y = 0.011111111111111111111111111111111111111111111111111111111111111111[01] thus underflow */
1079 1.1.1.3.4.1 christos mpfr_mul (z, y, x, MPFR_RNDN);
1080 1.1.1.3.4.1 christos MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -1) == 0);
1081 1.1.1.3.4.1 christos MPFR_ASSERTN(mpfr_underflow_p ());
1082 1.1.1.3.4.1 christos
1083 1.1.1.3.4.1 christos mpfr_clear (y);
1084 1.1.1.3.4.1 christos mpfr_clear (x);
1085 1.1.1.3.4.1 christos mpfr_clear (z);
1086 1.1.1.3.4.1 christos
1087 1.1.1.3.4.1 christos mpfr_set_emin (emin);
1088 1.1.1.3.4.1 christos }
1089 1.1.1.3.4.1 christos
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.4.1 christos 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.4.1 christos small_prec ();
1132 1.1 mrg
1133 1.1 mrg check_regression ();
1134 1.1.1.3.4.1 christos 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.4.1 christos bug20161209 ();
1141 1.1.1.3.4.1 christos bug20161209a ();
1142 1.1.1.3.4.1 christos bug20170602 ();
1143 1.1.1.3.4.1 christos 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