math/gen-auto-libm-tests.c

   1 /* Generate expected output for libm tests with MPFR and MPC.
   2    Copyright (C) 2013 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 /* Compile this program as:
  20
  21    gcc -std=gnu99 -O2 -Wall -Wextra gen-auto-libm-tests.c -lmpc -lmpfr -lgmp \
  22      -o gen-auto-libm-tests
  23
  24    (use of current MPC and MPFR versions recommended) and run it as:
  25
  26    gen-auto-libm-tests auto-libm-test-in auto-libm-test-out
  27
  28    The input file auto-libm-test-in contains three kinds of lines:
  29
  30    Lines beginning with "#" are comments, and are ignored, as are
  31    empty lines.
  32
  33    Other lines are test lines, of the form "function input1 input2
  34    ... [flag1 flag2 ...]".  Inputs are either finite real numbers or
  35    integers, depending on the function under test.  Real numbers may
  36    be in any form acceptable to mpfr_strtofr (base 0); integers in any
  37    form acceptable to mpz_set_str (base 0).  In addition, real numbers
  38    may be certain special strings such as "pi", as listed in the
  39    special_real_inputs array.
  40
  41    Each flag is a flag name possibly followed by a series of
  42    ":condition".  Conditions may be any of the names of floating-point
  43    formats in the floating_point_formats array, "long32" and "long64"
  44    to indicate the number of bits in the "long" type, or other strings
  45    for which libm-test.inc defines a TEST_COND_<condition> macro (with
  46    "-"- changed to "_" in the condition name) evaluating to nonzero
  47    when the condition is true and zero when the condition is false.
  48    The meaning is that the flag applies to the test if all the listed
  49    conditions are true.  "flag:cond1:cond2 flag:cond3:cond4" means the
  50    flag applies if ((cond1 && cond2) || (cond3 && cond4)).
  51
  52    A real number specified as an input is considered to represent the
  53    set of real numbers arising from rounding the given number in any
  54    direction for any supported floating-point format; any roundings
  55    that give infinity are ignored.  Each input on a test line has all
  56    the possible roundings considered independently.  Each resulting
  57    choice of the tuple of inputs to the function is ignored if the
  58    mathematical result of the function involves a NaN or an exact
  59    infinity, and is otherwise considered for each floating-point
  60    format for which all those inputs are exactly representable.  Thus
  61    tests may result in "overflow", "underflow" and "inexact"
  62    exceptions; "invalid" may arise only when the final result type is
  63    an integer type and it is the conversion of a mathematically
  64    defined finite result to integer type that results in that
  65    exception.
  66
  67    By default, it is assumed that "overflow" and "underflow"
  68    exceptions should be correct, but that "inexact" exceptions should
  69    only be correct for functions listed as exactly determined.  For
  70    such functions, "underflow" exceptions should respect whether the
  71    machine has before-rounding or after-rounding tininess detection.
  72    For other functions, it is considered that if the exact result is
  73    somewhere between the greatest magnitude subnormal of a given sign
  74    (exclusive) and the least magnitude normal of that sign
  75    (inclusive), underflow exceptions are permitted but optional on all
  76    machines, and they are also permitted but optional for smaller
  77    subnormal exact results for functions that are not exactly
  78    determined.  errno setting is expected for overflow to infinity and
  79    underflow to zero (for real functions), and for out-of-range
  80    conversion of a finite result to integer type, and is considered
  81    permitted but optional for all other cases where overflow
  82    exceptions occur, and where underflow exceptions occur or are
  83    permitted.  In other cases (where no overflow or underflow is
  84    permitted), errno is expected to be left unchanged.
  85
  86    The flag "no-test-inline" indicates a test is disabled for inline
  87    function testing; "xfail" indicates the test is disabled as
  88    expected to produce incorrect results, "xfail-rounding" indicates
  89    the test is disabled only in rounding modes other than
  90    round-to-nearest.  Otherwise, test flags are of the form
  91    "spurious-<exception>" and "missing-<exception>", for any exception
  92    ("overflow", "underflow", "inexact", "invalid", "divbyzero"),
  93    "spurious-errno" and "missing-errno", to indicate when tests are
  94    expected to deviate from the exception and errno settings
  95    corresponding to the mathematical results.  "xfail",
  96    "xfail-rounding", "spurious-" and "missing-" flags should be
  97    accompanied by a comment referring to an open bug in glibc
  98    Bugzilla.
  99
 100    The output file auto-libm-test-out contains the test lines from
 101    auto-libm-test-in, and, after the line for a given test, some
 102    number of output test lines.  An output test line is of the form "=
 103    function rounding-mode format input1 input2 ... : output1 output2
 104    ... : flags".  rounding-mode is "tonearest", "towardzero", "upward"
 105    or "downward".  format is a name from the floating_point_formats
 106    array, possibly followed by a sequence of ":flag" for flags from
 107    "long32", "long64", "before-rounding" and "after-rounding" (the
 108    last two indicating tests where expectations for underflow
 109    exceptions depend on how the architecture detects tininess).
 110    Inputs and outputs are specified as hex floats with the required
 111    suffix for the floating-point type, or plus_infty or minus_infty
 112    for infinite expected results, or as integer constant expressions
 113    (not necessarily with the right type) or IGNORE for integer inputs
 114    and outputs.  Flags are "no-test-inline", "xfail", "<exception>",
 115    "<exception>-ok", "errno-<value>", "errno-<value>-ok", where
 116    "<exception>" and "errno-<value>" are unconditional, indicating
 117    that a correct result means the given exception should be raised or
 118    errno should be set to the given value, and other settings may be
 119    conditional or unconditional; "-ok" means not to test for the given
 120    exception or errno value (whether because it was marked as possibly
 121    missing or spurious, or because the calculation of correct results
 122    indicated it was optional).  */
 123
 124 #define _GNU_SOURCE
 125
 126 #include <assert.h>
 127 #include <ctype.h>
 128 #include <errno.h>
 129 #include <error.h>
 130 #include <stdbool.h>
 131 #include <stdint.h>
 132 #include <stdio.h>
 133 #include <stdlib.h>
 134 #include <string.h>
 135
 136 #include <gmp.h>
 137 #include <mpfr.h>
 138 #include <mpc.h>
 139
 140 #define ARRAY_SIZE(A) (sizeof (A) / sizeof ((A)[0]))
 141
 142 /* The supported floating-point formats.  */
 143 typedef enum
 144   {
 145     fp_flt_32,
 146     fp_dbl_64,
 147     fp_ldbl_96_intel,
 148     fp_ldbl_96_m68k,
 149     fp_ldbl_128,
 150     fp_ldbl_128ibm,
 151     fp_num_formats,
 152     fp_first_format = 0
 153   } fp_format;
 154
 155 /* Structure describing a single floating-point format.  */
 156 typedef struct
 157 {
 158   /* The name of the format.  */
 159   const char *name;
 160   /* The suffix to use on floating-point constants with this
 161      format.  */
 162   const char *suffix;
 163   /* A string for the largest normal value, or NULL for IEEE formats
 164      where this can be determined automatically.  */
 165   const char *max_string;
 166   /* The number of mantissa bits.  */
 167   int mant_dig;
 168   /* The least N such that 2^N overflows.  */
 169   int max_exp;
 170   /* One more than the least N such that 2^N is normal.  */
 171   int min_exp;
 172   /* The largest normal value.  */
 173   mpfr_t max;
 174   /* The least positive normal value, 2^(MIN_EXP-1).  */
 175   mpfr_t min;
 176   /* The greatest positive subnormal value.  */
 177   mpfr_t subnorm_max;
 178   /* The least positive subnormal value, 2^(MIN_EXP-MANT_DIG).  */
 179   mpfr_t subnorm_min;
 180 } fp_format_desc;
 181
 182 /* List of floating-point formats, in the same order as the fp_format
 183    enumeration.  */
 184 static fp_format_desc fp_formats[fp_num_formats] =
 185   {
 186     { "flt-32", "f", NULL, 24, 128, -125, {}, {}, {}, {} },
 187     { "dbl-64", "", NULL, 53, 1024, -1021, {}, {}, {}, {} },
 188     { "ldbl-96-intel", "L", NULL, 64, 16384, -16381, {}, {}, {}, {} },
 189     { "ldbl-96-m68k", "L", NULL, 64, 16384, -16382, {}, {}, {}, {} },
 190     { "ldbl-128", "L", NULL, 113, 16384, -16381, {}, {}, {}, {} },
 191     { "ldbl-128ibm", "L", "0x1.fffffffffffff7ffffffffffff8p+1023",
 192       106, 1024, -968, {}, {}, {}, {} },
 193   };
 194
 195 /* The supported rounding modes.  */
 196 typedef enum
 197   {
 198     rm_downward,
 199     rm_tonearest,
 200     rm_towardzero,
 201     rm_upward,
 202     rm_num_modes,
 203     rm_first_mode = 0
 204   } rounding_mode;
 205
 206 /* Structure describing a single rounding mode.  */
 207 typedef struct
 208 {
 209   /* The name of the rounding mode.  */
 210   const char *name;
 211   /* The MPFR rounding mode.  */
 212   mpfr_rnd_t mpfr_mode;
 213 } rounding_mode_desc;
 214
 215 /* List of rounding modes, in the same order as the rounding_mode
 216    enumeration.  */
 217 static const rounding_mode_desc rounding_modes[rm_num_modes] =
 218   {
 219     { "downward", MPFR_RNDD },
 220     { "tonearest", MPFR_RNDN },
 221     { "towardzero", MPFR_RNDZ },
 222     { "upward", MPFR_RNDU },
 223   };
 224
 225 /* The supported exceptions.  */
 226 typedef enum
 227   {
 228     exc_divbyzero,
 229     exc_inexact,
 230     exc_invalid,
 231     exc_overflow,
 232     exc_underflow,
 233     exc_num_exceptions,
 234     exc_first_exception = 0
 235   } fp_exception;
 236
 237 /* List of exceptions, in the same order as the fp_exception
 238    enumeration.  */
 239 static const char *const exceptions[exc_num_exceptions] =
 240   {
 241     "divbyzero",
 242     "inexact",
 243     "invalid",
 244     "overflow",
 245     "underflow",
 246   };
 247
 248 /* The internal precision to use for most MPFR calculations, which
 249    must be at least 2 more than the greatest precision of any
 250    supported floating-point format.  */
 251 static int internal_precision;
 252
 253 /* A value that overflows all supported floating-point formats.  */
 254 static mpfr_t global_max;
 255
 256 /* A value that is at most half the least subnormal in any
 257    floating-point format and so is rounded the same way as all
 258    sufficiently small positive values.  */
 259 static mpfr_t global_min;
 260
 261 /* The maximum number of (real or integer) arguments to a function
 262    handled by this program (complex arguments count as two real
 263    arguments).  */
 264 #define MAX_NARGS 4
 265
 266 /* The maximum number of (real or integer) return values from a
 267    function handled by this program.  */
 268 #define MAX_NRET 2
 269
 270 /* A type of a function argument or return value.  */
 271 typedef enum
 272   {
 273     /* No type (not a valid argument or return value).  */
 274     type_none,
 275     /* A floating-point value with the type corresponding to that of
 276        the function.  */
 277     type_fp,
 278     /* An integer value of type int.  */
 279     type_int,
 280     /* An integer value of type long.  */
 281     type_long,
 282     /* An integer value of type long long.  */
 283     type_long_long,
 284   } arg_ret_type;
 285
 286 /* A type of a generic real or integer value.  */
 287 typedef enum
 288   {
 289     /* No type.  */
 290     gtype_none,
 291     /* Floating-point (represented with MPFR).  */
 292     gtype_fp,
 293     /* Integer (represented with GMP).  */
 294     gtype_int,
 295   } generic_value_type;
 296
 297 /* A generic value (argument or result).  */
 298 typedef struct
 299 {
 300   /* The type of this value.  */
 301   generic_value_type type;
 302   /* Its value.  */
 303   union
 304   {
 305     mpfr_t f;
 306     mpz_t i;
 307   } value;
 308 } generic_value;
 309
 310 /* A type of input flag.  */
 311 typedef enum
 312   {
 313     flag_no_test_inline,
 314     flag_xfail,
 315     flag_xfail_rounding,
 316     /* The "spurious" and "missing" flags must be in the same order as
 317        the fp_exception enumeration.  */
 318     flag_spurious_divbyzero,
 319     flag_spurious_inexact,
 320     flag_spurious_invalid,
 321     flag_spurious_overflow,
 322     flag_spurious_underflow,
 323     flag_spurious_errno,
 324     flag_missing_divbyzero,
 325     flag_missing_inexact,
 326     flag_missing_invalid,
 327     flag_missing_overflow,
 328     flag_missing_underflow,
 329     flag_missing_errno,
 330     num_input_flag_types,
 331     flag_first_flag = 0,
 332     flag_spurious_first = flag_spurious_divbyzero,
 333     flag_missing_first = flag_missing_divbyzero
 334   } input_flag_type;
 335
 336 /* List of flags, in the same order as the input_flag_type
 337    enumeration.  */
 338 static const char *const input_flags[num_input_flag_types] =
 339   {
 340     "no-test-inline",
 341     "xfail",
 342     "xfail-rounding",
 343     "spurious-divbyzero",
 344     "spurious-inexact",
 345     "spurious-invalid",
 346     "spurious-overflow",
 347     "spurious-underflow",
 348     "spurious-errno",
 349     "missing-divbyzero",
 350     "missing-inexact",
 351     "missing-invalid",
 352     "missing-overflow",
 353     "missing-underflow",
 354     "missing-errno",
 355   };
 356
 357 /* An input flag, possibly conditional.  */
 358 typedef struct
 359 {
 360   /* The type of this flag.  */
 361   input_flag_type type;
 362   /* The conditions on this flag, as a string ":cond1:cond2..." or
 363      NULL.  */
 364   const char *cond;
 365 } input_flag;
 366
 367 /* Structure describing a single test from the input file (which may
 368    expand into many tests in the output).  The choice of function,
 369    which implies the numbers and types of arguments and results, is
 370    implicit rather than stored in this structure (except as part of
 371    the source line).  */
 372 typedef struct
 373 {
 374   /* The text of the input line describing the test, including the
 375      trailing newline.  */
 376   const char *line;
 377   /* The number of combinations of interpretations of input values for
 378      different floating-point formats and rounding modes.  */
 379   size_t num_input_cases;
 380   /* The corresponding lists of inputs.  */
 381   generic_value **inputs;
 382   /* The number of flags for this test.  */
 383   size_t num_flags;
 384   /* The corresponding list of flags.  */
 385   input_flag *flags;
 386   /* The old output for this test.  */
 387   const char *old_output;
 388 } input_test;
 389
 390 /* Ways to calculate a function.  */
 391 typedef enum
 392   {
 393     /* MPFR function with a single argument and result.  */
 394     mpfr_f_f,
 395   } func_calc_method;
 396
 397 /* Description of how to calculate a function.  */
 398 typedef struct
 399 {
 400   /* Which method is used to calculate the function.  */
 401   func_calc_method method;
 402   /* The specific function called.  */
 403   union
 404   {
 405     int (*mpfr_f_f) (mpfr_t, const mpfr_t, mpfr_rnd_t);
 406   } func;
 407 } func_calc_desc;
 408
 409 /* Structure describing a function handled by this program.  */
 410 typedef struct
 411 {
 412   /* The name of the function.  */
 413   const char *name;
 414   /* The number of arguments.  */
 415   size_t num_args;
 416   /* The types of the arguments.  */
 417   arg_ret_type arg_types[MAX_NARGS];
 418   /* The number of return values.  */
 419   size_t num_ret;
 420   /* The types of the return values.  */
 421   arg_ret_type ret_types[MAX_NRET];
 422   /* Whether the function has exactly determined results and
 423      exceptions.  */
 424   bool exact;
 425   /* Whether the function is a complex function, so errno setting is
 426      optional.  */
 427   bool complex_fn;
 428   /* How to calculate this function.  */
 429   func_calc_desc calc;
 430   /* The number of tests allocated for this function.  */
 431   size_t num_tests_alloc;
 432   /* The number of tests for this function.  */
 433   size_t num_tests;
 434   /* The tests themselves.  */
 435   input_test *tests;
 436 } test_function;
 437
 438 #define FUNC_mpfr_f_f(NAME, MPFR_FUNC, EXACT)           \
 439   {                                                     \
 440     NAME, 1, { type_fp }, 1, { type_fp }, EXACT, false, \
 441     { mpfr_f_f, { .mpfr_f_f = MPFR_FUNC } }, 0, 0, NULL \
 442   }
 443
 444 /* List of functions handled by this program.  */
 445 static test_function test_functions[] =
 446   {
 447     FUNC_mpfr_f_f ("acos", mpfr_acos, false),
 448     FUNC_mpfr_f_f ("acosh", mpfr_acosh, false),
 449     FUNC_mpfr_f_f ("asin", mpfr_asin, false),
 450     FUNC_mpfr_f_f ("asinh", mpfr_asinh, false),
 451     FUNC_mpfr_f_f ("atan", mpfr_atan, false),
 452     FUNC_mpfr_f_f ("atanh", mpfr_atanh, false),
 453     FUNC_mpfr_f_f ("cbrt", mpfr_cbrt, false),
 454     FUNC_mpfr_f_f ("cos", mpfr_cos, false),
 455     FUNC_mpfr_f_f ("cosh", mpfr_cosh, false),
 456     FUNC_mpfr_f_f ("erf", mpfr_erf, false),
 457     FUNC_mpfr_f_f ("erfc", mpfr_erfc, false),
 458     FUNC_mpfr_f_f ("exp", mpfr_exp, false),
 459     FUNC_mpfr_f_f ("exp10", mpfr_exp10, false),
 460     FUNC_mpfr_f_f ("exp2", mpfr_exp2, false),
 461     FUNC_mpfr_f_f ("expm1", mpfr_expm1, false),
 462     FUNC_mpfr_f_f ("j0", mpfr_j0, false),
 463     FUNC_mpfr_f_f ("j1", mpfr_j1, false),
 464     FUNC_mpfr_f_f ("log", mpfr_log, false),
 465     FUNC_mpfr_f_f ("log10", mpfr_log10, false),
 466     FUNC_mpfr_f_f ("log1p", mpfr_log1p, false),
 467     FUNC_mpfr_f_f ("log2", mpfr_log2, false),
 468     FUNC_mpfr_f_f ("sin", mpfr_sin, false),
 469     FUNC_mpfr_f_f ("sinh", mpfr_sinh, false),
 470     FUNC_mpfr_f_f ("sqrt", mpfr_sqrt, true),
 471     FUNC_mpfr_f_f ("tan", mpfr_tan, false),
 472     FUNC_mpfr_f_f ("tanh", mpfr_tanh, false),
 473     FUNC_mpfr_f_f ("tgamma", mpfr_gamma, false),
 474     FUNC_mpfr_f_f ("y0", mpfr_y0, false),
 475     FUNC_mpfr_f_f ("y1", mpfr_y1, false),
 476   };
 477
 478 /* Allocate memory, with error checking.  */
 479
 480 static void *
 481 xmalloc (size_t n)
 482 {
 483   void *p = malloc (n);
 484   if (p == NULL)
 485     error (EXIT_FAILURE, errno, "xmalloc failed");
 486   return p;
 487 }
 488
 489 static void *
 490 xrealloc (void *p, size_t n)
 491 {
 492   p = realloc (p, n);
 493   if (p == NULL)
 494     error (EXIT_FAILURE, errno, "xrealloc failed");
 495   return p;
 496 }
 497
 498 static char *
 499 xstrdup (const char *s)
 500 {
 501   char *p = strdup (s);
 502   if (p == NULL)
 503     error (EXIT_FAILURE, errno, "xstrdup failed");
 504   return p;
 505 }
 506
 507 /* Assert that the result of an MPFR operation was exact; that is,
 508    that the returned ternary value was 0.  */
 509
 510 static void
 511 assert_exact (int i)
 512 {
 513   assert (i == 0);
 514 }
 515
 516 /* Return the generic type of an argument or return value type T.  */
 517
 518 static generic_value_type
 519 generic_arg_ret_type (arg_ret_type t)
 520 {
 521   switch (t)
 522     {
 523     case type_fp:
 524       return gtype_fp;
 525
 526     case type_int:
 527     case type_long:
 528     case type_long_long:
 529       return gtype_int;
 530
 531     default:
 532       abort ();
 533     }
 534 }
 535
 536 /* Free a generic_value *V.  */
 537
 538 static void
 539 generic_value_free (generic_value *v)
 540 {
 541   switch (v->type)
 542     {
 543     case gtype_fp:
 544       mpfr_clear (v->value.f);
 545       break;
 546
 547     case gtype_int:
 548       mpz_clear (v->value.i);
 549       break;
 550
 551     default:
 552       abort ();
 553     }
 554 }
 555
 556 /* Copy a generic_value *SRC to *DEST.  */
 557
 558 static void
 559 generic_value_copy (generic_value *dest, const generic_value *src)
 560 {
 561   dest->type = src->type;
 562   switch (src->type)
 563     {
 564     case gtype_fp:
 565       mpfr_init (dest->value.f);
 566       assert_exact (mpfr_set (dest->value.f, src->value.f, MPFR_RNDN));
 567       break;
 568
 569     case gtype_int:
 570       mpz_init (dest->value.i);
 571       mpz_set (dest->value.i, src->value.i);
 572       break;
 573
 574     default:
 575       abort ();
 576     }
 577 }
 578
 579 /* Initialize data for floating-point formats.  */
 580
 581 static void
 582 init_fp_formats ()
 583 {
 584   int global_max_exp = 0, global_min_subnorm_exp = 0;
 585   for (fp_format f = fp_first_format; f < fp_num_formats; f++)
 586     {
 587       if (fp_formats[f].mant_dig + 2 > internal_precision)
 588         internal_precision = fp_formats[f].mant_dig + 2;
 589       if (fp_formats[f].max_exp > global_max_exp)
 590         global_max_exp = fp_formats[f].max_exp;
 591       int min_subnorm_exp = fp_formats[f].min_exp - fp_formats[f].mant_dig;
 592       if (min_subnorm_exp < global_min_subnorm_exp)
 593         global_min_subnorm_exp = min_subnorm_exp;
 594       mpfr_init2 (fp_formats[f].max, fp_formats[f].mant_dig);
 595       if (fp_formats[f].max_string != NULL)
 596         {
 597           char *ep = NULL;
 598           assert_exact (mpfr_strtofr (fp_formats[f].max,
 599                                       fp_formats[f].max_string,
 600                                       &ep, 0, MPFR_RNDN));
 601           assert (*ep == 0);
 602         }
 603       else
 604         {
 605           assert_exact (mpfr_set_ui_2exp (fp_formats[f].max, 1,
 606                                           fp_formats[f].max_exp,
 607                                           MPFR_RNDN));
 608           mpfr_nextbelow (fp_formats[f].max);
 609         }
 610       mpfr_init2 (fp_formats[f].min, fp_formats[f].mant_dig);
 611       assert_exact (mpfr_set_ui_2exp (fp_formats[f].min, 1,
 612                                       fp_formats[f].min_exp - 1,
 613                                       MPFR_RNDN));
 614       mpfr_init2 (fp_formats[f].subnorm_max, fp_formats[f].mant_dig);
 615       assert_exact (mpfr_set (fp_formats[f].subnorm_max, fp_formats[f].min,
 616                               MPFR_RNDN));
 617       mpfr_nextbelow (fp_formats[f].subnorm_max);
 618       mpfr_nextbelow (fp_formats[f].subnorm_max);
 619       mpfr_init2 (fp_formats[f].subnorm_min, fp_formats[f].mant_dig);
 620       assert_exact (mpfr_set_ui_2exp (fp_formats[f].subnorm_min, 1,
 621                                       min_subnorm_exp, MPFR_RNDN));
 622     }
 623   mpfr_set_default_prec (internal_precision);
 624   mpfr_init (global_max);
 625   assert_exact (mpfr_set_ui_2exp (global_max, 1, global_max_exp, MPFR_RNDN));
 626   mpfr_init (global_min);
 627   assert_exact (mpfr_set_ui_2exp (global_min, 1, global_min_subnorm_exp - 1,
 628                                   MPFR_RNDN));
 629 }
 630
 631 /* Fill in mpfr_t values for special strings in input arguments.  */
 632
 633 static size_t
 634 special_fill_max (mpfr_t res0, mpfr_t res1 __attribute__ ((unused)),
 635                   fp_format format)
 636 {
 637   mpfr_init2 (res0, fp_formats[format].mant_dig);
 638   assert_exact (mpfr_set (res0, fp_formats[format].max, MPFR_RNDN));
 639   return 1;
 640 }
 641
 642 static size_t
 643 special_fill_minus_max (mpfr_t res0, mpfr_t res1 __attribute__ ((unused)),
 644                         fp_format format)
 645 {
 646   mpfr_init2 (res0, fp_formats[format].mant_dig);
 647   assert_exact (mpfr_neg (res0, fp_formats[format].max, MPFR_RNDN));
 648   return 1;
 649 }
 650
 651 static size_t
 652 special_fill_pi (mpfr_t res0, mpfr_t res1, fp_format format)
 653 {
 654   mpfr_init2 (res0, fp_formats[format].mant_dig);
 655   mpfr_const_pi (res0, MPFR_RNDU);
 656   mpfr_init2 (res1, fp_formats[format].mant_dig);
 657   mpfr_const_pi (res1, MPFR_RNDD);
 658   return 2;
 659 }
 660
 661 static size_t
 662 special_fill_minus_pi (mpfr_t res0, mpfr_t res1, fp_format format)
 663 {
 664   mpfr_init2 (res0, fp_formats[format].mant_dig);
 665   mpfr_const_pi (res0, MPFR_RNDU);
 666   assert_exact (mpfr_neg (res0, res0, MPFR_RNDN));
 667   mpfr_init2 (res1, fp_formats[format].mant_dig);
 668   mpfr_const_pi (res1, MPFR_RNDD);
 669   assert_exact (mpfr_neg (res1, res1, MPFR_RNDN));
 670   return 2;
 671 }
 672
 673 static size_t
 674 special_fill_pi_2 (mpfr_t res0, mpfr_t res1, fp_format format)
 675 {
 676   mpfr_init2 (res0, fp_formats[format].mant_dig);
 677   mpfr_const_pi (res0, MPFR_RNDU);
 678   assert_exact (mpfr_div_ui (res0, res0, 2, MPFR_RNDN));
 679   mpfr_init2 (res1, fp_formats[format].mant_dig);
 680   mpfr_const_pi (res1, MPFR_RNDD);
 681   assert_exact (mpfr_div_ui (res1, res1, 2, MPFR_RNDN));
 682   return 2;
 683 }
 684
 685 static size_t
 686 special_fill_minus_pi_2 (mpfr_t res0, mpfr_t res1, fp_format format)
 687 {
 688   mpfr_init2 (res0, fp_formats[format].mant_dig);
 689   mpfr_const_pi (res0, MPFR_RNDU);
 690   assert_exact (mpfr_div_ui (res0, res0, 2, MPFR_RNDN));
 691   assert_exact (mpfr_neg (res0, res0, MPFR_RNDN));
 692   mpfr_init2 (res1, fp_formats[format].mant_dig);
 693   mpfr_const_pi (res1, MPFR_RNDD);
 694   assert_exact (mpfr_div_ui (res1, res1, 2, MPFR_RNDN));
 695   assert_exact (mpfr_neg (res1, res1, MPFR_RNDN));
 696   return 2;
 697 }
 698
 699 static size_t
 700 special_fill_pi_4 (mpfr_t res0, mpfr_t res1, fp_format format)
 701 {
 702   mpfr_init2 (res0, fp_formats[format].mant_dig);
 703   assert_exact (mpfr_set_si (res0, 1, MPFR_RNDN));
 704   mpfr_atan (res0, res0, MPFR_RNDU);
 705   mpfr_init2 (res1, fp_formats[format].mant_dig);
 706   assert_exact (mpfr_set_si (res1, 1, MPFR_RNDN));
 707   mpfr_atan (res1, res1, MPFR_RNDD);
 708   return 2;
 709 }
 710
 711 static size_t
 712 special_fill_pi_6 (mpfr_t res0, mpfr_t res1, fp_format format)
 713 {
 714   mpfr_init2 (res0, fp_formats[format].mant_dig);
 715   assert_exact (mpfr_set_si_2exp (res0, 1, -1, MPFR_RNDN));
 716   mpfr_asin (res0, res0, MPFR_RNDU);
 717   mpfr_init2 (res1, fp_formats[format].mant_dig);
 718   assert_exact (mpfr_set_si_2exp (res1, 1, -1, MPFR_RNDN));
 719   mpfr_asin (res1, res1, MPFR_RNDD);
 720   return 2;
 721 }
 722
 723 static size_t
 724 special_fill_minus_pi_6 (mpfr_t res0, mpfr_t res1, fp_format format)
 725 {
 726   mpfr_init2 (res0, fp_formats[format].mant_dig);
 727   assert_exact (mpfr_set_si_2exp (res0, -1, -1, MPFR_RNDN));
 728   mpfr_asin (res0, res0, MPFR_RNDU);
 729   mpfr_init2 (res1, fp_formats[format].mant_dig);
 730   assert_exact (mpfr_set_si_2exp (res1, -1, -1, MPFR_RNDN));
 731   mpfr_asin (res1, res1, MPFR_RNDD);
 732   return 2;
 733 }
 734
 735 static size_t
 736 special_fill_pi_3 (mpfr_t res0, mpfr_t res1, fp_format format)
 737 {
 738   mpfr_init2 (res0, fp_formats[format].mant_dig);
 739   assert_exact (mpfr_set_si_2exp (res0, 1, -1, MPFR_RNDN));
 740   mpfr_acos (res0, res0, MPFR_RNDU);
 741   mpfr_init2 (res1, fp_formats[format].mant_dig);
 742   assert_exact (mpfr_set_si_2exp (res1, 1, -1, MPFR_RNDN));
 743   mpfr_acos (res1, res1, MPFR_RNDD);
 744   return 2;
 745 }
 746
 747 static size_t
 748 special_fill_2pi_3 (mpfr_t res0, mpfr_t res1, fp_format format)
 749 {
 750   mpfr_init2 (res0, fp_formats[format].mant_dig);
 751   assert_exact (mpfr_set_si_2exp (res0, -1, -1, MPFR_RNDN));
 752   mpfr_acos (res0, res0, MPFR_RNDU);
 753   mpfr_init2 (res1, fp_formats[format].mant_dig);
 754   assert_exact (mpfr_set_si_2exp (res1, -1, -1, MPFR_RNDN));
 755   mpfr_acos (res1, res1, MPFR_RNDD);
 756   return 2;
 757 }
 758
 759 static size_t
 760 special_fill_e (mpfr_t res0, mpfr_t res1, fp_format format)
 761 {
 762   mpfr_init2 (res0, fp_formats[format].mant_dig);
 763   assert_exact (mpfr_set_si (res0, 1, MPFR_RNDN));
 764   mpfr_exp (res0, res0, MPFR_RNDU);
 765   mpfr_init2 (res1, fp_formats[format].mant_dig);
 766   assert_exact (mpfr_set_si (res1, 1, MPFR_RNDN));
 767   mpfr_exp (res1, res1, MPFR_RNDD);
 768   return 2;
 769 }
 770
 771 static size_t
 772 special_fill_1_e (mpfr_t res0, mpfr_t res1, fp_format format)
 773 {
 774   mpfr_init2 (res0, fp_formats[format].mant_dig);
 775   assert_exact (mpfr_set_si (res0, -1, MPFR_RNDN));
 776   mpfr_exp (res0, res0, MPFR_RNDU);
 777   mpfr_init2 (res1, fp_formats[format].mant_dig);
 778   assert_exact (mpfr_set_si (res1, -1, MPFR_RNDN));
 779   mpfr_exp (res1, res1, MPFR_RNDD);
 780   return 2;
 781 }
 782
 783 static size_t
 784 special_fill_e_minus_1 (mpfr_t res0, mpfr_t res1, fp_format format)
 785 {
 786   mpfr_init2 (res0, fp_formats[format].mant_dig);
 787   assert_exact (mpfr_set_si (res0, 1, MPFR_RNDN));
 788   mpfr_expm1 (res0, res0, MPFR_RNDU);
 789   mpfr_init2 (res1, fp_formats[format].mant_dig);
 790   assert_exact (mpfr_set_si (res1, 1, MPFR_RNDN));
 791   mpfr_expm1 (res1, res1, MPFR_RNDD);
 792   return 2;
 793 }
 794
 795 /* A special string accepted in input arguments.  */
 796 typedef struct
 797 {
 798   /* The string.  */
 799   const char *str;
 800   /* The function that interprets it for a given floating-point
 801      format, filling in up to two mpfr_t values and returning the
 802      number of values filled.  */
 803   size_t (*func) (mpfr_t, mpfr_t, fp_format);
 804 } special_real_input;
 805
 806 /* List of special strings accepted in input arguments.  */
 807
 808 static const special_real_input special_real_inputs[] =
 809   {
 810     { "max", special_fill_max },
 811     { "-max", special_fill_minus_max },
 812     { "pi", special_fill_pi },
 813     { "-pi", special_fill_minus_pi },
 814     { "pi/2", special_fill_pi_2 },
 815     { "-pi/2", special_fill_minus_pi_2 },
 816     { "pi/4", special_fill_pi_4 },
 817     { "pi/6", special_fill_pi_6 },
 818     { "-pi/6", special_fill_minus_pi_6 },
 819     { "pi/3", special_fill_pi_3 },
 820     { "2pi/3", special_fill_2pi_3 },
 821     { "e", special_fill_e },
 822     { "1/e", special_fill_1_e },
 823     { "e-1", special_fill_e_minus_1 },
 824   };
 825
 826 /* Given a real number R computed in round-to-zero mode, set the
 827    lowest bit as a sticky bit if INEXACT, and saturate the exponent
 828    range for very large or small values.  */
 829
 830 static void
 831 adjust_real (mpfr_t r, bool inexact)
 832 {
 833   if (!inexact)
 834     return;
 835   /* NaNs are exact, as are infinities in round-to-zero mode.  */
 836   assert (mpfr_number_p (r));
 837   if (mpfr_cmpabs (r, global_min) < 0)
 838     assert_exact (mpfr_copysign (r, global_min, r, MPFR_RNDN));
 839   else if (mpfr_cmpabs (r, global_max) > 0)
 840     assert_exact (mpfr_copysign (r, global_max, r, MPFR_RNDN));
 841   else
 842     {
 843       mpz_t tmp;
 844       mpz_init (tmp);
 845       mpfr_exp_t e = mpfr_get_z_2exp (tmp, r);
 846       mpz_setbit (tmp, 0);
 847       assert_exact (mpfr_set_z_2exp (r, tmp, e, MPFR_RNDN));
 848       mpz_clear (tmp);
 849     }
 850 }
 851
 852 /* Given a finite real number R with sticky bit, compute the roundings
 853    to FORMAT in each rounding mode, storing the results in RES, the
 854    before-rounding exceptions in EXC_BEFORE and the after-rounding
 855    exceptions in EXC_AFTER.  */
 856
 857 static void
 858 round_real (mpfr_t res[rm_num_modes],
 859             unsigned int exc_before[rm_num_modes],
 860             unsigned int exc_after[rm_num_modes],
 861             mpfr_t r, fp_format format)
 862 {
 863   assert (mpfr_number_p (r));
 864   for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
 865     {
 866       mpfr_init2 (res[m], fp_formats[format].mant_dig);
 867       exc_before[m] = exc_after[m] = 0;
 868       bool inexact = mpfr_set (res[m], r, rounding_modes[m].mpfr_mode);
 869       if (mpfr_cmpabs (res[m], fp_formats[format].max) > 0)
 870         {
 871           inexact = true;
 872           exc_before[m] |= 1U << exc_overflow;
 873           exc_after[m] |= 1U << exc_overflow;
 874           bool overflow_inf;
 875           switch (m)
 876             {
 877             case rm_tonearest:
 878               overflow_inf = true;
 879               break;
 880             case rm_towardzero:
 881               overflow_inf = false;
 882               break;
 883             case rm_downward:
 884               overflow_inf = mpfr_signbit (res[m]);
 885               break;
 886             case rm_upward:
 887               overflow_inf = !mpfr_signbit (res[m]);
 888               break;
 889             default:
 890               abort ();
 891             }
 892           if (overflow_inf)
 893             mpfr_set_inf (res[m], mpfr_signbit (res[m]) ? -1 : 1);
 894           else
 895             assert_exact (mpfr_copysign (res[m], fp_formats[format].max,
 896                                          res[m], MPFR_RNDN));
 897         }
 898       if (mpfr_cmpabs (r, fp_formats[format].min) < 0)
 899         {
 900           /* Tiny before rounding; may or may not be tiny after
 901              rounding, and underflow applies only if also inexact
 902              around rounding to a possibly subnormal value.  */
 903           bool tiny_after_rounding
 904             = mpfr_cmpabs (res[m], fp_formats[format].min) < 0;
 905           /* To round to a possibly subnormal value, and determine
 906              inexactness as a subnormal in the process, scale up and
 907              round to integer, then scale back down.  */
 908           mpfr_t tmp;
 909           mpfr_init (tmp);
 910           assert_exact (mpfr_mul_2si (tmp, r, (fp_formats[format].mant_dig
 911                                                - fp_formats[format].min_exp),
 912                                       MPFR_RNDN));
 913           int rint_res = mpfr_rint (tmp, tmp, rounding_modes[m].mpfr_mode);
 914           /* The integer must be representable.  */
 915           assert (rint_res == 0 || rint_res == 2 || rint_res == -2);
 916           /* If rounding to full precision was inexact, so must
 917              rounding to subnormal precision be inexact.  */
 918           if (inexact)
 919             assert (rint_res != 0);
 920           else
 921             inexact = rint_res != 0;
 922           assert_exact (mpfr_mul_2si (res[m], tmp,
 923                                       (fp_formats[format].min_exp
 924                                        - fp_formats[format].mant_dig),
 925                                       MPFR_RNDN));
 926           mpfr_clear (tmp);
 927           if (inexact)
 928             {
 929               exc_before[m] |= 1U << exc_underflow;
 930               if (tiny_after_rounding)
 931                 exc_after[m] |= 1U << exc_underflow;
 932             }
 933         }
 934       if (inexact)
 935         {
 936           exc_before[m] |= 1U << exc_inexact;
 937           exc_after[m] |= 1U << exc_inexact;
 938         }
 939     }
 940 }
 941
 942 /* Handle the input argument at ARG (NUL-terminated), updating the
 943    lists of test inputs in IT accordingly.  NUM_PREV_ARGS arguments
 944    are already in those lists.  The argument, of type GTYPE, comes
 945    from file FILENAME, line LINENO.  */
 946
 947 static void
 948 handle_input_arg (const char *arg, input_test *it, size_t num_prev_args,
 949                   generic_value_type gtype,
 950                   const char *filename, unsigned int lineno)
 951 {
 952   size_t num_values = 0;
 953   generic_value values[2 * fp_num_formats];
 954   switch (gtype)
 955     {
 956     case gtype_fp:
 957       for (fp_format f = fp_first_format; f < fp_num_formats; f++)
 958         {
 959           mpfr_t extra_values[2];
 960           size_t num_extra_values = 0;
 961           for (size_t i = 0; i < ARRAY_SIZE (special_real_inputs); i++)
 962             {
 963               if (strcmp (arg, special_real_inputs[i].str) == 0)
 964                 {
 965                   num_extra_values
 966                     = special_real_inputs[i].func (extra_values[0],
 967                                                    extra_values[1], f);
 968                   assert (num_extra_values > 0
 969                           && num_extra_values <= ARRAY_SIZE (extra_values));
 970                   break;
 971                 }
 972             }
 973           if (num_extra_values == 0)
 974             {
 975               mpfr_t tmp;
 976               char *ep;
 977               mpfr_init (tmp);
 978               bool inexact = mpfr_strtofr (tmp, arg, &ep, 0, MPFR_RNDZ);
 979               if (*ep != 0 || !mpfr_number_p (tmp))
 980                 error_at_line (EXIT_FAILURE, 0, filename, lineno,
 981                                "bad floating-point argument: '%s'", arg);
 982               adjust_real (tmp, inexact);
 983               mpfr_t rounded[rm_num_modes];
 984               unsigned int exc_before[rm_num_modes];
 985               unsigned int exc_after[rm_num_modes];
 986               round_real (rounded, exc_before, exc_after, tmp, f);
 987               mpfr_clear (tmp);
 988               if (mpfr_number_p (rounded[rm_upward]))
 989                 {
 990                   mpfr_init2 (extra_values[num_extra_values],
 991                               fp_formats[f].mant_dig);
 992                   assert_exact (mpfr_set (extra_values[num_extra_values],
 993                                           rounded[rm_upward], MPFR_RNDN));
 994                   num_extra_values++;
 995                 }
 996               if (mpfr_number_p (rounded[rm_downward]))
 997                 {
 998                   mpfr_init2 (extra_values[num_extra_values],
 999                               fp_formats[f].mant_dig);
1000                   assert_exact (mpfr_set (extra_values[num_extra_values],
1001                                           rounded[rm_downward], MPFR_RNDN));
1002                   num_extra_values++;
1003                 }
1004               for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
1005                 mpfr_clear (rounded[m]);
1006             }
1007           for (size_t i = 0; i < num_extra_values; i++)
1008             {
1009               bool found = false;
1010               for (size_t j = 0; j < num_values; j++)
1011                 {
1012                   if (mpfr_equal_p (values[j].value.f, extra_values[i])
1013                       && ((mpfr_signbit (values[j].value.f) != 0)
1014                           == (mpfr_signbit (extra_values[i]) != 0)))
1015                     {
1016                       found = true;
1017                       break;
1018                     }
1019                 }
1020               if (!found)
1021                 {
1022                   assert (num_values < ARRAY_SIZE (values));
1023                   values[num_values].type = gtype_fp;
1024                   mpfr_init2 (values[num_values].value.f,
1025                               fp_formats[f].mant_dig);
1026                   assert_exact (mpfr_set (values[num_values].value.f,
1027                                           extra_values[i], MPFR_RNDN));
1028                   num_values++;
1029                 }
1030               mpfr_clear (extra_values[i]);
1031             }
1032         }
1033       break;
1034
1035     case gtype_int:
1036       num_values = 1;
1037       values[0].type = gtype_int;
1038       int ret = mpz_init_set_str (values[0].value.i, arg, 0);
1039       if (ret != 0)
1040         error_at_line (EXIT_FAILURE, 0, filename, lineno,
1041                        "bad integer argument: '%s'", arg);
1042       break;
1043
1044     default:
1045       abort ();
1046     }
1047   assert (num_values > 0 && num_values <= ARRAY_SIZE (values));
1048   if (it->num_input_cases >= SIZE_MAX / num_values)
1049     error_at_line (EXIT_FAILURE, 0, filename, lineno, "too many input cases");
1050   generic_value **old_inputs = it->inputs;
1051   size_t new_num_input_cases = it->num_input_cases * num_values;
1052   generic_value **new_inputs = xmalloc (new_num_input_cases
1053                                         * sizeof (new_inputs[0]));
1054   for (size_t i = 0; i < it->num_input_cases; i++)
1055     {
1056       for (size_t j = 0; j < num_values; j++)
1057         {
1058           size_t idx = i * num_values + j;
1059           new_inputs[idx] = xmalloc ((num_prev_args + 1)
1060                                      * sizeof (new_inputs[idx][0]));
1061           for (size_t k = 0; k < num_prev_args; k++)
1062             generic_value_copy (&new_inputs[idx][k], &old_inputs[i][k]);
1063           generic_value_copy (&new_inputs[idx][num_prev_args], &values[j]);
1064         }
1065       for (size_t j = 0; j < num_prev_args; j++)
1066         generic_value_free (&old_inputs[i][j]);
1067       free (old_inputs[i]);
1068     }
1069   free (old_inputs);
1070   for (size_t i = 0; i < num_values; i++)
1071     generic_value_free (&values[i]);
1072   it->inputs = new_inputs;
1073   it->num_input_cases = new_num_input_cases;
1074 }
1075
1076 /* Handle the input flag ARG (NUL-terminated), storing it in *FLAG.
1077    The flag comes from file FILENAME, line LINENO.  */
1078
1079 static void
1080 handle_input_flag (char *arg, input_flag *flag,
1081                    const char *filename, unsigned int lineno)
1082 {
1083   char *ep = strchr (arg, ':');
1084   if (ep == NULL)
1085     {
1086       ep = strchr (arg, 0);
1087       assert (ep != NULL);
1088     }
1089   char c = *ep;
1090   *ep = 0;
1091   bool found = false;
1092   for (input_flag_type i = flag_first_flag; i <= num_input_flag_types; i++)
1093     {
1094       if (strcmp (arg, input_flags[i]) == 0)
1095         {
1096           found = true;
1097           flag->type = i;
1098           break;
1099         }
1100     }
1101   if (!found)
1102     error_at_line (EXIT_FAILURE, 0, filename, lineno, "unknown flag: '%s'",
1103                    arg);
1104   *ep = c;
1105   if (c == 0)
1106     flag->cond = NULL;
1107   else
1108     flag->cond = xstrdup (ep);
1109 }
1110
1111 /* Add the test LINE (file FILENAME, line LINENO) to the test
1112    data.  */
1113
1114 static void
1115 add_test (char *line, const char *filename, unsigned int lineno)
1116 {
1117   size_t num_tokens = 1;
1118   char *p = line;
1119   while ((p = strchr (p, ' ')) != NULL)
1120     {
1121       num_tokens++;
1122       p++;
1123     }
1124   if (num_tokens < 2)
1125     error_at_line (EXIT_FAILURE, 0, filename, lineno,
1126                    "line too short: '%s'", line);
1127   p = strchr (line, ' ');
1128   size_t func_name_len = p - line;
1129   for (size_t i = 0; i < ARRAY_SIZE (test_functions); i++)
1130     {
1131       if (func_name_len == strlen (test_functions[i].name)
1132           && strncmp (line, test_functions[i].name, func_name_len) == 0)
1133         {
1134           test_function *tf = &test_functions[i];
1135           if (num_tokens < 1 + tf->num_args)
1136             error_at_line (EXIT_FAILURE, 0, filename, lineno,
1137                            "line too short: '%s'", line);
1138           if (tf->num_tests == tf->num_tests_alloc)
1139             {
1140               tf->num_tests_alloc = 2 * tf->num_tests_alloc + 16;
1141               tf->tests
1142                 = xrealloc (tf->tests,
1143                             tf->num_tests_alloc * sizeof (tf->tests[0]));
1144             }
1145           input_test *it = &tf->tests[tf->num_tests];
1146           it->line = line;
1147           it->num_input_cases = 1;
1148           it->inputs = xmalloc (sizeof (it->inputs[0]));
1149           it->inputs[0] = NULL;
1150           it->old_output = NULL;
1151           p++;
1152           for (size_t j = 0; j < tf->num_args; j++)
1153             {
1154               char *ep = strchr (p, ' ');
1155               if (ep == NULL)
1156                 {
1157                   ep = strchr (p, '\n');
1158                   assert (ep != NULL);
1159                 }
1160               if (ep == p)
1161                 error_at_line (EXIT_FAILURE, 0, filename, lineno,
1162                                "empty token in line: '%s'", line);
1163               for (char *t = p; t < ep; t++)
1164                 if (isspace ((unsigned char) *t))
1165                   error_at_line (EXIT_FAILURE, 0, filename, lineno,
1166                                  "whitespace in token in line: '%s'", line);
1167               char c = *ep;
1168               *ep = 0;
1169               handle_input_arg (p, it, j,
1170                                 generic_arg_ret_type (tf->arg_types[j]),
1171                                 filename, lineno);
1172               *ep = c;
1173               p = ep + 1;
1174             }
1175           it->num_flags = num_tokens - 1 - tf->num_args;
1176           it->flags = xmalloc (it->num_flags * sizeof (it->flags[0]));
1177           for (size_t j = 0; j < it->num_flags; j++)
1178             {
1179               char *ep = strchr (p, ' ');
1180               if (ep == NULL)
1181                 {
1182                   ep = strchr (p, '\n');
1183                   assert (ep != NULL);
1184                 }
1185               if (ep == p)
1186                 error_at_line (EXIT_FAILURE, 0, filename, lineno,
1187                                "empty token in line: '%s'", line);
1188               for (char *t = p; t < ep; t++)
1189                 if (isspace ((unsigned char) *t))
1190                   error_at_line (EXIT_FAILURE, 0, filename, lineno,
1191                                  "whitespace in token in line: '%s'", line);
1192               char c = *ep;
1193               *ep = 0;
1194               handle_input_flag (p, &it->flags[j], filename, lineno);
1195               *ep = c;
1196               p = ep + 1;
1197             }
1198           assert (*p == 0);
1199           tf->num_tests++;
1200           return;
1201         }
1202     }
1203   error_at_line (EXIT_FAILURE, 0, filename, lineno,
1204                  "unknown function in line: '%s'", line);
1205 }
1206
1207 /* Read in the test input data from FILENAME.  */
1208
1209 static void
1210 read_input (const char *filename)
1211 {
1212   FILE *fp = fopen (filename, "r");
1213   if (fp == NULL)
1214     error (EXIT_FAILURE, errno, "open '%s'", filename);
1215   unsigned int lineno = 0;
1216   for (;;)
1217     {
1218       size_t size = 0;
1219       char *line = NULL;
1220       ssize_t ret = getline (&line, &size, fp);
1221       if (ret == -1)
1222         break;
1223       lineno++;
1224       if (line[0] == '#' || line[0] == '\n')
1225         continue;
1226       add_test (line, filename, lineno);
1227     }
1228   if (ferror (fp))
1229     error (EXIT_FAILURE, errno, "read from '%s'", filename);
1230   if (fclose (fp) != 0)
1231     error (EXIT_FAILURE, errno, "close '%s'", filename);
1232 }
1233
1234 /* Calculate the generic results (round-to-zero with sticky bit) for
1235    the function described by CALC, with inputs INPUTS.  */
1236
1237 static void
1238 calc_generic_results (generic_value *outputs, generic_value *inputs,
1239                       const func_calc_desc *calc)
1240 {
1241   bool inexact;
1242   switch (calc->method)
1243     {
1244     case mpfr_f_f:
1245       assert (inputs[0].type == gtype_fp);
1246       outputs[0].type = gtype_fp;
1247       mpfr_init (outputs[0].value.f);
1248       inexact = calc->func.mpfr_f_f (outputs[0].value.f, inputs[0].value.f,
1249                                      MPFR_RNDZ);
1250       adjust_real (outputs[0].value.f, inexact);
1251       break;
1252
1253     default:
1254       abort ();
1255     }
1256 }
1257
1258 /* Return the number of bits for integer type TYPE, where "long" has
1259    LONG_BITS bits (32 or 64).  */
1260
1261 static int
1262 int_type_bits (arg_ret_type type, int long_bits)
1263 {
1264   assert (long_bits == 32 || long_bits == 64);
1265   switch (type)
1266     {
1267     case type_int:
1268       return 32;
1269       break;
1270
1271     case type_long:
1272       return long_bits;
1273       break;
1274
1275     case type_long_long:
1276       return 64;
1277       break;
1278
1279     default:
1280       abort ();
1281     }
1282 }
1283
1284 /* Check whether an integer Z fits a given type TYPE, where "long" has
1285    LONG_BITS bits (32 or 64).  */
1286
1287 static bool
1288 int_fits_type (mpz_t z, arg_ret_type type, int long_bits)
1289 {
1290   int bits = int_type_bits (type, long_bits);
1291   bool ret = true;
1292   mpz_t t;
1293   mpz_init (t);
1294   mpz_ui_pow_ui (t, 2, bits - 1);
1295   if (mpz_cmp (z, t) >= 0)
1296     ret = false;
1297   mpz_neg (t, t);
1298   if (mpz_cmp (z, t) < 0)
1299     ret = false;
1300   mpz_clear (t);
1301   return ret;
1302 }
1303
1304 /* Print a generic value V to FP (name FILENAME), preceded by a space,
1305    for type TYPE, floating-point format FORMAT, LONG_BITS bits per
1306    long, printing " IGNORE" instead if IGNORE.  */
1307
1308 static void
1309 output_generic_value (FILE *fp, const char *filename, const generic_value *v,
1310                       bool ignore, arg_ret_type type, fp_format format,
1311                       int long_bits)
1312 {
1313   if (ignore)
1314     {
1315       if (fputs (" IGNORE", fp) < 0)
1316         error (EXIT_FAILURE, errno, "write to '%s'", filename);
1317       return;
1318     }
1319   assert (v->type == generic_arg_ret_type (type));
1320   const char *suffix;
1321   switch (type)
1322     {
1323     case type_fp:
1324       suffix = fp_formats[format].suffix;
1325       break;
1326
1327     case type_int:
1328       suffix = "";
1329       break;
1330
1331     case type_long:
1332       suffix = "L";
1333       break;
1334
1335     case type_long_long:
1336       suffix = "LL";
1337       break;
1338
1339     default:
1340       abort ();
1341     }
1342   switch (v->type)
1343     {
1344     case gtype_fp:
1345       if (mpfr_inf_p (v->value.f))
1346         {
1347           if (fputs ((mpfr_signbit (v->value.f)
1348                       ? " minus_infty" : " plus_infty"), fp) < 0)
1349             error (EXIT_FAILURE, errno, "write to '%s'", filename);
1350         }
1351       else
1352         {
1353           assert (mpfr_number_p (v->value.f));
1354           if (mpfr_fprintf (fp, " %Ra%s", v->value.f, suffix) < 0)
1355             error (EXIT_FAILURE, errno, "mpfr_fprintf to '%s'", filename);
1356         }
1357       break;
1358
1359     case gtype_int: ;
1360       int bits = int_type_bits (type, long_bits);
1361       mpz_t tmp;
1362       mpz_init (tmp);
1363       mpz_ui_pow_ui (tmp, 2, bits - 1);
1364       mpz_neg (tmp, tmp);
1365       if (mpz_cmp (v->value.i, tmp) == 0)
1366         {
1367           mpz_add_ui (tmp, tmp, 1);
1368           if (mpfr_fprintf (fp, " (%Zd%s-1)", tmp, suffix) < 0)
1369             error (EXIT_FAILURE, errno, "mpfr_fprintf to '%s'", filename);
1370         }
1371       else
1372         {
1373           if (mpfr_fprintf (fp, " %Zd%s", v->value.i, suffix) < 0)
1374             error (EXIT_FAILURE, errno, "mpfr_fprintf to '%s'", filename);
1375         }
1376       mpz_clear (tmp);
1377       break;
1378
1379     default:
1380       abort ();
1381     }
1382 }
1383
1384 /* Generate test output to FP (name FILENAME) for test function TF,
1385    input test IT, choice of input values INPUTS.  */
1386
1387 static void
1388 output_for_one_input_case (FILE *fp, const char *filename, test_function *tf,
1389                            input_test *it, generic_value *inputs)
1390 {
1391   bool long_bits_matters = false;
1392   bool fits_long32 = true;
1393   for (size_t i = 0; i < tf->num_args; i++)
1394     {
1395       generic_value_type gtype = generic_arg_ret_type (tf->arg_types[i]);
1396       assert (inputs[i].type == gtype);
1397       if (gtype == gtype_int)
1398         {
1399           bool fits_64 = int_fits_type (inputs[i].value.i, tf->arg_types[i],
1400                                         64);
1401           if (!fits_64)
1402             return;
1403           if (tf->arg_types[i] == type_long
1404               && !int_fits_type (inputs[i].value.i, tf->arg_types[i], 32))
1405             {
1406               long_bits_matters = true;
1407               fits_long32 = false;
1408             }
1409         }
1410     }
1411   generic_value generic_outputs[MAX_NRET];
1412   calc_generic_results (generic_outputs, inputs, &tf->calc);
1413   bool ignore_output_long32[MAX_NRET] = { false };
1414   bool ignore_output_long64[MAX_NRET] = { false };
1415   for (size_t i = 0; i < tf->num_ret; i++)
1416     {
1417       assert (generic_outputs[i].type
1418               == generic_arg_ret_type (tf->ret_types[i]));
1419       switch (generic_outputs[i].type)
1420         {
1421         case gtype_fp:
1422           if (!mpfr_number_p (generic_outputs[i].value.f))
1423             goto out; /* Result is NaN or exact infinity.  */
1424           break;
1425
1426         case gtype_int:
1427           ignore_output_long32[i] = !int_fits_type (generic_outputs[i].value.i,
1428                                                     tf->ret_types[i], 32);
1429           ignore_output_long64[i] = !int_fits_type (generic_outputs[i].value.i,
1430                                                     tf->ret_types[i], 64);
1431           if (ignore_output_long32[i] != ignore_output_long64[i])
1432             long_bits_matters = true;
1433           break;
1434
1435         default:
1436           abort ();
1437         }
1438     }
1439   /* Iterate over relevant sizes of long and floating-point formats.  */
1440   for (int long_bits = 32; long_bits <= 64; long_bits += 32)
1441     {
1442       if (long_bits == 32 && !fits_long32)
1443         continue;
1444       if (long_bits == 64 && !long_bits_matters)
1445         continue;
1446       const char *long_cond;
1447       if (long_bits_matters)
1448         long_cond = (long_bits == 32 ? ":long32" : ":long64");
1449       else
1450         long_cond = "";
1451       bool *ignore_output = (long_bits == 32
1452                              ? ignore_output_long32
1453                              : ignore_output_long64);
1454       for (fp_format f = fp_first_format; f < fp_num_formats; f++)
1455         {
1456           bool fits = true;
1457           mpfr_t res[rm_num_modes];
1458           unsigned int exc_before[rm_num_modes];
1459           unsigned int exc_after[rm_num_modes];
1460           for (size_t i = 0; i < tf->num_args; i++)
1461             {
1462               if (inputs[i].type == gtype_fp)
1463                 round_real (res, exc_before, exc_after, inputs[i].value.f, f);
1464               if (!mpfr_equal_p (res[rm_tonearest], inputs[i].value.f))
1465                 fits = false;
1466               for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
1467                 mpfr_clear (res[m]);
1468               if (!fits)
1469                 break;
1470             }
1471           if (!fits)
1472             continue;
1473           /* The inputs fit this type, so compute the ideal outputs
1474              and exceptions.  */
1475           mpfr_t all_res[MAX_NRET][rm_num_modes];
1476           unsigned int all_exc_before[MAX_NRET][rm_num_modes];
1477           unsigned int all_exc_after[MAX_NRET][rm_num_modes];
1478           unsigned int merged_exc_before[rm_num_modes] = { 0 };
1479           unsigned int merged_exc_after[rm_num_modes] = { 0 };
1480           /* For functions not exactly determined, track whether
1481              underflow is required (some result is inexact, and
1482              magnitude does not exceed the greatest magnitude
1483              subnormal), and permitted (not an exact zero, and
1484              magnitude does not exceed the least magnitude
1485              normal).  */
1486           bool must_underflow = false;
1487           bool may_underflow = false;
1488           for (size_t i = 0; i < tf->num_ret; i++)
1489             {
1490               switch (generic_outputs[i].type)
1491                 {
1492                 case gtype_fp:
1493                   round_real (all_res[i], all_exc_before[i], all_exc_after[i],
1494                               generic_outputs[i].value.f, f);
1495                   for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
1496                     {
1497                       merged_exc_before[m] |= all_exc_before[i][m];
1498                       merged_exc_after[m] |= all_exc_after[i][m];
1499                       if (!tf->exact)
1500                         {
1501                           must_underflow
1502                             |= ((all_exc_before[i][m]
1503                                  & (1U << exc_inexact)) != 0
1504                                 && (mpfr_cmpabs (generic_outputs[i].value.f,
1505                                                 fp_formats[f].subnorm_max)
1506                                     <= 0));
1507                           may_underflow
1508                             |= (!mpfr_zero_p (generic_outputs[i].value.f)
1509                                 && mpfr_cmpabs (generic_outputs[i].value.f,
1510                                                 fp_formats[f].min) <= 0);
1511                         }
1512                     }
1513                   break;
1514
1515                 case gtype_int:
1516                   if (ignore_output[i])
1517                     for (rounding_mode m = rm_first_mode;
1518                          m < rm_num_modes;
1519                          m++)
1520                       {
1521                         merged_exc_before[m] |= 1U << exc_invalid;
1522                         merged_exc_after[m] |= 1U << exc_invalid;
1523                       }
1524                   break;
1525
1526                 default:
1527                   abort ();
1528                 }
1529             }
1530           assert (may_underflow || !must_underflow);
1531           for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
1532             {
1533               bool before_after_matters
1534                 = tf->exact && merged_exc_before[m] != merged_exc_after[m];
1535               for (int after = 0; after <= 1; after++)
1536                 {
1537                   if (after == 1 && !before_after_matters)
1538                     continue;
1539                   const char *after_cond;
1540                   if (before_after_matters)
1541                     after_cond = (after
1542                                   ? ":after-rounding"
1543                                   : ":before-rounding");
1544                   else
1545                     after_cond = "";
1546                   unsigned int merged_exc = (after
1547                                              ? merged_exc_after[m]
1548                                              : merged_exc_before[m]);
1549                   if (fprintf (fp, "= %s %s %s%s%s", tf->name,
1550                                rounding_modes[m].name, fp_formats[f].name,
1551                                long_cond, after_cond) < 0)
1552                     error (EXIT_FAILURE, errno, "write to '%s'", filename);
1553                   /* Print inputs.  */
1554                   for (size_t i = 0; i < tf->num_args; i++)
1555                     output_generic_value (fp, filename, &inputs[i], false,
1556                                           tf->arg_types[i], f, long_bits);
1557                   if (fputs (" :", fp) < 0)
1558                     error (EXIT_FAILURE, errno, "write to '%s'", filename);
1559                   /* Print outputs.  */
1560                   bool must_erange = false;
1561                   for (size_t i = 0; i < tf->num_ret; i++)
1562                     {
1563                       generic_value g;
1564                       g.type = generic_outputs[i].type;
1565                       switch (g.type)
1566                         {
1567                         case gtype_fp:
1568                           if (mpfr_inf_p (all_res[i][m])
1569                               && (all_exc_before[i][m]
1570                                   & (1U << exc_overflow)) != 0)
1571                             must_erange = true;
1572                           if (mpfr_zero_p (all_res[i][m])
1573                               && (tf->exact
1574                                   || mpfr_zero_p (all_res[i][rm_tonearest]))
1575                               && (all_exc_before[i][m]
1576                                   & (1U << exc_underflow)) != 0)
1577                             must_erange = true;
1578                           mpfr_init2 (g.value.f, fp_formats[f].mant_dig);
1579                           assert_exact (mpfr_set (g.value.f, all_res[i][m],
1580                                                   MPFR_RNDN));
1581                           break;
1582
1583                         case gtype_int:
1584                           mpz_init (g.value.i);
1585                           mpz_set (g.value.i, generic_outputs[i].value.i);
1586                           break;
1587
1588                         default:
1589                           abort ();
1590                         }
1591                       output_generic_value (fp, filename, &g, ignore_output[i],
1592                                             tf->ret_types[i], f, long_bits);
1593                       generic_value_free (&g);
1594                     }
1595                   if (fputs (" :", fp) < 0)
1596                     error (EXIT_FAILURE, errno, "write to '%s'", filename);
1597                   /* Print miscellaneous flags (passed through from
1598                      input).  */
1599                   for (size_t i = 0; i < it->num_flags; i++)
1600                     switch (it->flags[i].type)
1601                       {
1602                       case flag_no_test_inline:
1603                       case flag_xfail:
1604                         if (fprintf (fp, " %s%s",
1605                                      input_flags[it->flags[i].type],
1606                                      (it->flags[i].cond
1607                                       ? it->flags[i].cond
1608                                       : "")) < 0)
1609                           error (EXIT_FAILURE, errno, "write to '%s'",
1610                                  filename);
1611                         break;
1612                       case flag_xfail_rounding:
1613                         if (m != rm_tonearest)
1614                           if (fprintf (fp, " xfail%s",
1615                                        (it->flags[i].cond
1616                                         ? it->flags[i].cond
1617                                         : "")) < 0)
1618                             error (EXIT_FAILURE, errno, "write to '%s'",
1619                                    filename);
1620                         break;
1621                       default:
1622                         break;
1623                       }
1624                   /* Print exception flags and compute errno
1625                      expectations where not already computed.  */
1626                   bool may_edom = false;
1627                   bool must_edom = false;
1628                   bool may_erange = must_erange || may_underflow;
1629                   for (fp_exception e = exc_first_exception;
1630                        e < exc_num_exceptions;
1631                        e++)
1632                     {
1633                       bool expect_e = (merged_exc & (1U << e)) != 0;
1634                       bool e_optional = false;
1635                       switch (e)
1636                         {
1637                         case exc_divbyzero:
1638                           if (expect_e)
1639                             may_erange = must_erange = true;
1640                           break;
1641
1642                         case exc_inexact:
1643                           if (!tf->exact)
1644                             e_optional = true;
1645                           break;
1646
1647                         case exc_invalid:
1648                           if (expect_e)
1649                             may_edom = must_edom = true;
1650                           break;
1651
1652                         case exc_overflow:
1653                           if (expect_e)
1654                             may_erange = true;
1655                           break;
1656
1657                         case exc_underflow:
1658                           if (expect_e)
1659                             may_erange = true;
1660                           if (must_underflow)
1661                             assert (expect_e);
1662                           if (may_underflow && !must_underflow)
1663                             e_optional = true;
1664                           break;
1665
1666                         default:
1667                           abort ();
1668                         }
1669                       if (e_optional)
1670                         {
1671                           if (fprintf (fp, " %s-ok", exceptions[e]) < 0)
1672                             error (EXIT_FAILURE, errno, "write to '%s'",
1673                                    filename);
1674                         }
1675                       else
1676                         {
1677                           if (expect_e)
1678                             if (fprintf (fp, " %s", exceptions[e]) < 0)
1679                               error (EXIT_FAILURE, errno, "write to '%s'",
1680                                      filename);
1681                           input_flag_type okflag;
1682                           okflag = (expect_e
1683                                     ? flag_missing_first
1684                                     : flag_spurious_first) + e;
1685                           for (size_t i = 0; i < it->num_flags; i++)
1686                             if (it->flags[i].type == okflag)
1687                               if (fprintf (fp, " %s-ok%s",
1688                                            exceptions[e],
1689                                            (it->flags[i].cond
1690                                             ? it->flags[i].cond
1691                                             : "")) < 0)
1692                                 error (EXIT_FAILURE, errno, "write to '%s'",
1693                                        filename);
1694                         }
1695                     }
1696                   /* Print errno expectations.  */
1697                   if (tf->complex_fn)
1698                     {
1699                       must_edom = false;
1700                       must_erange = false;
1701                     }
1702                   if (may_edom && !must_edom)
1703                     {
1704                       if (fputs (" errno-edom-ok", fp) < 0)
1705                         error (EXIT_FAILURE, errno, "write to '%s'",
1706                                filename);
1707                     }
1708                   else
1709                     {
1710                       if (must_edom)
1711                         if (fputs (" errno-edom", fp) < 0)
1712                           error (EXIT_FAILURE, errno, "write to '%s'",
1713                                  filename);
1714                       input_flag_type okflag = (must_edom
1715                                                 ? flag_missing_errno
1716                                                 : flag_spurious_errno);
1717                       for (size_t i = 0; i < it->num_flags; i++)
1718                         if (it->flags[i].type == okflag)
1719                           if (fprintf (fp, " errno-edom-ok%s",
1720                                        (it->flags[i].cond
1721                                         ? it->flags[i].cond
1722                                         : "")) < 0)
1723                             error (EXIT_FAILURE, errno, "write to '%s'",
1724                                    filename);
1725                     }
1726                   if (may_erange && !must_erange)
1727                     {
1728                       if (fputs (" errno-erange-ok", fp) < 0)
1729                         error (EXIT_FAILURE, errno, "write to '%s'",
1730                                filename);
1731                     }
1732                   else
1733                     {
1734                       if (must_erange)
1735                         if (fputs (" errno-erange", fp) < 0)
1736                           error (EXIT_FAILURE, errno, "write to '%s'",
1737                                  filename);
1738                       input_flag_type okflag = (must_erange
1739                                                 ? flag_missing_errno
1740                                                 : flag_spurious_errno);
1741                       for (size_t i = 0; i < it->num_flags; i++)
1742                         if (it->flags[i].type == okflag)
1743                           if (fprintf (fp, " errno-erange-ok%s",
1744                                        (it->flags[i].cond
1745                                         ? it->flags[i].cond
1746                                         : "")) < 0)
1747                             error (EXIT_FAILURE, errno, "write to '%s'",
1748                                    filename);
1749                     }
1750                   if (putc ('\n', fp) < 0)
1751                     error (EXIT_FAILURE, errno, "write to '%s'", filename);
1752                 }
1753             }
1754           for (size_t i = 0; i < tf->num_ret; i++)
1755             {
1756               if (generic_outputs[i].type == gtype_fp)
1757                 for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
1758                   mpfr_clear (all_res[i][m]);
1759             }
1760         }
1761     }
1762  out:
1763   for (size_t i = 0; i < tf->num_ret; i++)
1764     generic_value_free (&generic_outputs[i]);
1765 }
1766
1767 /* Generate test output data to FILENAME.  */
1768
1769 static void
1770 generate_output (const char *filename)
1771 {
1772   FILE *fp = fopen (filename, "w");
1773   if (fp == NULL)
1774     error (EXIT_FAILURE, errno, "open '%s'", filename);
1775   for (size_t i = 0; i < ARRAY_SIZE (test_functions); i++)
1776     {
1777       test_function *tf = &test_functions[i];
1778       for (size_t j = 0; j < tf->num_tests; j++)
1779         {
1780           input_test *it = &tf->tests[j];
1781           if (fputs (it->line, fp) < 0)
1782             error (EXIT_FAILURE, errno, "write to '%s'", filename);
1783           for (size_t k = 0; k < it->num_input_cases; k++)
1784             output_for_one_input_case (fp, filename, tf, it, it->inputs[k]);
1785         }
1786     }
1787   if (fclose (fp) != 0)
1788     error (EXIT_FAILURE, errno, "close '%s'", filename);
1789 }
1790
1791 int
1792 main (int argc, char **argv)
1793 {
1794   if (argc != 3)
1795     error (EXIT_FAILURE, 0, "usage: gen-auto-libm-tests <input> <output>");
1796   const char *input_filename = argv[1];
1797   const char *output_filename = argv[2];
1798   init_fp_formats ();
1799   read_input (input_filename);
1800   generate_output (output_filename);
1801   exit (EXIT_SUCCESS);
1802 }