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