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