math/libm-test-support.c

   1 /* Support code for testing libm functions (compiled once per type).
   2    Copyright (C) 1997-2021 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4
   5    The GNU C Library is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU Lesser General Public
   7    License as published by the Free Software Foundation; either
   8    version 2.1 of the License, or (at your option) any later version.
   9
  10    The GNU C Library is distributed in the hope that it will be useful,
  11    but WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    Lesser General Public License for more details.
  14
  15    You should have received a copy of the GNU Lesser General Public
  16    License along with the GNU C Library; if not, see
  17    <https://www.gnu.org/licenses/>.  */
  18
  19 /* Part of testsuite for libm.
  20
  21    libm-test-support.c contains functions shared by tests of different
  22    libm functions and types; it is compiled once per type.
  23    libm-test-driver.c defines the main function, and various variables
  24    that are used to configure the code in libm-test-support.c for
  25    different types and for variants such as testing inline functions.
  26
  27    The tests of individual functions are in .inc files processed by
  28    gen-libm-test.py, with the resulting files included together with
  29    libm-test-driver.c.
  30
  31    The per-type headers included both before libm-test-support.c and
  32    for the tests of individual functions must define the following
  33    macros:
  34
  35    FUNC(function): Convert general function name (like cos) to name
  36    with correct suffix (e.g. cosl or cosf).
  37
  38    FLOAT: Floating-point type to test.
  39
  40    BUILD_COMPLEX(real, imag): Create a complex number by calling a
  41    macro such as CMPLX.
  42
  43    PREFIX: The prefix for <float.h> macros for the type (e.g. LDBL,
  44    DBL, or FLT).
  45
  46    TYPE_STR: The name of the type as used in ulps files, as a string.
  47
  48    ULP_IDX: The array indexes for ulps values for this function.
  49
  50    LIT: Append the correct suffix to a literal.
  51
  52    LITM: Append the correct suffix to an M_* macro name.
  53
  54    FTOSTR: A function similar in type to strfromf which converts a
  55    FLOAT to a string.
  56
  57    snan_value_MACRO: The macro such as SNAN for a signaling NaN for
  58    the type.
  59
  60 */
  61
  62 /* Parameter handling is primitive in the moment:
  63    --verbose=[0..3] for different levels of output:
  64    0: only error count
  65    1: basic report on failed tests (default)
  66    2: full report on all tests
  67    -v for full output (equals --verbose=3)
  68    -u for generation of an ULPs file
  69  */
  70
  71 /* "Philosophy":
  72
  73    This suite tests some aspects of the correct implementation of
  74    mathematical functions in libm.  Some simple, specific parameters
  75    are tested for correctness but there's no exhaustive
  76    testing.  Handling of specific inputs (e.g. infinity, not-a-number)
  77    is also tested.  Correct handling of exceptions is checked
  78    against.  These implemented tests should check all cases that are
  79    specified in ISO C99.
  80
  81    NaN values: The payload of NaNs is set in inputs for functions
  82    where it is significant, and is examined in the outputs of some
  83    functions.
  84
  85    Inline functions: Inlining functions should give an improvement in
  86    speed - but not in precission.  The inlined functions return
  87    reasonable values for a reasonable range of input values.  The
  88    result is not necessarily correct for all values and exceptions are
  89    not correctly raised in all cases.  Problematic input and return
  90    values are infinity, not-a-number and minus zero.  This suite
  91    therefore does not check these specific inputs and the exception
  92    handling for inlined mathematical functions - just the "reasonable"
  93    values are checked.
  94
  95    Beware: The tests might fail for any of the following reasons:
  96    - Tests are wrong
  97    - Functions are wrong
  98    - Floating Point Unit not working properly
  99    - Compiler has errors
 100
 101    With e.g. gcc 2.7.2.2 the test for cexp fails because of a compiler error.
 102
 103
 104    To Do: All parameter should be numbers that can be represented as
 105    exact floating point values.  Currently some values cannot be
 106    represented exactly and therefore the result is not the expected
 107    result.  For this we will use 36 digits so that numbers can be
 108    represented exactly.  */
 109
 110 #include "libm-test-support.h"
 111
 112 #include <argp.h>
 113 #include <errno.h>
 114 #include <string.h>
 115
 116 /* This header defines func_ulps, func_real_ulps and func_imag_ulps
 117    arrays.  */
 118 #include "libm-test-ulps.h"
 119
 120 /* Maximum character buffer to store a stringitized FLOAT value.  */
 121 #define FSTR_MAX (128)
 122
 123 #define ulps_file_name "ULPs"   /* Name of the ULPs file.  */
 124 static FILE *ulps_file;         /* File to document difference.  */
 125 static int output_ulps;         /* Should ulps printed?  */
 126 static char *output_dir;        /* Directory where generated files will be written.  */
 127
 128 static int noErrors;    /* number of errors */
 129 static int noTests;     /* number of tests (without testing exceptions) */
 130 static int noExcTests;  /* number of tests for exception flags */
 131 static int noErrnoTests;/* number of tests for errno values */
 132
 133 static int verbose;
 134 static int output_max_error;    /* Should the maximal errors printed?  */
 135 static int output_points;       /* Should the single function results printed?  */
 136 static int ignore_max_ulp;      /* Should we ignore max_ulp?  */
 137 static int test_ibm128;         /* Is argument or result IBM long double?  */
 138
 139 static FLOAT max_error, real_max_error, imag_max_error;
 140
 141 static FLOAT prev_max_error, prev_real_max_error, prev_imag_max_error;
 142
 143 static FLOAT max_valid_error;
 144
 145 /* Sufficient numbers of digits to represent any floating-point value
 146    unambiguously (for any choice of the number of bits in the first
 147    hex digit, in the case of TYPE_HEX_DIG).  When used with printf
 148    formats where the precision counts only digits after the point, 1
 149    is subtracted from these values. */
 150 #define TYPE_DECIMAL_DIG __CONCATX (PREFIX, _DECIMAL_DIG)
 151 #define TYPE_HEX_DIG ((MANT_DIG + 6) / 4)
 152
 153 /* Converts VALUE (a floating-point number) to string and writes it to DEST.
 154    PRECISION specifies the number of fractional digits that should be printed.
 155    CONVERSION is the conversion specifier, such as in printf, e.g. 'f' or 'a'.
 156    The output is prepended with an empty space if VALUE is non-negative.  */
 157 static void
 158 fmt_ftostr (char *dest, size_t size, int precision, const char *conversion,
 159             FLOAT value)
 160 {
 161   char format[64];
 162   char *ptr_format;
 163   int ret;
 164
 165   /* Generate the format string.  */
 166   ptr_format = stpcpy (format, "%.");
 167   ret = sprintf (ptr_format, "%d", precision);
 168   ptr_format += ret;
 169   ptr_format = stpcpy (ptr_format, conversion);
 170
 171   /* Add a space to the beginning of the output string, if the floating-point
 172      number is non-negative.  This mimics the behavior of the space (' ') flag
 173      in snprintf, which is not available on strfrom.  */
 174   if (! signbit (value))
 175     {
 176       *dest = ' ';
 177       dest++;
 178       size--;
 179     }
 180
 181   /* Call the float to string conversion function, e.g.: strfromd.  */
 182   FTOSTR (dest, size, format, value);
 183 }
 184
 185 /* Compare KEY (a string, with the name of a function) with ULP (a
 186    pointer to a struct ulp_data structure), returning a value less
 187    than, equal to or greater than zero for use in bsearch.  */
 188
 189 static int
 190 compare_ulp_data (const void *key, const void *ulp)
 191 {
 192   const char *keystr = key;
 193   const struct ulp_data *ulpdat = ulp;
 194   return strcmp (keystr, ulpdat->name);
 195 }
 196
 197 static const int ulp_idx = ULP_IDX;
 198
 199 /* Return the ulps for NAME in array DATA with NMEMB elements, or 0 if
 200    no ulps listed.  */
 201
 202 static FLOAT
 203 find_ulps (const char *name, const struct ulp_data *data, size_t nmemb)
 204 {
 205   const struct ulp_data *entry = bsearch (name, data, nmemb, sizeof (*data),
 206                                           compare_ulp_data);
 207   if (entry == NULL)
 208     return 0;
 209   else
 210     return entry->max_ulp[ulp_idx];
 211 }
 212
 213 void
 214 init_max_error (const char *name, int exact, int testing_ibm128)
 215 {
 216   max_error = 0;
 217   real_max_error = 0;
 218   imag_max_error = 0;
 219   test_ibm128 = testing_ibm128;
 220   prev_max_error = find_ulps (name, func_ulps,
 221                               sizeof (func_ulps) / sizeof (func_ulps[0]));
 222   prev_real_max_error = find_ulps (name, func_real_ulps,
 223                                    (sizeof (func_real_ulps)
 224                                     / sizeof (func_real_ulps[0])));
 225   prev_imag_max_error = find_ulps (name, func_imag_ulps,
 226                                    (sizeof (func_imag_ulps)
 227                                     / sizeof (func_imag_ulps[0])));
 228   if (testing_ibm128)
 229     /* The documented accuracy of IBM long double division is 3ulp
 230        (see libgcc/config/rs6000/ibm-ldouble-format), so do not
 231        require better accuracy for libm functions that are exactly
 232        defined for other formats.  */
 233     max_valid_error = exact ? 3 : 16;
 234   else
 235     max_valid_error = exact ? 0 : 9;
 236   prev_max_error = (prev_max_error <= max_valid_error
 237                     ? prev_max_error
 238                     : max_valid_error);
 239   prev_real_max_error = (prev_real_max_error <= max_valid_error
 240                          ? prev_real_max_error
 241                          : max_valid_error);
 242   prev_imag_max_error = (prev_imag_max_error <= max_valid_error
 243                          ? prev_imag_max_error
 244                          : max_valid_error);
 245   feclearexcept (FE_ALL_EXCEPT);
 246   errno = 0;
 247 }
 248
 249 static void
 250 set_max_error (FLOAT current, FLOAT *curr_max_error)
 251 {
 252   if (current > *curr_max_error && current <= max_valid_error)
 253     *curr_max_error = current;
 254 }
 255
 256
 257 /* Print a FLOAT.  */
 258 static void
 259 print_float (FLOAT f)
 260 {
 261   /* As printf doesn't differ between a sNaN and a qNaN, do this manually.  */
 262   if (issignaling (f))
 263     printf ("sNaN\n");
 264   else if (isnan (f))
 265     printf ("qNaN\n");
 266   else
 267     {
 268       char fstrn[FSTR_MAX], fstrx[FSTR_MAX];
 269       fmt_ftostr (fstrn, FSTR_MAX, TYPE_DECIMAL_DIG - 1, "e", f);
 270       fmt_ftostr (fstrx, FSTR_MAX, TYPE_HEX_DIG - 1, "a", f);
 271       printf ("%s  %s\n", fstrn, fstrx);
 272     }
 273 }
 274
 275 /* Should the message print to screen?  This depends on the verbose flag,
 276    and the test status.  */
 277 static int
 278 print_screen (int ok)
 279 {
 280   if (output_points
 281       && (verbose > 1
 282           || (verbose == 1 && ok == 0)))
 283     return 1;
 284   return 0;
 285 }
 286
 287
 288 /* Should the message print to screen?  This depends on the verbose flag,
 289    and the test status.  */
 290 static int
 291 print_screen_max_error (int ok)
 292 {
 293   if (output_max_error
 294       && (verbose > 1
 295           || ((verbose == 1) && (ok == 0))))
 296     return 1;
 297   return 0;
 298 }
 299
 300 /* Update statistic counters.  */
 301 static void
 302 update_stats (int ok)
 303 {
 304   ++noTests;
 305   if (!ok)
 306     ++noErrors;
 307 }
 308
 309 static void
 310 print_function_ulps (const char *function_name, FLOAT ulp)
 311 {
 312   if (output_ulps)
 313     {
 314       char ustrn[FSTR_MAX];
 315       FTOSTR (ustrn, FSTR_MAX, "%.0f", FUNC (ceil) (ulp));
 316       fprintf (ulps_file, "Function: \"%s\":\n", function_name);
 317       fprintf (ulps_file, "%s: %s\n", qtype_str, ustrn);
 318     }
 319 }
 320
 321
 322 static void
 323 print_complex_function_ulps (const char *function_name, FLOAT real_ulp,
 324                              FLOAT imag_ulp)
 325 {
 326   if (output_ulps)
 327     {
 328       char fstrn[FSTR_MAX];
 329       if (real_ulp != 0.0)
 330         {
 331           FTOSTR (fstrn, FSTR_MAX, "%.0f", FUNC (ceil) (real_ulp));
 332           fprintf (ulps_file, "Function: Real part of \"%s\":\n", function_name);
 333           fprintf (ulps_file, "%s: %s\n", qtype_str, fstrn);
 334         }
 335       if (imag_ulp != 0.0)
 336         {
 337           FTOSTR (fstrn, FSTR_MAX, "%.0f", FUNC (ceil) (imag_ulp));
 338           fprintf (ulps_file, "Function: Imaginary part of \"%s\":\n", function_name);
 339           fprintf (ulps_file, "%s: %s\n", qtype_str, fstrn);
 340         }
 341
 342
 343     }
 344 }
 345
 346
 347
 348 /* Test if Floating-Point stack hasn't changed */
 349 static void
 350 fpstack_test (const char *test_name)
 351 {
 352 #if defined (__i386__) || defined (__x86_64__)
 353   static int old_stack;
 354   int sw;
 355
 356   asm ("fnstsw" : "=a" (sw));
 357   sw >>= 11;
 358   sw &= 7;
 359
 360   if (sw != old_stack)
 361     {
 362       printf ("FP-Stack wrong after test %s (%d, should be %d)\n",
 363               test_name, sw, old_stack);
 364       ++noErrors;
 365       old_stack = sw;
 366     }
 367 #endif
 368 }
 369
 370
 371 void
 372 print_max_error (const char *func_name)
 373 {
 374   int ok = 0;
 375
 376   if (max_error == 0.0 || (max_error <= prev_max_error && !ignore_max_ulp))
 377     {
 378       ok = 1;
 379     }
 380
 381   if (!ok)
 382     print_function_ulps (func_name, max_error);
 383
 384
 385   if (print_screen_max_error (ok))
 386     {
 387       char mestr[FSTR_MAX], pmestr[FSTR_MAX];
 388       FTOSTR (mestr, FSTR_MAX, "%.0f", FUNC (ceil) (max_error));
 389       FTOSTR (pmestr, FSTR_MAX, "%.0f", FUNC (ceil) (prev_max_error));
 390       printf ("Maximal error of `%s'\n", func_name);
 391       printf (" is      : %s ulp\n", mestr);
 392       printf (" accepted: %s ulp\n", pmestr);
 393     }
 394
 395   update_stats (ok);
 396 }
 397
 398
 399 void
 400 print_complex_max_error (const char *func_name)
 401 {
 402   int real_ok = 0, imag_ok = 0, ok;
 403
 404   if (real_max_error == 0
 405       || (real_max_error <= prev_real_max_error && !ignore_max_ulp))
 406     {
 407       real_ok = 1;
 408     }
 409
 410   if (imag_max_error == 0
 411       || (imag_max_error <= prev_imag_max_error && !ignore_max_ulp))
 412     {
 413       imag_ok = 1;
 414     }
 415
 416   ok = real_ok && imag_ok;
 417
 418   if (!ok)
 419     print_complex_function_ulps (func_name,
 420                                  real_ok ? 0 : real_max_error,
 421                                  imag_ok ? 0 : imag_max_error);
 422
 423   if (print_screen_max_error (ok))
 424     {
 425       char rmestr[FSTR_MAX], prmestr[FSTR_MAX];
 426       char imestr[FSTR_MAX], pimestr[FSTR_MAX];
 427       FTOSTR (rmestr, FSTR_MAX, "%.0f", FUNC (ceil) (real_max_error));
 428       FTOSTR (prmestr, FSTR_MAX, "%.0f", FUNC (ceil) (prev_real_max_error));
 429       FTOSTR (imestr, FSTR_MAX, "%.0f", FUNC (ceil) (imag_max_error));
 430       FTOSTR (pimestr, FSTR_MAX, "%.0f", FUNC (ceil) (prev_imag_max_error));
 431       printf ("Maximal error of real part of: %s\n", func_name);
 432       printf (" is      : %s ulp\n", rmestr);
 433       printf (" accepted: %s ulp\n", prmestr);
 434       printf ("Maximal error of imaginary part of: %s\n", func_name);
 435       printf (" is      : %s ulp\n", imestr);
 436       printf (" accepted: %s ulp\n", pimestr);
 437     }
 438
 439   update_stats (ok);
 440 }
 441
 442
 443 #if FE_ALL_EXCEPT
 444 /* Test whether a given exception was raised.  */
 445 static void
 446 test_single_exception (const char *test_name,
 447                        int exception,
 448                        int exc_flag,
 449                        int fe_flag,
 450                        const char *flag_name)
 451 {
 452   int ok = 1;
 453   if (exception & exc_flag)
 454     {
 455       if (fetestexcept (fe_flag))
 456         {
 457           if (print_screen (1))
 458             printf ("Pass: %s: Exception \"%s\" set\n", test_name, flag_name);
 459         }
 460       else
 461         {
 462           ok = 0;
 463           if (print_screen (0))
 464             printf ("Failure: %s: Exception \"%s\" not set\n",
 465                     test_name, flag_name);
 466         }
 467     }
 468   else
 469     {
 470       if (fetestexcept (fe_flag))
 471         {
 472           ok = 0;
 473           if (print_screen (0))
 474             printf ("Failure: %s: Exception \"%s\" set\n",
 475                     test_name, flag_name);
 476         }
 477       else
 478         {
 479           if (print_screen (1))
 480             printf ("%s: Exception \"%s\" not set\n", test_name,
 481                     flag_name);
 482         }
 483     }
 484   if (!ok)
 485     ++noErrors;
 486 }
 487 #endif
 488
 489 /* Test whether exceptions given by EXCEPTION are raised.  Ignore thereby
 490    allowed but not required exceptions.
 491 */
 492 static void
 493 test_exceptions (const char *test_name, int exception)
 494 {
 495   if (flag_test_exceptions && EXCEPTION_TESTS (FLOAT))
 496     {
 497       ++noExcTests;
 498 #ifdef FE_DIVBYZERO
 499       if ((exception & DIVIDE_BY_ZERO_EXCEPTION_OK) == 0)
 500         test_single_exception (test_name, exception,
 501                                DIVIDE_BY_ZERO_EXCEPTION, FE_DIVBYZERO,
 502                                "Divide by zero");
 503 #endif
 504 #ifdef FE_INVALID
 505       if ((exception & INVALID_EXCEPTION_OK) == 0)
 506         test_single_exception (test_name, exception,
 507                                INVALID_EXCEPTION, FE_INVALID,
 508                                "Invalid operation");
 509 #endif
 510 #ifdef FE_OVERFLOW
 511       if ((exception & OVERFLOW_EXCEPTION_OK) == 0)
 512         test_single_exception (test_name, exception, OVERFLOW_EXCEPTION,
 513                                FE_OVERFLOW, "Overflow");
 514 #endif
 515       /* Spurious "underflow" and "inexact" exceptions are always
 516          allowed for IBM long double, in line with the underlying
 517          arithmetic.  */
 518 #ifdef FE_UNDERFLOW
 519       if ((exception & UNDERFLOW_EXCEPTION_OK) == 0
 520           && !(test_ibm128
 521                && (exception & UNDERFLOW_EXCEPTION) == 0))
 522         test_single_exception (test_name, exception, UNDERFLOW_EXCEPTION,
 523                                FE_UNDERFLOW, "Underflow");
 524 #endif
 525 #ifdef FE_INEXACT
 526       if ((exception & (INEXACT_EXCEPTION | NO_INEXACT_EXCEPTION)) != 0
 527           && !(test_ibm128
 528                && (exception & NO_INEXACT_EXCEPTION) != 0))
 529         test_single_exception (test_name, exception, INEXACT_EXCEPTION,
 530                                FE_INEXACT, "Inexact");
 531 #endif
 532     }
 533   feclearexcept (FE_ALL_EXCEPT);
 534 }
 535
 536 /* Test whether errno for TEST_NAME, set to ERRNO_VALUE, has value
 537    EXPECTED_VALUE (description EXPECTED_NAME).  */
 538 static void
 539 test_single_errno (const char *test_name, int errno_value,
 540                    int expected_value, const char *expected_name)
 541 {
 542   if (errno_value == expected_value)
 543     {
 544       if (print_screen (1))
 545         printf ("Pass: %s: errno set to %d (%s)\n", test_name, errno_value,
 546                 expected_name);
 547     }
 548   else
 549     {
 550       ++noErrors;
 551       if (print_screen (0))
 552         printf ("Failure: %s: errno set to %d, expected %d (%s)\n",
 553                 test_name, errno_value, expected_value, expected_name);
 554     }
 555 }
 556
 557 /* Test whether errno (value ERRNO_VALUE) has been for TEST_NAME set
 558    as required by EXCEPTIONS.  */
 559 static void
 560 test_errno (const char *test_name, int errno_value, int exceptions)
 561 {
 562   if (flag_test_errno)
 563     {
 564       ++noErrnoTests;
 565       if (exceptions & ERRNO_UNCHANGED)
 566         test_single_errno (test_name, errno_value, 0, "unchanged");
 567       if (exceptions & ERRNO_EDOM)
 568         test_single_errno (test_name, errno_value, EDOM, "EDOM");
 569       if (exceptions & ERRNO_ERANGE)
 570         test_single_errno (test_name, errno_value, ERANGE, "ERANGE");
 571     }
 572 }
 573
 574 /* Returns the number of ulps that GIVEN is away from EXPECTED.  */
 575 #define ULPDIFF(given, expected) \
 576         (FUNC(fabs) ((given) - (expected)) / ulp (expected))
 577
 578 /* Returns the size of an ulp for VALUE.  */
 579 static FLOAT
 580 ulp (FLOAT value)
 581 {
 582   FLOAT ulp;
 583
 584   switch (fpclassify (value))
 585     {
 586       case FP_ZERO:
 587         /* We compute the distance to the next FP which is the same as the
 588            value of the smallest subnormal number. Previously we used
 589            2^-(MANT_DIG - 1) which is too large a value to be useful. Note that we
 590            can't use ilogb(0), since that isn't a valid thing to do. As a point
 591            of comparison Java's ulp returns the next normal value e.g.
 592            2^(1 - MAX_EXP) for ulp(0), but that is not what we want for
 593            glibc.  */
 594         /* Fall through...  */
 595       case FP_SUBNORMAL:
 596         /* The next closest subnormal value is a constant distance away.  */
 597         ulp = FUNC(ldexp) (1.0, MIN_EXP - MANT_DIG);
 598         break;
 599
 600       case FP_NORMAL:
 601         ulp = FUNC(ldexp) (1.0, FUNC(ilogb) (value) - MANT_DIG + 1);
 602         break;
 603
 604       default:
 605         /* It should never happen. */
 606         abort ();
 607         break;
 608     }
 609   return ulp;
 610 }
 611
 612 static void
 613 check_float_internal (const char *test_name, FLOAT computed, FLOAT expected,
 614                       int exceptions,
 615                       FLOAT *curr_max_error, FLOAT max_ulp)
 616 {
 617   int ok = 0;
 618   int print_diff = 0;
 619   FLOAT diff = 0;
 620   FLOAT ulps = 0;
 621   int errno_value = errno;
 622
 623   test_exceptions (test_name, exceptions);
 624   test_errno (test_name, errno_value, exceptions);
 625   if (exceptions & IGNORE_RESULT)
 626     goto out;
 627   if (issignaling (computed) && issignaling (expected))
 628     {
 629       if ((exceptions & TEST_NAN_SIGN) != 0
 630           && signbit (computed) != signbit (expected))
 631         {
 632           ok = 0;
 633           printf ("signaling NaN has wrong sign.\n");
 634         }
 635       else if ((exceptions & TEST_NAN_PAYLOAD) != 0
 636                && (FUNC (getpayload) (&computed)
 637                    != FUNC (getpayload) (&expected)))
 638         {
 639           ok = 0;
 640           printf ("signaling NaN has wrong payload.\n");
 641         }
 642       else
 643         ok = 1;
 644     }
 645   else if (issignaling (computed) || issignaling (expected))
 646     ok = 0;
 647   else if (isnan (computed) && isnan (expected))
 648     {
 649       if ((exceptions & TEST_NAN_SIGN) != 0
 650           && signbit (computed) != signbit (expected))
 651         {
 652           ok = 0;
 653           printf ("quiet NaN has wrong sign.\n");
 654         }
 655       else if ((exceptions & TEST_NAN_PAYLOAD) != 0
 656                && (FUNC (getpayload) (&computed)
 657                    != FUNC (getpayload) (&expected)))
 658         {
 659           ok = 0;
 660           printf ("quiet NaN has wrong payload.\n");
 661         }
 662       else
 663         ok = 1;
 664     }
 665   else if (isinf (computed) && isinf (expected))
 666     {
 667       /* Test for sign of infinities.  */
 668       if ((exceptions & IGNORE_ZERO_INF_SIGN) == 0
 669           && signbit (computed) != signbit (expected))
 670         {
 671           ok = 0;
 672           printf ("infinity has wrong sign.\n");
 673         }
 674       else
 675         ok = 1;
 676     }
 677   /* Don't calculate ULPs for infinities or any kind of NaNs.  */
 678   else if (isinf (computed) || isnan (computed)
 679            || isinf (expected) || isnan (expected))
 680     ok = 0;
 681   else
 682     {
 683       diff = FUNC(fabs) (computed - expected);
 684       ulps = ULPDIFF (computed, expected);
 685       set_max_error (ulps, curr_max_error);
 686       print_diff = 1;
 687       if ((exceptions & IGNORE_ZERO_INF_SIGN) == 0
 688           && computed == 0.0 && expected == 0.0
 689           && signbit(computed) != signbit (expected))
 690         ok = 0;
 691       else if (ulps <= max_ulp && !ignore_max_ulp)
 692         ok = 1;
 693       else
 694         ok = 0;
 695     }
 696   if (print_screen (ok))
 697     {
 698       if (!ok)
 699         printf ("Failure: ");
 700       printf ("Test: %s\n", test_name);
 701       printf ("Result:\n");
 702       printf (" is:         ");
 703       print_float (computed);
 704       printf (" should be:  ");
 705       print_float (expected);
 706       if (print_diff)
 707         {
 708           char dstrn[FSTR_MAX], dstrx[FSTR_MAX];
 709           char ustrn[FSTR_MAX], mustrn[FSTR_MAX];
 710           fmt_ftostr (dstrn, FSTR_MAX, TYPE_DECIMAL_DIG - 1, "e", diff);
 711           fmt_ftostr (dstrx, FSTR_MAX, TYPE_HEX_DIG - 1, "a", diff);
 712           fmt_ftostr (ustrn, FSTR_MAX, 4, "f", ulps);
 713           fmt_ftostr (mustrn, FSTR_MAX, 4, "f", max_ulp);
 714           printf (" difference: %s  %s\n", dstrn, dstrx);
 715           printf (" ulp       : %s\n", ustrn);
 716           printf (" max.ulp   : %s\n", mustrn);
 717         }
 718     }
 719   update_stats (ok);
 720
 721  out:
 722   fpstack_test (test_name);
 723   feclearexcept (FE_ALL_EXCEPT);
 724   errno = 0;
 725 }
 726
 727
 728 void
 729 check_float (const char *test_name, FLOAT computed, FLOAT expected,
 730              int exceptions)
 731 {
 732   check_float_internal (test_name, computed, expected,
 733                         exceptions, &max_error, prev_max_error);
 734 }
 735
 736
 737 void
 738 check_complex (const char *test_name, CFLOAT computed,
 739                CFLOAT expected,
 740                int exception)
 741 {
 742   FLOAT part_comp, part_exp;
 743   char *str;
 744
 745   if (asprintf (&str, "Real part of: %s", test_name) == -1)
 746     abort ();
 747
 748   part_comp = __real__ computed;
 749   part_exp = __real__ expected;
 750
 751   check_float_internal (str, part_comp, part_exp,
 752                         exception, &real_max_error, prev_real_max_error);
 753   free (str);
 754
 755   if (asprintf (&str, "Imaginary part of: %s", test_name) == -1)
 756     abort ();
 757
 758   part_comp = __imag__ computed;
 759   part_exp = __imag__ expected;
 760
 761   /* Don't check again for exceptions or errno, just pass through the
 762      other relevant flags.  */
 763   check_float_internal (str, part_comp, part_exp,
 764                         exception & (IGNORE_ZERO_INF_SIGN
 765                                      | TEST_NAN_SIGN
 766                                      | IGNORE_RESULT),
 767                         &imag_max_error, prev_imag_max_error);
 768   free (str);
 769 }
 770
 771
 772 /* Check that computed and expected values are equal (int values).  */
 773 void
 774 check_int (const char *test_name, int computed, int expected,
 775            int exceptions)
 776 {
 777   int ok = 0;
 778   int errno_value = errno;
 779
 780   test_exceptions (test_name, exceptions);
 781   test_errno (test_name, errno_value, exceptions);
 782   if (exceptions & IGNORE_RESULT)
 783     goto out;
 784   noTests++;
 785   if (computed == expected)
 786     ok = 1;
 787
 788   if (print_screen (ok))
 789     {
 790       if (!ok)
 791         printf ("Failure: ");
 792       printf ("Test: %s\n", test_name);
 793       printf ("Result:\n");
 794       printf (" is:         %d\n", computed);
 795       printf (" should be:  %d\n", expected);
 796     }
 797
 798   update_stats (ok);
 799  out:
 800   fpstack_test (test_name);
 801   feclearexcept (FE_ALL_EXCEPT);
 802   errno = 0;
 803 }
 804
 805
 806 /* Check that computed and expected values are equal (long int values).  */
 807 void
 808 check_long (const char *test_name, long int computed, long int expected,
 809             int exceptions)
 810 {
 811   int ok = 0;
 812   int errno_value = errno;
 813
 814   test_exceptions (test_name, exceptions);
 815   test_errno (test_name, errno_value, exceptions);
 816   if (exceptions & IGNORE_RESULT)
 817     goto out;
 818   noTests++;
 819   if (computed == expected)
 820     ok = 1;
 821
 822   if (print_screen (ok))
 823     {
 824       if (!ok)
 825         printf ("Failure: ");
 826       printf ("Test: %s\n", test_name);
 827       printf ("Result:\n");
 828       printf (" is:         %ld\n", computed);
 829       printf (" should be:  %ld\n", expected);
 830     }
 831
 832   update_stats (ok);
 833  out:
 834   fpstack_test (test_name);
 835   feclearexcept (FE_ALL_EXCEPT);
 836   errno = 0;
 837 }
 838
 839
 840 /* Check that computed value is true/false.  */
 841 void
 842 check_bool (const char *test_name, int computed, int expected,
 843             int exceptions)
 844 {
 845   int ok = 0;
 846   int errno_value = errno;
 847
 848   test_exceptions (test_name, exceptions);
 849   test_errno (test_name, errno_value, exceptions);
 850   if (exceptions & IGNORE_RESULT)
 851     goto out;
 852   noTests++;
 853   if ((computed == 0) == (expected == 0))
 854     ok = 1;
 855
 856   if (print_screen (ok))
 857     {
 858       if (!ok)
 859         printf ("Failure: ");
 860       printf ("Test: %s\n", test_name);
 861       printf ("Result:\n");
 862       printf (" is:         %d\n", computed);
 863       printf (" should be:  %d\n", expected);
 864     }
 865
 866   update_stats (ok);
 867  out:
 868   fpstack_test (test_name);
 869   feclearexcept (FE_ALL_EXCEPT);
 870   errno = 0;
 871 }
 872
 873
 874 /* check that computed and expected values are equal (long int values) */
 875 void
 876 check_longlong (const char *test_name, long long int computed,
 877                 long long int expected,
 878                 int exceptions)
 879 {
 880   int ok = 0;
 881   int errno_value = errno;
 882
 883   test_exceptions (test_name, exceptions);
 884   test_errno (test_name, errno_value, exceptions);
 885   if (exceptions & IGNORE_RESULT)
 886     goto out;
 887   noTests++;
 888   if (computed == expected)
 889     ok = 1;
 890
 891   if (print_screen (ok))
 892     {
 893       if (!ok)
 894         printf ("Failure:");
 895       printf ("Test: %s\n", test_name);
 896       printf ("Result:\n");
 897       printf (" is:         %lld\n", computed);
 898       printf (" should be:  %lld\n", expected);
 899     }
 900
 901   update_stats (ok);
 902  out:
 903   fpstack_test (test_name);
 904   feclearexcept (FE_ALL_EXCEPT);
 905   errno = 0;
 906 }
 907
 908
 909 /* Check that computed and expected values are equal (intmax_t values).  */
 910 void
 911 check_intmax_t (const char *test_name, intmax_t computed,
 912                 intmax_t expected, int exceptions)
 913 {
 914   int ok = 0;
 915   int errno_value = errno;
 916
 917   test_exceptions (test_name, exceptions);
 918   test_errno (test_name, errno_value, exceptions);
 919   if (exceptions & IGNORE_RESULT)
 920     goto out;
 921   noTests++;
 922   if (computed == expected)
 923     ok = 1;
 924
 925   if (print_screen (ok))
 926     {
 927       if (!ok)
 928         printf ("Failure:");
 929       printf ("Test: %s\n", test_name);
 930       printf ("Result:\n");
 931       printf (" is:         %jd\n", computed);
 932       printf (" should be:  %jd\n", expected);
 933     }
 934
 935   update_stats (ok);
 936  out:
 937   fpstack_test (test_name);
 938   feclearexcept (FE_ALL_EXCEPT);
 939   errno = 0;
 940 }
 941
 942
 943 /* Check that computed and expected values are equal (uintmax_t values).  */
 944 void
 945 check_uintmax_t (const char *test_name, uintmax_t computed,
 946                  uintmax_t expected, int exceptions)
 947 {
 948   int ok = 0;
 949   int errno_value = errno;
 950
 951   test_exceptions (test_name, exceptions);
 952   test_errno (test_name, errno_value, exceptions);
 953   if (exceptions & IGNORE_RESULT)
 954     goto out;
 955   noTests++;
 956   if (computed == expected)
 957     ok = 1;
 958
 959   if (print_screen (ok))
 960     {
 961       if (!ok)
 962         printf ("Failure:");
 963       printf ("Test: %s\n", test_name);
 964       printf ("Result:\n");
 965       printf (" is:         %ju\n", computed);
 966       printf (" should be:  %ju\n", expected);
 967     }
 968
 969   update_stats (ok);
 970  out:
 971   fpstack_test (test_name);
 972   feclearexcept (FE_ALL_EXCEPT);
 973   errno = 0;
 974 }
 975
 976 /* Return whether a test with flags EXCEPTIONS should be run.  */
 977 int
 978 enable_test (int exceptions)
 979 {
 980   if (exceptions & XFAIL_TEST)
 981     return 0;
 982   if ((!SNAN_TESTS (FLOAT) || !snan_tests_arg)
 983       && (exceptions & TEST_SNAN) != 0)
 984     return 0;
 985   if (flag_test_mathvec && (exceptions & NO_TEST_MATHVEC) != 0)
 986     return 0;
 987
 988   return 1;
 989 }
 990
 991 static void
 992 initialize (void)
 993 {
 994   fpstack_test ("start *init*");
 995
 996   /* Clear all exceptions.  From now on we must not get random exceptions.  */
 997   feclearexcept (FE_ALL_EXCEPT);
 998   errno = 0;
 999
1000   /* Test to make sure we start correctly.  */
1001   fpstack_test ("end *init*");
1002 }
1003
1004 /* Definitions of arguments for argp functions.  */
1005 static const struct argp_option options[] =
1006 {
1007   { "verbose", 'v', "NUMBER", 0, "Level of verbosity (0..3)"},
1008   { "ulps-file", 'u', NULL, 0, "Output ulps to file ULPs"},
1009   { "no-max-error", 'f', NULL, 0,
1010     "Don't output maximal errors of functions"},
1011   { "no-points", 'p', NULL, 0,
1012     "Don't output results of functions invocations"},
1013   { "ignore-max-ulp", 'i', "yes/no", 0,
1014     "Ignore given maximal errors"},
1015   { "output-dir", 'o', "DIR", 0,
1016     "Directory where generated files will be placed"},
1017   { NULL, 0, NULL, 0, NULL }
1018 };
1019
1020 /* Prototype for option handler.  */
1021 static error_t parse_opt (int key, char *arg, struct argp_state *state);
1022
1023 /* Data structure to communicate with argp functions.  */
1024 static struct argp argp =
1025 {
1026   options, parse_opt, NULL, doc,
1027 };
1028
1029
1030 /* Handle program arguments.  */
1031 static error_t
1032 parse_opt (int key, char *arg, struct argp_state *state)
1033 {
1034   switch (key)
1035     {
1036     case 'f':
1037       output_max_error = 0;
1038       break;
1039     case 'i':
1040       if (strcmp (arg, "yes") == 0)
1041         ignore_max_ulp = 1;
1042       else if (strcmp (arg, "no") == 0)
1043         ignore_max_ulp = 0;
1044       break;
1045     case 'o':
1046       output_dir = (char *) malloc (strlen (arg) + 1);
1047       if (output_dir != NULL)
1048         strcpy (output_dir, arg);
1049       else
1050         return errno;
1051       break;
1052     case 'p':
1053       output_points = 0;
1054       break;
1055     case 'u':
1056       output_ulps = 1;
1057       break;
1058     case 'v':
1059       if (optarg)
1060         verbose = (unsigned int) strtoul (optarg, NULL, 0);
1061       else
1062         verbose = 3;
1063       break;
1064     default:
1065       return ARGP_ERR_UNKNOWN;
1066     }
1067   return 0;
1068 }
1069
1070 /* Verify that our ulp () implementation is behaving as expected
1071    or abort.  */
1072 static void
1073 check_ulp (void)
1074 {
1075    FLOAT ulps, ulpx, value;
1076    int i;
1077    /* Check ulp of zero is a subnormal value...  */
1078    ulps = ulp (0x0.0p0);
1079    if (fpclassify (ulps) != FP_SUBNORMAL)
1080      {
1081        fprintf (stderr, "ulp (0x0.0p0) is not FP_SUBNORMAL!\n");
1082        exit (EXIT_FAILURE);
1083      }
1084    /* Check that the ulp of one is a normal value... */
1085    ulps = ulp (LIT(1.0));
1086    if (fpclassify (ulps) != FP_NORMAL)
1087      {
1088        fprintf (stderr, "ulp (1.0L) is not FP_NORMAL\n");
1089        exit (EXIT_FAILURE);
1090      }
1091
1092    /* Compute the next subnormal value using nextafter to validate ulp.
1093       We allow +/- 1 ulp around the represented value.  */
1094    value = FUNC(nextafter) (0, 1);
1095    ulps = ULPDIFF (value, 0);
1096    ulpx = ulp (LIT(1.0));
1097    if (ulps < (LIT(1.0) - ulpx) || ulps > (LIT(1.0) + ulpx))
1098      {
1099        fprintf (stderr, "Value outside of 1 +/- 1ulp.\n");
1100        exit (EXIT_FAILURE);
1101      }
1102    /* Compute the nearest representable number from 10 towards 20.
1103       The result is 10 + 1ulp.  We use this to check the ulp function.
1104       We allow +/- 1 ulp around the represented value.  */
1105    value = FUNC(nextafter) (10, 20);
1106    ulps = ULPDIFF (value, 10);
1107    ulpx = ulp (LIT(1.0));
1108    if (ulps < (LIT(1.0) - ulpx) || ulps > (LIT(1.0) + ulpx))
1109      {
1110        fprintf (stderr, "Value outside of 1 +/- 1ulp.\n");
1111        exit (EXIT_FAILURE);
1112      }
1113    /* This gives one more ulp.  */
1114    value = FUNC(nextafter) (value, 20);
1115    ulps = ULPDIFF (value, 10);
1116    ulpx = ulp (LIT(2.0));
1117    if (ulps < (LIT(2.0) - ulpx) || ulps > (LIT(2.0) + ulpx))
1118      {
1119        fprintf (stderr, "Value outside of 2 +/- 1ulp.\n");
1120        exit (EXIT_FAILURE);
1121      }
1122    /* And now calculate 100 ulp.  */
1123    for (i = 2; i < 100; i++)
1124      value = FUNC(nextafter) (value, 20);
1125    ulps = ULPDIFF (value, 10);
1126    ulpx = ulp (LIT(100.0));
1127    if (ulps < (LIT(100.0) - ulpx) || ulps > (LIT(100.0) + ulpx))
1128      {
1129        fprintf (stderr, "Value outside of 100 +/- 1ulp.\n");
1130        exit (EXIT_FAILURE);
1131      }
1132 }
1133
1134 /* Do all initialization for a test run with arguments given by ARGC
1135    and ARGV.  */
1136 void
1137 libm_test_init (int argc, char **argv)
1138 {
1139   int remaining;
1140   char *ulps_file_path;
1141   size_t dir_len = 0;
1142
1143   verbose = 1;
1144   output_ulps = 0;
1145   output_max_error = 1;
1146   output_points = 1;
1147   output_dir = NULL;
1148   /* XXX set to 0 for releases.  */
1149   ignore_max_ulp = 0;
1150
1151   /* Parse and process arguments.  */
1152   argp_parse (&argp, argc, argv, 0, &remaining, NULL);
1153
1154   if (remaining != argc)
1155     {
1156       fprintf (stderr, "wrong number of arguments");
1157       argp_help (&argp, stdout, ARGP_HELP_SEE, program_invocation_short_name);
1158       exit (EXIT_FAILURE);
1159     }
1160
1161   if (output_ulps)
1162     {
1163       if (output_dir != NULL)
1164         dir_len = strlen (output_dir);
1165       ulps_file_path = (char *) malloc (dir_len + strlen (ulps_file_name) + 1);
1166       if (ulps_file_path == NULL)
1167         {
1168           perror ("can't allocate path for `ULPs' file: ");
1169           exit (1);
1170         }
1171       sprintf (ulps_file_path, "%s%s", output_dir == NULL ? "" : output_dir, ulps_file_name);
1172       ulps_file = fopen (ulps_file_path, "a");
1173       if (ulps_file == NULL)
1174         {
1175           perror ("can't open file `ULPs' for writing: ");
1176           exit (1);
1177         }
1178     }
1179
1180
1181   initialize ();
1182   fputs (test_msg, stdout);
1183
1184   check_ulp ();
1185 }
1186
1187 /* Process the test results, returning the exit status.  */
1188 int
1189 libm_test_finish (void)
1190 {
1191   if (output_ulps)
1192     fclose (ulps_file);
1193
1194   printf ("\nTest suite completed:\n");
1195   printf ("  %d test cases plus %d tests for exception flags and\n"
1196           "    %d tests for errno executed.\n",
1197           noTests, noExcTests, noErrnoTests);
1198   if (noErrors)
1199     {
1200       printf ("  %d errors occurred.\n", noErrors);
1201       return 1;
1202     }
1203   printf ("  All tests passed successfully.\n");
1204
1205   return 0;
1206 }