1 /* Test file for mpfr_exp.
2
3 Copyright 1999, 2001-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #include "mpfr-test.h"
24
25 #ifdef CHECK_EXTERNAL
26 static int
test_exp(mpfr_ptr a,mpfr_srcptr b,mpfr_rnd_t rnd_mode)27 test_exp (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode)
28 {
29 int res;
30 int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_get_prec (a)>=53;
31 if (ok)
32 {
33 mpfr_print_raw (b);
34 }
35 res = mpfr_exp (a, b, rnd_mode);
36 if (ok)
37 {
38 printf (" ");
39 mpfr_print_raw (a);
40 printf ("\n");
41 }
42 return res;
43 }
44 #else
45 #define test_exp mpfr_exp
46 #endif
47
48 /* returns the number of ulp of error */
49 static void
check3(const char * op,mpfr_rnd_t rnd,const char * res)50 check3 (const char *op, mpfr_rnd_t rnd, const char *res)
51 {
52 mpfr_t x, y;
53
54 mpfr_inits2 (53, x, y, (mpfr_ptr) 0);
55 /* y negative. If we forget to set the sign in mpfr_exp, we'll see it. */
56 mpfr_set_si (y, -1, MPFR_RNDN);
57 mpfr_set_str1 (x, op);
58 test_exp (y, x, rnd);
59 if (mpfr_cmp_str1 (y, res) )
60 {
61 printf ("mpfr_exp failed for x=%s, rnd=%s\n",
62 op, mpfr_print_rnd_mode (rnd));
63 printf ("expected result is %s, got ", res);
64 mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN);
65 putchar('\n');
66 exit (1);
67 }
68 mpfr_clears (x, y, (mpfr_ptr) 0);
69 }
70
71 /* expx is the value of exp(X) rounded toward -infinity */
72 static void
check_worst_case(const char * Xs,const char * expxs)73 check_worst_case (const char *Xs, const char *expxs)
74 {
75 mpfr_t x, y;
76
77 mpfr_inits2 (53, x, y, (mpfr_ptr) 0);
78 mpfr_set_str1(x, Xs);
79 test_exp(y, x, MPFR_RNDD);
80 if (mpfr_cmp_str1 (y, expxs))
81 {
82 printf ("exp(x) rounded toward -infinity is wrong\n");
83 exit(1);
84 }
85 mpfr_set_str1(x, Xs);
86 test_exp(x, x, MPFR_RNDU);
87 mpfr_nexttoinf (y);
88 if (mpfr_cmp(x,y))
89 {
90 printf ("exp(x) rounded toward +infinity is wrong\n");
91 exit(1);
92 }
93 mpfr_clears (x, y, (mpfr_ptr) 0);
94 }
95
96 /* worst cases communicated by Jean-Michel Muller and Vincent Lefevre */
97 static int
check_worst_cases(void)98 check_worst_cases (void)
99 {
100 mpfr_t x; mpfr_t y;
101
102 mpfr_init(x);
103 mpfr_set_prec (x, 53);
104
105 check_worst_case("4.44089209850062517562e-16", "1.00000000000000022204");
106 check_worst_case("6.39488462184069720009e-14", "1.00000000000006372680");
107 check_worst_case("1.84741111297455401935e-12", "1.00000000000184718907");
108 check_worst_case("1.76177628026265550074e-10", "1.00000000017617751702");
109 check3("1.76177628026265550074e-10", MPFR_RNDN, "1.00000000017617773906");
110 check_worst_case("7.54175277499595900852e-10", "1.00000000075417516676");
111 check3("7.54175277499595900852e-10", MPFR_RNDN, "1.00000000075417538881");
112 /* bug found by Vincent Lefe`vre on December 8, 1999 */
113 check3("-5.42410311287441459172e+02", MPFR_RNDN, "2.7176584868845723e-236");
114 /* further cases communicated by Vincent Lefe`vre on January 27, 2000 */
115 check3("-1.32920285897904911589e-10", MPFR_RNDN, "0.999999999867079769622");
116 check3("-1.44037948245738330735e-10", MPFR_RNDN, "0.9999999998559621072757");
117 check3("-1.66795910430705305937e-10", MPFR_RNDZ, "0.9999999998332040895832");
118 check3("-1.64310953745426656203e-10", MPFR_RNDN, "0.9999999998356891017792");
119 check3("-1.38323574826034659172e-10", MPFR_RNDZ, "0.9999999998616764251835");
120 check3("-1.23621668465115401498e-10", MPFR_RNDZ, "0.9999999998763783315425");
121
122 mpfr_set_prec (x, 601);
123 mpfr_set_str (x, "0.88b6ba510e10450edc258748bc9dfdd466f21b47ed264cdf24aa8f64af1f3fad9ec2301d43c0743f534b5aa20091ff6d352df458ef1ba519811ef6f5b11853534fd8fa32764a0a6d2d0dd20@0", 16, MPFR_RNDZ);
124 mpfr_init2 (y, 601);
125 mpfr_exp_2 (y, x, MPFR_RNDD);
126 mpfr_exp_3 (x, x, MPFR_RNDD);
127 if (mpfr_cmp (x, y))
128 {
129 printf ("mpfr_exp_2 and mpfr_exp_3 differ for prec=601\n");
130 printf ("mpfr_exp_2 gives ");
131 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
132 printf ("\nmpfr_exp_3 gives ");
133 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
134 printf ("\n");
135 exit (1);
136 }
137
138 mpfr_set_prec (x, 13001);
139 mpfr_set_prec (y, 13001);
140 mpfr_urandomb (x, RANDS);
141 mpfr_exp_3 (y, x, MPFR_RNDN);
142 mpfr_exp_2 (x, x, MPFR_RNDN);
143 if (mpfr_cmp (x, y))
144 {
145 printf ("mpfr_exp_2 and mpfr_exp_3 differ for prec=13001\n");
146 exit (1);
147 }
148
149 mpfr_set_prec (x, 118);
150 mpfr_set_str_binary (x, "0.1110010100011101010000111110011000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E-86");
151 mpfr_set_prec (y, 118);
152 mpfr_exp_2 (y, x, MPFR_RNDU);
153 mpfr_exp_3 (x, x, MPFR_RNDU);
154 if (mpfr_cmp (x, y))
155 {
156 printf ("mpfr_exp_2 and mpfr_exp_3 differ for prec=118\n");
157 printf ("mpfr_exp_2 gives ");
158 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
159 printf ("\nmpfr_exp_3 gives ");
160 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
161 printf ("\n");
162 exit (1);
163 }
164
165 mpfr_clear (x);
166 mpfr_clear (y);
167 return 0;
168 }
169
170 static void
compare_exp2_exp3(mpfr_prec_t p0,mpfr_prec_t p1)171 compare_exp2_exp3 (mpfr_prec_t p0, mpfr_prec_t p1)
172 {
173 mpfr_t x, y, z;
174 mpfr_prec_t prec;
175 mpfr_rnd_t rnd;
176
177 mpfr_init (x);
178 mpfr_init (y);
179 mpfr_init (z);
180 for (prec = p0; prec <= p1; prec ++)
181 {
182 mpfr_set_prec (x, prec);
183 mpfr_set_prec (y, prec);
184 mpfr_set_prec (z, prec);
185 do
186 mpfr_urandomb (x, RANDS);
187 while (MPFR_IS_ZERO (x)); /* 0 is handled by mpfr_exp only */
188 rnd = RND_RAND ();
189 mpfr_exp_2 (y, x, rnd);
190 mpfr_exp_3 (z, x, rnd);
191 if (mpfr_cmp (y,z))
192 {
193 printf ("mpfr_exp_2 and mpfr_exp_3 disagree for rnd=%s and\nx=",
194 mpfr_print_rnd_mode (rnd));
195 mpfr_dump (x);
196 printf ("mpfr_exp_2 gives ");
197 mpfr_dump (y);
198 printf ("mpfr_exp_3 gives ");
199 mpfr_dump (z);
200 exit (1);
201 }
202 }
203
204 mpfr_clear (x);
205 mpfr_clear (y);
206 mpfr_clear (z);
207 }
208
209 static void
check_large(void)210 check_large (void)
211 {
212 mpfr_t x, z;
213 mpfr_prec_t prec;
214
215 /* bug found by Patrick Pe'lissier on 7 Jun 2004 */
216 prec = 203780;
217 mpfr_init2 (x, prec);
218 mpfr_init2 (z, prec);
219 mpfr_set_ui (x, 3, MPFR_RNDN);
220 mpfr_sqrt (x, x, MPFR_RNDN);
221 mpfr_sub_ui (x, x, 1, MPFR_RNDN);
222 mpfr_exp_3 (z, x, MPFR_RNDN);
223 mpfr_clear (x);
224 mpfr_clear (z);
225 }
226
227 #define TEST_FUNCTION test_exp
228 #define TEST_RANDOM_EMIN -36
229 #define TEST_RANDOM_EMAX 36
230 #include "tgeneric.c"
231
232 static void
check_special(void)233 check_special (void)
234 {
235 mpfr_t x, y, z;
236 mpfr_exp_t emin, emax;
237
238 emin = mpfr_get_emin ();
239 emax = mpfr_get_emax ();
240
241 mpfr_init (x);
242 mpfr_init (y);
243 mpfr_init (z);
244
245 /* check exp(NaN) = NaN */
246 mpfr_set_nan (x);
247 test_exp (y, x, MPFR_RNDN);
248 if (!mpfr_nan_p (y))
249 {
250 printf ("Error for exp(NaN)\n");
251 exit (1);
252 }
253
254 /* check exp(+inf) = +inf */
255 mpfr_set_inf (x, 1);
256 test_exp (y, x, MPFR_RNDN);
257 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
258 {
259 printf ("Error for exp(+inf)\n");
260 exit (1);
261 }
262
263 /* check exp(-inf) = +0 */
264 mpfr_set_inf (x, -1);
265 test_exp (y, x, MPFR_RNDN);
266 if (MPFR_NOTZERO (y) || MPFR_IS_NEG (y))
267 {
268 printf ("Error for exp(-inf)\n");
269 exit (1);
270 }
271
272 /* Check overflow. Corner case of mpfr_exp_2 */
273 mpfr_set_prec (x, 64);
274 set_emax (MPFR_EMAX_DEFAULT);
275 set_emin (MPFR_EMIN_DEFAULT);
276 mpfr_set_str (x,
277 "0.1011000101110010000101111111010100001100000001110001100111001101E30",
278 2, MPFR_RNDN);
279 mpfr_exp (x, x, MPFR_RNDD);
280 if (mpfr_cmp_str (x,
281 ".1111111111111111111111111111111111111111111111111111111111111111E1073741823",
282 2, MPFR_RNDN) != 0)
283 {
284 printf ("Wrong overflow detection in mpfr_exp\n");
285 mpfr_dump (x);
286 exit (1);
287 }
288 /* Check underflow. Corner case of mpfr_exp_2 */
289 mpfr_set_str (x,
290 "-0.1011000101110010000101111111011111010001110011110111100110101100E30",
291 2, MPFR_RNDN);
292 mpfr_exp (x, x, MPFR_RNDN);
293 if (mpfr_cmp_str (x, "0.1E-1073741823", 2, MPFR_RNDN) != 0)
294 {
295 printf ("Wrong underflow (1) detection in mpfr_exp\n");
296 mpfr_dump (x);
297 exit (1);
298 }
299 mpfr_set_str (x,
300 "-0.1011001101110010000101111111011111010001110011110111100110111101E30",
301 2, MPFR_RNDN);
302 mpfr_exp (x, x, MPFR_RNDN);
303 if (mpfr_cmp_ui (x, 0) != 0)
304 {
305 printf ("Wrong underflow (2) detection in mpfr_exp\n");
306 mpfr_dump (x);
307 exit (1);
308 }
309 /* Check overflow. Corner case of mpfr_exp_3 */
310 if (MPFR_PREC_MAX >= MPFR_EXP_THRESHOLD + 10 && MPFR_PREC_MAX >= 64)
311 {
312 /* this ensures that for small MPFR_EXP_THRESHOLD, the following
313 mpfr_set_str conversion is exact */
314 mpfr_set_prec (x, (MPFR_EXP_THRESHOLD + 10 > 64)
315 ? MPFR_EXP_THRESHOLD + 10 : 64);
316 mpfr_set_str (x,
317 "0.1011000101110010000101111111010100001100000001110001100111001101E30",
318 2, MPFR_RNDN);
319 mpfr_clear_overflow ();
320 mpfr_exp (x, x, MPFR_RNDD);
321 if (!mpfr_overflow_p ())
322 {
323 printf ("Wrong overflow detection in mpfr_exp_3\n");
324 mpfr_dump (x);
325 exit (1);
326 }
327 /* Check underflow. Corner case of mpfr_exp_3 */
328 mpfr_set_str (x,
329 "-0.1011000101110010000101111111011111010001110011110111100110101100E30",
330 2, MPFR_RNDN);
331 mpfr_clear_underflow ();
332 mpfr_exp (x, x, MPFR_RNDN);
333 if (!mpfr_underflow_p ())
334 {
335 printf ("Wrong underflow detection in mpfr_exp_3\n");
336 mpfr_dump (x);
337 exit (1);
338 }
339 mpfr_set_prec (x, 53);
340 }
341
342 /* check overflow */
343 set_emax (10);
344 mpfr_set_ui (x, 7, MPFR_RNDN);
345 test_exp (y, x, MPFR_RNDN);
346 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
347 {
348 printf ("Error for exp(7) for emax=10\n");
349 exit (1);
350 }
351 set_emax (emax);
352
353 /* check underflow */
354 set_emin (-10);
355 mpfr_set_si (x, -9, MPFR_RNDN);
356 test_exp (y, x, MPFR_RNDN);
357 if (MPFR_NOTZERO (y) || MPFR_IS_NEG (y))
358 {
359 printf ("Error for exp(-9) for emin=-10\n");
360 printf ("Expected +0\n");
361 printf ("Got "); mpfr_dump (y);
362 exit (1);
363 }
364 set_emin (emin);
365
366 /* check case EXP(x) < -precy */
367 mpfr_set_prec (y, 2);
368 mpfr_set_str_binary (x, "-0.1E-3");
369 test_exp (y, x, MPFR_RNDD);
370 if (mpfr_cmp_ui_2exp (y, 3, -2))
371 {
372 printf ("Error for exp(-1/16), prec=2, RNDD\n");
373 printf ("expected 0.11, got ");
374 mpfr_dump (y);
375 exit (1);
376 }
377 test_exp (y, x, MPFR_RNDZ);
378 if (mpfr_cmp_ui_2exp (y, 3, -2))
379 {
380 printf ("Error for exp(-1/16), prec=2, RNDZ\n");
381 printf ("expected 0.11, got ");
382 mpfr_dump (y);
383 exit (1);
384 }
385 mpfr_set_str_binary (x, "0.1E-3");
386 test_exp (y, x, MPFR_RNDN);
387 if (mpfr_cmp_ui (y, 1))
388 {
389 printf ("Error for exp(1/16), prec=2, RNDN\n");
390 exit (1);
391 }
392 test_exp (y, x, MPFR_RNDU);
393 if (mpfr_cmp_ui_2exp (y, 3, -1))
394 {
395 printf ("Error for exp(1/16), prec=2, RNDU\n");
396 exit (1);
397 }
398
399 /* bug reported by Franky Backeljauw, 28 Mar 2003 */
400 mpfr_set_prec (x, 53);
401 mpfr_set_prec (y, 53);
402 mpfr_set_str_binary (x, "1.1101011000111101011110000111010010101001101001110111e28");
403 test_exp (y, x, MPFR_RNDN);
404
405 mpfr_set_prec (x, 153);
406 mpfr_set_prec (z, 153);
407 mpfr_set_str_binary (x, "1.1101011000111101011110000111010010101001101001110111e28");
408 test_exp (z, x, MPFR_RNDN);
409 mpfr_prec_round (z, 53, MPFR_RNDN);
410
411 if (mpfr_cmp (y, z))
412 {
413 printf ("Error in mpfr_exp for large argument\n");
414 exit (1);
415 }
416
417 /* corner cases in mpfr_exp_3 */
418 mpfr_set_prec (x, 2);
419 mpfr_set_ui (x, 1, MPFR_RNDN);
420 mpfr_set_prec (y, 2);
421 mpfr_exp_3 (y, x, MPFR_RNDN);
422
423 /* Check some little things about overflow detection */
424 set_emin (-125);
425 set_emax (128);
426 mpfr_set_prec (x, 107);
427 mpfr_set_prec (y, 107);
428 mpfr_set_str_binary (x, "0.11110000000000000000000000000000000000000000000"
429 "0000000000000000000000000000000000000000000000000000"
430 "00000000E4");
431 test_exp (y, x, MPFR_RNDN);
432 if (mpfr_cmp_str (y, "0.11000111100001100110010101111101011010010101010000"
433 "1101110111100010111001011111111000110111001011001101010"
434 "01E22", 2, MPFR_RNDN))
435 {
436 printf ("Special overflow error (1)\n");
437 mpfr_dump (y);
438 exit (1);
439 }
440
441 set_emin (emin);
442 set_emax (emax);
443
444 /* Check for overflow producing a segfault with HUGE exponent */
445 mpfr_set_ui (x, 3, MPFR_RNDN);
446 mpfr_mul_2ui (x, x, 32, MPFR_RNDN);
447 test_exp (y, x, MPFR_RNDN); /* Can't test return value: May overflow or not*/
448
449 /* Bug due to wrong approximation of (x)/log2 */
450 mpfr_set_prec (x, 163);
451
452 mpfr_set_str (x, "-4.28ac8fceeadcda06bb56359017b1c81b85b392e7", 16,
453 MPFR_RNDN);
454 mpfr_exp (x, x, MPFR_RNDN);
455 if (mpfr_cmp_str (x, "3.fffffffffffffffffffffffffffffffffffffffe8@-2",
456 16, MPFR_RNDN))
457 {
458 printf ("Error for x= -4.28ac8fceeadcda06bb56359017b1c81b85b392e7");
459 printf ("expected 3.fffffffffffffffffffffffffffffffffffffffe8@-2");
460 printf ("Got ");
461 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN); putchar ('\n');
462 }
463
464 /* bug found by Guillaume Melquiond, 13 Sep 2005 */
465 mpfr_set_prec (x, 53);
466 mpfr_set_str_binary (x, "-1E-400");
467 mpfr_exp (x, x, MPFR_RNDZ);
468 if (mpfr_cmp_ui (x, 1) == 0)
469 {
470 printf ("Error for exp(-2^(-400))\n");
471 exit (1);
472 }
473
474 mpfr_clear (x);
475 mpfr_clear (y);
476 mpfr_clear (z);
477 }
478
479 /* check sign of inexact flag */
480 static void
check_inexact(void)481 check_inexact (void)
482 {
483 mpfr_t x, y;
484 int inexact;
485
486 mpfr_init2 (x, 53);
487 mpfr_init2 (y, 53);
488
489 mpfr_set_str_binary (x,
490 "1.0000000000001001000110100100101000001101101011100101e2");
491 inexact = test_exp (y, x, MPFR_RNDN);
492 if (inexact <= 0)
493 {
494 printf ("Wrong inexact flag (Got %d instead of 1)\n", inexact);
495 exit (1);
496 }
497
498 mpfr_clear (x);
499 mpfr_clear (y);
500 }
501
502 static void
check_exp10(void)503 check_exp10(void)
504 {
505 mpfr_t x;
506 int inexact;
507
508 mpfr_init2 (x, 200);
509 mpfr_set_ui(x, 4, MPFR_RNDN);
510
511 inexact = mpfr_exp10 (x, x, MPFR_RNDN);
512 if (mpfr_cmp_ui(x, 10*10*10*10))
513 {
514 printf ("exp10: Wrong returned value\n");
515 exit (1);
516 }
517 if (inexact != 0)
518 {
519 printf ("exp10: Wrong inexact flag\n");
520 exit (1);
521 }
522
523 mpfr_clear (x);
524 }
525
526 static void
overflowed_exp0(void)527 overflowed_exp0 (void)
528 {
529 mpfr_t x, y;
530 int emax, i, inex, rnd, err = 0;
531 mpfr_exp_t old_emax;
532
533 old_emax = mpfr_get_emax ();
534
535 mpfr_init2 (x, 8);
536 mpfr_init2 (y, 8);
537
538 for (emax = -1; emax <= 0; emax++)
539 {
540 mpfr_set_ui_2exp (y, 1, emax, MPFR_RNDN);
541 mpfr_nextbelow (y);
542 set_emax (emax); /* 1 is not representable. */
543 /* and if emax < 0, 1 - eps is not representable either. */
544 for (i = -1; i <= 1; i++)
545 RND_LOOP (rnd)
546 {
547 mpfr_set_si_2exp (x, i, -512 * ABS (i), MPFR_RNDN);
548 mpfr_clear_flags ();
549 inex = mpfr_exp (x, x, (mpfr_rnd_t) rnd);
550 if ((i >= 0 || emax < 0 || rnd == MPFR_RNDN || rnd == MPFR_RNDU) &&
551 ! mpfr_overflow_p ())
552 {
553 printf ("Error in overflowed_exp0 (i = %d, rnd = %s):\n"
554 " The overflow flag is not set.\n",
555 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
556 err = 1;
557 }
558 if (rnd == MPFR_RNDZ || rnd == MPFR_RNDD)
559 {
560 if (inex >= 0)
561 {
562 printf ("Error in overflowed_exp0 (i = %d, rnd = %s):\n"
563 " The inexact value must be negative.\n",
564 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
565 err = 1;
566 }
567 if (! mpfr_equal_p (x, y))
568 {
569 printf ("Error in overflowed_exp0 (i = %d, rnd = %s):\n"
570 " Got ", i,
571 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
572 mpfr_dump (x);
573 printf (" instead of 0.11111111E%d.\n", emax);
574 err = 1;
575 }
576 }
577 else if (rnd != MPFR_RNDF)
578 {
579 if (inex <= 0)
580 {
581 printf ("Error in overflowed_exp0 (i = %d, rnd = %s):\n"
582 " The inexact value must be positive.\n",
583 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
584 err = 1;
585 }
586 if (! (mpfr_inf_p (x) && MPFR_IS_POS (x)))
587 {
588 printf ("Error in overflowed_exp0 (i = %d, rnd = %s):\n"
589 " Got ", i,
590 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
591 mpfr_dump (x);
592 printf (" instead of +Inf.\n");
593 err = 1;
594 }
595 }
596 }
597 set_emax (old_emax);
598 }
599
600 if (err)
601 exit (1);
602 mpfr_clear (x);
603 mpfr_clear (y);
604 }
605
606 /* This bug occurs in mpfr_exp_2 on a Linux-64 machine, r5475. */
607 static void
bug20080731(void)608 bug20080731 (void)
609 {
610 mpfr_exp_t emin;
611 mpfr_t x, y1, y2;
612 mpfr_prec_t prec = 64;
613
614 emin = mpfr_get_emin ();
615 set_emin (MPFR_EMIN_MIN);
616
617 mpfr_init2 (x, 200);
618 mpfr_set_str (x, "-2.c5c85fdf473de6af278ece700fcbdabd03cd0cb9ca62d8b62c@7",
619 16, MPFR_RNDN);
620
621 /* exp(x) is just below 0xf.fffffffffffffffp-1073741828 */
622
623 mpfr_init2 (y1, prec);
624 mpfr_exp (y1, x, MPFR_RNDU);
625
626 /* Compute the result with a higher internal precision. */
627 mpfr_init2 (y2, 300);
628 mpfr_exp (y2, x, MPFR_RNDU);
629 mpfr_prec_round (y2, prec, MPFR_RNDU);
630
631 if (mpfr_cmp0 (y1, y2) != 0)
632 {
633 printf ("Error in bug20080731\nExpected ");
634 mpfr_dump (y2);
635 printf ("\nGot ");
636 mpfr_dump (y1);
637 printf ("\n");
638 exit (1);
639 }
640
641 mpfr_clears (x, y1, y2, (mpfr_ptr) 0);
642 set_emin (emin);
643 }
644
645 /* Emulate mpfr_exp with mpfr_exp_3 in the general case. */
646 static int
exp_3(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)647 exp_3 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
648 {
649 int inexact;
650
651 inexact = mpfr_exp_3 (y, x, rnd_mode);
652 return mpfr_check_range (y, inexact, rnd_mode);
653 }
654
655 static void
underflow_up(int extended_emin)656 underflow_up (int extended_emin)
657 {
658 mpfr_t minpos, x, y, t, t2;
659 int precx, precy;
660 int inex;
661 int rnd;
662 int e3;
663 int i, j;
664
665 mpfr_init2 (minpos, 2);
666 mpfr_set_ui (minpos, 0, MPFR_RNDN);
667 mpfr_nextabove (minpos);
668
669 /* Let's test values near the underflow boundary.
670 *
671 * Minimum representable positive number: minpos = 2^(emin - 1).
672 * Let's choose an MPFR number x = log(minpos) + eps, with |eps| small
673 * (note: eps cannot be 0, and cannot be a rational number either).
674 * Then exp(x) = minpos * exp(eps) ~= minpos * (1 + eps + eps^2).
675 * We will compute y = rnd(exp(x)) in some rounding mode, precision p.
676 * 1. If eps > 0, then in any rounding mode:
677 * rnd(exp(x)) >= minpos and no underflow.
678 * So, let's take x1 = rndu(log(minpos)) in some precision.
679 * 2. If eps < 0, then exp(x) < minpos and the result will be either 0
680 * or minpos. An underflow always occurs in MPFR_RNDZ and MPFR_RNDD,
681 * but not necessarily in MPFR_RNDN and MPFR_RNDU (this is underflow
682 * after rounding in an unbounded exponent range). If -a < eps < -b,
683 * minpos * (1 - a) < exp(x) < minpos * (1 - b + b^2).
684 * - If eps > -2^(-p), no underflow in MPFR_RNDU.
685 * - If eps > -2^(-p-1), no underflow in MPFR_RNDN.
686 * - If eps < - (2^(-p-1) + 2^(-2p-1)), underflow in MPFR_RNDN.
687 * - If eps < - (2^(-p) + 2^(-2p+1)), underflow in MPFR_RNDU.
688 * - In MPFR_RNDN, result is minpos iff exp(eps) > 1/2, i.e.
689 * - log(2) < eps < ...
690 *
691 * Moreover, since precy < MPFR_EXP_THRESHOLD (to avoid tests that take
692 * too much time), mpfr_exp() always selects mpfr_exp_2(); so, we need
693 * to test mpfr_exp_3() too. This will be done via the e3 variable:
694 * e3 = 0: mpfr_exp(), thus mpfr_exp_2().
695 * e3 = 1: mpfr_exp_3(), via the exp_3() wrapper.
696 * i.e.: inex = e3 ? exp_3 (y, x, rnd) : mpfr_exp (y, x, rnd);
697 */
698
699 /* Case eps > 0. In revision 5461 (trunk) on a 64-bit Linux machine:
700 * Incorrect flags in underflow_up, eps > 0, MPFR_RNDN and extended emin
701 * for precx = 96, precy = 16, mpfr_exp_3
702 * Got 9 instead of 8.
703 * Note: testing this case in several precisions for x and y introduces
704 * some useful random. Indeed, the bug is not always triggered.
705 * Fixed in r5469.
706 */
707 for (precx = 16; precx <= 128; precx += 16)
708 {
709 mpfr_init2 (x, precx);
710 mpfr_log (x, minpos, MPFR_RNDU);
711 for (precy = 16; precy <= 128; precy += 16)
712 {
713 mpfr_init2 (y, precy);
714
715 for (e3 = 0; e3 <= 1; e3++)
716 {
717 RND_LOOP (rnd)
718 {
719 int err = 0;
720
721 mpfr_clear_flags ();
722 inex = e3 ? exp_3 (y, x, (mpfr_rnd_t) rnd)
723 : mpfr_exp (y, x, (mpfr_rnd_t) rnd);
724 /* for MPFR_RNDF, the inexact flag is undefined */
725 if (__gmpfr_flags != MPFR_FLAGS_INEXACT && rnd != MPFR_RNDF)
726 {
727 printf ("Incorrect flags in underflow_up, eps > 0, %s",
728 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
729 if (extended_emin)
730 printf (" and extended emin");
731 printf ("\nfor precx = %d, precy = %d, %s\n",
732 precx, precy, e3 ? "mpfr_exp_3" : "mpfr_exp");
733 printf ("x="); mpfr_dump (x);
734 printf ("y="); mpfr_dump (y);
735 printf ("Got %u instead of %u.\n",
736 (unsigned int) __gmpfr_flags,
737 (unsigned int) MPFR_FLAGS_INEXACT);
738 err = 1;
739 }
740 if (mpfr_cmp0 (y, minpos) < 0)
741 {
742 printf ("Incorrect result in underflow_up, eps > 0, %s",
743 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
744 if (extended_emin)
745 printf (" and extended emin");
746 printf ("\nfor precx = %d, precy = %d, %s\n",
747 precx, precy, e3 ? "mpfr_exp_3" : "mpfr_exp");
748 mpfr_dump (y);
749 err = 1;
750 }
751 /* for MPFR_RNDF, the ternary value is undefined */
752 if (rnd != MPFR_RNDF)
753 MPFR_ASSERTN (inex != 0);
754 if (rnd == MPFR_RNDD || rnd == MPFR_RNDZ)
755 MPFR_ASSERTN (inex < 0);
756 if (rnd == MPFR_RNDU)
757 MPFR_ASSERTN (inex > 0);
758 if (err)
759 exit (1);
760 }
761 }
762
763 mpfr_clear (y);
764 }
765 mpfr_clear (x);
766 }
767
768 /* Case - log(2) < eps < 0 in MPFR_RNDN, starting with small-precision x;
769 * only check the result and the ternary value.
770 * Previous to r5453 (trunk), on 32-bit and 64-bit machines, this fails
771 * for precx = 65 and precy = 16, e.g.:
772 * exp_2.c:264: assertion failed: ...
773 * because mpfr_sub (r, x, r, MPFR_RNDU); yields a null value. This is
774 * fixed in r5453 by going to next Ziv's iteration.
775 */
776 for (precx = sizeof(mpfr_exp_t) * CHAR_BIT + 1; precx <= 81; precx += 8)
777 {
778 mpfr_init2 (x, precx);
779 mpfr_log (x, minpos, MPFR_RNDD); /* |ulp| <= 1/2 */
780 for (precy = 16; precy <= 128; precy += 16)
781 {
782 mpfr_init2 (y, precy);
783 inex = mpfr_exp (y, x, MPFR_RNDN);
784 if (inex <= 0 || mpfr_cmp0 (y, minpos) != 0)
785 {
786 printf ("Error in underflow_up, - log(2) < eps < 0");
787 if (extended_emin)
788 printf (" and extended emin");
789 printf (" for prec = %d\nExpected ", precy);
790 mpfr_out_str (stdout, 16, 0, minpos, MPFR_RNDN);
791 printf (" (minimum positive MPFR number) and inex > 0\nGot ");
792 mpfr_out_str (stdout, 16, 0, y, MPFR_RNDN);
793 printf ("\nwith inex = %d\n", inex);
794 exit (1);
795 }
796 mpfr_clear (y);
797 }
798 mpfr_clear (x);
799 }
800
801 /* Cases eps ~ -2^(-p) and eps ~ -2^(-p-1). More precisely,
802 * _ for j = 0, eps > -2^(-(p+i)),
803 * _ for j = 1, eps < - (2^(-(p+i)) + 2^(1-2(p+i))),
804 * where i = 0 or 1.
805 */
806 mpfr_inits2 (2, t, t2, (mpfr_ptr) 0);
807 for (precy = 16; precy <= 128; precy += 16)
808 {
809 mpfr_set_ui_2exp (t, 1, - precy, MPFR_RNDN); /* 2^(-p) */
810 mpfr_set_ui_2exp (t2, 1, 1 - 2 * precy, MPFR_RNDN); /* 2^(-2p+1) */
811 precx = sizeof(mpfr_exp_t) * CHAR_BIT + 2 * precy + 8;
812 mpfr_init2 (x, precx);
813 mpfr_init2 (y, precy);
814 for (i = 0; i <= 1; i++)
815 {
816 for (j = 0; j <= 1; j++)
817 {
818 if (j == 0)
819 {
820 /* Case eps > -2^(-(p+i)). */
821 mpfr_log (x, minpos, MPFR_RNDU);
822 }
823 else /* j == 1 */
824 {
825 /* Case eps < - (2^(-(p+i)) + 2^(1-2(p+i))). */
826 mpfr_log (x, minpos, MPFR_RNDD);
827 inex = mpfr_sub (x, x, t2, MPFR_RNDN);
828 MPFR_ASSERTN (inex == 0);
829 }
830 inex = mpfr_sub (x, x, t, MPFR_RNDN);
831 MPFR_ASSERTN (inex == 0);
832
833 RND_LOOP (rnd)
834 for (e3 = 0; e3 <= 1; e3++)
835 {
836 int err = 0;
837 unsigned int flags;
838
839 flags = MPFR_FLAGS_INEXACT |
840 (((rnd == MPFR_RNDU || rnd == MPFR_RNDA)
841 && (i == 1 || j == 0)) ||
842 (rnd == MPFR_RNDN && (i == 1 && j == 0)) ?
843 0 : MPFR_FLAGS_UNDERFLOW);
844 mpfr_clear_flags ();
845 inex = e3 ? exp_3 (y, x, (mpfr_rnd_t) rnd)
846 : mpfr_exp (y, x, (mpfr_rnd_t) rnd);
847 if (__gmpfr_flags != flags)
848 {
849 printf ("Incorrect flags in underflow_up, %s",
850 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
851 if (extended_emin)
852 printf (" and extended emin");
853 printf ("\nfor precx = %d, precy = %d, ",
854 precx, precy);
855 if (j == 0)
856 printf ("eps >~ -2^(-%d)", precy + i);
857 else
858 printf ("eps <~ - (2^(-%d) + 2^(%d))",
859 precy + i, 1 - 2 * (precy + i));
860 printf (", %s\n", e3 ? "mpfr_exp_3" : "mpfr_exp");
861 printf ("Got %u instead of %u.\n",
862 (unsigned int) __gmpfr_flags, flags);
863 err = 1;
864 }
865 if (rnd == MPFR_RNDF)
866 continue; /* the test below makes no sense, since RNDF
867 does not give a deterministic result */
868 if (rnd == MPFR_RNDU || rnd == MPFR_RNDA || rnd == MPFR_RNDN ?
869 mpfr_cmp0 (y, minpos) != 0 : MPFR_NOTZERO (y))
870 {
871 printf ("Incorrect result in underflow_up, %s",
872 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
873 if (extended_emin)
874 printf (" and extended emin");
875 printf ("\nfor precx = %d, precy = %d, ",
876 precx, precy);
877 if (j == 0)
878 printf ("eps >~ -2^(-%d)", precy + i);
879 else
880 printf ("eps <~ - (2^(-%d) + 2^(%d))",
881 precy + i, 1 - 2 * (precy + i));
882 printf (", %s\n", e3 ? "mpfr_exp_3" : "mpfr_exp");
883 mpfr_dump (y);
884 err = 1;
885 }
886 if (err)
887 exit (1);
888 } /* for (e3 ...) */
889 } /* for (j ...) */
890 mpfr_div_2si (t, t, 1, MPFR_RNDN);
891 mpfr_div_2si (t2, t2, 2, MPFR_RNDN);
892 } /* for (i ...) */
893 mpfr_clears (x, y, (mpfr_ptr) 0);
894 } /* for (precy ...) */
895 mpfr_clears (t, t2, (mpfr_ptr) 0);
896
897 /* Case exp(eps) ~= 1/2, i.e. eps ~= - log(2).
898 * We choose x0 and x1 with high enough precision such that:
899 * x0 = rndd(rndd(log(minpos)) - rndu(log(2)))
900 * x1 = rndu(rndu(log(minpos)) - rndd(log(2)))
901 * In revision 5507 (trunk) on a 64-bit Linux machine, this fails:
902 * Error in underflow_up, eps >~ - log(2) and extended emin
903 * for precy = 16, mpfr_exp
904 * Expected 1.0@-1152921504606846976 (minimum positive MPFR number),
905 * inex > 0 and flags = 9
906 * Got 0
907 * with inex = -1 and flags = 9
908 * due to a double-rounding problem in mpfr_mul_2si when rescaling
909 * the result.
910 */
911 mpfr_inits2 (sizeof(mpfr_exp_t) * CHAR_BIT + 64, x, t, (mpfr_ptr) 0);
912 for (i = 0; i <= 1; i++)
913 {
914 mpfr_log (x, minpos, i ? MPFR_RNDU : MPFR_RNDD);
915 mpfr_const_log2 (t, i ? MPFR_RNDD : MPFR_RNDU);
916 mpfr_sub (x, x, t, i ? MPFR_RNDU : MPFR_RNDD);
917 for (precy = 16; precy <= 128; precy += 16)
918 {
919 mpfr_init2 (y, precy);
920 for (e3 = 0; e3 <= 1; e3++)
921 {
922 unsigned int flags, uflags =
923 MPFR_FLAGS_INEXACT | MPFR_FLAGS_UNDERFLOW;
924
925 mpfr_clear_flags ();
926 inex = e3 ? exp_3 (y, x, MPFR_RNDN) : mpfr_exp (y, x, MPFR_RNDN);
927 flags = __gmpfr_flags;
928 if (flags != uflags ||
929 (i ? (inex <= 0 || mpfr_cmp0 (y, minpos) != 0)
930 : (inex >= 0 || MPFR_NOTZERO (y))))
931 {
932 printf ("Error in underflow_up, eps %c~ - log(2)",
933 i ? '>' : '<');
934 if (extended_emin)
935 printf (" and extended emin");
936 printf ("\nfor precy = %d, %s\nExpected ", precy,
937 e3 ? "mpfr_exp_3" : "mpfr_exp");
938 if (i)
939 {
940 mpfr_out_str (stdout, 16, 0, minpos, MPFR_RNDN);
941 printf (" (minimum positive MPFR number),\ninex >");
942 }
943 else
944 {
945 printf ("+0, inex <");
946 }
947 printf (" 0 and flags = %u\nGot ", uflags);
948 mpfr_out_str (stdout, 16, 0, y, MPFR_RNDN);
949 printf ("\nwith inex = %d and flags = %u\n", inex, flags);
950 exit (1);
951 }
952 }
953 mpfr_clear (y);
954 }
955 }
956 mpfr_clears (x, t, (mpfr_ptr) 0);
957
958 mpfr_clear (minpos);
959 }
960
961 static void
underflow(void)962 underflow (void)
963 {
964 mpfr_exp_t emin;
965
966 underflow_up (0);
967
968 emin = mpfr_get_emin ();
969 set_emin (MPFR_EMIN_MIN);
970 if (mpfr_get_emin () != emin)
971 {
972 underflow_up (1);
973 set_emin (emin);
974 }
975 }
976
977 /* bug found with GMP_CHECK_RANDOMIZE=1514290185 */
978 static void
bug20171223(void)979 bug20171223 (void)
980 {
981 mpfr_t x, y;
982 int inex;
983
984 mpfr_init2 (x, 372);
985 mpfr_init2 (y, 2);
986 mpfr_set_str (x, "-6.9314716128384587678466323621915206417385796077947874471662159283492445979241549847386366371775938082803907383582e-01", 10, MPFR_RNDN);
987 /* exp(x) = 0.500000009638..., should be rounded to 0.5 */
988 inex = mpfr_exp (y, x, MPFR_RNDD);
989 MPFR_ASSERTN(mpfr_cmp_ui_2exp (y, 1, -1) == 0);
990 MPFR_ASSERTN(inex < 0);
991 mpfr_clear (x);
992 mpfr_clear (y);
993 }
994
995 int
main(int argc,char * argv[])996 main (int argc, char *argv[])
997 {
998 tests_start_mpfr ();
999
1000 if (argc > 1)
1001 check_large ();
1002
1003 bug20171223 ();
1004
1005 check_inexact ();
1006 check_special ();
1007
1008 test_generic (MPFR_PREC_MIN, 100, 100);
1009
1010 compare_exp2_exp3 (20, 1000);
1011 check_worst_cases();
1012 check3("0.0", MPFR_RNDU, "1.0");
1013 check3("-1e-170", MPFR_RNDU, "1.0");
1014 check3("-1e-170", MPFR_RNDN, "1.0");
1015 check3("-8.88024741073346941839e-17", MPFR_RNDU, "1.0");
1016 check3("8.70772839244701057915e-01", MPFR_RNDN, "2.38875626491680437269");
1017 check3("1.0", MPFR_RNDN, "2.71828182845904509080");
1018 check3("-3.42135637628104173534e-07", MPFR_RNDZ, "0.999999657864420798958");
1019 /* worst case for argument reduction, very near from 5*log(2),
1020 thanks to Jean-Michel Muller */
1021 check3("3.4657359027997265421", MPFR_RNDN, "32.0");
1022 check3("3.4657359027997265421", MPFR_RNDU, "32.0");
1023 check3("3.4657359027997265421", MPFR_RNDD, "31.999999999999996447");
1024 check3("2.26523754332090625496e+01", MPFR_RNDD, "6.8833785261699581146e9");
1025 check3("1.31478962104089092122e+01", MPFR_RNDZ, "5.12930793917860137299e+05");
1026 check3("4.25637507920002378103e-01", MPFR_RNDU, "1.53056585656161181497e+00");
1027 check3("6.26551618962329307459e-16", MPFR_RNDU, "1.00000000000000066613e+00");
1028 check3("-3.35589513871216568383e-03",MPFR_RNDD, "9.96649729583626853291e-01");
1029 check3("1.95151388850007272424e+01", MPFR_RNDU, "2.98756340674767792225e+08");
1030 check3("2.45045953503350730784e+01", MPFR_RNDN, "4.38743344916128387451e+10");
1031 check3("2.58165606081678085104e+01", MPFR_RNDD, "1.62925781879432281494e+11");
1032 check3("-2.36539020084338638128e+01",MPFR_RNDZ, "5.33630792749924762447e-11");
1033 check3("2.39211946135858077866e+01", MPFR_RNDU, "2.44817704330214385986e+10");
1034 check3("-2.78190533055889162029e+01",MPFR_RNDZ, "8.2858803483596879512e-13");
1035 check3("2.64028186174889789584e+01", MPFR_RNDD, "2.9281844652878973388e11");
1036 check3("2.92086338843268329413e+01", MPFR_RNDZ, "4.8433797301907177734e12");
1037 check3("-2.46355324071459982349e+01",MPFR_RNDZ, "1.9995129297760994791e-11");
1038 check3("-2.23509444608605427618e+01",MPFR_RNDZ, "1.9638492867489702307e-10");
1039 check3("-2.41175390197331687148e+01",MPFR_RNDD, "3.3564940885530624592e-11");
1040 check3("2.46363885231578088053e+01", MPFR_RNDU, "5.0055014282693267822e10");
1041 check3("111.1263531080090984914932050742208957672119140625", MPFR_RNDN, "1.8262572323517295459e48");
1042 check3("-3.56196340354684821250e+02",MPFR_RNDN, "2.0225297096141478156e-155");
1043 check3("6.59678273772710895173e+02", MPFR_RNDU, "3.1234469273830195529e286");
1044 check3("5.13772529701934331570e+02", MPFR_RNDD, "1.3445427121297197752e223");
1045 check3("3.57430211008718345056e+02", MPFR_RNDD, "1.6981197246857298443e155");
1046 check3("3.82001814471465536371e+02", MPFR_RNDU, "7.9667300591087367805e165");
1047 check3("5.92396038219384422518e+02", MPFR_RNDD, "1.880747529554661989e257");
1048 check3("-5.02678550462488090034e+02",MPFR_RNDU, "4.8919201895446217839e-219");
1049 check3("5.30015757134837031117e+02", MPFR_RNDD, "1.5237672861171573939e230");
1050 check3("5.16239362447650933063e+02", MPFR_RNDZ, "1.5845518406744492105e224");
1051 check3("6.00812634798592370977e-01", MPFR_RNDN, "1.823600119339019443");
1052 check_exp10 ();
1053
1054 bug20080731 ();
1055
1056 overflowed_exp0 ();
1057 underflow ();
1058
1059 data_check ("data/exp", mpfr_exp, "mpfr_exp");
1060 bad_cases (mpfr_exp, mpfr_log, "mpfr_exp", 0, -256, 255, 4, 128, 800, 50);
1061
1062 tests_end_mpfr ();
1063 return 0;
1064 }
1065