math/libm-test-support.c

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