tests/tpow_all.c

   1 /* Test file for the various power functions
   2
   3 Copyright 2008, 2009, 2010, 2011, 2012 Free Software Foundation, Inc.
   4 Contributed by the AriC and Caramel 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 http://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 /* Note: some tests of the other tpow* test files could be moved there.
  24    The main goal of this test file is to test _all_ the power functions
  25    on special values, to make sure that they are consistent and give the
  26    expected result, in particular because such special cases are handled
  27    in different ways in each function. */
  28
  29 /* Execute with at least an argument to report all the errors found by
  30    comparisons. */
  31
  32 #include <stdlib.h>
  33
  34 #include "mpfr-test.h"
  35
  36 /* Behavior of cmpres (called by test_others):
  37  *   0: stop as soon as an error is found.
  38  *   1: report all errors found by test_others.
  39  *  -1: the 1 is changed to this value as soon as an error has been found.
  40  */
  41 static int all_cmpres_errors;
  42
  43 /* Non-zero when extended exponent range */
  44 static int ext = 0;
  45
  46 static const char *val[] =
  47   { "min", "min+", "max", "@NaN@", "-@Inf@", "-4", "-3", "-2", "-1.5",
  48     "-1", "-0.5", "-0", "0", "0.5", "1", "1.5", "2", "3", "4", "@Inf@" };
  49
  50 static void
  51 err (const char *s, int i, int j, int rnd, mpfr_srcptr z, int inex)
  52 {
  53   puts (s);
  54   if (ext)
  55     puts ("extended exponent range");
  56   printf ("x = %s, y = %s, %s\n", val[i], val[j],
  57           mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
  58   printf ("z = ");
  59   mpfr_out_str (stdout, 10, 0, z, MPFR_RNDN);
  60   printf ("\ninex = %d\n", inex);
  61   exit (1);
  62 }
  63
  64 /* Arguments:
  65  *   spx: non-zero if px is a stringm zero if px is a MPFR number.
  66  *   px: value of x (string or MPFR number).
  67  *   sy: value of y (string).
  68  *   rnd: rounding mode.
  69  *   z1: expected result (null pointer if unknown pure FP value).
  70  *   inex1: expected ternary value (if z1 is not a null pointer).
  71  *   z2: computed result.
  72  *   inex2: computed ternary value.
  73  *   flags1: expected flags (computed flags in __gmpfr_flags).
  74  *   s1, s2: strings about the context.
  75  */
  76 static void
  77 cmpres (int spx, const void *px, const char *sy, mpfr_rnd_t rnd,
  78         mpfr_srcptr z1, int inex1, mpfr_srcptr z2, int inex2,
  79         unsigned int flags1, const char *s1, const char *s2)
  80 {
  81   unsigned int flags2 = __gmpfr_flags;
  82
  83   if (flags1 == flags2)
  84     {
  85       /* Note: the test on the sign of z1 and z2 is needed
  86          in case they are both zeros. */
  87       if (z1 == NULL)
  88         {
  89           if (MPFR_IS_PURE_FP (z2))
  90             return;
  91         }
  92       else if (SAME_SIGN (inex1, inex2) &&
  93                ((MPFR_IS_NAN (z1) && MPFR_IS_NAN (z2)) ||
  94                 ((MPFR_IS_NEG (z1) ^ MPFR_IS_NEG (z2)) == 0 &&
  95                  mpfr_equal_p (z1, z2))))
  96         return;
  97     }
  98
  99   printf ("Error in %s\nwith %s%s\nx = ", s1, s2,
 100           ext ? ", extended exponent range" : "");
 101   if (spx)
 102     printf ("%s, ", (char *) px);
 103   else
 104     {
 105       mpfr_out_str (stdout, 16, 0, (mpfr_ptr) px, MPFR_RNDN);
 106       puts (",");
 107     }
 108   printf ("y = %s, %s\n", sy, mpfr_print_rnd_mode (rnd));
 109   printf ("Expected ");
 110   if (z1 == NULL)
 111     {
 112       printf ("pure FP value, flags = %u\n", flags1);
 113     }
 114   else
 115     {
 116       mpfr_out_str (stdout, 16, 0, z1, MPFR_RNDN);
 117       printf (", inex = %d, flags = %u\n", SIGN (inex1), flags1);
 118     }
 119   printf ("Got      ");
 120   mpfr_out_str (stdout, 16, 0, z2, MPFR_RNDN);
 121   printf (", inex = %d, flags = %u\n", SIGN (inex2), flags2);
 122   if (all_cmpres_errors != 0)
 123     all_cmpres_errors = -1;
 124   else
 125     exit (1);
 126 }
 127
 128 static int
 129 is_odd (mpfr_srcptr x)
 130 {
 131   /* works only with the values from val[] */
 132   return mpfr_integer_p (x) && mpfr_fits_slong_p (x, MPFR_RNDN) &&
 133     (mpfr_get_si (x, MPFR_RNDN) & 1);
 134 }
 135
 136 /* Compare the result (z1,inex1) of mpfr_pow with all flags cleared
 137    with those of mpfr_pow with all flags set and of the other power
 138    functions. Arguments x and y are the input values; sx and sy are
 139    their string representations (sx may be null); rnd contains the
 140    rounding mode; s is a string containing the function that called
 141    test_others. */
 142 static void
 143 test_others (const void *sx, const char *sy, mpfr_rnd_t rnd,
 144              mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z1,
 145              int inex1, unsigned int flags, const char *s)
 146 {
 147   mpfr_t z2;
 148   int inex2;
 149   int spx = sx != NULL;
 150
 151   if (!spx)
 152     sx = x;
 153
 154   mpfr_init2 (z2, mpfr_get_prec (z1));
 155
 156   __gmpfr_flags = MPFR_FLAGS_ALL;
 157   inex2 = mpfr_pow (z2, x, y, rnd);
 158   cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
 159           s, "mpfr_pow, flags set");
 160
 161   /* If y is an integer that fits in an unsigned long and is not -0,
 162      we can test mpfr_pow_ui. */
 163   if (MPFR_IS_POS (y) && mpfr_integer_p (y) &&
 164       mpfr_fits_ulong_p (y, MPFR_RNDN))
 165     {
 166       unsigned long yy = mpfr_get_ui (y, MPFR_RNDN);
 167
 168       mpfr_clear_flags ();
 169       inex2 = mpfr_pow_ui (z2, x, yy, rnd);
 170       cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
 171               s, "mpfr_pow_ui, flags cleared");
 172       __gmpfr_flags = MPFR_FLAGS_ALL;
 173       inex2 = mpfr_pow_ui (z2, x, yy, rnd);
 174       cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
 175               s, "mpfr_pow_ui, flags set");
 176
 177       /* If x is an integer that fits in an unsigned long and is not -0,
 178          we can also test mpfr_ui_pow_ui. */
 179       if (MPFR_IS_POS (x) && mpfr_integer_p (x) &&
 180           mpfr_fits_ulong_p (x, MPFR_RNDN))
 181         {
 182           unsigned long xx = mpfr_get_ui (x, MPFR_RNDN);
 183
 184           mpfr_clear_flags ();
 185           inex2 = mpfr_ui_pow_ui (z2, xx, yy, rnd);
 186           cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
 187                   s, "mpfr_ui_pow_ui, flags cleared");
 188           __gmpfr_flags = MPFR_FLAGS_ALL;
 189           inex2 = mpfr_ui_pow_ui (z2, xx, yy, rnd);
 190           cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
 191                   s, "mpfr_ui_pow_ui, flags set");
 192         }
 193     }
 194
 195   /* If y is an integer but not -0 and not huge, we can test mpfr_pow_z,
 196      and possibly mpfr_pow_si (and possibly mpfr_ui_div). */
 197   if (MPFR_IS_ZERO (y) ? MPFR_IS_POS (y) :
 198       (mpfr_integer_p (y) && MPFR_GET_EXP (y) < 256))
 199     {
 200       mpz_t yyy;
 201
 202       /* If y fits in a long, we can test mpfr_pow_si. */
 203       if (mpfr_fits_slong_p (y, MPFR_RNDN))
 204         {
 205           long yy = mpfr_get_si (y, MPFR_RNDN);
 206
 207           mpfr_clear_flags ();
 208           inex2 = mpfr_pow_si (z2, x, yy, rnd);
 209           cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
 210                   s, "mpfr_pow_si, flags cleared");
 211           __gmpfr_flags = MPFR_FLAGS_ALL;
 212           inex2 = mpfr_pow_si (z2, x, yy, rnd);
 213           cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
 214                   s, "mpfr_pow_si, flags set");
 215
 216           /* If y = -1, we can test mpfr_ui_div. */
 217           if (yy == -1)
 218             {
 219               mpfr_clear_flags ();
 220               inex2 = mpfr_ui_div (z2, 1, x, rnd);
 221               cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
 222                       s, "mpfr_ui_div, flags cleared");
 223               __gmpfr_flags = MPFR_FLAGS_ALL;
 224               inex2 = mpfr_ui_div (z2, 1, x, rnd);
 225               cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
 226                       s, "mpfr_ui_div, flags set");
 227             }
 228
 229           /* If y = 2, we can test mpfr_sqr. */
 230           if (yy == 2)
 231             {
 232               mpfr_clear_flags ();
 233               inex2 = mpfr_sqr (z2, x, rnd);
 234               cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
 235                       s, "mpfr_sqr, flags cleared");
 236               __gmpfr_flags = MPFR_FLAGS_ALL;
 237               inex2 = mpfr_sqr (z2, x, rnd);
 238               cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
 239                       s, "mpfr_sqr, flags set");
 240             }
 241         }
 242
 243       /* Test mpfr_pow_z. */
 244       mpz_init (yyy);
 245       mpfr_get_z (yyy, y, MPFR_RNDN);
 246       mpfr_clear_flags ();
 247       inex2 = mpfr_pow_z (z2, x, yyy, rnd);
 248       cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
 249               s, "mpfr_pow_z, flags cleared");
 250       __gmpfr_flags = MPFR_FLAGS_ALL;
 251       inex2 = mpfr_pow_z (z2, x, yyy, rnd);
 252       cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
 253               s, "mpfr_pow_z, flags set");
 254       mpz_clear (yyy);
 255     }
 256
 257   /* If y = 0.5, we can test mpfr_sqrt, except if x is -0 or -Inf (because
 258      the rule for mpfr_pow on these special values is different). */
 259   if (MPFR_IS_PURE_FP (y) && mpfr_cmp_str1 (y, "0.5") == 0 &&
 260       ! ((MPFR_IS_ZERO (x) || MPFR_IS_INF (x)) && MPFR_IS_NEG (x)))
 261     {
 262       mpfr_clear_flags ();
 263       inex2 = mpfr_sqrt (z2, x, rnd);
 264       cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
 265               s, "mpfr_sqrt, flags cleared");
 266       __gmpfr_flags = MPFR_FLAGS_ALL;
 267       inex2 = mpfr_sqrt (z2, x, rnd);
 268       cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
 269               s, "mpfr_sqrt, flags set");
 270     }
 271
 272 #if MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)
 273   /* If y = -0.5, we can test mpfr_rec_sqrt, except if x = -Inf
 274      (because the rule for mpfr_pow on -Inf is different). */
 275   if (MPFR_IS_PURE_FP (y) && mpfr_cmp_str1 (y, "-0.5") == 0 &&
 276       ! (MPFR_IS_INF (x) && MPFR_IS_NEG (x)))
 277     {
 278       mpfr_clear_flags ();
 279       inex2 = mpfr_rec_sqrt (z2, x, rnd);
 280       cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
 281               s, "mpfr_rec_sqrt, flags cleared");
 282       __gmpfr_flags = MPFR_FLAGS_ALL;
 283       inex2 = mpfr_rec_sqrt (z2, x, rnd);
 284       cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
 285               s, "mpfr_rec_sqrt, flags set");
 286     }
 287 #endif
 288
 289   /* If x is an integer that fits in an unsigned long and is not -0,
 290      we can test mpfr_ui_pow. */
 291   if (MPFR_IS_POS (x) && mpfr_integer_p (x) &&
 292       mpfr_fits_ulong_p (x, MPFR_RNDN))
 293     {
 294       unsigned long xx = mpfr_get_ui (x, MPFR_RNDN);
 295
 296       mpfr_clear_flags ();
 297       inex2 = mpfr_ui_pow (z2, xx, y, rnd);
 298       cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
 299               s, "mpfr_ui_pow, flags cleared");
 300       __gmpfr_flags = MPFR_FLAGS_ALL;
 301       inex2 = mpfr_ui_pow (z2, xx, y, rnd);
 302       cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
 303               s, "mpfr_ui_pow, flags set");
 304
 305       /* If x = 2, we can test mpfr_exp2. */
 306       if (xx == 2)
 307         {
 308           mpfr_clear_flags ();
 309           inex2 = mpfr_exp2 (z2, y, rnd);
 310           cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
 311                   s, "mpfr_exp2, flags cleared");
 312           __gmpfr_flags = MPFR_FLAGS_ALL;
 313           inex2 = mpfr_exp2 (z2, y, rnd);
 314           cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
 315                   s, "mpfr_exp2, flags set");
 316         }
 317
 318       /* If x = 10, we can test mpfr_exp10. */
 319       if (xx == 10)
 320         {
 321           mpfr_clear_flags ();
 322           inex2 = mpfr_exp10 (z2, y, rnd);
 323           cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
 324                   s, "mpfr_exp10, flags cleared");
 325           __gmpfr_flags = MPFR_FLAGS_ALL;
 326           inex2 = mpfr_exp10 (z2, y, rnd);
 327           cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
 328                   s, "mpfr_exp10, flags set");
 329         }
 330     }
 331
 332   mpfr_clear (z2);
 333 }
 334
 335 static int
 336 my_setstr (mpfr_ptr t, const char *s)
 337 {
 338   if (strcmp (s, "min") == 0)
 339     {
 340       mpfr_setmin (t, mpfr_get_emin ());
 341       MPFR_SET_POS (t);
 342       return 0;
 343     }
 344   if (strcmp (s, "min+") == 0)
 345     {
 346       mpfr_setmin (t, mpfr_get_emin ());
 347       MPFR_SET_POS (t);
 348       mpfr_nextabove (t);
 349       return 0;
 350     }
 351   if (strcmp (s, "max") == 0)
 352     {
 353       mpfr_setmax (t, mpfr_get_emax ());
 354       MPFR_SET_POS (t);
 355       return 0;
 356     }
 357   return mpfr_set_str (t, s, 10, MPFR_RNDN);
 358 }
 359
 360 static void
 361 tst (void)
 362 {
 363   int sv = sizeof (val) / sizeof (*val);
 364   int i, j;
 365   int rnd;
 366   mpfr_t x, y, z, tmp;
 367
 368   mpfr_inits2 (53, x, y, z, tmp, (mpfr_ptr) 0);
 369
 370   for (i = 0; i < sv; i++)
 371     for (j = 0; j < sv; j++)
 372       RND_LOOP (rnd)
 373         {
 374           int exact, inex;
 375           unsigned int flags;
 376
 377           if (my_setstr (x, val[i]) || my_setstr (y, val[j]))
 378             {
 379               printf ("internal error for (%d,%d,%d)\n", i, j, rnd);
 380               exit (1);
 381             }
 382           mpfr_clear_flags ();
 383           inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd);
 384           flags = __gmpfr_flags;
 385           if (! MPFR_IS_NAN (z) && mpfr_nanflag_p ())
 386             err ("got NaN flag without NaN value", i, j, rnd, z, inex);
 387           if (MPFR_IS_NAN (z) && ! mpfr_nanflag_p ())
 388             err ("got NaN value without NaN flag", i, j, rnd, z, inex);
 389           if (inex != 0 && ! mpfr_inexflag_p ())
 390             err ("got non-zero ternary value without inexact flag",
 391                  i, j, rnd, z, inex);
 392           if (inex == 0 && mpfr_inexflag_p ())
 393             err ("got null ternary value with inexact flag",
 394                  i, j, rnd, z, inex);
 395           if (i >= 3 && j >= 3)
 396             {
 397               if (mpfr_underflow_p ())
 398                 err ("got underflow", i, j, rnd, z, inex);
 399               if (mpfr_overflow_p ())
 400                 err ("got overflow", i, j, rnd, z, inex);
 401               exact = MPFR_IS_SINGULAR (z) ||
 402                 (mpfr_mul_2ui (tmp, z, 16, MPFR_RNDN), mpfr_integer_p (tmp));
 403               if (exact && inex != 0)
 404                 err ("got exact value with ternary flag different from 0",
 405                      i, j, rnd, z, inex);
 406               if (! exact && inex == 0)
 407                 err ("got inexact value with ternary flag equal to 0",
 408                      i, j, rnd, z, inex);
 409             }
 410           if (MPFR_IS_ZERO (x) && ! MPFR_IS_NAN (y) && MPFR_NOTZERO (y))
 411             {
 412               if (MPFR_IS_NEG (y) && ! MPFR_IS_INF (z))
 413                 err ("expected an infinity", i, j, rnd, z, inex);
 414               if (MPFR_IS_POS (y) && ! MPFR_IS_ZERO (z))
 415                 err ("expected a zero", i, j, rnd, z, inex);
 416               if ((MPFR_IS_NEG (x) && is_odd (y)) ^ MPFR_IS_NEG (z))
 417                 err ("wrong sign", i, j, rnd, z, inex);
 418             }
 419           if (! MPFR_IS_NAN (x) && mpfr_cmp_si (x, -1) == 0)
 420             {
 421               /* x = -1 */
 422               if (! (MPFR_IS_INF (y) || mpfr_integer_p (y)) &&
 423                   ! MPFR_IS_NAN (z))
 424                 err ("expected NaN", i, j, rnd, z, inex);
 425               if ((MPFR_IS_INF (y) || (mpfr_integer_p (y) && ! is_odd (y)))
 426                   && ! mpfr_equal_p (z, __gmpfr_one))
 427                 err ("expected 1", i, j, rnd, z, inex);
 428               if (is_odd (y) &&
 429                   (MPFR_IS_NAN (z) || mpfr_cmp_si (z, -1) != 0))
 430                 err ("expected -1", i, j, rnd, z, inex);
 431             }
 432           if ((mpfr_equal_p (x, __gmpfr_one) || MPFR_IS_ZERO (y)) &&
 433               ! mpfr_equal_p (z, __gmpfr_one))
 434             err ("expected 1", i, j, rnd, z, inex);
 435           if (MPFR_IS_PURE_FP (x) && MPFR_IS_NEG (x) &&
 436               MPFR_IS_FP (y) && ! mpfr_integer_p (y) &&
 437               ! MPFR_IS_NAN (z))
 438             err ("expected NaN", i, j, rnd, z, inex);
 439           if (MPFR_IS_INF (y) && MPFR_NOTZERO (x))
 440             {
 441               int cmpabs1 = mpfr_cmpabs (x, __gmpfr_one);
 442
 443               if ((MPFR_IS_NEG (y) ? (cmpabs1 < 0) : (cmpabs1 > 0)) &&
 444                   ! (MPFR_IS_POS (z) && MPFR_IS_INF (z)))
 445                 err ("expected +Inf", i, j, rnd, z, inex);
 446               if ((MPFR_IS_NEG (y) ? (cmpabs1 > 0) : (cmpabs1 < 0)) &&
 447                   ! (MPFR_IS_POS (z) && MPFR_IS_ZERO (z)))
 448                 err ("expected +0", i, j, rnd, z, inex);
 449             }
 450           if (MPFR_IS_INF (x) && ! MPFR_IS_NAN (y) && MPFR_NOTZERO (y))
 451             {
 452               if (MPFR_IS_POS (y) && ! MPFR_IS_INF (z))
 453                 err ("expected an infinity", i, j, rnd, z, inex);
 454               if (MPFR_IS_NEG (y) && ! MPFR_IS_ZERO (z))
 455                 err ("expected a zero", i, j, rnd, z, inex);
 456               if ((MPFR_IS_NEG (x) && is_odd (y)) ^ MPFR_IS_NEG (z))
 457                 err ("wrong sign", i, j, rnd, z, inex);
 458             }
 459           test_others (val[i], val[j], (mpfr_rnd_t) rnd, x, y, z, inex, flags,
 460                        "tst");
 461         }
 462   mpfr_clears (x, y, z, tmp, (mpfr_ptr) 0);
 463 }
 464
 465 static void
 466 underflow_up1 (void)
 467 {
 468   mpfr_t delta, x, y, z, z0;
 469   mpfr_exp_t n;
 470   int inex;
 471   int rnd;
 472   int i;
 473
 474   n = mpfr_get_emin ();
 475   if (n < LONG_MIN)
 476     return;
 477
 478   mpfr_init2 (delta, 2);
 479   inex = mpfr_set_ui_2exp (delta, 1, -2, MPFR_RNDN);
 480   MPFR_ASSERTN (inex == 0);
 481
 482   mpfr_init2 (x, 8);
 483   inex = mpfr_set_ui (x, 2, MPFR_RNDN);
 484   MPFR_ASSERTN (inex == 0);
 485
 486   mpfr_init2 (y, sizeof (long) * CHAR_BIT + 4);
 487   inex = mpfr_set_si (y, n, MPFR_RNDN);
 488   MPFR_ASSERTN (inex == 0);
 489
 490   mpfr_init2 (z0, 2);
 491   mpfr_set_ui (z0, 0, MPFR_RNDN);
 492
 493   mpfr_init2 (z, 32);
 494
 495   for (i = 0; i <= 12; i++)
 496     {
 497       unsigned int flags = 0;
 498       char sy[16];
 499
 500       /* Test 2^(emin - i/4).
 501        * --> Underflow iff i > 4.
 502        * --> Zero in MPFR_RNDN iff i >= 8.
 503        */
 504
 505       if (i != 0 && i != 4)
 506         flags |= MPFR_FLAGS_INEXACT;
 507       if (i > 4)
 508         flags |= MPFR_FLAGS_UNDERFLOW;
 509
 510       sprintf (sy, "emin - %d/4", i);
 511
 512       RND_LOOP (rnd)
 513         {
 514           int zero;
 515
 516           zero = (i > 4 && (rnd == MPFR_RNDZ || rnd == MPFR_RNDD)) ||
 517             (i >= 8 && rnd == MPFR_RNDN);
 518
 519           mpfr_clear_flags ();
 520           inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd);
 521           cmpres (1, "2", sy, (mpfr_rnd_t) rnd, zero ? z0 : (mpfr_ptr) NULL,
 522                   -1, z, inex, flags, "underflow_up1", "mpfr_pow");
 523           test_others ("2", sy, (mpfr_rnd_t) rnd, x, y, z, inex, flags,
 524                        "underflow_up1");
 525         }
 526
 527       inex = mpfr_sub (y, y, delta, MPFR_RNDN);
 528       MPFR_ASSERTN (inex == 0);
 529     }
 530
 531   mpfr_clears (delta, x, y, z, z0, (mpfr_ptr) 0);
 532 }
 533
 534 /* With pow.c r5497, the following test fails on a 64-bit Linux machine
 535  * due to a double-rounding problem when rescaling the result:
 536  *   Error with underflow_up2 and extended exponent range
 537  *   x = 7.fffffffffffffff0@-1,
 538  *   y = 4611686018427387904, MPFR_RNDN
 539  *   Expected 1.0000000000000000@-1152921504606846976, inex = 1, flags = 9
 540  *   Got      0, inex = -1, flags = 9
 541  * With pow_ui.c r5423, the following test fails on a 64-bit Linux machine
 542  * as underflows and overflows are not handled correctly (the approximation
 543  * error is ignored):
 544  *   Error with mpfr_pow_ui, flags cleared
 545  *   x = 7.fffffffffffffff0@-1,
 546  *   y = 4611686018427387904, MPFR_RNDN
 547  *   Expected 1.0000000000000000@-1152921504606846976, inex = 1, flags = 9
 548  *   Got      0, inex = -1, flags = 9
 549  */
 550 static void
 551 underflow_up2 (void)
 552 {
 553   mpfr_t x, y, z, z0, eps;
 554   mpfr_exp_t n;
 555   int inex;
 556   int rnd;
 557
 558   n = 1 - mpfr_get_emin ();
 559   MPFR_ASSERTN (n > 1);
 560   if (n > ULONG_MAX)
 561     return;
 562
 563   mpfr_init2 (eps, 2);
 564   mpfr_set_ui_2exp (eps, 1, -1, MPFR_RNDN);  /* 1/2 */
 565   mpfr_div_ui (eps, eps, n, MPFR_RNDZ);      /* 1/(2n) rounded toward zero */
 566
 567   mpfr_init2 (x, sizeof (unsigned long) * CHAR_BIT + 1);
 568   inex = mpfr_ui_sub (x, 1, eps, MPFR_RNDN);
 569   MPFR_ASSERTN (inex == 0);  /* since n < 2^(size_of_long_in_bits) */
 570   inex = mpfr_div_2ui (x, x, 1, MPFR_RNDN);  /* 1/2 - eps/2 exactly */
 571   MPFR_ASSERTN (inex == 0);
 572
 573   mpfr_init2 (y, sizeof (unsigned long) * CHAR_BIT);
 574   inex = mpfr_set_ui (y, n, MPFR_RNDN);
 575   MPFR_ASSERTN (inex == 0);
 576
 577   /* 0 < eps < 1 / (2n), thus (1 - eps)^n > 1/2,
 578      and 1/2 (1/2)^n < (1/2 - eps/2)^n < (1/2)^n. */
 579   mpfr_inits2 (64, z, z0, (mpfr_ptr) 0);
 580   RND_LOOP (rnd)
 581     {
 582       unsigned int ufinex = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
 583       int expected_inex;
 584       char sy[256];
 585
 586       mpfr_set_ui (z0, 0, MPFR_RNDN);
 587       expected_inex = rnd == MPFR_RNDN || rnd == MPFR_RNDU || rnd == MPFR_RNDA ?
 588         (mpfr_nextabove (z0), 1) : -1;
 589       sprintf (sy, "%lu", (unsigned long) n);
 590
 591       mpfr_clear_flags ();
 592       inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd);
 593       cmpres (0, x, sy, (mpfr_rnd_t) rnd, z0, expected_inex, z, inex, ufinex,
 594               "underflow_up2", "mpfr_pow");
 595       test_others (NULL, sy, (mpfr_rnd_t) rnd, x, y, z, inex, ufinex,
 596                    "underflow_up2");
 597     }
 598
 599   mpfr_clears (x, y, z, z0, eps, (mpfr_ptr) 0);
 600 }
 601
 602 static void
 603 underflow_up3 (void)
 604 {
 605   mpfr_t x, y, z, z0;
 606   int inex;
 607   int rnd;
 608   int i;
 609
 610   mpfr_init2 (x, 64);
 611   mpfr_init2 (y, sizeof (mpfr_exp_t) * CHAR_BIT);
 612   mpfr_init2 (z, 32);
 613   mpfr_init2 (z0, 2);
 614
 615   inex = mpfr_set_exp_t (y, mpfr_get_emin () - 2, MPFR_RNDN);
 616   MPFR_ASSERTN (inex == 0);
 617   for (i = -1; i <= 1; i++)
 618     RND_LOOP (rnd)
 619       {
 620         unsigned int ufinex = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
 621         int expected_inex;
 622
 623         mpfr_set_ui (x, 2, MPFR_RNDN);
 624         if (i < 0)
 625           mpfr_nextbelow (x);
 626         if (i > 0)
 627           mpfr_nextabove (x);
 628         /* x = 2 + i * eps, y = emin - 2, x^y ~= 2^(emin - 2) */
 629
 630         expected_inex = rnd == MPFR_RNDU || rnd == MPFR_RNDA
 631           || (rnd == MPFR_RNDN && i < 0) ? 1 : -1;
 632
 633         mpfr_set_ui (z0, 0, MPFR_RNDN);
 634         if (expected_inex > 0)
 635           mpfr_nextabove (z0);
 636
 637         mpfr_clear_flags ();
 638         inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd);
 639         cmpres (0, x, "emin - 2", (mpfr_rnd_t) rnd, z0, expected_inex, z, inex,
 640                 ufinex, "underflow_up3", "mpfr_pow");
 641         test_others (NULL, "emin - 2", (mpfr_rnd_t) rnd, x, y, z, inex, ufinex,
 642                      "underflow_up3");
 643       }
 644
 645   mpfr_clears (x, y, z, z0, (mpfr_ptr) 0);
 646 }
 647
 648 static void
 649 underflow_up (void)
 650 {
 651   underflow_up1 ();
 652   underflow_up2 ();
 653   underflow_up3 ();
 654 }
 655
 656 static void
 657 overflow_inv (void)
 658 {
 659   mpfr_t x, y, z;
 660   int precx;
 661   int s, t;
 662   int inex;
 663   int rnd;
 664
 665   mpfr_init2 (y, 2);
 666   mpfr_init2 (z, 8);
 667
 668   mpfr_set_si (y, -1, MPFR_RNDN);
 669   for (precx = 10; precx <= 100; precx += 90)
 670     {
 671       const char *sp = precx == 10 ?
 672         "overflow_inv (precx = 10)" : "overflow_inv (precx = 100)";
 673
 674       mpfr_init2 (x, precx);
 675       for (s = -1; s <= 1; s += 2)
 676         {
 677           inex = mpfr_set_si_2exp (x, s, - mpfr_get_emax (), MPFR_RNDN);
 678           MPFR_ASSERTN (inex == 0);
 679           for (t = 0; t <= 5; t++)
 680             {
 681               /* If precx = 10:
 682                * x = s * 2^(-emax) * (1 + t * 2^(-9)), so that
 683                * 1/x = s * 2^emax * (1 - t * 2^(-9) + eps) with eps > 0.
 684                * Values of (1/x) / 2^emax and overflow condition for x > 0:
 685                * t = 0: 1                           o: always
 686                * t = 1: 0.11111111 100000000011...  o: MPFR_RNDN and MPFR_RNDU
 687                * t = 2: 0.11111111 000000001111...  o: MPFR_RNDU
 688                * t = 3: 0.11111110 100000100011...  o: never
 689                *
 690                * If precx = 100:
 691                * t = 0: always overflow
 692                * t > 0: overflow for MPFR_RNDN and MPFR_RNDU.
 693                */
 694               RND_LOOP (rnd)
 695                 {
 696                   int inf, overflow;
 697                   mpfr_rnd_t rnd2;
 698
 699                   if (rnd == MPFR_RNDA)
 700                     rnd2 = s < 0 ? MPFR_RNDD : MPFR_RNDU;
 701                   else
 702                     rnd2 = (mpfr_rnd_t) rnd;
 703
 704                   overflow = t == 0 ||
 705                     ((mpfr_rnd_t) rnd == MPFR_RNDN &&
 706                      (precx > 10 || t == 1)) ||
 707                     (rnd2 == (s < 0 ? MPFR_RNDD : MPFR_RNDU) &&
 708                      (precx > 10 || t <= 2));
 709                   inf = overflow &&
 710                     ((mpfr_rnd_t) rnd == MPFR_RNDN ||
 711                      rnd2 == (s < 0 ? MPFR_RNDD : MPFR_RNDU));
 712                   mpfr_clear_flags ();
 713                   inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd);
 714                   if (overflow ^ !! mpfr_overflow_p ())
 715                     {
 716                       printf ("Bad overflow flag in %s\nfor mpfr_pow%s\n"
 717                               "s = %d, t = %d, %s\n", sp,
 718                               ext ? ", extended exponent range" : "",
 719                               s, t, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
 720                       exit (1);
 721                     }
 722                   if (overflow && (inf ^ !! MPFR_IS_INF (z)))
 723                     {
 724                       printf ("Bad value in %s\nfor mpfr_pow%s\n"
 725                               "s = %d, t = %d, %s\nGot ", sp,
 726                               ext ? ", extended exponent range" : "",
 727                               s, t, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
 728                       mpfr_out_str (stdout, 16, 0, z, MPFR_RNDN);
 729                       printf (" instead of %s value.\n",
 730                               inf ? "infinite" : "finite");
 731                       exit (1);
 732                     }
 733                   test_others (NULL, "-1", (mpfr_rnd_t) rnd, x, y, z,
 734                                inex, __gmpfr_flags, sp);
 735                 }  /* RND_LOOP */
 736               mpfr_nexttoinf (x);
 737             }  /* for (t = ...) */
 738         }  /* for (s = ...) */
 739       mpfr_clear (x);
 740     }  /* for (precx = ...) */
 741
 742   mpfr_clears (y, z, (mpfr_ptr) 0);
 743 }
 744
 745 static void
 746 alltst (void)
 747 {
 748   mpfr_exp_t emin, emax;
 749
 750   ext = 0;
 751   tst ();
 752   underflow_up ();
 753   overflow_inv ();
 754
 755   emin = mpfr_get_emin ();
 756   emax = mpfr_get_emax ();
 757   set_emin (MPFR_EMIN_MIN);
 758   set_emax (MPFR_EMAX_MAX);
 759   if (mpfr_get_emin () != emin || mpfr_get_emax () != emax)
 760     {
 761       ext = 1;
 762       tst ();
 763       underflow_up ();
 764       overflow_inv ();
 765       set_emin (emin);
 766       set_emax (emax);
 767     }
 768 }
 769
 770 int
 771 main (int argc, char *argv[])
 772 {
 773   tests_start_mpfr ();
 774   all_cmpres_errors = argc > 1;
 775   alltst ();
 776   tests_end_mpfr ();
 777   return all_cmpres_errors < 0;
 778 }