1/* Test file for mpfr_exp10.
2
3Copyright 2007-2023 Free Software Foundation, Inc.
4Contributed by the AriC and Caramba projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 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
30static void
31special_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
62static void
63emax_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
90static void
91exp_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
113static void
114overfl_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*/
209static void
210bug20230213 (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*/
285static void
286bug20230427 (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
374int
375main (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}
503