tests/tgeneric.c

   1 /* Generic test file for functions with one or two arguments (the second being
   2    either mpfr_t or double).
   3
   4 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Free Software Foundation, Inc.
   5 Contributed by the AriC and Caramel projects, INRIA.
   6
   7 This file is part of the GNU MPFR Library.
   8
   9 The GNU MPFR Library is free software; you can redistribute it and/or modify
  10 it under the terms of the GNU Lesser General Public License as published by
  11 the Free Software Foundation; either version 3 of the License, or (at your
  12 option) any later version.
  13
  14 The GNU MPFR Library is distributed in the hope that it will be useful, but
  15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  17 License for more details.
  18
  19 You should have received a copy of the GNU Lesser General Public License
  20 along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
  21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  23
  24 /* define TWO_ARGS for two-argument functions like mpfr_pow
  25    define DOUBLE_ARG1 or DOUBLE_ARG2 for function with a double operand in
  26    first or second place like sub_d or d_sub */
  27
  28 #ifndef TEST_RANDOM_POS
  29 /* For the random function: one number on two is negative. */
  30 #define TEST_RANDOM_POS 256
  31 #endif
  32
  33 #ifndef TEST_RANDOM_POS2
  34 /* For the random function: one number on two is negative. */
  35 #define TEST_RANDOM_POS2 256
  36 #endif
  37
  38 #ifndef TEST_RANDOM_EMIN
  39 #define TEST_RANDOM_EMIN -256
  40 #endif
  41
  42 #ifndef TEST_RANDOM_EMAX
  43 #define TEST_RANDOM_EMAX 255
  44 #endif
  45
  46 /* If the MPFR_SUSPICIOUS_OVERFLOW test fails but this is not a bug,
  47    then define TGENERIC_SO_TEST with an adequate test (possibly 0) to
  48    omit this particular case. */
  49 #ifndef TGENERIC_SO_TEST
  50 #define TGENERIC_SO_TEST 1
  51 #endif
  52
  53 #define STR(F) #F
  54 #define MAKE_STR(S) STR(S)
  55
  56 /* The (void *) below is needed to avoid a warning with gcc 4.2+ and functions
  57  * with 2 arguments. See <http://gcc.gnu.org/bugzilla/show_bug.cgi?id=36299>.
  58  */
  59 #define TGENERIC_FAIL(S, X, U)                                          \
  60   do                                                                    \
  61     {                                                                   \
  62       printf ("%s\nx = ", (S));                                         \
  63       mpfr_out_str (stdout, 2, 0, (X), MPFR_RNDN);                      \
  64       printf ("\n");                                                    \
  65       if ((void *) U != 0)                                              \
  66         {                                                               \
  67           printf ("u = ");                                              \
  68           mpfr_out_str (stdout, 2, 0, (U), MPFR_RNDN);                  \
  69           printf ("\n");                                                \
  70         }                                                               \
  71       printf ("yprec = %u, rnd_mode = %s, inexact = %d, flags = %u\n",  \
  72               (unsigned int) yprec, mpfr_print_rnd_mode (rnd), compare, \
  73               (unsigned int) __gmpfr_flags);                            \
  74       exit (1);                                                         \
  75     }                                                                   \
  76   while (0)
  77
  78 #define TGENERIC_CHECK_AUX(S, EXPR, U)                        \
  79   do                                                          \
  80     if (!(EXPR))                                              \
  81       TGENERIC_FAIL (S " in " MAKE_STR(TEST_FUNCTION), x, U); \
  82   while (0)
  83
  84 #undef TGENERIC_CHECK
  85 #if defined(TWO_ARGS) || defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
  86 #define TGENERIC_CHECK(S, EXPR) TGENERIC_CHECK_AUX(S, EXPR, u)
  87 #else
  88 #define TGENERIC_CHECK(S, EXPR) TGENERIC_CHECK_AUX(S, EXPR, 0)
  89 #endif
  90
  91 #ifdef DEBUG_TGENERIC
  92 #define TGENERIC_IAUX(F,P,X,U)                                          \
  93   do                                                                    \
  94     {                                                                   \
  95       printf ("tgeneric: testing function " STR(F)                      \
  96               ", %s, target prec = %lu\nx = ",                          \
  97               mpfr_print_rnd_mode (rnd), (unsigned long) (P));          \
  98       mpfr_out_str (stdout, 2, 0, (X), MPFR_RNDN);                      \
  99       printf ("\n");                                                    \
 100       if (U)                                                            \
 101         {                                                               \
 102           printf ("u = ");                                              \
 103           mpfr_out_str (stdout, 2, 0, (U), MPFR_RNDN);                  \
 104           printf ("\n");                                                \
 105         }                                                               \
 106     }                                                                   \
 107   while (0)
 108 #undef TGENERIC_INFO
 109 #if defined(TWO_ARGS) || defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
 110 #define TGENERIC_INFO(F,P) TGENERIC_IAUX(F,P,x,u)
 111 #else
 112 #define TGENERIC_INFO(F,P) TGENERIC_IAUX(F,P,x,0)
 113 #endif
 114 #endif
 115
 116 /* For some functions (for example cos), the argument reduction is too
 117    expensive when using mpfr_get_emax(). Then simply define REDUCE_EMAX
 118    to some reasonable value before including tgeneric.c. */
 119 #ifndef REDUCE_EMAX
 120 #define REDUCE_EMAX mpfr_get_emax ()
 121 #endif
 122
 123 static void
 124 test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N)
 125 {
 126   mpfr_prec_t prec, xprec, yprec;
 127   mpfr_t x, y, z, t, w;
 128 #ifdef TWO_ARGS
 129   mpfr_t u;
 130 #elif defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
 131   mpfr_t u;
 132   double d;
 133 #endif
 134   mpfr_rnd_t rnd;
 135   int inexact, compare, compare2;
 136   unsigned int n;
 137   unsigned long ctrt = 0, ctrn = 0;
 138   mpfr_exp_t old_emin, old_emax;
 139
 140   old_emin = mpfr_get_emin ();
 141   old_emax = mpfr_get_emax ();
 142
 143   mpfr_inits2 (MPFR_PREC_MIN, x, y, z, t, w, (mpfr_ptr) 0);
 144 #if defined(TWO_ARGS) || defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
 145   mpfr_init2 (u, MPFR_PREC_MIN);
 146 #endif
 147
 148   /* generic test */
 149   for (prec = p0; prec <= p1; prec++)
 150     {
 151       mpfr_set_prec (z, prec);
 152       mpfr_set_prec (t, prec);
 153       yprec = prec + 10;
 154       mpfr_set_prec (y, yprec);
 155       mpfr_set_prec (w, yprec);
 156
 157       /* Note: in precision p1, we test 4 special cases. */
 158       for (n = 0; n < (prec == p1 ? N + 4 : N); n++)
 159         {
 160           int infinite_input = 0;
 161
 162           xprec = prec;
 163           if (randlimb () & 1)
 164             {
 165               xprec *= (double) randlimb () / MP_LIMB_T_MAX;
 166               if (xprec < MPFR_PREC_MIN)
 167                 xprec = MPFR_PREC_MIN;
 168             }
 169           mpfr_set_prec (x, xprec);
 170 #ifdef TWO_ARGS
 171           mpfr_set_prec (u, xprec);
 172 #elif defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
 173           mpfr_set_prec (u, IEEE_DBL_MANT_DIG);
 174 #endif
 175
 176           if (n > 3 || prec < p1)
 177             {
 178 #if defined(RAND_FUNCTION)
 179               RAND_FUNCTION (x);
 180 #if defined(TWO_ARGS) || defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
 181               RAND_FUNCTION (u);
 182 #endif
 183 #else
 184               tests_default_random (x, TEST_RANDOM_POS,
 185                                     TEST_RANDOM_EMIN, TEST_RANDOM_EMAX);
 186 #if defined(TWO_ARGS) || defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
 187               tests_default_random (u, TEST_RANDOM_POS2,
 188                                     TEST_RANDOM_EMIN, TEST_RANDOM_EMAX);
 189 #endif
 190 #endif
 191             }
 192           else
 193             {
 194               /* Special cases tested in precision p1 if n <= 3. They are
 195                  useful really in the extended exponent range. */
 196 #if (defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)) && defined(MPFR_ERRDIVZERO)
 197               goto next_n;
 198 #endif
 199               set_emin (MPFR_EMIN_MIN);
 200               set_emax (MPFR_EMAX_MAX);
 201               if (n <= 1)
 202                 {
 203                   mpfr_set_si (x, n == 0 ? 1 : -1, MPFR_RNDN);
 204                   mpfr_set_exp (x, mpfr_get_emin ());
 205 #if defined(TWO_ARGS) || defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
 206                   mpfr_set_si (u, randlimb () % 2 == 0 ? 1 : -1, MPFR_RNDN);
 207                   mpfr_set_exp (u, mpfr_get_emin ());
 208 #endif
 209                 }
 210               else  /* 2 <= n <= 3 */
 211                 {
 212                   if (getenv ("MPFR_CHECK_MAX") == NULL)
 213                     goto next_n;
 214                   mpfr_set_si (x, n == 0 ? 1 : -1, MPFR_RNDN);
 215                   mpfr_setmax (x, REDUCE_EMAX);
 216 #if defined(TWO_ARGS) || defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
 217                   mpfr_set_si (u, randlimb () % 2 == 0 ? 1 : -1, MPFR_RNDN);
 218                   mpfr_setmax (u, mpfr_get_emax ());
 219 #endif
 220                 }
 221             }
 222
 223           rnd = RND_RAND ();
 224           mpfr_clear_flags ();
 225 #ifdef DEBUG_TGENERIC
 226           TGENERIC_INFO (TEST_FUNCTION, MPFR_PREC (y));
 227 #endif
 228 #if defined(TWO_ARGS)
 229           compare = TEST_FUNCTION (y, x, u, rnd);
 230 #elif defined(DOUBLE_ARG1)
 231           d = mpfr_get_d (u, rnd);
 232           compare = TEST_FUNCTION (y, d, x, rnd);
 233           /* d can be infinite due to overflow in mpfr_get_d */
 234           infinite_input |= DOUBLE_ISINF (d);
 235 #elif defined(DOUBLE_ARG2)
 236           d = mpfr_get_d (u, rnd);
 237           compare = TEST_FUNCTION (y, x, d, rnd);
 238           /* d can be infinite due to overflow in mpfr_get_d */
 239           infinite_input |= DOUBLE_ISINF (d);
 240 #else
 241           compare = TEST_FUNCTION (y, x, rnd);
 242 #endif
 243           TGENERIC_CHECK ("Bad inexact flag",
 244                           (compare != 0) ^ (mpfr_inexflag_p () == 0));
 245           ctrt++;
 246           /* Consistency test in a reduced exponent range. Doing it
 247              for the first 10 samples and for prec == p1 (which has
 248              some special cases) should be sufficient. */
 249           if (ctrt <= 10 || prec == p1)
 250             {
 251               unsigned int flags, oldflags = __gmpfr_flags;
 252               mpfr_exp_t e, emin, emax, oemin, oemax;
 253
 254               /* Determine the smallest exponent range containing the
 255                  exponents of the mpfr_t inputs (x, and u if TWO_ARGS)
 256                  and output (y). */
 257               emin = MPFR_EMAX_MAX;
 258               emax = MPFR_EMIN_MIN;
 259               if (MPFR_IS_PURE_FP (x))
 260                 {
 261                   e = MPFR_GET_EXP (x);
 262                   if (e < emin)
 263                     emin = e;
 264                   if (e > emax)
 265                     emax = e;
 266                 }
 267               if (MPFR_IS_PURE_FP (y))
 268                 {
 269                   e = MPFR_GET_EXP (y);
 270                   if (e < emin)
 271                     emin = e;
 272                   if (e > emax)
 273                     emax = e;
 274                 }
 275 #if defined(TWO_ARGS)
 276               if (MPFR_IS_PURE_FP (u))
 277                 {
 278                   e = MPFR_GET_EXP (u);
 279                   if (e < emin)
 280                     emin = e;
 281                   if (e > emax)
 282                     emax = e;
 283                 }
 284 #endif
 285               if (emin > emax)
 286                 emin = emax;  /* case where all values are singular */
 287               oemin = mpfr_get_emin ();
 288               oemax = mpfr_get_emax ();
 289               mpfr_set_emin (emin);
 290               mpfr_set_emax (emax);
 291 #ifdef DEBUG_TGENERIC
 292               /* Useful information in case of assertion failure. */
 293               printf ("tgeneric: reduced exponent range [%"
 294                       MPFR_EXP_FSPEC "d,%" MPFR_EXP_FSPEC "d]\n",
 295                       (mpfr_eexp_t) emin, (mpfr_eexp_t) emax);
 296 #endif
 297               mpfr_clear_flags ();
 298 #if defined(TWO_ARGS)
 299               inexact = TEST_FUNCTION (w, x, u, rnd);
 300 #elif defined(DOUBLE_ARG1)
 301               inexact = TEST_FUNCTION (w, d, x, rnd);
 302 #elif defined(DOUBLE_ARG2)
 303               inexact = TEST_FUNCTION (w, x, d, rnd);
 304 #else
 305               inexact = TEST_FUNCTION (w, x, rnd);
 306 #endif
 307               flags = __gmpfr_flags;
 308               mpfr_set_emin (oemin);
 309               mpfr_set_emax (oemax);
 310               if (! (SAME_VAL (w, y) &&
 311                      SAME_SIGN (inexact, compare) &&
 312                      flags == oldflags))
 313                 {
 314                   printf ("Error in " MAKE_STR(TEST_FUNCTION)
 315                           ", reduced exponent range [%"
 316                           MPFR_EXP_FSPEC "d,%" MPFR_EXP_FSPEC "d] on:\n",
 317                           (mpfr_eexp_t) emin, (mpfr_eexp_t) emax);
 318                   printf ("x = ");
 319                   mpfr_dump (x);
 320 #if defined(TWO_ARGS) || defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
 321                   printf ("u = ");
 322                   mpfr_dump (u);
 323 #endif
 324                   printf ("yprec = %u, rnd_mode = %s\n",
 325                           (unsigned int) yprec, mpfr_print_rnd_mode (rnd));
 326                   printf ("Expected:\n  y = ");
 327                   mpfr_dump (y);
 328                   printf ("  inex = %d, flags = %u\n",
 329                           SIGN (compare), oldflags);
 330                   printf ("Got:\n  w = ");
 331                   mpfr_dump (w);
 332                   printf ("  inex = %d, flags = %u\n",
 333                           SIGN (inexact), flags);
 334                   exit (1);
 335                 }
 336             }
 337           if (MPFR_IS_SINGULAR (y))
 338             {
 339               if (MPFR_IS_NAN (y) || mpfr_nanflag_p ())
 340                 TGENERIC_CHECK ("Bad NaN flag",
 341                                 MPFR_IS_NAN (y) && mpfr_nanflag_p ());
 342               else if (MPFR_IS_INF (y))
 343                 {
 344                   TGENERIC_CHECK ("Bad overflow flag",
 345                                   (compare != 0) ^ (mpfr_overflow_p () == 0));
 346                   TGENERIC_CHECK ("Bad divide-by-zero flag",
 347                                   (compare == 0 && !infinite_input) ^
 348                                   (mpfr_divby0_p () == 0));
 349                 }
 350               else if (MPFR_IS_ZERO (y))
 351                 TGENERIC_CHECK ("Bad underflow flag",
 352                                 (compare != 0) ^ (mpfr_underflow_p () == 0));
 353             }
 354           else if (mpfr_divby0_p ())
 355             {
 356               TGENERIC_CHECK ("Both overflow and divide-by-zero",
 357                               ! mpfr_overflow_p ());
 358               TGENERIC_CHECK ("Both underflow and divide-by-zero",
 359                               ! mpfr_underflow_p ());
 360               TGENERIC_CHECK ("Bad compare value (divide-by-zero)",
 361                               compare == 0);
 362             }
 363           else if (mpfr_overflow_p ())
 364             {
 365               TGENERIC_CHECK ("Both underflow and overflow",
 366                               ! mpfr_underflow_p ());
 367               TGENERIC_CHECK ("Bad compare value (overflow)", compare != 0);
 368               mpfr_nexttoinf (y);
 369               TGENERIC_CHECK ("Should have been max MPFR number",
 370                               MPFR_IS_INF (y));
 371             }
 372           else if (mpfr_underflow_p ())
 373             {
 374               TGENERIC_CHECK ("Bad compare value (underflow)", compare != 0);
 375               mpfr_nexttozero (y);
 376               TGENERIC_CHECK ("Should have been min MPFR number",
 377                               MPFR_IS_ZERO (y));
 378             }
 379           else if (mpfr_can_round (y, yprec, rnd, rnd, prec))
 380             {
 381               ctrn++;
 382               mpfr_set (t, y, rnd);
 383               /* Risk of failures are known when some flags are already set
 384                  before the function call. Do not set the erange flag, as
 385                  it will remain set after the function call and no checks
 386                  are performed in such a case (see the mpfr_erangeflag_p
 387                  test below). */
 388               if (randlimb () & 1)
 389                 __gmpfr_flags = MPFR_FLAGS_ALL ^ MPFR_FLAGS_ERANGE;
 390 #ifdef DEBUG_TGENERIC
 391               TGENERIC_INFO (TEST_FUNCTION, MPFR_PREC (z));
 392 #endif
 393               /* Let's increase the precision of the inputs in a random way.
 394                  In most cases, this doesn't make any difference, but for
 395                  the mpfr_fmod bug fixed in r6230, this triggers the bug. */
 396               mpfr_prec_round (x, mpfr_get_prec (x) + (randlimb () & 15),
 397                                MPFR_RNDN);
 398 #if defined(TWO_ARGS)
 399               mpfr_prec_round (u, mpfr_get_prec (u) + (randlimb () & 15),
 400                                MPFR_RNDN);
 401               inexact = TEST_FUNCTION (z, x, u, rnd);
 402 #elif defined(DOUBLE_ARG1)
 403               inexact = TEST_FUNCTION (z, d, x, rnd);
 404 #elif defined(DOUBLE_ARG2)
 405               inexact = TEST_FUNCTION (z, x, d, rnd);
 406 #else
 407               inexact = TEST_FUNCTION (z, x, rnd);
 408 #endif
 409               if (mpfr_erangeflag_p ())
 410                 goto next_n;
 411               if (mpfr_nan_p (z) || mpfr_cmp (t, z) != 0)
 412                 {
 413                   printf ("results differ for x=");
 414                   mpfr_out_str (stdout, 2, xprec, x, MPFR_RNDN);
 415 #ifdef TWO_ARGS
 416                   printf ("\nu=");
 417                   mpfr_out_str (stdout, 2, xprec, u, MPFR_RNDN);
 418 #elif defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
 419                   printf ("\nu=");
 420                   mpfr_out_str (stdout, 2, IEEE_DBL_MANT_DIG, u, MPFR_RNDN);
 421 #endif
 422                   printf (" prec=%u rnd_mode=%s\n", (unsigned) prec,
 423                           mpfr_print_rnd_mode (rnd));
 424                   printf ("got      ");
 425                   mpfr_out_str (stdout, 2, prec, z, MPFR_RNDN);
 426                   puts ("");
 427                   printf ("expected ");
 428                   mpfr_out_str (stdout, 2, prec, t, MPFR_RNDN);
 429                   puts ("");
 430                   printf ("approx   ");
 431                   mpfr_print_binary (y);
 432                   puts ("");
 433                   exit (1);
 434                 }
 435               compare2 = mpfr_cmp (t, y);
 436               /* if rounding to nearest, cannot know the sign of t - f(x)
 437                  because of composed rounding: y = o(f(x)) and t = o(y) */
 438               if (compare * compare2 >= 0)
 439                 compare = compare + compare2;
 440               else
 441                 compare = inexact; /* cannot determine sign(t-f(x)) */
 442               if (((inexact == 0) && (compare != 0)) ||
 443                   ((inexact > 0) && (compare <= 0)) ||
 444                   ((inexact < 0) && (compare >= 0)))
 445                 {
 446                   printf ("Wrong inexact flag for rnd=%s: expected %d, got %d"
 447                           "\n", mpfr_print_rnd_mode (rnd), compare, inexact);
 448                   printf ("x="); mpfr_print_binary (x); puts ("");
 449 #if defined(TWO_ARGS) || defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
 450                   printf ("u="); mpfr_print_binary (u); puts ("");
 451 #endif
 452                   printf ("y="); mpfr_print_binary (y); puts ("");
 453                   printf ("t="); mpfr_print_binary (t); puts ("");
 454                   exit (1);
 455                 }
 456             }
 457           else if (getenv ("MPFR_SUSPICIOUS_OVERFLOW") != NULL)
 458             {
 459               /* For developers only! */
 460               MPFR_ASSERTN (MPFR_IS_PURE_FP (y));
 461               mpfr_nexttoinf (y);
 462               if (MPFR_IS_INF (y) && MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG (y))
 463                   && !mpfr_overflow_p () && TGENERIC_SO_TEST)
 464                 {
 465                   printf ("Possible bug! |y| is the maximum finite number "
 466                           "and has been obtained when\nrounding toward zero"
 467                           " (%s). Thus there is a very probable overflow,\n"
 468                           "but the overflow flag is not set!\n",
 469                           mpfr_print_rnd_mode (rnd));
 470                   printf ("x="); mpfr_print_binary (x); puts ("");
 471 #if defined(TWO_ARGS) || defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
 472                   printf ("u="); mpfr_print_binary (u); puts ("");
 473 #endif
 474                   exit (1);
 475                 }
 476             }
 477
 478         next_n:
 479           /* In case the exponent range has been changed by
 480              tests_default_random() or for special values... */
 481           mpfr_set_emin (old_emin);
 482           mpfr_set_emax (old_emax);
 483         }
 484     }
 485
 486 #ifndef TGENERIC_NOWARNING
 487   if (3 * ctrn < 2 * ctrt)
 488     printf ("Warning! Too few normal cases in generic tests (%lu / %lu)\n",
 489             ctrn, ctrt);
 490 #endif
 491
 492   mpfr_clears (x, y, z, t, w, (mpfr_ptr) 0);
 493 #if defined(TWO_ARGS) || defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
 494   mpfr_clear (u);
 495 #endif
 496 }
 497
 498 #undef TEST_RANDOM_POS
 499 #undef TEST_RANDOM_POS2
 500 #undef TEST_RANDOM_EMIN
 501 #undef TEST_RANDOM_EMAX
 502 #undef RAND_FUNCTION
 503 #undef TWO_ARGS
 504 #undef TWO_ARGS_UI
 505 #undef TEST_FUNCTION
 506 #undef test_generic