math/gen-auto-libm-tests.c

   1 /* Generate expected output for libm tests with MPFR and MPC.
   2    Copyright (C) 2013-2014 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; "ignore-zero-inf-sign" indicates the the signs of
  88    zero and infinite results should be ignored; "xfail" indicates the
  89    test is disabled as expected to produce incorrect results,
  90    "xfail-rounding" indicates the test is disabled only in rounding
  91    modes other than round-to-nearest.  Otherwise, test flags are of
  92    the form "spurious-<exception>" and "missing-<exception>", for any
  93    exception ("overflow", "underflow", "inexact", "invalid",
  94    "divbyzero"), "spurious-errno" and "missing-errno", to indicate
  95    when tests are expected to deviate from the exception and errno
  96    settings corresponding to the mathematical results.  "xfail",
  97    "xfail-rounding", "spurious-" and "missing-" flags should be
  98    accompanied by a comment referring to an open bug in glibc
  99    Bugzilla.
 100
 101    The output file auto-libm-test-out contains the test lines from
 102    auto-libm-test-in, and, after the line for a given test, some
 103    number of output test lines.  An output test line is of the form "=
 104    function rounding-mode format input1 input2 ... : output1 output2
 105    ... : flags".  rounding-mode is "tonearest", "towardzero", "upward"
 106    or "downward".  format is a name from the floating_point_formats
 107    array, possibly followed by a sequence of ":flag" for flags from
 108    "long32" and "long64".  Inputs and outputs are specified as hex
 109    floats with the required suffix for the floating-point type, or
 110    plus_infty or minus_infty for infinite expected results, or as
 111    integer constant expressions (not necessarily with the right type)
 112    or IGNORE for integer inputs and outputs.  Flags are
 113    "no-test-inline", "ignore-zero-info-sign", "xfail", "<exception>",
 114    "<exception>-ok", "errno-<value>", "errno-<value>-ok", which may be
 115    unconditional or conditional.  "<exception>" indicates that a
 116    correct result means the given exception should be raised.
 117    "errno-<value>" indicates that a correct result means errno should
 118    be set to the given value.  "-ok" means not to test for the given
 119    exception or errno value (whether because it was marked as possibly
 120    missing or spurious, or because the calculation of correct results
 121    indicated it was optional).  Conditions "before-rounding" and
 122    "after-rounding" indicate tests where expectations for underflow
 123    exceptions depend on how the architecture detects tininess.  */
 124
 125 #define _GNU_SOURCE
 126
 127 #include <assert.h>
 128 #include <ctype.h>
 129 #include <errno.h>
 130 #include <error.h>
 131 #include <stdbool.h>
 132 #include <stdint.h>
 133 #include <stdio.h>
 134 #include <stdlib.h>
 135 #include <string.h>
 136
 137 #include <gmp.h>
 138 #include <mpfr.h>
 139 #include <mpc.h>
 140
 141 #define ARRAY_SIZE(A) (sizeof (A) / sizeof ((A)[0]))
 142
 143 /* The supported floating-point formats.  */
 144 typedef enum
 145   {
 146     fp_flt_32,
 147     fp_dbl_64,
 148     fp_ldbl_96_intel,
 149     fp_ldbl_96_m68k,
 150     fp_ldbl_128,
 151     fp_ldbl_128ibm,
 152     fp_num_formats,
 153     fp_first_format = 0
 154   } fp_format;
 155
 156 /* Structure describing a single floating-point format.  */
 157 typedef struct
 158 {
 159   /* The name of the format.  */
 160   const char *name;
 161   /* The suffix to use on floating-point constants with this
 162      format.  */
 163   const char *suffix;
 164   /* A string for the largest normal value, or NULL for IEEE formats
 165      where this can be determined automatically.  */
 166   const char *max_string;
 167   /* The number of mantissa bits.  */
 168   int mant_dig;
 169   /* The least N such that 2^N overflows.  */
 170   int max_exp;
 171   /* One more than the least N such that 2^N is normal.  */
 172   int min_exp;
 173   /* The largest normal value.  */
 174   mpfr_t max;
 175   /* The value 0.5ulp above the least positive normal value.  */
 176   mpfr_t min_plus_half;
 177   /* The least positive normal value, 2^(MIN_EXP-1).  */
 178   mpfr_t min;
 179   /* The greatest positive subnormal value.  */
 180   mpfr_t subnorm_max;
 181   /* The least positive subnormal value, 2^(MIN_EXP-MANT_DIG).  */
 182   mpfr_t subnorm_min;
 183 } fp_format_desc;
 184
 185 /* List of floating-point formats, in the same order as the fp_format
 186    enumeration.  */
 187 static fp_format_desc fp_formats[fp_num_formats] =
 188   {
 189     { "flt-32", "f", NULL, 24, 128, -125, {}, {}, {}, {}, {} },
 190     { "dbl-64", "", NULL, 53, 1024, -1021, {}, {}, {}, {}, {} },
 191     { "ldbl-96-intel", "L", NULL, 64, 16384, -16381, {}, {}, {}, {}, {} },
 192     { "ldbl-96-m68k", "L", NULL, 64, 16384, -16382, {}, {}, {}, {}, {} },
 193     { "ldbl-128", "L", NULL, 113, 16384, -16381, {}, {}, {}, {}, {} },
 194     { "ldbl-128ibm", "L", "0x1.fffffffffffff7ffffffffffff8p+1023",
 195       106, 1024, -968, {}, {}, {}, {}, {} },
 196   };
 197
 198 /* The supported rounding modes.  */
 199 typedef enum
 200   {
 201     rm_downward,
 202     rm_tonearest,
 203     rm_towardzero,
 204     rm_upward,
 205     rm_num_modes,
 206     rm_first_mode = 0
 207   } rounding_mode;
 208
 209 /* Structure describing a single rounding mode.  */
 210 typedef struct
 211 {
 212   /* The name of the rounding mode.  */
 213   const char *name;
 214   /* The MPFR rounding mode.  */
 215   mpfr_rnd_t mpfr_mode;
 216   /* The MPC rounding mode.  */
 217   mpc_rnd_t mpc_mode;
 218 } rounding_mode_desc;
 219
 220 /* List of rounding modes, in the same order as the rounding_mode
 221    enumeration.  */
 222 static const rounding_mode_desc rounding_modes[rm_num_modes] =
 223   {
 224     { "downward", MPFR_RNDD, MPC_RNDDD },
 225     { "tonearest", MPFR_RNDN, MPC_RNDNN },
 226     { "towardzero", MPFR_RNDZ, MPC_RNDZZ },
 227     { "upward", MPFR_RNDU, MPC_RNDUU },
 228   };
 229
 230 /* The supported exceptions.  */
 231 typedef enum
 232   {
 233     exc_divbyzero,
 234     exc_inexact,
 235     exc_invalid,
 236     exc_overflow,
 237     exc_underflow,
 238     exc_num_exceptions,
 239     exc_first_exception = 0
 240   } fp_exception;
 241
 242 /* List of exceptions, in the same order as the fp_exception
 243    enumeration.  */
 244 static const char *const exceptions[exc_num_exceptions] =
 245   {
 246     "divbyzero",
 247     "inexact",
 248     "invalid",
 249     "overflow",
 250     "underflow",
 251   };
 252
 253 /* The internal precision to use for most MPFR calculations, which
 254    must be at least 2 more than the greatest precision of any
 255    supported floating-point format.  */
 256 static int internal_precision;
 257
 258 /* A value that overflows all supported floating-point formats.  */
 259 static mpfr_t global_max;
 260
 261 /* A value that is at most half the least subnormal in any
 262    floating-point format and so is rounded the same way as all
 263    sufficiently small positive values.  */
 264 static mpfr_t global_min;
 265
 266 /* The maximum number of (real or integer) arguments to a function
 267    handled by this program (complex arguments count as two real
 268    arguments).  */
 269 #define MAX_NARGS 4
 270
 271 /* The maximum number of (real or integer) return values from a
 272    function handled by this program.  */
 273 #define MAX_NRET 2
 274
 275 /* A type of a function argument or return value.  */
 276 typedef enum
 277   {
 278     /* No type (not a valid argument or return value).  */
 279     type_none,
 280     /* A floating-point value with the type corresponding to that of
 281        the function.  */
 282     type_fp,
 283     /* An integer value of type int.  */
 284     type_int,
 285     /* An integer value of type long.  */
 286     type_long,
 287     /* An integer value of type long long.  */
 288     type_long_long,
 289   } arg_ret_type;
 290
 291 /* A type of a generic real or integer value.  */
 292 typedef enum
 293   {
 294     /* No type.  */
 295     gtype_none,
 296     /* Floating-point (represented with MPFR).  */
 297     gtype_fp,
 298     /* Integer (represented with GMP).  */
 299     gtype_int,
 300   } generic_value_type;
 301
 302 /* A generic value (argument or result).  */
 303 typedef struct
 304 {
 305   /* The type of this value.  */
 306   generic_value_type type;
 307   /* Its value.  */
 308   union
 309   {
 310     mpfr_t f;
 311     mpz_t i;
 312   } value;
 313 } generic_value;
 314
 315 /* A type of input flag.  */
 316 typedef enum
 317   {
 318     flag_no_test_inline,
 319     flag_ignore_zero_inf_sign,
 320     flag_xfail,
 321     flag_xfail_rounding,
 322     /* The "spurious" and "missing" flags must be in the same order as
 323        the fp_exception enumeration.  */
 324     flag_spurious_divbyzero,
 325     flag_spurious_inexact,
 326     flag_spurious_invalid,
 327     flag_spurious_overflow,
 328     flag_spurious_underflow,
 329     flag_spurious_errno,
 330     flag_missing_divbyzero,
 331     flag_missing_inexact,
 332     flag_missing_invalid,
 333     flag_missing_overflow,
 334     flag_missing_underflow,
 335     flag_missing_errno,
 336     num_input_flag_types,
 337     flag_first_flag = 0,
 338     flag_spurious_first = flag_spurious_divbyzero,
 339     flag_missing_first = flag_missing_divbyzero
 340   } input_flag_type;
 341
 342 /* List of flags, in the same order as the input_flag_type
 343    enumeration.  */
 344 static const char *const input_flags[num_input_flag_types] =
 345   {
 346     "no-test-inline",
 347     "ignore-zero-inf-sign",
 348     "xfail",
 349     "xfail-rounding",
 350     "spurious-divbyzero",
 351     "spurious-inexact",
 352     "spurious-invalid",
 353     "spurious-overflow",
 354     "spurious-underflow",
 355     "spurious-errno",
 356     "missing-divbyzero",
 357     "missing-inexact",
 358     "missing-invalid",
 359     "missing-overflow",
 360     "missing-underflow",
 361     "missing-errno",
 362   };
 363
 364 /* An input flag, possibly conditional.  */
 365 typedef struct
 366 {
 367   /* The type of this flag.  */
 368   input_flag_type type;
 369   /* The conditions on this flag, as a string ":cond1:cond2..." or
 370      NULL.  */
 371   const char *cond;
 372 } input_flag;
 373
 374 /* Structure describing a single test from the input file (which may
 375    expand into many tests in the output).  The choice of function,
 376    which implies the numbers and types of arguments and results, is
 377    implicit rather than stored in this structure (except as part of
 378    the source line).  */
 379 typedef struct
 380 {
 381   /* The text of the input line describing the test, including the
 382      trailing newline.  */
 383   const char *line;
 384   /* The number of combinations of interpretations of input values for
 385      different floating-point formats and rounding modes.  */
 386   size_t num_input_cases;
 387   /* The corresponding lists of inputs.  */
 388   generic_value **inputs;
 389   /* The number of flags for this test.  */
 390   size_t num_flags;
 391   /* The corresponding list of flags.  */
 392   input_flag *flags;
 393   /* The old output for this test.  */
 394   const char *old_output;
 395 } input_test;
 396
 397 /* Ways to calculate a function.  */
 398 typedef enum
 399   {
 400     /* MPFR function with a single argument and result.  */
 401     mpfr_f_f,
 402     /* MPFR function with two arguments and one result.  */
 403     mpfr_ff_f,
 404     /* MPFR function with three arguments and one result.  */
 405     mpfr_fff_f,
 406     /* MPFR function with a single argument and floating-point and
 407        integer results.  */
 408     mpfr_f_f1,
 409     /* MPFR function with integer and floating-point arguments and one
 410        result.  */
 411     mpfr_if_f,
 412     /* MPFR function with a single argument and two floating-point
 413        results.  */
 414     mpfr_f_11,
 415     /* MPC function with a single complex argument and one real
 416        result.  */
 417     mpc_c_f,
 418     /* MPC function with a single complex argument and one complex
 419        result.  */
 420     mpc_c_c,
 421     /* MPC function with two complex arguments and one complex
 422        result.  */
 423     mpc_cc_c,
 424   } func_calc_method;
 425
 426 /* Description of how to calculate a function.  */
 427 typedef struct
 428 {
 429   /* Which method is used to calculate the function.  */
 430   func_calc_method method;
 431   /* The specific function called.  */
 432   union
 433   {
 434     int (*mpfr_f_f) (mpfr_t, const mpfr_t, mpfr_rnd_t);
 435     int (*mpfr_ff_f) (mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t);
 436     int (*mpfr_fff_f) (mpfr_t, const mpfr_t, const mpfr_t, const mpfr_t,
 437                        mpfr_rnd_t);
 438     int (*mpfr_f_f1) (mpfr_t, int *, const mpfr_t, mpfr_rnd_t);
 439     int (*mpfr_if_f) (mpfr_t, long, const mpfr_t, mpfr_rnd_t);
 440     int (*mpfr_f_11) (mpfr_t, mpfr_t, const mpfr_t, mpfr_rnd_t);
 441     int (*mpc_c_f) (mpfr_t, const mpc_t, mpfr_rnd_t);
 442     int (*mpc_c_c) (mpc_t, const mpc_t, mpc_rnd_t);
 443     int (*mpc_cc_c) (mpc_t, const mpc_t, const mpc_t, mpc_rnd_t);
 444   } func;
 445 } func_calc_desc;
 446
 447 /* Structure describing a function handled by this program.  */
 448 typedef struct
 449 {
 450   /* The name of the function.  */
 451   const char *name;
 452   /* The number of arguments.  */
 453   size_t num_args;
 454   /* The types of the arguments.  */
 455   arg_ret_type arg_types[MAX_NARGS];
 456   /* The number of return values.  */
 457   size_t num_ret;
 458   /* The types of the return values.  */
 459   arg_ret_type ret_types[MAX_NRET];
 460   /* Whether the function has exactly determined results and
 461      exceptions.  */
 462   bool exact;
 463   /* Whether the function is a complex function, so errno setting is
 464      optional.  */
 465   bool complex_fn;
 466   /* Whether to treat arguments given as floating-point constants as
 467      exact only, rather than rounding them up and down to all
 468      formats.  */
 469   bool exact_args;
 470   /* How to calculate this function.  */
 471   func_calc_desc calc;
 472   /* The number of tests allocated for this function.  */
 473   size_t num_tests_alloc;
 474   /* The number of tests for this function.  */
 475   size_t num_tests;
 476   /* The tests themselves.  */
 477   input_test *tests;
 478 } test_function;
 479
 480 #define ARGS1(T1) 1, { T1 }
 481 #define ARGS2(T1, T2) 2, { T1, T2 }
 482 #define ARGS3(T1, T2, T3) 3, { T1, T2, T3 }
 483 #define ARGS4(T1, T2, T3, T4) 4, { T1, T2, T3, T4 }
 484 #define RET1(T1) 1, { T1 }
 485 #define RET2(T1, T2) 2, { T1, T2 }
 486 #define CALC(TYPE, FN) { TYPE, { .TYPE = FN } }
 487 #define FUNC(NAME, ARGS, RET, EXACT, COMPLEX_FN, EXACT_ARGS, CALC)      \
 488   {                                                                     \
 489     NAME, ARGS, RET, EXACT, COMPLEX_FN, EXACT_ARGS, CALC, 0, 0, NULL    \
 490   }
 491
 492 #define FUNC_mpfr_f_f(NAME, MPFR_FUNC, EXACT)                           \
 493   FUNC (NAME, ARGS1 (type_fp), RET1 (type_fp), EXACT, false, false,     \
 494         CALC (mpfr_f_f, MPFR_FUNC))
 495 #define FUNC_mpfr_ff_f(NAME, MPFR_FUNC, EXACT)                          \
 496   FUNC (NAME, ARGS2 (type_fp, type_fp), RET1 (type_fp), EXACT, false,   \
 497         false, CALC (mpfr_ff_f, MPFR_FUNC))
 498 #define FUNC_mpfr_if_f(NAME, MPFR_FUNC, EXACT)                          \
 499   FUNC (NAME, ARGS2 (type_int, type_fp), RET1 (type_fp), EXACT, false,  \
 500         false, CALC (mpfr_if_f, MPFR_FUNC))
 501 #define FUNC_mpc_c_f(NAME, MPFR_FUNC, EXACT)                            \
 502   FUNC (NAME, ARGS2 (type_fp, type_fp), RET1 (type_fp), EXACT, true,    \
 503         false, CALC (mpc_c_f, MPFR_FUNC))
 504 #define FUNC_mpc_c_c(NAME, MPFR_FUNC, EXACT)                            \
 505   FUNC (NAME, ARGS2 (type_fp, type_fp), RET2 (type_fp, type_fp), EXACT, \
 506         true, false, CALC (mpc_c_c, MPFR_FUNC))
 507
 508 /* List of functions handled by this program.  */
 509 static test_function test_functions[] =
 510   {
 511     FUNC_mpfr_f_f ("acos", mpfr_acos, false),
 512     FUNC_mpfr_f_f ("acosh", mpfr_acosh, false),
 513     FUNC_mpfr_f_f ("asin", mpfr_asin, false),
 514     FUNC_mpfr_f_f ("asinh", mpfr_asinh, false),
 515     FUNC_mpfr_f_f ("atan", mpfr_atan, false),
 516     FUNC_mpfr_ff_f ("atan2", mpfr_atan2, false),
 517     FUNC_mpfr_f_f ("atanh", mpfr_atanh, false),
 518     FUNC_mpc_c_f ("cabs", mpc_abs, false),
 519     FUNC_mpc_c_c ("cacos", mpc_acos, false),
 520     FUNC_mpc_c_c ("cacosh", mpc_acosh, false),
 521     FUNC_mpc_c_f ("carg", mpc_arg, false),
 522     FUNC_mpc_c_c ("casin", mpc_asin, false),
 523     FUNC_mpc_c_c ("casinh", mpc_asinh, false),
 524     FUNC_mpc_c_c ("catan", mpc_atan, false),
 525     FUNC_mpc_c_c ("catanh", mpc_atanh, false),
 526     FUNC_mpfr_f_f ("cbrt", mpfr_cbrt, false),
 527     FUNC_mpc_c_c ("ccos", mpc_cos, false),
 528     FUNC_mpc_c_c ("ccosh", mpc_cosh, false),
 529     FUNC_mpc_c_c ("cexp", mpc_exp, false),
 530     FUNC_mpc_c_c ("clog", mpc_log, false),
 531     FUNC_mpc_c_c ("clog10", mpc_log10, false),
 532     FUNC_mpfr_f_f ("cos", mpfr_cos, false),
 533     FUNC_mpfr_f_f ("cosh", mpfr_cosh, false),
 534     FUNC ("cpow", ARGS4 (type_fp, type_fp, type_fp, type_fp),
 535           RET2 (type_fp, type_fp), false, true, false,
 536           CALC (mpc_cc_c, mpc_pow)),
 537     FUNC_mpc_c_c ("csin", mpc_sin, false),
 538     FUNC_mpc_c_c ("csinh", mpc_sinh, false),
 539     FUNC_mpc_c_c ("csqrt", mpc_sqrt, false),
 540     FUNC_mpc_c_c ("ctan", mpc_tan, false),
 541     FUNC_mpc_c_c ("ctanh", mpc_tanh, false),
 542     FUNC_mpfr_f_f ("erf", mpfr_erf, false),
 543     FUNC_mpfr_f_f ("erfc", mpfr_erfc, false),
 544     FUNC_mpfr_f_f ("exp", mpfr_exp, false),
 545     FUNC_mpfr_f_f ("exp10", mpfr_exp10, false),
 546     FUNC_mpfr_f_f ("exp2", mpfr_exp2, false),
 547     FUNC_mpfr_f_f ("expm1", mpfr_expm1, false),
 548     FUNC ("fma", ARGS3 (type_fp, type_fp, type_fp), RET1 (type_fp),
 549           true, false, true, CALC (mpfr_fff_f, mpfr_fma)),
 550     FUNC_mpfr_ff_f ("hypot", mpfr_hypot, false),
 551     FUNC_mpfr_f_f ("j0", mpfr_j0, false),
 552     FUNC_mpfr_f_f ("j1", mpfr_j1, false),
 553     FUNC_mpfr_if_f ("jn", mpfr_jn, false),
 554     FUNC ("lgamma", ARGS1 (type_fp), RET2 (type_fp, type_int), false, false,
 555           false, CALC (mpfr_f_f1, mpfr_lgamma)),
 556     FUNC_mpfr_f_f ("log", mpfr_log, false),
 557     FUNC_mpfr_f_f ("log10", mpfr_log10, false),
 558     FUNC_mpfr_f_f ("log1p", mpfr_log1p, false),
 559     FUNC_mpfr_f_f ("log2", mpfr_log2, false),
 560     FUNC_mpfr_ff_f ("pow", mpfr_pow, false),
 561     FUNC_mpfr_f_f ("sin", mpfr_sin, false),
 562     FUNC ("sincos", ARGS1 (type_fp), RET2 (type_fp, type_fp), false, false,
 563           false, CALC (mpfr_f_11, mpfr_sin_cos)),
 564     FUNC_mpfr_f_f ("sinh", mpfr_sinh, false),
 565     FUNC_mpfr_f_f ("sqrt", mpfr_sqrt, true),
 566     FUNC_mpfr_f_f ("tan", mpfr_tan, false),
 567     FUNC_mpfr_f_f ("tanh", mpfr_tanh, false),
 568     FUNC_mpfr_f_f ("tgamma", mpfr_gamma, false),
 569     FUNC_mpfr_f_f ("y0", mpfr_y0, false),
 570     FUNC_mpfr_f_f ("y1", mpfr_y1, false),
 571     FUNC_mpfr_if_f ("yn", mpfr_yn, false),
 572   };
 573
 574 /* Allocate memory, with error checking.  */
 575
 576 static void *
 577 xmalloc (size_t n)
 578 {
 579   void *p = malloc (n);
 580   if (p == NULL)
 581     error (EXIT_FAILURE, errno, "xmalloc failed");
 582   return p;
 583 }
 584
 585 static void *
 586 xrealloc (void *p, size_t n)
 587 {
 588   p = realloc (p, n);
 589   if (p == NULL)
 590     error (EXIT_FAILURE, errno, "xrealloc failed");
 591   return p;
 592 }
 593
 594 static char *
 595 xstrdup (const char *s)
 596 {
 597   char *p = strdup (s);
 598   if (p == NULL)
 599     error (EXIT_FAILURE, errno, "xstrdup failed");
 600   return p;
 601 }
 602
 603 /* Assert that the result of an MPFR operation was exact; that is,
 604    that the returned ternary value was 0.  */
 605
 606 static void
 607 assert_exact (int i)
 608 {
 609   assert (i == 0);
 610 }
 611
 612 /* Return the generic type of an argument or return value type T.  */
 613
 614 static generic_value_type
 615 generic_arg_ret_type (arg_ret_type t)
 616 {
 617   switch (t)
 618     {
 619     case type_fp:
 620       return gtype_fp;
 621
 622     case type_int:
 623     case type_long:
 624     case type_long_long:
 625       return gtype_int;
 626
 627     default:
 628       abort ();
 629     }
 630 }
 631
 632 /* Free a generic_value *V.  */
 633
 634 static void
 635 generic_value_free (generic_value *v)
 636 {
 637   switch (v->type)
 638     {
 639     case gtype_fp:
 640       mpfr_clear (v->value.f);
 641       break;
 642
 643     case gtype_int:
 644       mpz_clear (v->value.i);
 645       break;
 646
 647     default:
 648       abort ();
 649     }
 650 }
 651
 652 /* Copy a generic_value *SRC to *DEST.  */
 653
 654 static void
 655 generic_value_copy (generic_value *dest, const generic_value *src)
 656 {
 657   dest->type = src->type;
 658   switch (src->type)
 659     {
 660     case gtype_fp:
 661       mpfr_init (dest->value.f);
 662       assert_exact (mpfr_set (dest->value.f, src->value.f, MPFR_RNDN));
 663       break;
 664
 665     case gtype_int:
 666       mpz_init (dest->value.i);
 667       mpz_set (dest->value.i, src->value.i);
 668       break;
 669
 670     default:
 671       abort ();
 672     }
 673 }
 674
 675 /* Initialize data for floating-point formats.  */
 676
 677 static void
 678 init_fp_formats ()
 679 {
 680   int global_max_exp = 0, global_min_subnorm_exp = 0;
 681   for (fp_format f = fp_first_format; f < fp_num_formats; f++)
 682     {
 683       if (fp_formats[f].mant_dig + 2 > internal_precision)
 684         internal_precision = fp_formats[f].mant_dig + 2;
 685       if (fp_formats[f].max_exp > global_max_exp)
 686         global_max_exp = fp_formats[f].max_exp;
 687       int min_subnorm_exp = fp_formats[f].min_exp - fp_formats[f].mant_dig;
 688       if (min_subnorm_exp < global_min_subnorm_exp)
 689         global_min_subnorm_exp = min_subnorm_exp;
 690       mpfr_init2 (fp_formats[f].max, fp_formats[f].mant_dig);
 691       if (fp_formats[f].max_string != NULL)
 692         {
 693           char *ep = NULL;
 694           assert_exact (mpfr_strtofr (fp_formats[f].max,
 695                                       fp_formats[f].max_string,
 696                                       &ep, 0, MPFR_RNDN));
 697           assert (*ep == 0);
 698         }
 699       else
 700         {
 701           assert_exact (mpfr_set_ui_2exp (fp_formats[f].max, 1,
 702                                           fp_formats[f].max_exp,
 703                                           MPFR_RNDN));
 704           mpfr_nextbelow (fp_formats[f].max);
 705         }
 706       mpfr_init2 (fp_formats[f].min, fp_formats[f].mant_dig);
 707       assert_exact (mpfr_set_ui_2exp (fp_formats[f].min, 1,
 708                                       fp_formats[f].min_exp - 1,
 709                                       MPFR_RNDN));
 710       mpfr_init2 (fp_formats[f].min_plus_half, fp_formats[f].mant_dig + 1);
 711       assert_exact (mpfr_set (fp_formats[f].min_plus_half,
 712                               fp_formats[f].min, MPFR_RNDN));
 713       mpfr_nextabove (fp_formats[f].min_plus_half);
 714       mpfr_init2 (fp_formats[f].subnorm_max, fp_formats[f].mant_dig);
 715       assert_exact (mpfr_set (fp_formats[f].subnorm_max, fp_formats[f].min,
 716                               MPFR_RNDN));
 717       mpfr_nextbelow (fp_formats[f].subnorm_max);
 718       mpfr_nextbelow (fp_formats[f].subnorm_max);
 719       mpfr_init2 (fp_formats[f].subnorm_min, fp_formats[f].mant_dig);
 720       assert_exact (mpfr_set_ui_2exp (fp_formats[f].subnorm_min, 1,
 721                                       min_subnorm_exp, MPFR_RNDN));
 722     }
 723   mpfr_set_default_prec (internal_precision);
 724   mpfr_init (global_max);
 725   assert_exact (mpfr_set_ui_2exp (global_max, 1, global_max_exp, MPFR_RNDN));
 726   mpfr_init (global_min);
 727   assert_exact (mpfr_set_ui_2exp (global_min, 1, global_min_subnorm_exp - 1,
 728                                   MPFR_RNDN));
 729 }
 730
 731 /* Fill in mpfr_t values for special strings in input arguments.  */
 732
 733 static size_t
 734 special_fill_max (mpfr_t res0, mpfr_t res1 __attribute__ ((unused)),
 735                   fp_format format)
 736 {
 737   mpfr_init2 (res0, fp_formats[format].mant_dig);
 738   assert_exact (mpfr_set (res0, fp_formats[format].max, MPFR_RNDN));
 739   return 1;
 740 }
 741
 742 static size_t
 743 special_fill_minus_max (mpfr_t res0, mpfr_t res1 __attribute__ ((unused)),
 744                         fp_format format)
 745 {
 746   mpfr_init2 (res0, fp_formats[format].mant_dig);
 747   assert_exact (mpfr_neg (res0, fp_formats[format].max, MPFR_RNDN));
 748   return 1;
 749 }
 750
 751 static size_t
 752 special_fill_min (mpfr_t res0, mpfr_t res1 __attribute__ ((unused)),
 753                   fp_format format)
 754 {
 755   mpfr_init2 (res0, fp_formats[format].mant_dig);
 756   assert_exact (mpfr_set (res0, fp_formats[format].min, MPFR_RNDN));
 757   return 1;
 758 }
 759
 760 static size_t
 761 special_fill_minus_min (mpfr_t res0, mpfr_t res1 __attribute__ ((unused)),
 762                         fp_format format)
 763 {
 764   mpfr_init2 (res0, fp_formats[format].mant_dig);
 765   assert_exact (mpfr_neg (res0, fp_formats[format].min, MPFR_RNDN));
 766   return 1;
 767 }
 768
 769 static size_t
 770 special_fill_min_subnorm (mpfr_t res0, mpfr_t res1 __attribute__ ((unused)),
 771                           fp_format format)
 772 {
 773   mpfr_init2 (res0, fp_formats[format].mant_dig);
 774   assert_exact (mpfr_set (res0, fp_formats[format].subnorm_min, MPFR_RNDN));
 775   return 1;
 776 }
 777
 778 static size_t
 779 special_fill_minus_min_subnorm (mpfr_t res0,
 780                                 mpfr_t res1 __attribute__ ((unused)),
 781                                 fp_format format)
 782 {
 783   mpfr_init2 (res0, fp_formats[format].mant_dig);
 784   assert_exact (mpfr_neg (res0, fp_formats[format].subnorm_min, MPFR_RNDN));
 785   return 1;
 786 }
 787
 788 static size_t
 789 special_fill_min_subnorm_p120 (mpfr_t res0,
 790                                mpfr_t res1 __attribute__ ((unused)),
 791                                fp_format format)
 792 {
 793   mpfr_init2 (res0, fp_formats[format].mant_dig);
 794   assert_exact (mpfr_mul_2ui (res0, fp_formats[format].subnorm_min,
 795                               120, MPFR_RNDN));
 796   return 1;
 797 }
 798
 799 static size_t
 800 special_fill_pi (mpfr_t res0, mpfr_t res1, fp_format format)
 801 {
 802   mpfr_init2 (res0, fp_formats[format].mant_dig);
 803   mpfr_const_pi (res0, MPFR_RNDU);
 804   mpfr_init2 (res1, fp_formats[format].mant_dig);
 805   mpfr_const_pi (res1, MPFR_RNDD);
 806   return 2;
 807 }
 808
 809 static size_t
 810 special_fill_minus_pi (mpfr_t res0, mpfr_t res1, fp_format format)
 811 {
 812   mpfr_init2 (res0, fp_formats[format].mant_dig);
 813   mpfr_const_pi (res0, MPFR_RNDU);
 814   assert_exact (mpfr_neg (res0, res0, MPFR_RNDN));
 815   mpfr_init2 (res1, fp_formats[format].mant_dig);
 816   mpfr_const_pi (res1, MPFR_RNDD);
 817   assert_exact (mpfr_neg (res1, res1, MPFR_RNDN));
 818   return 2;
 819 }
 820
 821 static size_t
 822 special_fill_pi_2 (mpfr_t res0, mpfr_t res1, fp_format format)
 823 {
 824   mpfr_init2 (res0, fp_formats[format].mant_dig);
 825   mpfr_const_pi (res0, MPFR_RNDU);
 826   assert_exact (mpfr_div_ui (res0, res0, 2, MPFR_RNDN));
 827   mpfr_init2 (res1, fp_formats[format].mant_dig);
 828   mpfr_const_pi (res1, MPFR_RNDD);
 829   assert_exact (mpfr_div_ui (res1, res1, 2, MPFR_RNDN));
 830   return 2;
 831 }
 832
 833 static size_t
 834 special_fill_minus_pi_2 (mpfr_t res0, mpfr_t res1, fp_format format)
 835 {
 836   mpfr_init2 (res0, fp_formats[format].mant_dig);
 837   mpfr_const_pi (res0, MPFR_RNDU);
 838   assert_exact (mpfr_div_ui (res0, res0, 2, MPFR_RNDN));
 839   assert_exact (mpfr_neg (res0, res0, MPFR_RNDN));
 840   mpfr_init2 (res1, fp_formats[format].mant_dig);
 841   mpfr_const_pi (res1, MPFR_RNDD);
 842   assert_exact (mpfr_div_ui (res1, res1, 2, MPFR_RNDN));
 843   assert_exact (mpfr_neg (res1, res1, MPFR_RNDN));
 844   return 2;
 845 }
 846
 847 static size_t
 848 special_fill_pi_4 (mpfr_t res0, mpfr_t res1, fp_format format)
 849 {
 850   mpfr_init2 (res0, fp_formats[format].mant_dig);
 851   assert_exact (mpfr_set_si (res0, 1, MPFR_RNDN));
 852   mpfr_atan (res0, res0, MPFR_RNDU);
 853   mpfr_init2 (res1, fp_formats[format].mant_dig);
 854   assert_exact (mpfr_set_si (res1, 1, MPFR_RNDN));
 855   mpfr_atan (res1, res1, MPFR_RNDD);
 856   return 2;
 857 }
 858
 859 static size_t
 860 special_fill_pi_6 (mpfr_t res0, mpfr_t res1, fp_format format)
 861 {
 862   mpfr_init2 (res0, fp_formats[format].mant_dig);
 863   assert_exact (mpfr_set_si_2exp (res0, 1, -1, MPFR_RNDN));
 864   mpfr_asin (res0, res0, MPFR_RNDU);
 865   mpfr_init2 (res1, fp_formats[format].mant_dig);
 866   assert_exact (mpfr_set_si_2exp (res1, 1, -1, MPFR_RNDN));
 867   mpfr_asin (res1, res1, MPFR_RNDD);
 868   return 2;
 869 }
 870
 871 static size_t
 872 special_fill_minus_pi_6 (mpfr_t res0, mpfr_t res1, fp_format format)
 873 {
 874   mpfr_init2 (res0, fp_formats[format].mant_dig);
 875   assert_exact (mpfr_set_si_2exp (res0, -1, -1, MPFR_RNDN));
 876   mpfr_asin (res0, res0, MPFR_RNDU);
 877   mpfr_init2 (res1, fp_formats[format].mant_dig);
 878   assert_exact (mpfr_set_si_2exp (res1, -1, -1, MPFR_RNDN));
 879   mpfr_asin (res1, res1, MPFR_RNDD);
 880   return 2;
 881 }
 882
 883 static size_t
 884 special_fill_pi_3 (mpfr_t res0, mpfr_t res1, fp_format format)
 885 {
 886   mpfr_init2 (res0, fp_formats[format].mant_dig);
 887   assert_exact (mpfr_set_si_2exp (res0, 1, -1, MPFR_RNDN));
 888   mpfr_acos (res0, res0, MPFR_RNDU);
 889   mpfr_init2 (res1, fp_formats[format].mant_dig);
 890   assert_exact (mpfr_set_si_2exp (res1, 1, -1, MPFR_RNDN));
 891   mpfr_acos (res1, res1, MPFR_RNDD);
 892   return 2;
 893 }
 894
 895 static size_t
 896 special_fill_2pi_3 (mpfr_t res0, mpfr_t res1, fp_format format)
 897 {
 898   mpfr_init2 (res0, fp_formats[format].mant_dig);
 899   assert_exact (mpfr_set_si_2exp (res0, -1, -1, MPFR_RNDN));
 900   mpfr_acos (res0, res0, MPFR_RNDU);
 901   mpfr_init2 (res1, fp_formats[format].mant_dig);
 902   assert_exact (mpfr_set_si_2exp (res1, -1, -1, MPFR_RNDN));
 903   mpfr_acos (res1, res1, MPFR_RNDD);
 904   return 2;
 905 }
 906
 907 static size_t
 908 special_fill_2pi (mpfr_t res0, mpfr_t res1, fp_format format)
 909 {
 910   mpfr_init2 (res0, fp_formats[format].mant_dig);
 911   mpfr_const_pi (res0, MPFR_RNDU);
 912   assert_exact (mpfr_mul_ui (res0, res0, 2, MPFR_RNDN));
 913   mpfr_init2 (res1, fp_formats[format].mant_dig);
 914   mpfr_const_pi (res1, MPFR_RNDD);
 915   assert_exact (mpfr_mul_ui (res1, res1, 2, MPFR_RNDN));
 916   return 2;
 917 }
 918
 919 static size_t
 920 special_fill_e (mpfr_t res0, mpfr_t res1, fp_format format)
 921 {
 922   mpfr_init2 (res0, fp_formats[format].mant_dig);
 923   assert_exact (mpfr_set_si (res0, 1, MPFR_RNDN));
 924   mpfr_exp (res0, res0, MPFR_RNDU);
 925   mpfr_init2 (res1, fp_formats[format].mant_dig);
 926   assert_exact (mpfr_set_si (res1, 1, MPFR_RNDN));
 927   mpfr_exp (res1, res1, MPFR_RNDD);
 928   return 2;
 929 }
 930
 931 static size_t
 932 special_fill_1_e (mpfr_t res0, mpfr_t res1, fp_format format)
 933 {
 934   mpfr_init2 (res0, fp_formats[format].mant_dig);
 935   assert_exact (mpfr_set_si (res0, -1, MPFR_RNDN));
 936   mpfr_exp (res0, res0, MPFR_RNDU);
 937   mpfr_init2 (res1, fp_formats[format].mant_dig);
 938   assert_exact (mpfr_set_si (res1, -1, MPFR_RNDN));
 939   mpfr_exp (res1, res1, MPFR_RNDD);
 940   return 2;
 941 }
 942
 943 static size_t
 944 special_fill_e_minus_1 (mpfr_t res0, mpfr_t res1, fp_format format)
 945 {
 946   mpfr_init2 (res0, fp_formats[format].mant_dig);
 947   assert_exact (mpfr_set_si (res0, 1, MPFR_RNDN));
 948   mpfr_expm1 (res0, res0, MPFR_RNDU);
 949   mpfr_init2 (res1, fp_formats[format].mant_dig);
 950   assert_exact (mpfr_set_si (res1, 1, MPFR_RNDN));
 951   mpfr_expm1 (res1, res1, MPFR_RNDD);
 952   return 2;
 953 }
 954
 955 /* A special string accepted in input arguments.  */
 956 typedef struct
 957 {
 958   /* The string.  */
 959   const char *str;
 960   /* The function that interprets it for a given floating-point
 961      format, filling in up to two mpfr_t values and returning the
 962      number of values filled.  */
 963   size_t (*func) (mpfr_t, mpfr_t, fp_format);
 964 } special_real_input;
 965
 966 /* List of special strings accepted in input arguments.  */
 967
 968 static const special_real_input special_real_inputs[] =
 969   {
 970     { "max", special_fill_max },
 971     { "-max", special_fill_minus_max },
 972     { "min", special_fill_min },
 973     { "-min", special_fill_minus_min },
 974     { "min_subnorm", special_fill_min_subnorm },
 975     { "-min_subnorm", special_fill_minus_min_subnorm },
 976     { "min_subnorm_p120", special_fill_min_subnorm_p120 },
 977     { "pi", special_fill_pi },
 978     { "-pi", special_fill_minus_pi },
 979     { "pi/2", special_fill_pi_2 },
 980     { "-pi/2", special_fill_minus_pi_2 },
 981     { "pi/4", special_fill_pi_4 },
 982     { "pi/6", special_fill_pi_6 },
 983     { "-pi/6", special_fill_minus_pi_6 },
 984     { "pi/3", special_fill_pi_3 },
 985     { "2pi/3", special_fill_2pi_3 },
 986     { "2pi", special_fill_2pi },
 987     { "e", special_fill_e },
 988     { "1/e", special_fill_1_e },
 989     { "e-1", special_fill_e_minus_1 },
 990   };
 991
 992 /* Given a real number R computed in round-to-zero mode, set the
 993    lowest bit as a sticky bit if INEXACT, and saturate the exponent
 994    range for very large or small values.  */
 995
 996 static void
 997 adjust_real (mpfr_t r, bool inexact)
 998 {
 999   if (!inexact)
1000     return;
1001   /* NaNs are exact, as are infinities in round-to-zero mode.  */
1002   assert (mpfr_number_p (r));
1003   if (mpfr_cmpabs (r, global_min) < 0)
1004     assert_exact (mpfr_copysign (r, global_min, r, MPFR_RNDN));
1005   else if (mpfr_cmpabs (r, global_max) > 0)
1006     assert_exact (mpfr_copysign (r, global_max, r, MPFR_RNDN));
1007   else
1008     {
1009       mpz_t tmp;
1010       mpz_init (tmp);
1011       mpfr_exp_t e = mpfr_get_z_2exp (tmp, r);
1012       if (mpz_sgn (tmp) < 0)
1013         {
1014           mpz_neg (tmp, tmp);
1015           mpz_setbit (tmp, 0);
1016           mpz_neg (tmp, tmp);
1017         }
1018       else
1019         mpz_setbit (tmp, 0);
1020       assert_exact (mpfr_set_z_2exp (r, tmp, e, MPFR_RNDN));
1021       mpz_clear (tmp);
1022     }
1023 }
1024
1025 /* Given a finite real number R with sticky bit, compute the roundings
1026    to FORMAT in each rounding mode, storing the results in RES, the
1027    before-rounding exceptions in EXC_BEFORE and the after-rounding
1028    exceptions in EXC_AFTER.  */
1029
1030 static void
1031 round_real (mpfr_t res[rm_num_modes],
1032             unsigned int exc_before[rm_num_modes],
1033             unsigned int exc_after[rm_num_modes],
1034             mpfr_t r, fp_format format)
1035 {
1036   assert (mpfr_number_p (r));
1037   for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
1038     {
1039       mpfr_init2 (res[m], fp_formats[format].mant_dig);
1040       exc_before[m] = exc_after[m] = 0;
1041       bool inexact = mpfr_set (res[m], r, rounding_modes[m].mpfr_mode);
1042       if (mpfr_cmpabs (res[m], fp_formats[format].max) > 0)
1043         {
1044           inexact = true;
1045           exc_before[m] |= 1U << exc_overflow;
1046           exc_after[m] |= 1U << exc_overflow;
1047           bool overflow_inf;
1048           switch (m)
1049             {
1050             case rm_tonearest:
1051               overflow_inf = true;
1052               break;
1053             case rm_towardzero:
1054               overflow_inf = false;
1055               break;
1056             case rm_downward:
1057               overflow_inf = mpfr_signbit (res[m]);
1058               break;
1059             case rm_upward:
1060               overflow_inf = !mpfr_signbit (res[m]);
1061               break;
1062             default:
1063               abort ();
1064             }
1065           if (overflow_inf)
1066             mpfr_set_inf (res[m], mpfr_signbit (res[m]) ? -1 : 1);
1067           else
1068             assert_exact (mpfr_copysign (res[m], fp_formats[format].max,
1069                                          res[m], MPFR_RNDN));
1070         }
1071       if (mpfr_cmpabs (r, fp_formats[format].min) < 0)
1072         {
1073           /* Tiny before rounding; may or may not be tiny after
1074              rounding, and underflow applies only if also inexact
1075              around rounding to a possibly subnormal value.  */
1076           bool tiny_after_rounding
1077             = mpfr_cmpabs (res[m], fp_formats[format].min) < 0;
1078           /* To round to a possibly subnormal value, and determine
1079              inexactness as a subnormal in the process, scale up and
1080              round to integer, then scale back down.  */
1081           mpfr_t tmp;
1082           mpfr_init (tmp);
1083           assert_exact (mpfr_mul_2si (tmp, r, (fp_formats[format].mant_dig
1084                                                - fp_formats[format].min_exp),
1085                                       MPFR_RNDN));
1086           int rint_res = mpfr_rint (tmp, tmp, rounding_modes[m].mpfr_mode);
1087           /* The integer must be representable.  */
1088           assert (rint_res == 0 || rint_res == 2 || rint_res == -2);
1089           /* If rounding to full precision was inexact, so must
1090              rounding to subnormal precision be inexact.  */
1091           if (inexact)
1092             assert (rint_res != 0);
1093           else
1094             inexact = rint_res != 0;
1095           assert_exact (mpfr_mul_2si (res[m], tmp,
1096                                       (fp_formats[format].min_exp
1097                                        - fp_formats[format].mant_dig),
1098                                       MPFR_RNDN));
1099           mpfr_clear (tmp);
1100           if (inexact)
1101             {
1102               exc_before[m] |= 1U << exc_underflow;
1103               if (tiny_after_rounding)
1104                 exc_after[m] |= 1U << exc_underflow;
1105             }
1106         }
1107       if (inexact)
1108         {
1109           exc_before[m] |= 1U << exc_inexact;
1110           exc_after[m] |= 1U << exc_inexact;
1111         }
1112     }
1113 }
1114
1115 /* Handle the input argument at ARG (NUL-terminated), updating the
1116    lists of test inputs in IT accordingly.  NUM_PREV_ARGS arguments
1117    are already in those lists.  If EXACT_ARGS, interpret a value given
1118    as a floating-point constant exactly (it must be exact for some
1119    supported format) rather than rounding up and down.  The argument,
1120    of type GTYPE, comes from file FILENAME, line LINENO.  */
1121
1122 static void
1123 handle_input_arg (const char *arg, input_test *it, size_t num_prev_args,
1124                   generic_value_type gtype, bool exact_args,
1125                   const char *filename, unsigned int lineno)
1126 {
1127   size_t num_values = 0;
1128   generic_value values[2 * fp_num_formats];
1129   bool check_empty_list = false;
1130   switch (gtype)
1131     {
1132     case gtype_fp:
1133       for (fp_format f = fp_first_format; f < fp_num_formats; f++)
1134         {
1135           mpfr_t extra_values[2];
1136           size_t num_extra_values = 0;
1137           for (size_t i = 0; i < ARRAY_SIZE (special_real_inputs); i++)
1138             {
1139               if (strcmp (arg, special_real_inputs[i].str) == 0)
1140                 {
1141                   num_extra_values
1142                     = special_real_inputs[i].func (extra_values[0],
1143                                                    extra_values[1], f);
1144                   assert (num_extra_values > 0
1145                           && num_extra_values <= ARRAY_SIZE (extra_values));
1146                   break;
1147                 }
1148             }
1149           if (num_extra_values == 0)
1150             {
1151               mpfr_t tmp;
1152               char *ep;
1153               if (exact_args)
1154                 check_empty_list = true;
1155               mpfr_init (tmp);
1156               bool inexact = mpfr_strtofr (tmp, arg, &ep, 0, MPFR_RNDZ);
1157               if (*ep != 0 || !mpfr_number_p (tmp))
1158                 error_at_line (EXIT_FAILURE, 0, filename, lineno,
1159                                "bad floating-point argument: '%s'", arg);
1160               adjust_real (tmp, inexact);
1161               mpfr_t rounded[rm_num_modes];
1162               unsigned int exc_before[rm_num_modes];
1163               unsigned int exc_after[rm_num_modes];
1164               round_real (rounded, exc_before, exc_after, tmp, f);
1165               mpfr_clear (tmp);
1166               if (mpfr_number_p (rounded[rm_upward])
1167                   && (!exact_args || mpfr_equal_p (rounded[rm_upward],
1168                                                    rounded[rm_downward])))
1169                 {
1170                   mpfr_init2 (extra_values[num_extra_values],
1171                               fp_formats[f].mant_dig);
1172                   assert_exact (mpfr_set (extra_values[num_extra_values],
1173                                           rounded[rm_upward], MPFR_RNDN));
1174                   num_extra_values++;
1175                 }
1176               if (mpfr_number_p (rounded[rm_downward]) && !exact_args)
1177                 {
1178                   mpfr_init2 (extra_values[num_extra_values],
1179                               fp_formats[f].mant_dig);
1180                   assert_exact (mpfr_set (extra_values[num_extra_values],
1181                                           rounded[rm_downward], MPFR_RNDN));
1182                   num_extra_values++;
1183                 }
1184               for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
1185                 mpfr_clear (rounded[m]);
1186             }
1187           for (size_t i = 0; i < num_extra_values; i++)
1188             {
1189               bool found = false;
1190               for (size_t j = 0; j < num_values; j++)
1191                 {
1192                   if (mpfr_equal_p (values[j].value.f, extra_values[i])
1193                       && ((mpfr_signbit (values[j].value.f) != 0)
1194                           == (mpfr_signbit (extra_values[i]) != 0)))
1195                     {
1196                       found = true;
1197                       break;
1198                     }
1199                 }
1200               if (!found)
1201                 {
1202                   assert (num_values < ARRAY_SIZE (values));
1203                   values[num_values].type = gtype_fp;
1204                   mpfr_init2 (values[num_values].value.f,
1205                               fp_formats[f].mant_dig);
1206                   assert_exact (mpfr_set (values[num_values].value.f,
1207                                           extra_values[i], MPFR_RNDN));
1208                   num_values++;
1209                 }
1210               mpfr_clear (extra_values[i]);
1211             }
1212         }
1213       break;
1214
1215     case gtype_int:
1216       num_values = 1;
1217       values[0].type = gtype_int;
1218       int ret = mpz_init_set_str (values[0].value.i, arg, 0);
1219       if (ret != 0)
1220         error_at_line (EXIT_FAILURE, 0, filename, lineno,
1221                        "bad integer argument: '%s'", arg);
1222       break;
1223
1224     default:
1225       abort ();
1226     }
1227   if (check_empty_list && num_values == 0)
1228     error_at_line (EXIT_FAILURE, 0, filename, lineno,
1229                    "floating-point argument not exact for any format: '%s'",
1230                    arg);
1231   assert (num_values > 0 && num_values <= ARRAY_SIZE (values));
1232   if (it->num_input_cases >= SIZE_MAX / num_values)
1233     error_at_line (EXIT_FAILURE, 0, filename, lineno, "too many input cases");
1234   generic_value **old_inputs = it->inputs;
1235   size_t new_num_input_cases = it->num_input_cases * num_values;
1236   generic_value **new_inputs = xmalloc (new_num_input_cases
1237                                         * sizeof (new_inputs[0]));
1238   for (size_t i = 0; i < it->num_input_cases; i++)
1239     {
1240       for (size_t j = 0; j < num_values; j++)
1241         {
1242           size_t idx = i * num_values + j;
1243           new_inputs[idx] = xmalloc ((num_prev_args + 1)
1244                                      * sizeof (new_inputs[idx][0]));
1245           for (size_t k = 0; k < num_prev_args; k++)
1246             generic_value_copy (&new_inputs[idx][k], &old_inputs[i][k]);
1247           generic_value_copy (&new_inputs[idx][num_prev_args], &values[j]);
1248         }
1249       for (size_t j = 0; j < num_prev_args; j++)
1250         generic_value_free (&old_inputs[i][j]);
1251       free (old_inputs[i]);
1252     }
1253   free (old_inputs);
1254   for (size_t i = 0; i < num_values; i++)
1255     generic_value_free (&values[i]);
1256   it->inputs = new_inputs;
1257   it->num_input_cases = new_num_input_cases;
1258 }
1259
1260 /* Handle the input flag ARG (NUL-terminated), storing it in *FLAG.
1261    The flag comes from file FILENAME, line LINENO.  */
1262
1263 static void
1264 handle_input_flag (char *arg, input_flag *flag,
1265                    const char *filename, unsigned int lineno)
1266 {
1267   char *ep = strchr (arg, ':');
1268   if (ep == NULL)
1269     {
1270       ep = strchr (arg, 0);
1271       assert (ep != NULL);
1272     }
1273   char c = *ep;
1274   *ep = 0;
1275   bool found = false;
1276   for (input_flag_type i = flag_first_flag; i <= num_input_flag_types; i++)
1277     {
1278       if (strcmp (arg, input_flags[i]) == 0)
1279         {
1280           found = true;
1281           flag->type = i;
1282           break;
1283         }
1284     }
1285   if (!found)
1286     error_at_line (EXIT_FAILURE, 0, filename, lineno, "unknown flag: '%s'",
1287                    arg);
1288   *ep = c;
1289   if (c == 0)
1290     flag->cond = NULL;
1291   else
1292     flag->cond = xstrdup (ep);
1293 }
1294
1295 /* Add the test LINE (file FILENAME, line LINENO) to the test
1296    data.  */
1297
1298 static void
1299 add_test (char *line, const char *filename, unsigned int lineno)
1300 {
1301   size_t num_tokens = 1;
1302   char *p = line;
1303   while ((p = strchr (p, ' ')) != NULL)
1304     {
1305       num_tokens++;
1306       p++;
1307     }
1308   if (num_tokens < 2)
1309     error_at_line (EXIT_FAILURE, 0, filename, lineno,
1310                    "line too short: '%s'", line);
1311   p = strchr (line, ' ');
1312   size_t func_name_len = p - line;
1313   for (size_t i = 0; i < ARRAY_SIZE (test_functions); i++)
1314     {
1315       if (func_name_len == strlen (test_functions[i].name)
1316           && strncmp (line, test_functions[i].name, func_name_len) == 0)
1317         {
1318           test_function *tf = &test_functions[i];
1319           if (num_tokens < 1 + tf->num_args)
1320             error_at_line (EXIT_FAILURE, 0, filename, lineno,
1321                            "line too short: '%s'", line);
1322           if (tf->num_tests == tf->num_tests_alloc)
1323             {
1324               tf->num_tests_alloc = 2 * tf->num_tests_alloc + 16;
1325               tf->tests
1326                 = xrealloc (tf->tests,
1327                             tf->num_tests_alloc * sizeof (tf->tests[0]));
1328             }
1329           input_test *it = &tf->tests[tf->num_tests];
1330           it->line = line;
1331           it->num_input_cases = 1;
1332           it->inputs = xmalloc (sizeof (it->inputs[0]));
1333           it->inputs[0] = NULL;
1334           it->old_output = NULL;
1335           p++;
1336           for (size_t j = 0; j < tf->num_args; j++)
1337             {
1338               char *ep = strchr (p, ' ');
1339               if (ep == NULL)
1340                 {
1341                   ep = strchr (p, '\n');
1342                   assert (ep != NULL);
1343                 }
1344               if (ep == p)
1345                 error_at_line (EXIT_FAILURE, 0, filename, lineno,
1346                                "empty token in line: '%s'", line);
1347               for (char *t = p; t < ep; t++)
1348                 if (isspace ((unsigned char) *t))
1349                   error_at_line (EXIT_FAILURE, 0, filename, lineno,
1350                                  "whitespace in token in line: '%s'", line);
1351               char c = *ep;
1352               *ep = 0;
1353               handle_input_arg (p, it, j,
1354                                 generic_arg_ret_type (tf->arg_types[j]),
1355                                 tf->exact_args, filename, lineno);
1356               *ep = c;
1357               p = ep + 1;
1358             }
1359           it->num_flags = num_tokens - 1 - tf->num_args;
1360           it->flags = xmalloc (it->num_flags * sizeof (it->flags[0]));
1361           for (size_t j = 0; j < it->num_flags; j++)
1362             {
1363               char *ep = strchr (p, ' ');
1364               if (ep == NULL)
1365                 {
1366                   ep = strchr (p, '\n');
1367                   assert (ep != NULL);
1368                 }
1369               if (ep == p)
1370                 error_at_line (EXIT_FAILURE, 0, filename, lineno,
1371                                "empty token in line: '%s'", line);
1372               for (char *t = p; t < ep; t++)
1373                 if (isspace ((unsigned char) *t))
1374                   error_at_line (EXIT_FAILURE, 0, filename, lineno,
1375                                  "whitespace in token in line: '%s'", line);
1376               char c = *ep;
1377               *ep = 0;
1378               handle_input_flag (p, &it->flags[j], filename, lineno);
1379               *ep = c;
1380               p = ep + 1;
1381             }
1382           assert (*p == 0);
1383           tf->num_tests++;
1384           return;
1385         }
1386     }
1387   error_at_line (EXIT_FAILURE, 0, filename, lineno,
1388                  "unknown function in line: '%s'", line);
1389 }
1390
1391 /* Read in the test input data from FILENAME.  */
1392
1393 static void
1394 read_input (const char *filename)
1395 {
1396   FILE *fp = fopen (filename, "r");
1397   if (fp == NULL)
1398     error (EXIT_FAILURE, errno, "open '%s'", filename);
1399   unsigned int lineno = 0;
1400   for (;;)
1401     {
1402       size_t size = 0;
1403       char *line = NULL;
1404       ssize_t ret = getline (&line, &size, fp);
1405       if (ret == -1)
1406         break;
1407       lineno++;
1408       if (line[0] == '#' || line[0] == '\n')
1409         continue;
1410       add_test (line, filename, lineno);
1411     }
1412   if (ferror (fp))
1413     error (EXIT_FAILURE, errno, "read from '%s'", filename);
1414   if (fclose (fp) != 0)
1415     error (EXIT_FAILURE, errno, "close '%s'", filename);
1416 }
1417
1418 /* Calculate the generic results (round-to-zero with sticky bit) for
1419    the function described by CALC, with inputs INPUTS, if MODE is
1420    rm_towardzero; for other modes, calculate results in that mode,
1421    which must be exact zero results.  */
1422
1423 static void
1424 calc_generic_results (generic_value *outputs, generic_value *inputs,
1425                       const func_calc_desc *calc, rounding_mode mode)
1426 {
1427   bool inexact;
1428   int mpc_ternary;
1429   mpc_t ci1, ci2, co;
1430   mpfr_rnd_t mode_mpfr = rounding_modes[mode].mpfr_mode;
1431   mpc_rnd_t mode_mpc = rounding_modes[mode].mpc_mode;
1432
1433   switch (calc->method)
1434     {
1435     case mpfr_f_f:
1436       assert (inputs[0].type == gtype_fp);
1437       outputs[0].type = gtype_fp;
1438       mpfr_init (outputs[0].value.f);
1439       inexact = calc->func.mpfr_f_f (outputs[0].value.f, inputs[0].value.f,
1440                                      mode_mpfr);
1441       if (mode != rm_towardzero)
1442         assert (!inexact && mpfr_zero_p (outputs[0].value.f));
1443       adjust_real (outputs[0].value.f, inexact);
1444       break;
1445
1446     case mpfr_ff_f:
1447       assert (inputs[0].type == gtype_fp);
1448       assert (inputs[1].type == gtype_fp);
1449       outputs[0].type = gtype_fp;
1450       mpfr_init (outputs[0].value.f);
1451       inexact = calc->func.mpfr_ff_f (outputs[0].value.f, inputs[0].value.f,
1452                                       inputs[1].value.f, mode_mpfr);
1453       if (mode != rm_towardzero)
1454         assert (!inexact && mpfr_zero_p (outputs[0].value.f));
1455       adjust_real (outputs[0].value.f, inexact);
1456       break;
1457
1458     case mpfr_fff_f:
1459       assert (inputs[0].type == gtype_fp);
1460       assert (inputs[1].type == gtype_fp);
1461       assert (inputs[2].type == gtype_fp);
1462       outputs[0].type = gtype_fp;
1463       mpfr_init (outputs[0].value.f);
1464       inexact = calc->func.mpfr_fff_f (outputs[0].value.f, inputs[0].value.f,
1465                                        inputs[1].value.f, inputs[2].value.f,
1466                                        mode_mpfr);
1467       if (mode != rm_towardzero)
1468         assert (!inexact && mpfr_zero_p (outputs[0].value.f));
1469       adjust_real (outputs[0].value.f, inexact);
1470       break;
1471
1472     case mpfr_f_f1:
1473       assert (inputs[0].type == gtype_fp);
1474       outputs[0].type = gtype_fp;
1475       outputs[1].type = gtype_int;
1476       mpfr_init (outputs[0].value.f);
1477       int i = 0;
1478       inexact = calc->func.mpfr_f_f1 (outputs[0].value.f, &i,
1479                                       inputs[0].value.f, mode_mpfr);
1480       if (mode != rm_towardzero)
1481         assert (!inexact && mpfr_zero_p (outputs[0].value.f));
1482       adjust_real (outputs[0].value.f, inexact);
1483       mpz_init_set_si (outputs[1].value.i, i);
1484       break;
1485
1486     case mpfr_if_f:
1487       assert (inputs[0].type == gtype_int);
1488       assert (inputs[1].type == gtype_fp);
1489       outputs[0].type = gtype_fp;
1490       mpfr_init (outputs[0].value.f);
1491       assert (mpz_fits_slong_p (inputs[0].value.i));
1492       long l = mpz_get_si (inputs[0].value.i);
1493       inexact = calc->func.mpfr_if_f (outputs[0].value.f, l,
1494                                       inputs[1].value.f, mode_mpfr);
1495       if (mode != rm_towardzero)
1496         assert (!inexact && mpfr_zero_p (outputs[0].value.f));
1497       adjust_real (outputs[0].value.f, inexact);
1498       break;
1499
1500     case mpfr_f_11:
1501       assert (inputs[0].type == gtype_fp);
1502       outputs[0].type = gtype_fp;
1503       mpfr_init (outputs[0].value.f);
1504       outputs[1].type = gtype_fp;
1505       mpfr_init (outputs[1].value.f);
1506       int comb_ternary = calc->func.mpfr_f_11 (outputs[0].value.f,
1507                                                outputs[1].value.f,
1508                                                inputs[0].value.f,
1509                                                mode_mpfr);
1510       if (mode != rm_towardzero)
1511         assert (((comb_ternary & 0x3) == 0
1512                  && mpfr_zero_p (outputs[0].value.f))
1513                 || ((comb_ternary & 0xc) == 0
1514                     && mpfr_zero_p (outputs[1].value.f)));
1515       adjust_real (outputs[0].value.f, (comb_ternary & 0x3) != 0);
1516       adjust_real (outputs[1].value.f, (comb_ternary & 0xc) != 0);
1517       break;
1518
1519     case mpc_c_f:
1520       assert (inputs[0].type == gtype_fp);
1521       assert (inputs[1].type == gtype_fp);
1522       outputs[0].type = gtype_fp;
1523       mpfr_init (outputs[0].value.f);
1524       mpc_init2 (ci1, internal_precision);
1525       assert_exact (mpc_set_fr_fr (ci1, inputs[0].value.f, inputs[1].value.f,
1526                                    MPC_RNDNN));
1527       inexact = calc->func.mpc_c_f (outputs[0].value.f, ci1, mode_mpfr);
1528       if (mode != rm_towardzero)
1529         assert (!inexact && mpfr_zero_p (outputs[0].value.f));
1530       adjust_real (outputs[0].value.f, inexact);
1531       mpc_clear (ci1);
1532       break;
1533
1534     case mpc_c_c:
1535       assert (inputs[0].type == gtype_fp);
1536       assert (inputs[1].type == gtype_fp);
1537       outputs[0].type = gtype_fp;
1538       mpfr_init (outputs[0].value.f);
1539       outputs[1].type = gtype_fp;
1540       mpfr_init (outputs[1].value.f);
1541       mpc_init2 (ci1, internal_precision);
1542       mpc_init2 (co, internal_precision);
1543       assert_exact (mpc_set_fr_fr (ci1, inputs[0].value.f, inputs[1].value.f,
1544                                    MPC_RNDNN));
1545       mpc_ternary = calc->func.mpc_c_c (co, ci1, mode_mpc);
1546       if (mode != rm_towardzero)
1547         assert ((!MPC_INEX_RE (mpc_ternary)
1548                  && mpfr_zero_p (mpc_realref (co)))
1549                 || (!MPC_INEX_IM (mpc_ternary)
1550                     && mpfr_zero_p (mpc_imagref (co))));
1551       assert_exact (mpfr_set (outputs[0].value.f, mpc_realref (co),
1552                               MPFR_RNDN));
1553       assert_exact (mpfr_set (outputs[1].value.f, mpc_imagref (co),
1554                               MPFR_RNDN));
1555       adjust_real (outputs[0].value.f, MPC_INEX_RE (mpc_ternary));
1556       adjust_real (outputs[1].value.f, MPC_INEX_IM (mpc_ternary));
1557       mpc_clear (ci1);
1558       mpc_clear (co);
1559       break;
1560
1561     case mpc_cc_c:
1562       assert (inputs[0].type == gtype_fp);
1563       assert (inputs[1].type == gtype_fp);
1564       assert (inputs[2].type == gtype_fp);
1565       assert (inputs[3].type == gtype_fp);
1566       outputs[0].type = gtype_fp;
1567       mpfr_init (outputs[0].value.f);
1568       outputs[1].type = gtype_fp;
1569       mpfr_init (outputs[1].value.f);
1570       mpc_init2 (ci1, internal_precision);
1571       mpc_init2 (ci2, internal_precision);
1572       mpc_init2 (co, internal_precision);
1573       assert_exact (mpc_set_fr_fr (ci1, inputs[0].value.f, inputs[1].value.f,
1574                                    MPC_RNDNN));
1575       assert_exact (mpc_set_fr_fr (ci2, inputs[2].value.f, inputs[3].value.f,
1576                                    MPC_RNDNN));
1577       mpc_ternary = calc->func.mpc_cc_c (co, ci1, ci2, mode_mpc);
1578       if (mode != rm_towardzero)
1579         assert ((!MPC_INEX_RE (mpc_ternary)
1580                  && mpfr_zero_p (mpc_realref (co)))
1581                 || (!MPC_INEX_IM (mpc_ternary)
1582                     && mpfr_zero_p (mpc_imagref (co))));
1583       assert_exact (mpfr_set (outputs[0].value.f, mpc_realref (co),
1584                               MPFR_RNDN));
1585       assert_exact (mpfr_set (outputs[1].value.f, mpc_imagref (co),
1586                               MPFR_RNDN));
1587       adjust_real (outputs[0].value.f, MPC_INEX_RE (mpc_ternary));
1588       adjust_real (outputs[1].value.f, MPC_INEX_IM (mpc_ternary));
1589       mpc_clear (ci1);
1590       mpc_clear (ci2);
1591       mpc_clear (co);
1592       break;
1593
1594     default:
1595       abort ();
1596     }
1597 }
1598
1599 /* Return the number of bits for integer type TYPE, where "long" has
1600    LONG_BITS bits (32 or 64).  */
1601
1602 static int
1603 int_type_bits (arg_ret_type type, int long_bits)
1604 {
1605   assert (long_bits == 32 || long_bits == 64);
1606   switch (type)
1607     {
1608     case type_int:
1609       return 32;
1610       break;
1611
1612     case type_long:
1613       return long_bits;
1614       break;
1615
1616     case type_long_long:
1617       return 64;
1618       break;
1619
1620     default:
1621       abort ();
1622     }
1623 }
1624
1625 /* Check whether an integer Z fits a given type TYPE, where "long" has
1626    LONG_BITS bits (32 or 64).  */
1627
1628 static bool
1629 int_fits_type (mpz_t z, arg_ret_type type, int long_bits)
1630 {
1631   int bits = int_type_bits (type, long_bits);
1632   bool ret = true;
1633   mpz_t t;
1634   mpz_init (t);
1635   mpz_ui_pow_ui (t, 2, bits - 1);
1636   if (mpz_cmp (z, t) >= 0)
1637     ret = false;
1638   mpz_neg (t, t);
1639   if (mpz_cmp (z, t) < 0)
1640     ret = false;
1641   mpz_clear (t);
1642   return ret;
1643 }
1644
1645 /* Print a generic value V to FP (name FILENAME), preceded by a space,
1646    for type TYPE, floating-point format FORMAT, LONG_BITS bits per
1647    long, printing " IGNORE" instead if IGNORE.  */
1648
1649 static void
1650 output_generic_value (FILE *fp, const char *filename, const generic_value *v,
1651                       bool ignore, arg_ret_type type, fp_format format,
1652                       int long_bits)
1653 {
1654   if (ignore)
1655     {
1656       if (fputs (" IGNORE", fp) < 0)
1657         error (EXIT_FAILURE, errno, "write to '%s'", filename);
1658       return;
1659     }
1660   assert (v->type == generic_arg_ret_type (type));
1661   const char *suffix;
1662   switch (type)
1663     {
1664     case type_fp:
1665       suffix = fp_formats[format].suffix;
1666       break;
1667
1668     case type_int:
1669       suffix = "";
1670       break;
1671
1672     case type_long:
1673       suffix = "L";
1674       break;
1675
1676     case type_long_long:
1677       suffix = "LL";
1678       break;
1679
1680     default:
1681       abort ();
1682     }
1683   switch (v->type)
1684     {
1685     case gtype_fp:
1686       if (mpfr_inf_p (v->value.f))
1687         {
1688           if (fputs ((mpfr_signbit (v->value.f)
1689                       ? " minus_infty" : " plus_infty"), fp) < 0)
1690             error (EXIT_FAILURE, errno, "write to '%s'", filename);
1691         }
1692       else
1693         {
1694           assert (mpfr_number_p (v->value.f));
1695           if (mpfr_fprintf (fp, " %Ra%s", v->value.f, suffix) < 0)
1696             error (EXIT_FAILURE, errno, "mpfr_fprintf to '%s'", filename);
1697         }
1698       break;
1699
1700     case gtype_int: ;
1701       int bits = int_type_bits (type, long_bits);
1702       mpz_t tmp;
1703       mpz_init (tmp);
1704       mpz_ui_pow_ui (tmp, 2, bits - 1);
1705       mpz_neg (tmp, tmp);
1706       if (mpz_cmp (v->value.i, tmp) == 0)
1707         {
1708           mpz_add_ui (tmp, tmp, 1);
1709           if (mpfr_fprintf (fp, " (%Zd%s-1)", tmp, suffix) < 0)
1710             error (EXIT_FAILURE, errno, "mpfr_fprintf to '%s'", filename);
1711         }
1712       else
1713         {
1714           if (mpfr_fprintf (fp, " %Zd%s", v->value.i, suffix) < 0)
1715             error (EXIT_FAILURE, errno, "mpfr_fprintf to '%s'", filename);
1716         }
1717       mpz_clear (tmp);
1718       break;
1719
1720     default:
1721       abort ();
1722     }
1723 }
1724
1725 /* Generate test output to FP (name FILENAME) for test function TF,
1726    input test IT, choice of input values INPUTS.  */
1727
1728 static void
1729 output_for_one_input_case (FILE *fp, const char *filename, test_function *tf,
1730                            input_test *it, generic_value *inputs)
1731 {
1732   bool long_bits_matters = false;
1733   bool fits_long32 = true;
1734   for (size_t i = 0; i < tf->num_args; i++)
1735     {
1736       generic_value_type gtype = generic_arg_ret_type (tf->arg_types[i]);
1737       assert (inputs[i].type == gtype);
1738       if (gtype == gtype_int)
1739         {
1740           bool fits_64 = int_fits_type (inputs[i].value.i, tf->arg_types[i],
1741                                         64);
1742           if (!fits_64)
1743             return;
1744           if (tf->arg_types[i] == type_long
1745               && !int_fits_type (inputs[i].value.i, tf->arg_types[i], 32))
1746             {
1747               long_bits_matters = true;
1748               fits_long32 = false;
1749             }
1750         }
1751     }
1752   generic_value generic_outputs[MAX_NRET];
1753   calc_generic_results (generic_outputs, inputs, &tf->calc, rm_towardzero);
1754   bool ignore_output_long32[MAX_NRET] = { false };
1755   bool ignore_output_long64[MAX_NRET] = { false };
1756   for (size_t i = 0; i < tf->num_ret; i++)
1757     {
1758       assert (generic_outputs[i].type
1759               == generic_arg_ret_type (tf->ret_types[i]));
1760       switch (generic_outputs[i].type)
1761         {
1762         case gtype_fp:
1763           if (!mpfr_number_p (generic_outputs[i].value.f))
1764             goto out; /* Result is NaN or exact infinity.  */
1765           break;
1766
1767         case gtype_int:
1768           ignore_output_long32[i] = !int_fits_type (generic_outputs[i].value.i,
1769                                                     tf->ret_types[i], 32);
1770           ignore_output_long64[i] = !int_fits_type (generic_outputs[i].value.i,
1771                                                     tf->ret_types[i], 64);
1772           if (ignore_output_long32[i] != ignore_output_long64[i])
1773             long_bits_matters = true;
1774           break;
1775
1776         default:
1777           abort ();
1778         }
1779     }
1780   /* Iterate over relevant sizes of long and floating-point formats.  */
1781   for (int long_bits = 32; long_bits <= 64; long_bits += 32)
1782     {
1783       if (long_bits == 32 && !fits_long32)
1784         continue;
1785       if (long_bits == 64 && !long_bits_matters)
1786         continue;
1787       const char *long_cond;
1788       if (long_bits_matters)
1789         long_cond = (long_bits == 32 ? ":long32" : ":long64");
1790       else
1791         long_cond = "";
1792       bool *ignore_output = (long_bits == 32
1793                              ? ignore_output_long32
1794                              : ignore_output_long64);
1795       for (fp_format f = fp_first_format; f < fp_num_formats; f++)
1796         {
1797           bool fits = true;
1798           mpfr_t res[rm_num_modes];
1799           unsigned int exc_before[rm_num_modes];
1800           unsigned int exc_after[rm_num_modes];
1801           for (size_t i = 0; i < tf->num_args; i++)
1802             {
1803               if (inputs[i].type == gtype_fp)
1804                 {
1805                   round_real (res, exc_before, exc_after, inputs[i].value.f,
1806                               f);
1807                   if (!mpfr_equal_p (res[rm_tonearest], inputs[i].value.f))
1808                     fits = false;
1809                   for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
1810                     mpfr_clear (res[m]);
1811                   if (!fits)
1812                     break;
1813                 }
1814             }
1815           if (!fits)
1816             continue;
1817           /* The inputs fit this type, so compute the ideal outputs
1818              and exceptions.  */
1819           mpfr_t all_res[MAX_NRET][rm_num_modes];
1820           unsigned int all_exc_before[MAX_NRET][rm_num_modes];
1821           unsigned int all_exc_after[MAX_NRET][rm_num_modes];
1822           unsigned int merged_exc_before[rm_num_modes] = { 0 };
1823           unsigned int merged_exc_after[rm_num_modes] = { 0 };
1824           /* For functions not exactly determined, track whether
1825              underflow is required (some result is inexact, and
1826              magnitude does not exceed the greatest magnitude
1827              subnormal), and permitted (not an exact zero, and
1828              magnitude does not exceed the least magnitude
1829              normal).  */
1830           bool must_underflow = false;
1831           bool may_underflow = false;
1832           for (size_t i = 0; i < tf->num_ret; i++)
1833             {
1834               switch (generic_outputs[i].type)
1835                 {
1836                 case gtype_fp:
1837                   round_real (all_res[i], all_exc_before[i], all_exc_after[i],
1838                               generic_outputs[i].value.f, f);
1839                   for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
1840                     {
1841                       merged_exc_before[m] |= all_exc_before[i][m];
1842                       merged_exc_after[m] |= all_exc_after[i][m];
1843                       if (!tf->exact)
1844                         {
1845                           must_underflow
1846                             |= ((all_exc_before[i][m]
1847                                  & (1U << exc_inexact)) != 0
1848                                 && (mpfr_cmpabs (generic_outputs[i].value.f,
1849                                                 fp_formats[f].subnorm_max)
1850                                     <= 0));
1851                           may_underflow
1852                             |= (!mpfr_zero_p (generic_outputs[i].value.f)
1853                                 && (mpfr_cmpabs (generic_outputs[i].value.f,
1854                                                  fp_formats[f].min_plus_half)
1855                                     <= 0));
1856                         }
1857                       /* If the result is an exact zero, the sign may
1858                          depend on the rounding mode, so recompute it
1859                          directly in that mode.  */
1860                       if (mpfr_zero_p (all_res[i][m])
1861                           && (all_exc_before[i][m] & (1U << exc_inexact)) == 0)
1862                         {
1863                           generic_value outputs_rm[MAX_NRET];
1864                           calc_generic_results (outputs_rm, inputs,
1865                                                 &tf->calc, m);
1866                           assert_exact (mpfr_set (all_res[i][m],
1867                                                   outputs_rm[i].value.f,
1868                                                   MPFR_RNDN));
1869                           for (size_t j = 0; j < tf->num_ret; j++)
1870                             generic_value_free (&outputs_rm[j]);
1871                         }
1872                     }
1873                   break;
1874
1875                 case gtype_int:
1876                   if (ignore_output[i])
1877                     for (rounding_mode m = rm_first_mode;
1878                          m < rm_num_modes;
1879                          m++)
1880                       {
1881                         merged_exc_before[m] |= 1U << exc_invalid;
1882                         merged_exc_after[m] |= 1U << exc_invalid;
1883                       }
1884                   break;
1885
1886                 default:
1887                   abort ();
1888                 }
1889             }
1890           assert (may_underflow || !must_underflow);
1891           for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
1892             {
1893               bool before_after_matters
1894                 = tf->exact && merged_exc_before[m] != merged_exc_after[m];
1895               if (before_after_matters)
1896                 {
1897                   assert ((merged_exc_before[m] ^ merged_exc_after[m])
1898                           == (1U << exc_underflow));
1899                   assert ((merged_exc_before[m] & (1U << exc_underflow)) != 0);
1900                 }
1901               unsigned int merged_exc = merged_exc_before[m];
1902               if (fprintf (fp, "= %s %s %s%s", tf->name,
1903                            rounding_modes[m].name, fp_formats[f].name,
1904                            long_cond) < 0)
1905                 error (EXIT_FAILURE, errno, "write to '%s'", filename);
1906               /* Print inputs.  */
1907               for (size_t i = 0; i < tf->num_args; i++)
1908                 output_generic_value (fp, filename, &inputs[i], false,
1909                                       tf->arg_types[i], f, long_bits);
1910               if (fputs (" :", fp) < 0)
1911                 error (EXIT_FAILURE, errno, "write to '%s'", filename);
1912               /* Print outputs.  */
1913               bool must_erange = false;
1914               for (size_t i = 0; i < tf->num_ret; i++)
1915                 {
1916                   generic_value g;
1917                   g.type = generic_outputs[i].type;
1918                   switch (g.type)
1919                     {
1920                     case gtype_fp:
1921                       if (mpfr_inf_p (all_res[i][m])
1922                           && (all_exc_before[i][m]
1923                               & (1U << exc_overflow)) != 0)
1924                         must_erange = true;
1925                       if (mpfr_zero_p (all_res[i][m])
1926                           && (tf->exact
1927                               || mpfr_zero_p (all_res[i][rm_tonearest]))
1928                           && (all_exc_before[i][m]
1929                               & (1U << exc_underflow)) != 0)
1930                         must_erange = true;
1931                       mpfr_init2 (g.value.f, fp_formats[f].mant_dig);
1932                       assert_exact (mpfr_set (g.value.f, all_res[i][m],
1933                                               MPFR_RNDN));
1934                       break;
1935
1936                     case gtype_int:
1937                       mpz_init (g.value.i);
1938                       mpz_set (g.value.i, generic_outputs[i].value.i);
1939                       break;
1940
1941                     default:
1942                       abort ();
1943                     }
1944                   output_generic_value (fp, filename, &g, ignore_output[i],
1945                                         tf->ret_types[i], f, long_bits);
1946                   generic_value_free (&g);
1947                 }
1948               if (fputs (" :", fp) < 0)
1949                 error (EXIT_FAILURE, errno, "write to '%s'", filename);
1950               /* Print miscellaneous flags (passed through from
1951                  input).  */
1952               for (size_t i = 0; i < it->num_flags; i++)
1953                 switch (it->flags[i].type)
1954                   {
1955                   case flag_no_test_inline:
1956                   case flag_ignore_zero_inf_sign:
1957                   case flag_xfail:
1958                     if (fprintf (fp, " %s%s",
1959                                  input_flags[it->flags[i].type],
1960                                  (it->flags[i].cond
1961                                   ? it->flags[i].cond
1962                                   : "")) < 0)
1963                       error (EXIT_FAILURE, errno, "write to '%s'",
1964                              filename);
1965                     break;
1966                   case flag_xfail_rounding:
1967                     if (m != rm_tonearest)
1968                       if (fprintf (fp, " xfail%s",
1969                                    (it->flags[i].cond
1970                                     ? it->flags[i].cond
1971                                     : "")) < 0)
1972                         error (EXIT_FAILURE, errno, "write to '%s'",
1973                                filename);
1974                     break;
1975                   default:
1976                     break;
1977                   }
1978               /* Print exception flags and compute errno
1979                  expectations where not already computed.  */
1980               bool may_edom = false;
1981               bool must_edom = false;
1982               bool may_erange = must_erange || may_underflow;
1983               for (fp_exception e = exc_first_exception;
1984                    e < exc_num_exceptions;
1985                    e++)
1986                 {
1987                   bool expect_e = (merged_exc & (1U << e)) != 0;
1988                   bool e_optional = false;
1989                   switch (e)
1990                     {
1991                     case exc_divbyzero:
1992                       if (expect_e)
1993                         may_erange = must_erange = true;
1994                       break;
1995
1996                     case exc_inexact:
1997                       if (!tf->exact)
1998                         e_optional = true;
1999                       break;
2000
2001                     case exc_invalid:
2002                       if (expect_e)
2003                         may_edom = must_edom = true;
2004                       break;
2005
2006                     case exc_overflow:
2007                       if (expect_e)
2008                         may_erange = true;
2009                       break;
2010
2011                     case exc_underflow:
2012                       if (expect_e)
2013                         may_erange = true;
2014                       if (must_underflow)
2015                         assert (expect_e);
2016                       if (may_underflow && !must_underflow)
2017                         e_optional = true;
2018                       break;
2019
2020                     default:
2021                       abort ();
2022                     }
2023                   if (e_optional)
2024                     {
2025                       assert (!before_after_matters);
2026                       if (fprintf (fp, " %s-ok", exceptions[e]) < 0)
2027                         error (EXIT_FAILURE, errno, "write to '%s'",
2028                                filename);
2029                     }
2030                   else
2031                     {
2032                       if (expect_e)
2033                         if (fprintf (fp, " %s", exceptions[e]) < 0)
2034                           error (EXIT_FAILURE, errno, "write to '%s'",
2035                                  filename);
2036                       if (before_after_matters && e == exc_underflow)
2037                         if (fputs (":before-rounding", fp) < 0)
2038                           error (EXIT_FAILURE, errno, "write to '%s'",
2039                                  filename);
2040                       for (int after = 0; after <= 1; after++)
2041                         {
2042                           bool expect_e_here = expect_e;
2043                           if (after == 1 && (!before_after_matters
2044                                              || e != exc_underflow))
2045                             continue;
2046                           const char *after_cond;
2047                           if (before_after_matters && e == exc_underflow)
2048                             {
2049                               after_cond = (after
2050                                             ? ":after-rounding"
2051                                             : ":before-rounding");
2052                               expect_e_here = !after;
2053                             }
2054                           else
2055                             after_cond = "";
2056                           input_flag_type okflag;
2057                           okflag = (expect_e_here
2058                                     ? flag_missing_first
2059                                     : flag_spurious_first) + e;
2060                           for (size_t i = 0; i < it->num_flags; i++)
2061                             if (it->flags[i].type == okflag)
2062                               if (fprintf (fp, " %s-ok%s%s",
2063                                            exceptions[e],
2064                                            (it->flags[i].cond
2065                                             ? it->flags[i].cond
2066                                             : ""), after_cond) < 0)
2067                                 error (EXIT_FAILURE, errno, "write to '%s'",
2068                                        filename);
2069                         }
2070                     }
2071                 }
2072               /* Print errno expectations.  */
2073               if (tf->complex_fn)
2074                 {
2075                   must_edom = false;
2076                   must_erange = false;
2077                 }
2078               if (may_edom && !must_edom)
2079                 {
2080                   if (fputs (" errno-edom-ok", fp) < 0)
2081                     error (EXIT_FAILURE, errno, "write to '%s'",
2082                            filename);
2083                 }
2084               else
2085                 {
2086                   if (must_edom)
2087                     if (fputs (" errno-edom", fp) < 0)
2088                       error (EXIT_FAILURE, errno, "write to '%s'",
2089                              filename);
2090                   input_flag_type okflag = (must_edom
2091                                             ? flag_missing_errno
2092                                             : flag_spurious_errno);
2093                   for (size_t i = 0; i < it->num_flags; i++)
2094                     if (it->flags[i].type == okflag)
2095                       if (fprintf (fp, " errno-edom-ok%s",
2096                                    (it->flags[i].cond
2097                                     ? it->flags[i].cond
2098                                     : "")) < 0)
2099                         error (EXIT_FAILURE, errno, "write to '%s'",
2100                                filename);
2101                 }
2102               if (before_after_matters)
2103                 assert (may_erange && !must_erange);
2104               if (may_erange && !must_erange)
2105                 {
2106                   if (fprintf (fp, " errno-erange-ok%s",
2107                                (before_after_matters
2108                                 ? ":before-rounding"
2109                                 : "")) < 0)
2110                     error (EXIT_FAILURE, errno, "write to '%s'",
2111                            filename);
2112                 }
2113               if (before_after_matters || !(may_erange && !must_erange))
2114                 {
2115                   if (must_erange)
2116                     if (fputs (" errno-erange", fp) < 0)
2117                       error (EXIT_FAILURE, errno, "write to '%s'",
2118                              filename);
2119                   input_flag_type okflag = (must_erange
2120                                             ? flag_missing_errno
2121                                             : flag_spurious_errno);
2122                   for (size_t i = 0; i < it->num_flags; i++)
2123                     if (it->flags[i].type == okflag)
2124                       if (fprintf (fp, " errno-erange-ok%s%s",
2125                                    (it->flags[i].cond
2126                                     ? it->flags[i].cond
2127                                     : ""),
2128                                    (before_after_matters
2129                                     ? ":after-rounding"
2130                                     : "")) < 0)
2131                         error (EXIT_FAILURE, errno, "write to '%s'",
2132                                filename);
2133                 }
2134               if (putc ('\n', fp) < 0)
2135                 error (EXIT_FAILURE, errno, "write to '%s'", filename);
2136             }
2137           for (size_t i = 0; i < tf->num_ret; i++)
2138             {
2139               if (generic_outputs[i].type == gtype_fp)
2140                 for (rounding_mode m = rm_first_mode; m < rm_num_modes; m++)
2141                   mpfr_clear (all_res[i][m]);
2142             }
2143         }
2144     }
2145  out:
2146   for (size_t i = 0; i < tf->num_ret; i++)
2147     generic_value_free (&generic_outputs[i]);
2148 }
2149
2150 /* Generate test output data to FILENAME.  */
2151
2152 static void
2153 generate_output (const char *filename)
2154 {
2155   FILE *fp = fopen (filename, "w");
2156   if (fp == NULL)
2157     error (EXIT_FAILURE, errno, "open '%s'", filename);
2158   for (size_t i = 0; i < ARRAY_SIZE (test_functions); i++)
2159     {
2160       test_function *tf = &test_functions[i];
2161       for (size_t j = 0; j < tf->num_tests; j++)
2162         {
2163           input_test *it = &tf->tests[j];
2164           if (fputs (it->line, fp) < 0)
2165             error (EXIT_FAILURE, errno, "write to '%s'", filename);
2166           for (size_t k = 0; k < it->num_input_cases; k++)
2167             output_for_one_input_case (fp, filename, tf, it, it->inputs[k]);
2168         }
2169     }
2170   if (fclose (fp) != 0)
2171     error (EXIT_FAILURE, errno, "close '%s'", filename);
2172 }
2173
2174 int
2175 main (int argc, char **argv)
2176 {
2177   if (argc != 3)
2178     error (EXIT_FAILURE, 0, "usage: gen-auto-libm-tests <input> <output>");
2179   const char *input_filename = argv[1];
2180   const char *output_filename = argv[2];
2181   init_fp_formats ();
2182   read_input (input_filename);
2183   generate_output (output_filename);
2184   exit (EXIT_SUCCESS);
2185 }