1 /* Test file for mpfr_exp10.
2
3 Copyright 2007-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 #define TEST_FUNCTION mpfr_exp10
26 #define TEST_RANDOM_EMIN -36
27 #define TEST_RANDOM_EMAX 36
28 #include "tgeneric.c"
29
30 static void
31 special_overflow (void)
32 {
33 mpfr_t x, y;
34 int inex;
35 mpfr_exp_t emin, emax;
36
37 emin = mpfr_get_emin ();
38 emax = mpfr_get_emax ();
39
40 set_emin (-125);
41 set_emax (128);
42
43 mpfr_init2 (x, 24);
44 mpfr_init2 (y, 24);
45
46 mpfr_set_str_binary (x, "0.101100100000000000110100E15");
47 inex = mpfr_exp10 (y, x, MPFR_RNDN);
48 if (!mpfr_inf_p (y) || inex <= 0)
49 {
50 printf ("Overflow error.\n");
51 mpfr_dump (y);
52 printf ("inex = %d\n", inex);
53 exit (1);
54 }
55
56 mpfr_clear (y);
57 mpfr_clear (x);
58 set_emin (emin);
59 set_emax (emax);
60 }
61
62 static void
63 emax_m_eps (void)
64 {
65 if (mpfr_get_emax () <= LONG_MAX)
66 {
67 mpfr_t x, y;
68 int inex, ov;
69
70 mpfr_init2 (x, sizeof(mpfr_exp_t) * CHAR_BIT * 4);
71 mpfr_init2 (y, 8);
72 mpfr_set_si (x, mpfr_get_emax (), MPFR_RNDN);
73
74 mpfr_clear_flags ();
75 inex = mpfr_exp10 (y, x, MPFR_RNDN);
76 ov = mpfr_overflow_p ();
77 if (!ov || !mpfr_inf_p (y) || inex <= 0)
78 {
79 printf ("Overflow error for x = emax, MPFR_RNDN.\n");
80 mpfr_dump (y);
81 printf ("inex = %d, %soverflow\n", inex, ov ? "" : "no ");
82 exit (1);
83 }
84
85 mpfr_clear (x);
86 mpfr_clear (y);
87 }
88 }
89
90 static void
91 exp_range (void)
92 {
93 mpfr_t x;
94 mpfr_exp_t emin;
95
96 emin = mpfr_get_emin ();
97 set_emin (3);
98 mpfr_init2 (x, 16);
99 mpfr_set_ui (x, 4, MPFR_RNDN);
100 mpfr_exp10 (x, x, MPFR_RNDN);
101 set_emin (emin);
102 if (mpfr_nan_p (x) || mpfr_cmp_ui (x, 10000) != 0)
103 {
104 printf ("Error in mpfr_exp10 for x = 4, with emin = 3\n");
105 printf ("Expected 10000, got ");
106 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
107 printf ("\n");
108 exit (1);
109 }
110 mpfr_clear (x);
111 }
112
113 static void
114 overfl_exp10_0 (void)
115 {
116 mpfr_t x, y;
117 int emax, i, inex, rnd, err = 0;
118 mpfr_exp_t old_emax;
119
120 old_emax = mpfr_get_emax ();
121
122 mpfr_init2 (x, 8);
123 mpfr_init2 (y, 8);
124
125 for (emax = -1; emax <= 0; emax++)
126 {
127 mpfr_set_ui_2exp (y, 1, emax, MPFR_RNDN);
128 mpfr_nextbelow (y);
129 set_emax (emax); /* 1 is not representable. */
130 /* and if emax < 0, 1 - eps is not representable either. */
131 for (i = -1; i <= 1; i++)
132 RND_LOOP (rnd)
133 {
134 mpfr_set_si_2exp (x, i, -512 * ABS (i), MPFR_RNDN);
135 mpfr_clear_flags ();
136 inex = mpfr_exp10 (x, x, (mpfr_rnd_t) rnd);
137 if ((i >= 0 || emax < 0 || rnd == MPFR_RNDN || rnd == MPFR_RNDU) &&
138 ! mpfr_overflow_p ())
139 {
140 printf ("Error in overfl_exp10_0 (i = %d, rnd = %s):\n"
141 " The overflow flag is not set.\n",
142 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
143 err = 1;
144 }
145 if (rnd == MPFR_RNDZ || rnd == MPFR_RNDD)
146 {
147 if (inex >= 0)
148 {
149 printf ("Error in overfl_exp10_0 (i = %d, rnd = %s):\n"
150 " The inexact value must be negative.\n",
151 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
152 err = 1;
153 }
154 if (! mpfr_equal_p (x, y))
155 {
156 printf ("Error in overfl_exp10_0 (i = %d, rnd = %s):\n"
157 " Got ", i,
158 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
159 mpfr_dump (x);
160 printf (" instead of 0.11111111E%d.\n", emax);
161 err = 1;
162 }
163 }
164 else if (rnd != MPFR_RNDF)
165 {
166 if (inex <= 0)
167 {
168 printf ("Error in overfl_exp10_0 (i = %d, rnd = %s):\n"
169 " The inexact value must be positive.\n",
170 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
171 err = 1;
172 }
173 if (! (mpfr_inf_p (x) && MPFR_IS_POS (x)))
174 {
175 printf ("Error in overfl_exp10_0 (i = %d, rnd = %s):\n"
176 " Got ", i,
177 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
178 mpfr_dump (x);
179 printf (" instead of +Inf.\n");
180 err = 1;
181 }
182 }
183 }
184 set_emax (old_emax);
185 }
186
187 if (err)
188 exit (1);
189 mpfr_clear (x);
190 mpfr_clear (y);
191 }
192
193 /* Bug in mpfr_pow_general found by ofuf_thresholds (on 2023-02-13 for
194 a 32-bit exponent, changed on 2023-03-06 for a 64-bit exponent too),
195 fixed in commit b62966df913f73f08b3c5252e1d0c702bc20442f.
196 With a 32-bit exponent, failure for i=0.
197 expected 0.1111E1073741823
198 got @Inf@
199 expected flags = inexact (8)
200 got flags = overflow inexact (10)
201 With a 64-bit exponent, failure for i=1.
202 expected 0.11111111111111111111111E4611686018427387903
203 got @Inf@
204 expected flags = inexact (8)
205 got flags = overflow inexact (10)
206 Note: ofuf_thresholds was added to the master branch, but for the
207 time being, there are issues with these tests.
208 */
209 static void
210 bug20230213 (void)
211 {
212 const char *s[2] = {
213 "0x1.34413504b3ccdbd5dd8p+28",
214 "0x1.34413509f79fef2c4e0dd14a7ae0ecfbacdbp+60"
215 };
216 mpfr_t x1, x2, y1, y2;
217 mpfr_prec_t px[2] = { 74, 147 };
218 mpfr_prec_t py[2] = { 4, 23 };
219 mpfr_exp_t old_emax, emax;
220 mpfr_flags_t flags1, flags2;
221 int i;
222
223 old_emax = mpfr_get_emax ();
224
225 for (i = 0; i < 2; i++)
226 {
227 if (i != 0)
228 set_emax (MPFR_EMAX_MAX);
229
230 emax = mpfr_get_emax ();
231
232 mpfr_inits2 (px[i], x1, x2, (mpfr_ptr) 0);
233 mpfr_inits2 (py[i], y1, y2, (mpfr_ptr) 0);
234
235 mpfr_setmax (y1, emax);
236 mpfr_log10 (x1, y1, MPFR_RNDD);
237 mpfr_set_str (x2, s[i], 0, MPFR_RNDN);
238 /* For i == 0, emax == 2^30, so that the value can be checked.
239 For i != 0, check the value for the case emax == 2^62.
240 The "0UL" ensures that the shifts are valid. */
241 if (i == 0 || (((0UL + MPFR_EMAX_MAX) >> 31) >> 30) == 1)
242 {
243 /* printf ("Checking x1 for i=%d\n", i); */
244 MPFR_ASSERTN (mpfr_equal_p (x1, x2));
245 }
246
247 /* Let MAXF be the maximum finite value (y1 above).
248 Since x1 < log10(MAXF), one should have exp10(x1) < MAXF, and
249 therefore, y2 = RU(exp10(x1)) <= RU(MAXF) = MAXF (no overflow). */
250 flags1 = MPFR_FLAGS_INEXACT;
251 mpfr_clear_flags ();
252 mpfr_exp10 (y2, x1, MPFR_RNDU);
253 flags2 = __gmpfr_flags;
254
255 if (! (mpfr_lessequal_p (y2, y1) && flags2 == flags1))
256 {
257 printf ("Error in bug20230213 for i=%d\n", i);
258 printf ("emax = %" MPFR_EXP_FSPEC "d\n", (mpfr_eexp_t) emax);
259 printf ("expected "); mpfr_dump (y1);
260 printf ("got "); mpfr_dump (y2);
261 printf ("expected flags =");
262 flags_out (flags1);
263 printf ("got flags =");
264 flags_out (flags2);
265 exit (1);
266 }
267
268 mpfr_clears (x1, x2, y1, y2, (mpfr_ptr) 0);
269 }
270
271 set_emax (old_emax);
272 }
273
274 /* Bug in mpfr_pow_general in precision 1 in the particular case of
275 rounding to nearest, z * 2^k = 2^(emin - 2) and real result larger
276 than this value; fixed in ff5012b61d5e5fee5156c57b8aa8fc1739c2a771
277 (which is simplified in 4f5de980be290687ac1409aa02873e9e0dd1a030);
278 initially found by ofuf_thresholds (though the test was incorrect).
279 With a 32-bit exponent, failure for i=0.
280 With a 64-bit exponent, failure for i=1.
281 The result was correct, but the underflow flag was missing.
282 Note: ofuf_thresholds was added to the master branch, but for the
283 time being, there are issues with these tests.
284 */
285 static void
286 bug20230427 (void)
287 {
288 const char *s[2] = {
289 "-0.1001101000100000100110101000011E29",
290 "-0.100110100010000010011010100001001111101111001111111101111001101E61"
291 };
292 mpfr_t x, y, z, t1, t2;
293 mpfr_exp_t old_emin;
294 mpfr_flags_t flags, ex_flags;
295 int i, inex;
296
297 old_emin = mpfr_get_emin ();
298
299 mpfr_init2 (x, 63);
300 mpfr_inits2 (1, y, z, (mpfr_ptr) 0);
301 mpfr_inits2 (128, t1, t2, (mpfr_ptr) 0);
302
303 for (i = 0; i < 2; i++)
304 {
305 if (i == 0)
306 {
307 /* Basic check: the default emin should be -2^30 (exactly). */
308 if (mpfr_get_emin () != -1073741823)
309 abort ();
310 }
311 else
312 {
313 /* This test assumes that MPFR_EMIN_MIN = -2^62 (exactly).
314 The "0UL" ensures that the shifts are valid. */
315 if ((((0UL - MPFR_EMIN_MIN) >> 31) >> 30) != 1)
316 break;
317
318 set_emin (MPFR_EMIN_MIN);
319 }
320
321 mpfr_set_str_binary (x, s[i]);
322
323 /* We will test 10^x rounded to nearest in precision 1.
324 Check that 2^(emin - 2) < 10^x < (3/2) * 2^(emin - 2).
325 This is approximate, but by outputting the values, one can check
326 that one is not too close to the boundaries:
327 emin - 2 = -4611686018427387905
328 log2(10^x) ~= -4611686018427387904.598
329 emin - 2 + log2(3/2) ~= -4611686018427387904.415
330 Thus the result should be the smallest positive number 2^(emin - 1)
331 because 10^x is closer to this number than to 0, the midpoint being
332 2^(emin - 2). And there should be an underflow in precision 1 because
333 the result rounded to nearest in an unbounded exponent range should
334 have been 2^(emin - 2), the midpoint being (3/2) * 2^(emin - 2).
335 */
336 mpfr_set_ui (t1, 10, MPFR_RNDN);
337 mpfr_log2 (t2, t1, MPFR_RNDN);
338 mpfr_mul (t1, t2, x, MPFR_RNDN);
339 inex = mpfr_set_exp_t (t2, mpfr_get_emin () - 2, MPFR_RNDN);
340 MPFR_ASSERTN (inex == 0);
341 MPFR_ASSERTN (mpfr_greater_p (t1, t2)); /* log2(10^x) > emin - 2 */
342 inex = mpfr_sub (t1, t1, t2, MPFR_RNDN);
343 MPFR_ASSERTN (inex == 0);
344 mpfr_set_ui (t2, 3, MPFR_RNDN);
345 mpfr_log2 (t2, t2, MPFR_RNDN);
346 mpfr_sub_ui (t2, t2, 1, MPFR_RNDN); /* log2(3/2) */
347 MPFR_ASSERTN (mpfr_less_p (t1, t2));
348
349 mpfr_clear_flags ();
350 mpfr_exp10 (y, x, MPFR_RNDN);
351 flags = __gmpfr_flags;
352 ex_flags = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
353
354 mpfr_setmin (z, mpfr_get_emin ()); /* z = 0.1@emin */
355 if (! (mpfr_equal_p (y, z) && flags == ex_flags))
356 {
357 printf ("Error in bug20230427 for i=%d\n", i);
358 printf ("expected "); mpfr_dump (z);
359 printf ("got "); mpfr_dump (y);
360 printf ("emin = %" MPFR_EXP_FSPEC "d\n",
361 (mpfr_eexp_t) mpfr_get_emin ());
362 printf ("expected flags =");
363 flags_out (ex_flags);
364 printf ("got flags =");
365 flags_out (flags);
366 exit (1);
367 }
368 }
369
370 mpfr_clears (x, y, z, t1, t2, (mpfr_ptr) 0);
371 set_emin (old_emin);
372 }
373
374 int
375 main (int argc, char *argv[])
376 {
377 mpfr_t x, y;
378 mpfr_exp_t emin, emax;
379 int inex, ov;
380
381 tests_start_mpfr ();
382
383 bug20230213 ();
384 bug20230427 ();
385
386 special_overflow ();
387 emax_m_eps ();
388 exp_range ();
389
390 mpfr_init (x);
391 mpfr_init (y);
392
393 mpfr_set_ui (x, 4, MPFR_RNDN);
394 mpfr_exp10 (y, x, MPFR_RNDN);
395 if (mpfr_cmp_ui (y, 10000) != 0)
396 {
397 printf ("Error for 10^4, MPFR_RNDN\n");
398 exit (1);
399 }
400 mpfr_exp10 (y, x, MPFR_RNDD);
401 if (mpfr_cmp_ui (y, 10000) != 0)
402 {
403 printf ("Error for 10^4, MPFR_RNDD\n");
404 exit (1);
405 }
406 mpfr_exp10 (y, x, MPFR_RNDU);
407 if (mpfr_cmp_ui (y, 10000) != 0)
408 {
409 printf ("Error for 10^4, MPFR_RNDU\n");
410 exit (1);
411 }
412
413 mpfr_set_prec (x, 10);
414 mpfr_set_prec (y, 10);
415 /* save emin */
416 emin = mpfr_get_emin ();
417 set_emin (-11);
418 mpfr_set_si (x, -4, MPFR_RNDN);
419 mpfr_exp10 (y, x, MPFR_RNDN);
420 if (MPFR_NOTZERO (y) || MPFR_IS_NEG (y))
421 {
422 printf ("Error for emin = -11, x = -4, RNDN\n");
423 printf ("Expected +0\n");
424 printf ("Got "); mpfr_dump (y);
425 exit (1);
426 }
427 /* restore emin */
428 set_emin (emin);
429
430 /* save emax */
431 emax = mpfr_get_emax ();
432 set_emax (13);
433 mpfr_set_ui (x, 4, MPFR_RNDN);
434 mpfr_exp10 (y, x, MPFR_RNDN);
435 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
436 {
437 printf ("Error for emax = 13, x = 4, RNDN\n");
438 printf ("Expected +inf\n");
439 printf ("Got "); mpfr_dump (y);
440 exit (1);
441 }
442 /* restore emax */
443 set_emax (emax);
444
445 MPFR_SET_INF (x);
446 MPFR_SET_POS (x);
447 mpfr_exp10 (y, x, MPFR_RNDN);
448 if (!MPFR_IS_INF (y))
449 {
450 printf ("evaluation of function in INF does not return INF\n");
451 exit (1);
452 }
453
454 MPFR_CHANGE_SIGN (x);
455 mpfr_exp10 (y, x, MPFR_RNDN);
456 if (!MPFR_IS_ZERO (y))
457 {
458 printf ("evaluation of function in -INF does not return 0\n");
459 exit (1);
460 }
461
462 MPFR_SET_NAN (x);
463 mpfr_exp10 (y, x, MPFR_RNDN);
464 if (!MPFR_IS_NAN (y))
465 {
466 printf ("evaluation of function in NaN does not return NaN\n");
467 exit (1);
468 }
469
470 if ((mpfr_uexp_t) 8 << 31 != 0 ||
471 mpfr_get_emax () <= (mpfr_uexp_t) 100000 * 100000)
472 {
473 /* emax <= 10000000000 */
474 mpfr_set_prec (x, 40);
475 mpfr_set_prec (y, 40);
476 mpfr_set_str (x, "3010299957", 10, MPFR_RNDN);
477 mpfr_clear_flags ();
478 inex = mpfr_exp10 (y, x, MPFR_RNDN);
479 ov = mpfr_overflow_p ();
480 if (!(MPFR_IS_INF (y) && MPFR_IS_POS (y) && ov))
481 {
482 printf ("Overflow error for x = 3010299957, MPFR_RNDN.\n");
483 mpfr_dump (y);
484 printf ("inex = %d, %soverflow\n", inex, ov ? "" : "no ");
485 exit (1);
486 }
487 }
488
489 test_generic (MPFR_PREC_MIN, 100, 100);
490
491 mpfr_clear (x);
492 mpfr_clear (y);
493
494 overfl_exp10_0 ();
495
496 data_check ("data/exp10", mpfr_exp10, "mpfr_exp10");
497 bad_cases (mpfr_exp10, mpfr_log10, "mpfr_exp10",
498 0, -256, 255, 4, 128, 800, 50);
499
500 tests_end_mpfr ();
501 return 0;
502 }