1	/*
2	* dotest.c - actually generate mathlib test cases
3	*
4	* Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
5	* See https://llvm.org/LICENSE.txt for license information.
6	* SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
7	*/
8
9	#include <stdio.h>
10	#include <string.h>
11	#include <stdlib.h>
12	#include <stdint.h>
13	#include <assert.h>
14	#include <limits.h>
15
16	#include "semi.h"
17	#include "intern.h"
18	#include "random.h"
19
20	#define MPFR_PREC 96 /* good enough for float or double + a few extra bits */
21
22	extern int lib_fo, lib_no_arith, ntests;
23
24	/*
25	* Prototypes.
26	*/
27	static void cases_biased(uint32 *, uint32, uint32);
28	static void cases_biased_positive(uint32 *, uint32, uint32);
29	static void cases_biased_float(uint32 *, uint32, uint32);
30	static void cases_uniform(uint32 *, uint32, uint32);
31	static void cases_uniform_positive(uint32 *, uint32, uint32);
32	static void cases_uniform_float(uint32 *, uint32, uint32);
33	static void cases_uniform_float_positive(uint32 *, uint32, uint32);
34	static void log_cases(uint32 *, uint32, uint32);
35	static void log_cases_float(uint32 *, uint32, uint32);
36	static void log1p_cases(uint32 *, uint32, uint32);
37	static void log1p_cases_float(uint32 *, uint32, uint32);
38	static void minmax_cases(uint32 *, uint32, uint32);
39	static void minmax_cases_float(uint32 *, uint32, uint32);
40	static void atan2_cases(uint32 *, uint32, uint32);
41	static void atan2_cases_float(uint32 *, uint32, uint32);
42	static void pow_cases(uint32 *, uint32, uint32);
43	static void pow_cases_float(uint32 *, uint32, uint32);
44	static void rred_cases(uint32 *, uint32, uint32);
45	static void rred_cases_float(uint32 *, uint32, uint32);
46	static void cases_semi1(uint32 *, uint32, uint32);
47	static void cases_semi1_float(uint32 *, uint32, uint32);
48	static void cases_semi2(uint32 *, uint32, uint32);
49	static void cases_semi2_float(uint32 *, uint32, uint32);
50	static void cases_ldexp(uint32 *, uint32, uint32);
51	static void cases_ldexp_float(uint32 *, uint32, uint32);
52
53	static void complex_cases_uniform(uint32 *, uint32, uint32);
54	static void complex_cases_uniform_float(uint32 *, uint32, uint32);
55	static void complex_cases_biased(uint32 *, uint32, uint32);
56	static void complex_cases_biased_float(uint32 *, uint32, uint32);
57	static void complex_log_cases(uint32 *, uint32, uint32);
58	static void complex_log_cases_float(uint32 *, uint32, uint32);
59	static void complex_pow_cases(uint32 *, uint32, uint32);
60	static void complex_pow_cases_float(uint32 *, uint32, uint32);
61	static void complex_arithmetic_cases(uint32 *, uint32, uint32);
62	static void complex_arithmetic_cases_float(uint32 *, uint32, uint32);
63
64	static uint32 doubletop(int x, int scale);
65	static uint32 floatval(int x, int scale);
66
67	/*
68	* Convert back and forth between IEEE bit patterns and the
69	* mpfr_t/mpc_t types.
70	*/
71	static void set_mpfr_d(mpfr_t x, uint32 h, uint32 l)
72	{
73	uint64_t hl = ((uint64_t)h << 32) | l;
74	uint32 exp = (hl >> 52) & 0x7ff;
75	int64_t mantissa = hl & (((uint64_t)1 << 52) - 1);
76	int sign = (hl >> 63) ? -1 : +1;
77	if (exp == 0x7ff) {
78	if (mantissa == 0)
79	mpfr_set_inf(x, sign);
80	else
81	mpfr_set_nan(x);
82	} else if (exp == 0 && mantissa == 0) {
83	mpfr_set_ui(x, 0, GMP_RNDN);
84	mpfr_setsign(x, x, sign < 0, GMP_RNDN);
85	} else {
86	if (exp != 0)
87	mantissa |= ((uint64_t)1 << 52);
88	else
89	exp++;
90	mpfr_set_sj_2exp(x, mantissa * sign, (int)exp - 0x3ff - 52, GMP_RNDN);
91	}
92	}
93	static void set_mpfr_f(mpfr_t x, uint32 f)
94	{
95	uint32 exp = (f >> 23) & 0xff;
96	int32 mantissa = f & ((1 << 23) - 1);
97	int sign = (f >> 31) ? -1 : +1;
98	if (exp == 0xff) {
99	if (mantissa == 0)
100	mpfr_set_inf(x, sign);
101	else
102	mpfr_set_nan(x);
103	} else if (exp == 0 && mantissa == 0) {
104	mpfr_set_ui(x, 0, GMP_RNDN);
105	mpfr_setsign(x, x, sign < 0, GMP_RNDN);
106	} else {
107	if (exp != 0)
108	mantissa |= (1 << 23);
109	else
110	exp++;
111	mpfr_set_sj_2exp(x, mantissa * sign, (int)exp - 0x7f - 23, GMP_RNDN);
112	}
113	}
114	static void set_mpc_d(mpc_t z, uint32 rh, uint32 rl, uint32 ih, uint32 il)
115	{
116	mpfr_t x, y;
117	mpfr_init2(x, MPFR_PREC);
118	mpfr_init2(y, MPFR_PREC);
119	set_mpfr_d(x, h: rh, l: rl);
120	set_mpfr_d(x: y, h: ih, l: il);
121	mpc_set_fr_fr(z, x, y, MPC_RNDNN);
122	mpfr_clear(x);
123	mpfr_clear(y);
124	}
125	static void set_mpc_f(mpc_t z, uint32 r, uint32 i)
126	{
127	mpfr_t x, y;
128	mpfr_init2(x, MPFR_PREC);
129	mpfr_init2(y, MPFR_PREC);
130	set_mpfr_f(x, f: r);
131	set_mpfr_f(x: y, f: i);
132	mpc_set_fr_fr(z, x, y, MPC_RNDNN);
133	mpfr_clear(x);
134	mpfr_clear(y);
135	}
136	static void get_mpfr_d(const mpfr_t x, uint32 *h, uint32 *l, uint32 *extra)
137	{
138	uint32_t sign, expfield, mantfield;
139	mpfr_t significand;
140	int exp;
141
142	if (mpfr_nan_p(x)) {
143	*h = 0x7ff80000;
144	*l = 0;
145	*extra = 0;
146	return;
147	}
148
149	sign = mpfr_signbit(x) ? 0x80000000U : 0;
150
151	if (mpfr_inf_p(x)) {
152	*h = 0x7ff00000 | sign;
153	*l = 0;
154	*extra = 0;
155	return;
156	}
157
158	if (mpfr_zero_p(x)) {
159	*h = 0x00000000 | sign;
160	*l = 0;
161	*extra = 0;
162	return;
163	}
164
165	mpfr_init2(significand, MPFR_PREC);
166	mpfr_set(significand, x, GMP_RNDN);
167	exp = mpfr_get_exp(significand);
168	mpfr_set_exp(significand, 0);
169
170	/* Now significand is in [1/2,1), and significand * 2^exp == x.
171	* So the IEEE exponent corresponding to exp==0 is 0x3fe. */
172	if (exp > 0x400) {
173	/* overflow to infinity anyway */
174	*h = 0x7ff00000 | sign;
175	*l = 0;
176	*extra = 0;
177	mpfr_clear(significand);
178	return;
179	}
180
181	if (exp <= -0x3fe || mpfr_zero_p(x))
182	exp = -0x3fd; /* denormalise */
183	expfield = exp + 0x3fd; /* offset to cancel leading mantissa bit */
184
185	mpfr_div_2si(significand, x, exp - 21, GMP_RNDN);
186	mpfr_abs(significand, significand, GMP_RNDN);
187	mantfield = mpfr_get_ui(significand, GMP_RNDZ);
188	*h = sign + ((uint64_t)expfield << 20) + mantfield;
189	mpfr_sub_ui(significand, significand, mantfield, GMP_RNDN);
190	mpfr_mul_2ui(significand, significand, 32, GMP_RNDN);
191	mantfield = mpfr_get_ui(significand, GMP_RNDZ);
192	*l = mantfield;
193	mpfr_sub_ui(significand, significand, mantfield, GMP_RNDN);
194	mpfr_mul_2ui(significand, significand, 32, GMP_RNDN);
195	mantfield = mpfr_get_ui(significand, GMP_RNDZ);
196	*extra = mantfield;
197
198	mpfr_clear(significand);
199	}
200	static void get_mpfr_f(const mpfr_t x, uint32 *f, uint32 *extra)
201	{
202	uint32_t sign, expfield, mantfield;
203	mpfr_t significand;
204	int exp;
205
206	if (mpfr_nan_p(x)) {
207	*f = 0x7fc00000;
208	*extra = 0;
209	return;
210	}
211
212	sign = mpfr_signbit(x) ? 0x80000000U : 0;
213
214	if (mpfr_inf_p(x)) {
215	*f = 0x7f800000 | sign;
216	*extra = 0;
217	return;
218	}
219
220	if (mpfr_zero_p(x)) {
221	*f = 0x00000000 | sign;
222	*extra = 0;
223	return;
224	}
225
226	mpfr_init2(significand, MPFR_PREC);
227	mpfr_set(significand, x, GMP_RNDN);
228	exp = mpfr_get_exp(significand);
229	mpfr_set_exp(significand, 0);
230
231	/* Now significand is in [1/2,1), and significand * 2^exp == x.
232	* So the IEEE exponent corresponding to exp==0 is 0x7e. */
233	if (exp > 0x80) {
234	/* overflow to infinity anyway */
235	*f = 0x7f800000 | sign;
236	*extra = 0;
237	mpfr_clear(significand);
238	return;
239	}
240
241	if (exp <= -0x7e || mpfr_zero_p(x))
242	exp = -0x7d; /* denormalise */
243	expfield = exp + 0x7d; /* offset to cancel leading mantissa bit */
244
245	mpfr_div_2si(significand, x, exp - 24, GMP_RNDN);
246	mpfr_abs(significand, significand, GMP_RNDN);
247	mantfield = mpfr_get_ui(significand, GMP_RNDZ);
248	*f = sign + ((uint64_t)expfield << 23) + mantfield;
249	mpfr_sub_ui(significand, significand, mantfield, GMP_RNDN);
250	mpfr_mul_2ui(significand, significand, 32, GMP_RNDN);
251	mantfield = mpfr_get_ui(significand, GMP_RNDZ);
252	*extra = mantfield;
253
254	mpfr_clear(significand);
255	}
256	static void get_mpc_d(const mpc_t z,
257	uint32 *rh, uint32 *rl, uint32 *rextra,
258	uint32 *ih, uint32 *il, uint32 *iextra)
259	{
260	mpfr_t x, y;
261	mpfr_init2(x, MPFR_PREC);
262	mpfr_init2(y, MPFR_PREC);
263	mpc_real(x, z, GMP_RNDN);
264	mpc_imag(y, z, GMP_RNDN);
265	get_mpfr_d(x, h: rh, l: rl, extra: rextra);
266	get_mpfr_d(x: y, h: ih, l: il, extra: iextra);
267	mpfr_clear(x);
268	mpfr_clear(y);
269	}
270	static void get_mpc_f(const mpc_t z,
271	uint32 *r, uint32 *rextra,
272	uint32 *i, uint32 *iextra)
273	{
274	mpfr_t x, y;
275	mpfr_init2(x, MPFR_PREC);
276	mpfr_init2(y, MPFR_PREC);
277	mpc_real(x, z, GMP_RNDN);
278	mpc_imag(y, z, GMP_RNDN);
279	get_mpfr_f(x, f: r, extra: rextra);
280	get_mpfr_f(x: y, f: i, extra: iextra);
281	mpfr_clear(x);
282	mpfr_clear(y);
283	}
284
285	/*
286	* Implementation of mathlib functions that aren't trivially
287	* implementable using an existing mpfr or mpc function.
288	*/
289	int test_rred(mpfr_t ret, const mpfr_t x, int *quadrant)
290	{
291	mpfr_t halfpi;
292	long quo;
293	int status;
294
295	/*
296	* In the worst case of range reduction, we get an input of size
297	* around 2^1024, and must find its remainder mod pi, which means
298	* we need 1024 bits of pi at least. Plus, the remainder might
299	* happen to come out very very small if we're unlucky. How
300	* unlucky can we be? Well, conveniently, I once went through and
301	* actually worked that out using Paxson's modular minimisation
302	* algorithm, and it turns out that the smallest exponent you can
303	* get out of a nontrivial[1] double precision range reduction is
304	* 0x3c2, i.e. of the order of 2^-61. So we need 1024 bits of pi
305	* to get us down to the units digit, another 61 or so bits (say
306	* 64) to get down to the highest set bit of the output, and then
307	* some bits to make the actual mantissa big enough.
308	*
309	* [1] of course the output of range reduction can have an
310	* arbitrarily small exponent in the trivial case, where the
311	* input is so small that it's the identity function. That
312	* doesn't count.
313	*/
314	mpfr_init2(halfpi, MPFR_PREC + 1024 + 64);
315	mpfr_const_pi(halfpi, GMP_RNDN);
316	mpfr_div_ui(halfpi, halfpi, 2, GMP_RNDN);
317
318	status = mpfr_remquo(ret, &quo, x, halfpi, GMP_RNDN);
319	*quadrant = quo & 3;
320
321	mpfr_clear(halfpi);
322
323	return status;
324	}
325	int test_lgamma(mpfr_t ret, const mpfr_t x, mpfr_rnd_t rnd)
326	{
327	/*
328	* mpfr_lgamma takes an extra int * parameter to hold the output
329	* sign. We don't bother testing that, so this wrapper throws away
330	* the sign and hence fits into the same function prototype as all
331	* the other real->real mpfr functions.
332	*
333	* There is also mpfr_lngamma which has no sign output and hence
334	* has the right prototype already, but unfortunately it returns
335	* NaN in cases where gamma(x) < 0, so it's no use to us.
336	*/
337	int sign;
338	return mpfr_lgamma(ret, &sign, x, rnd);
339	}
340	int test_cpow(mpc_t ret, const mpc_t x, const mpc_t y, mpc_rnd_t rnd)
341	{
342	/*
343	* For complex pow, we must bump up the precision by a huge amount
344	* if we want it to get the really difficult cases right. (Not
345	* that we expect the library under test to be getting those cases
346	* right itself, but we'd at least like the test suite to report
347	* them as wrong for the _right reason_.)
348	*
349	* This works around a bug in mpc_pow(), fixed by r1455 in the MPC
350	* svn repository (2014-10-14) and expected to be in any MPC
351	* release after 1.0.2 (which was the latest release already made
352	* at the time of the fix). So as and when we update to an MPC
353	* with the fix in it, we could remove this workaround.
354	*
355	* For the reasons for choosing this amount of extra precision,
356	* see analysis in complex/cpownotes.txt for the rationale for the
357	* amount.
358	*/
359	mpc_t xbig, ybig, retbig;
360	int status;
361
362	mpc_init2(xbig, 1034 + 53 + 60 + MPFR_PREC);
363	mpc_init2(ybig, 1034 + 53 + 60 + MPFR_PREC);
364	mpc_init2(retbig, 1034 + 53 + 60 + MPFR_PREC);
365
366	mpc_set(xbig, x, MPC_RNDNN);
367	mpc_set(ybig, y, MPC_RNDNN);
368	status = mpc_pow(retbig, xbig, ybig, rnd);
369	mpc_set(ret, retbig, rnd);
370
371	mpc_clear(xbig);
372	mpc_clear(ybig);
373	mpc_clear(retbig);
374
375	return status;
376	}
377
378	/*
379	* Identify 'hard' values (NaN, Inf, nonzero denormal) for deciding
380	* whether microlib will decline to run a test.
381	*/
382	#define is_shard(in) ( \
383	(((in)[0] & 0x7F800000) == 0x7F800000 || \
384	(((in)[0] & 0x7F800000) == 0 && ((in)[0]&0x7FFFFFFF) != 0)))
385
386	#define is_dhard(in) ( \
387	(((in)[0] & 0x7FF00000) == 0x7FF00000 || \
388	(((in)[0] & 0x7FF00000) == 0 && (((in)[0] & 0xFFFFF) | (in)[1]) != 0)))
389
390	/*
391	* Identify integers.
392	*/
393	int is_dinteger(uint32 *in)
394	{
395	uint32 out[3];
396	if ((0x7FF00000 & ~in[0]) == 0)
397	return 0; /* not finite, hence not integer */
398	test_ceil(in, out);
399	return in[0] == out[0] && in[1] == out[1];
400	}
401	int is_sinteger(uint32 *in)
402	{
403	uint32 out[3];
404	if ((0x7F800000 & ~in[0]) == 0)
405	return 0; /* not finite, hence not integer */
406	test_ceilf(in, out);
407	return in[0] == out[0];
408	}
409
410	/*
411	* Identify signalling NaNs.
412	*/
413	int is_dsnan(const uint32 *in)
414	{
415	if ((in[0] & 0x7FF00000) != 0x7FF00000)
416	return 0; /* not the inf/nan exponent */
417	if ((in[0] << 12) == 0 && in[1] == 0)
418	return 0; /* inf */
419	if (in[0] & 0x00080000)
420	return 0; /* qnan */
421	return 1;
422	}
423	int is_ssnan(const uint32 *in)
424	{
425	if ((in[0] & 0x7F800000) != 0x7F800000)
426	return 0; /* not the inf/nan exponent */
427	if ((in[0] << 9) == 0)
428	return 0; /* inf */
429	if (in[0] & 0x00400000)
430	return 0; /* qnan */
431	return 1;
432	}
433	int is_snan(const uint32 *in, int size)
434	{
435	return size == 2 ? is_dsnan(in) : is_ssnan(in);
436	}
437
438	/*
439	* Wrapper functions called to fix up unusual results after the main
440	* test function has run.
441	*/
442	void universal_wrapper(wrapperctx *ctx)
443	{
444	/*
445	* Any SNaN input gives rise to a QNaN output.
446	*/
447	int op;
448	for (op = 0; op < wrapper_get_nops(ctx); op++) {
449	int size = wrapper_get_size(ctx, op);
450
451	if (!wrapper_is_complex(ctx, op) &&
452	is_snan(in: wrapper_get_ieee(ctx, op), size)) {
453	wrapper_set_nan(ctx);
454	}
455	}
456	}
457
458	Testable functions[] = {
459	/*
460	* Trig functions: sin, cos, tan. We test the core function
461	* between -16 and +16: we assume that range reduction exists
462	* and will be used for larger arguments, and we'll test that
463	* separately. Also we only go down to 2^-27 in magnitude,
464	* because below that sin(x)=tan(x)=x and cos(x)=1 as far as
465	* double precision can tell, which is boring.
466	*/
467	{"sin", (funcptr)mpfr_sin, args1, {NULL},
468	cases_uniform, 0x3e400000, 0x40300000},
469	{"sinf", (funcptr)mpfr_sin, args1f, {NULL},
470	cases_uniform_float, 0x39800000, 0x41800000},
471	{"cos", (funcptr)mpfr_cos, args1, {NULL},
472	cases_uniform, 0x3e400000, 0x40300000},
473	{"cosf", (funcptr)mpfr_cos, args1f, {NULL},
474	cases_uniform_float, 0x39800000, 0x41800000},
475	{"tan", (funcptr)mpfr_tan, args1, {NULL},
476	cases_uniform, 0x3e400000, 0x40300000},
477	{"tanf", (funcptr)mpfr_tan, args1f, {NULL},
478	cases_uniform_float, 0x39800000, 0x41800000},
479	{"sincosf_sinf", (funcptr)mpfr_sin, args1f, {NULL},
480	cases_uniform_float, 0x39800000, 0x41800000},
481	{"sincosf_cosf", (funcptr)mpfr_cos, args1f, {NULL},
482	cases_uniform_float, 0x39800000, 0x41800000},
483	/*
484	* Inverse trig: asin, acos. Between 1 and -1, of course. acos
485	* goes down to 2^-54, asin to 2^-27.
486	*/
487	{"asin", (funcptr)mpfr_asin, args1, {NULL},
488	cases_uniform, 0x3e400000, 0x3fefffff},
489	{"asinf", (funcptr)mpfr_asin, args1f, {NULL},
490	cases_uniform_float, 0x39800000, 0x3f7fffff},
491	{"acos", (funcptr)mpfr_acos, args1, {NULL},
492	cases_uniform, 0x3c900000, 0x3fefffff},
493	{"acosf", (funcptr)mpfr_acos, args1f, {NULL},
494	cases_uniform_float, 0x33800000, 0x3f7fffff},
495	/*
496	* Inverse trig: atan. atan is stable (in double prec) with
497	* argument magnitude past 2^53, so we'll test up to there.
498	* atan(x) is boringly just x below 2^-27.
499	*/
500	{"atan", (funcptr)mpfr_atan, args1, {NULL},
501	cases_uniform, 0x3e400000, 0x43400000},
502	{"atanf", (funcptr)mpfr_atan, args1f, {NULL},
503	cases_uniform_float, 0x39800000, 0x4b800000},
504	/*
505	* atan2. Interesting cases arise when the exponents of the
506	* arguments differ by at most about 50.
507	*/
508	{"atan2", (funcptr)mpfr_atan2, args2, {NULL},
509	atan2_cases, 0},
510	{"atan2f", (funcptr)mpfr_atan2, args2f, {NULL},
511	atan2_cases_float, 0},
512	/*
513	* The exponentials: exp, sinh, cosh. They overflow at around
514	* 710. exp and sinh are boring below 2^-54, cosh below 2^-27.
515	*/
516	{"exp", (funcptr)mpfr_exp, args1, {NULL},
517	cases_uniform, 0x3c900000, 0x40878000},
518	{"expf", (funcptr)mpfr_exp, args1f, {NULL},
519	cases_uniform_float, 0x33800000, 0x42dc0000},
520	{"sinh", (funcptr)mpfr_sinh, args1, {NULL},
521	cases_uniform, 0x3c900000, 0x40878000},
522	{"sinhf", (funcptr)mpfr_sinh, args1f, {NULL},
523	cases_uniform_float, 0x33800000, 0x42dc0000},
524	{"cosh", (funcptr)mpfr_cosh, args1, {NULL},
525	cases_uniform, 0x3e400000, 0x40878000},
526	{"coshf", (funcptr)mpfr_cosh, args1f, {NULL},
527	cases_uniform_float, 0x39800000, 0x42dc0000},
528	/*
529	* tanh is stable past around 20. It's boring below 2^-27.
530	*/
531	{"tanh", (funcptr)mpfr_tanh, args1, {NULL},
532	cases_uniform, 0x3e400000, 0x40340000},
533	{"tanhf", (funcptr)mpfr_tanh, args1f, {NULL},
534	cases_uniform, 0x39800000, 0x41100000},
535	/*
536	* log must be tested only on positive numbers, but can cover
537	* the whole range of positive nonzero finite numbers. It never
538	* gets boring.
539	*/
540	{"log", (funcptr)mpfr_log, args1, {NULL}, log_cases, 0},
541	{"logf", (funcptr)mpfr_log, args1f, {NULL}, log_cases_float, 0},
542	{"log10", (funcptr)mpfr_log10, args1, {NULL}, log_cases, 0},
543	{"log10f", (funcptr)mpfr_log10, args1f, {NULL}, log_cases_float, 0},
544	/*
545	* pow.
546	*/
547	{"pow", (funcptr)mpfr_pow, args2, {NULL}, pow_cases, 0},
548	{"powf", (funcptr)mpfr_pow, args2f, {NULL}, pow_cases_float, 0},
549	/*
550	* Trig range reduction. We are able to test this for all
551	* finite values, but will only bother for things between 2^-3
552	* and 2^+52.
553	*/
554	{"rred", (funcptr)test_rred, rred, {NULL}, rred_cases, 0},
555	{"rredf", (funcptr)test_rred, rredf, {NULL}, rred_cases_float, 0},
556	/*
557	* Square and cube root.
558	*/
559	{"sqrt", (funcptr)mpfr_sqrt, args1, {NULL}, log_cases, 0},
560	{"sqrtf", (funcptr)mpfr_sqrt, args1f, {NULL}, log_cases_float, 0},
561	{"cbrt", (funcptr)mpfr_cbrt, args1, {NULL}, log_cases, 0},
562	{"cbrtf", (funcptr)mpfr_cbrt, args1f, {NULL}, log_cases_float, 0},
563	{"hypot", (funcptr)mpfr_hypot, args2, {NULL}, atan2_cases, 0},
564	{"hypotf", (funcptr)mpfr_hypot, args2f, {NULL}, atan2_cases_float, 0},
565	/*
566	* Seminumerical functions.
567	*/
568	{"ceil", (funcptr)test_ceil, semi1, {NULL}, cases_semi1},
569	{"ceilf", (funcptr)test_ceilf, semi1f, {NULL}, cases_semi1_float},
570	{"floor", (funcptr)test_floor, semi1, {NULL}, cases_semi1},
571	{"floorf", (funcptr)test_floorf, semi1f, {NULL}, cases_semi1_float},
572	{"fmod", (funcptr)test_fmod, semi2, {NULL}, cases_semi2},
573	{"fmodf", (funcptr)test_fmodf, semi2f, {NULL}, cases_semi2_float},
574	{"ldexp", (funcptr)test_ldexp, t_ldexp, {NULL}, cases_ldexp},
575	{"ldexpf", (funcptr)test_ldexpf, t_ldexpf, {NULL}, cases_ldexp_float},
576	{"frexp", (funcptr)test_frexp, t_frexp, {NULL}, cases_semi1},
577	{"frexpf", (funcptr)test_frexpf, t_frexpf, {NULL}, cases_semi1_float},
578	{"modf", (funcptr)test_modf, t_modf, {NULL}, cases_semi1},
579	{"modff", (funcptr)test_modff, t_modff, {NULL}, cases_semi1_float},
580
581	/*
582	* Classification and more semi-numericals
583	*/
584	{"copysign", (funcptr)test_copysign, semi2, {NULL}, cases_semi2},
585	{"copysignf", (funcptr)test_copysignf, semi2f, {NULL}, cases_semi2_float},
586	{"isfinite", (funcptr)test_isfinite, classify, {NULL}, cases_uniform, 0, 0x7fffffff},
587	{"isfinitef", (funcptr)test_isfinitef, classifyf, {NULL}, cases_uniform_float, 0, 0x7fffffff},
588	{"isinf", (funcptr)test_isinf, classify, {NULL}, cases_uniform, 0, 0x7fffffff},
589	{"isinff", (funcptr)test_isinff, classifyf, {NULL}, cases_uniform_float, 0, 0x7fffffff},
590	{"isnan", (funcptr)test_isnan, classify, {NULL}, cases_uniform, 0, 0x7fffffff},
591	{"isnanf", (funcptr)test_isnanf, classifyf, {NULL}, cases_uniform_float, 0, 0x7fffffff},
592	{"isnormal", (funcptr)test_isnormal, classify, {NULL}, cases_uniform, 0, 0x7fffffff},
593	{"isnormalf", (funcptr)test_isnormalf, classifyf, {NULL}, cases_uniform_float, 0, 0x7fffffff},
594	{"signbit", (funcptr)test_signbit, classify, {NULL}, cases_uniform, 0, 0x7fffffff},
595	{"signbitf", (funcptr)test_signbitf, classifyf, {NULL}, cases_uniform_float, 0, 0x7fffffff},
596	{"fpclassify", (funcptr)test_fpclassify, classify, {NULL}, cases_uniform, 0, 0x7fffffff},
597	{"fpclassifyf", (funcptr)test_fpclassifyf, classifyf, {NULL}, cases_uniform_float, 0, 0x7fffffff},
598	/*
599	* Comparisons
600	*/
601	{"isgreater", (funcptr)test_isgreater, compare, {NULL}, cases_uniform, 0, 0x7fffffff},
602	{"isgreaterequal", (funcptr)test_isgreaterequal, compare, {NULL}, cases_uniform, 0, 0x7fffffff},
603	{"isless", (funcptr)test_isless, compare, {NULL}, cases_uniform, 0, 0x7fffffff},
604	{"islessequal", (funcptr)test_islessequal, compare, {NULL}, cases_uniform, 0, 0x7fffffff},
605	{"islessgreater", (funcptr)test_islessgreater, compare, {NULL}, cases_uniform, 0, 0x7fffffff},
606	{"isunordered", (funcptr)test_isunordered, compare, {NULL}, cases_uniform, 0, 0x7fffffff},
607
608	{"isgreaterf", (funcptr)test_isgreaterf, comparef, {NULL}, cases_uniform_float, 0, 0x7fffffff},
609	{"isgreaterequalf", (funcptr)test_isgreaterequalf, comparef, {NULL}, cases_uniform_float, 0, 0x7fffffff},
610	{"islessf", (funcptr)test_islessf, comparef, {NULL}, cases_uniform_float, 0, 0x7fffffff},
611	{"islessequalf", (funcptr)test_islessequalf, comparef, {NULL}, cases_uniform_float, 0, 0x7fffffff},
612	{"islessgreaterf", (funcptr)test_islessgreaterf, comparef, {NULL}, cases_uniform_float, 0, 0x7fffffff},
613	{"isunorderedf", (funcptr)test_isunorderedf, comparef, {NULL}, cases_uniform_float, 0, 0x7fffffff},
614
615	/*
616	* Inverse Hyperbolic functions
617	*/
618	{"atanh", (funcptr)mpfr_atanh, args1, {NULL}, cases_uniform, 0x3e400000, 0x3fefffff},
619	{"asinh", (funcptr)mpfr_asinh, args1, {NULL}, cases_uniform, 0x3e400000, 0x3fefffff},
620	{"acosh", (funcptr)mpfr_acosh, args1, {NULL}, cases_uniform_positive, 0x3ff00000, 0x7fefffff},
621
622	{"atanhf", (funcptr)mpfr_atanh, args1f, {NULL}, cases_uniform_float, 0x32000000, 0x3f7fffff},
623	{"asinhf", (funcptr)mpfr_asinh, args1f, {NULL}, cases_uniform_float, 0x32000000, 0x3f7fffff},
624	{"acoshf", (funcptr)mpfr_acosh, args1f, {NULL}, cases_uniform_float_positive, 0x3f800000, 0x7f800000},
625
626	/*
627	* Everything else (sitting in a section down here at the bottom
628	* because historically they were not tested because we didn't
629	* have reference implementations for them)
630	*/
631	{"csin", (funcptr)mpc_sin, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000},
632	{"csinf", (funcptr)mpc_sin, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000},
633	{"ccos", (funcptr)mpc_cos, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000},
634	{"ccosf", (funcptr)mpc_cos, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000},
635	{"ctan", (funcptr)mpc_tan, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000},
636	{"ctanf", (funcptr)mpc_tan, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000},
637
638	{"casin", (funcptr)mpc_asin, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000},
639	{"casinf", (funcptr)mpc_asin, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000},
640	{"cacos", (funcptr)mpc_acos, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000},
641	{"cacosf", (funcptr)mpc_acos, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000},
642	{"catan", (funcptr)mpc_atan, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000},
643	{"catanf", (funcptr)mpc_atan, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000},
644
645	{"csinh", (funcptr)mpc_sinh, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000},
646	{"csinhf", (funcptr)mpc_sinh, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000},
647	{"ccosh", (funcptr)mpc_cosh, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000},
648	{"ccoshf", (funcptr)mpc_cosh, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000},
649	{"ctanh", (funcptr)mpc_tanh, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000},
650	{"ctanhf", (funcptr)mpc_tanh, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000},
651
652	{"casinh", (funcptr)mpc_asinh, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000},
653	{"casinhf", (funcptr)mpc_asinh, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000},
654	{"cacosh", (funcptr)mpc_acosh, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000},
655	{"cacoshf", (funcptr)mpc_acosh, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000},
656	{"catanh", (funcptr)mpc_atanh, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000},
657	{"catanhf", (funcptr)mpc_atanh, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000},
658
659	{"cexp", (funcptr)mpc_exp, args1c, {NULL}, complex_cases_uniform, 0x3c900000, 0x40862000},
660	{"cpow", (funcptr)test_cpow, args2c, {NULL}, complex_pow_cases, 0x3fc00000, 0x40000000},
661	{"clog", (funcptr)mpc_log, args1c, {NULL}, complex_log_cases, 0, 0},
662	{"csqrt", (funcptr)mpc_sqrt, args1c, {NULL}, complex_log_cases, 0, 0},
663
664	{"cexpf", (funcptr)mpc_exp, args1fc, {NULL}, complex_cases_uniform_float, 0x24800000, 0x42b00000},
665	{"cpowf", (funcptr)test_cpow, args2fc, {NULL}, complex_pow_cases_float, 0x3e000000, 0x41000000},
666	{"clogf", (funcptr)mpc_log, args1fc, {NULL}, complex_log_cases_float, 0, 0},
667	{"csqrtf", (funcptr)mpc_sqrt, args1fc, {NULL}, complex_log_cases_float, 0, 0},
668
669	{"cdiv", (funcptr)mpc_div, args2c, {NULL}, complex_arithmetic_cases, 0, 0},
670	{"cmul", (funcptr)mpc_mul, args2c, {NULL}, complex_arithmetic_cases, 0, 0},
671	{"cadd", (funcptr)mpc_add, args2c, {NULL}, complex_arithmetic_cases, 0, 0},
672	{"csub", (funcptr)mpc_sub, args2c, {NULL}, complex_arithmetic_cases, 0, 0},
673
674	{"cdivf", (funcptr)mpc_div, args2fc, {NULL}, complex_arithmetic_cases_float, 0, 0},
675	{"cmulf", (funcptr)mpc_mul, args2fc, {NULL}, complex_arithmetic_cases_float, 0, 0},
676	{"caddf", (funcptr)mpc_add, args2fc, {NULL}, complex_arithmetic_cases_float, 0, 0},
677	{"csubf", (funcptr)mpc_sub, args2fc, {NULL}, complex_arithmetic_cases_float, 0, 0},
678
679	{"cabsf", (funcptr)mpc_abs, args1fcr, {NULL}, complex_arithmetic_cases_float, 0, 0},
680	{"cabs", (funcptr)mpc_abs, args1cr, {NULL}, complex_arithmetic_cases, 0, 0},
681	{"cargf", (funcptr)mpc_arg, args1fcr, {NULL}, complex_arithmetic_cases_float, 0, 0},
682	{"carg", (funcptr)mpc_arg, args1cr, {NULL}, complex_arithmetic_cases, 0, 0},
683	{"cimagf", (funcptr)mpc_imag, args1fcr, {NULL}, complex_arithmetic_cases_float, 0, 0},
684	{"cimag", (funcptr)mpc_imag, args1cr, {NULL}, complex_arithmetic_cases, 0, 0},
685	{"conjf", (funcptr)mpc_conj, args1fc, {NULL}, complex_arithmetic_cases_float, 0, 0},
686	{"conj", (funcptr)mpc_conj, args1c, {NULL}, complex_arithmetic_cases, 0, 0},
687	{"cprojf", (funcptr)mpc_proj, args1fc, {NULL}, complex_arithmetic_cases_float, 0, 0},
688	{"cproj", (funcptr)mpc_proj, args1c, {NULL}, complex_arithmetic_cases, 0, 0},
689	{"crealf", (funcptr)mpc_real, args1fcr, {NULL}, complex_arithmetic_cases_float, 0, 0},
690	{"creal", (funcptr)mpc_real, args1cr, {NULL}, complex_arithmetic_cases, 0, 0},
691	{"erfcf", (funcptr)mpfr_erfc, args1f, {NULL}, cases_biased_float, 0x1e800000, 0x41000000},
692	{"erfc", (funcptr)mpfr_erfc, args1, {NULL}, cases_biased, 0x3bd00000, 0x403c0000},
693	{"erff", (funcptr)mpfr_erf, args1f, {NULL}, cases_biased_float, 0x03800000, 0x40700000},
694	{"erf", (funcptr)mpfr_erf, args1, {NULL}, cases_biased, 0x00800000, 0x40200000},
695	{"exp2f", (funcptr)mpfr_exp2, args1f, {NULL}, cases_uniform_float, 0x33800000, 0x43c00000},
696	{"exp2", (funcptr)mpfr_exp2, args1, {NULL}, cases_uniform, 0x3ca00000, 0x40a00000},
697	{"expm1f", (funcptr)mpfr_expm1, args1f, {NULL}, cases_uniform_float, 0x33000000, 0x43800000},
698	{"expm1", (funcptr)mpfr_expm1, args1, {NULL}, cases_uniform, 0x3c900000, 0x409c0000},
699	{"fmaxf", (funcptr)mpfr_max, args2f, {NULL}, minmax_cases_float, 0, 0x7f7fffff},
700	{"fmax", (funcptr)mpfr_max, args2, {NULL}, minmax_cases, 0, 0x7fefffff},
701	{"fminf", (funcptr)mpfr_min, args2f, {NULL}, minmax_cases_float, 0, 0x7f7fffff},
702	{"fmin", (funcptr)mpfr_min, args2, {NULL}, minmax_cases, 0, 0x7fefffff},
703	{"lgammaf", (funcptr)test_lgamma, args1f, {NULL}, cases_uniform_float, 0x01800000, 0x7f800000},
704	{"lgamma", (funcptr)test_lgamma, args1, {NULL}, cases_uniform, 0x00100000, 0x7ff00000},
705	{"log1pf", (funcptr)mpfr_log1p, args1f, {NULL}, log1p_cases_float, 0, 0},
706	{"log1p", (funcptr)mpfr_log1p, args1, {NULL}, log1p_cases, 0, 0},
707	{"log2f", (funcptr)mpfr_log2, args1f, {NULL}, log_cases_float, 0, 0},
708	{"log2", (funcptr)mpfr_log2, args1, {NULL}, log_cases, 0, 0},
709	{"tgammaf", (funcptr)mpfr_gamma, args1f, {NULL}, cases_uniform_float, 0x2f800000, 0x43000000},
710	{"tgamma", (funcptr)mpfr_gamma, args1, {NULL}, cases_uniform, 0x3c000000, 0x40800000},
711	};
712
713	const int nfunctions = ( sizeof(functions)/sizeof(*functions) );
714
715	#define random_sign ( random_upto(1) ? 0x80000000 : 0 )
716
717	static int iszero(uint32 *x) {
718	return !((x[0] & 0x7FFFFFFF) || x[1]);
719	}
720
721
722	static void complex_log_cases(uint32 *out, uint32 param1,
723	uint32 param2) {
724	cases_uniform(out,0x00100000,0x7fefffff);
725	cases_uniform(out+2,0x00100000,0x7fefffff);
726	}
727
728
729	static void complex_log_cases_float(uint32 *out, uint32 param1,
730	uint32 param2) {
731	cases_uniform_float(out,0x00800000,0x7f7fffff);
732	cases_uniform_float(out+2,0x00800000,0x7f7fffff);
733	}
734
735	static void complex_cases_biased(uint32 *out, uint32 lowbound,
736	uint32 highbound) {
737	cases_biased(out,lowbound,highbound);
738	cases_biased(out+2,lowbound,highbound);
739	}
740
741	static void complex_cases_biased_float(uint32 *out, uint32 lowbound,
742	uint32 highbound) {
743	cases_biased_float(out,lowbound,highbound);
744	cases_biased_float(out+2,lowbound,highbound);
745	}
746
747	static void complex_cases_uniform(uint32 *out, uint32 lowbound,
748	uint32 highbound) {
749	cases_uniform(out,lowbound,highbound);
750	cases_uniform(out+2,lowbound,highbound);
751	}
752
753	static void complex_cases_uniform_float(uint32 *out, uint32 lowbound,
754	uint32 highbound) {
755	cases_uniform_float(out,lowbound,highbound);
756	cases_uniform(out+2,lowbound,highbound);
757	}
758
759	static void complex_pow_cases(uint32 *out, uint32 lowbound,
760	uint32 highbound) {
761	/*
762	* Generating non-overflowing cases for complex pow:
763	*
764	* Our base has both parts within the range [1/2,2], and hence
765	* its magnitude is within [1/2,2*sqrt(2)]. The magnitude of its
766	* logarithm in base 2 is therefore at most the magnitude of
767	* (log2(2*sqrt(2)) + i*pi/log(2)), or in other words
768	* hypot(3/2,pi/log(2)) = 4.77. So the magnitude of the exponent
769	* input must be at most our output magnitude limit (as a power
770	* of two) divided by that.
771	*
772	* I also set the output magnitude limit a bit low, because we
773	* don't guarantee (and neither does glibc) to prevent internal
774	* overflow in cases where the output _magnitude_ overflows but
775	* scaling it back down by cos and sin of the argument brings it
776	* back in range.
777	*/
778	cases_uniform(out,0x3fe00000, 0x40000000);
779	cases_uniform(out+2,0x3fe00000, 0x40000000);
780	cases_uniform(out+4,0x3f800000, 0x40600000);
781	cases_uniform(out+6,0x3f800000, 0x40600000);
782	}
783
784	static void complex_pow_cases_float(uint32 *out, uint32 lowbound,
785	uint32 highbound) {
786	/*
787	* Reasoning as above, though of course the detailed numbers are
788	* all different.
789	*/
790	cases_uniform_float(out,0x3f000000, 0x40000000);
791	cases_uniform_float(out+2,0x3f000000, 0x40000000);
792	cases_uniform_float(out+4,0x3d600000, 0x41900000);
793	cases_uniform_float(out+6,0x3d600000, 0x41900000);
794	}
795
796	static void complex_arithmetic_cases(uint32 *out, uint32 lowbound,
797	uint32 highbound) {
798	cases_uniform(out,0,0x7fefffff);
799	cases_uniform(out+2,0,0x7fefffff);
800	cases_uniform(out+4,0,0x7fefffff);
801	cases_uniform(out+6,0,0x7fefffff);
802	}
803
804	static void complex_arithmetic_cases_float(uint32 *out, uint32 lowbound,
805	uint32 highbound) {
806	cases_uniform_float(out,0,0x7f7fffff);
807	cases_uniform_float(out+2,0,0x7f7fffff);
808	cases_uniform_float(out+4,0,0x7f7fffff);
809	cases_uniform_float(out+6,0,0x7f7fffff);
810	}
811
812	/*
813	* Included from fplib test suite, in a compact self-contained
814	* form.
815	*/
816
817	void float32_case(uint32 *ret) {
818	int n, bits;
819	uint32 f;
820	static int premax, preptr;
821	static uint32 *specifics = NULL;
822
823	if (!ret) {
824	if (specifics)
825	free(ptr: specifics);
826	specifics = NULL;
827	premax = preptr = 0;
828	return;
829	}
830
831	if (!specifics) {
832	int exps[] = {
833	-127, -126, -125, -24, -4, -3, -2, -1, 0, 1, 2, 3, 4,
834	24, 29, 30, 31, 32, 61, 62, 63, 64, 126, 127, 128
835	};
836	int sign, eptr;
837	uint32 se, j;
838	/*
839	* We want a cross product of:
840	* - each of two sign bits (2)
841	* - each of the above (unbiased) exponents (25)
842	* - the following list of fraction parts:
843	* * zero (1)
844	* * all bits (1)
845	* * one-bit-set (23)
846	* * one-bit-clear (23)
847	* * one-bit-and-above (20: 3 are duplicates)
848	* * one-bit-and-below (20: 3 are duplicates)
849	* (total 88)
850	* (total 4400)
851	*/
852	specifics = malloc(size: 4400 * sizeof(*specifics));
853	preptr = 0;
854	for (sign = 0; sign <= 1; sign++) {
855	for (eptr = 0; eptr < sizeof(exps)/sizeof(*exps); eptr++) {
856	se = (sign ? 0x80000000 : 0) | ((exps[eptr]+127) << 23);
857	/*
858	* Zero.
859	*/
860	specifics[preptr++] = se | 0;
861	/*
862	* All bits.
863	*/
864	specifics[preptr++] = se | 0x7FFFFF;
865	/*
866	* One-bit-set.
867	*/
868	for (j = 1; j && j <= 0x400000; j <<= 1)
869	specifics[preptr++] = se | j;
870	/*
871	* One-bit-clear.
872	*/
873	for (j = 1; j && j <= 0x400000; j <<= 1)
874	specifics[preptr++] = se | (0x7FFFFF ^ j);
875	/*
876	* One-bit-and-everything-below.
877	*/
878	for (j = 2; j && j <= 0x100000; j <<= 1)
879	specifics[preptr++] = se | (2*j-1);
880	/*
881	* One-bit-and-everything-above.
882	*/
883	for (j = 4; j && j <= 0x200000; j <<= 1)
884	specifics[preptr++] = se | (0x7FFFFF ^ (j-1));
885	/*
886	* Done.
887	*/
888	}
889	}
890	assert(preptr == 4400);
891	premax = preptr;
892	}
893
894	/*
895	* Decide whether to return a pre or a random case.
896	*/
897	n = random32() % (premax+1);
898	if (n < preptr) {
899	/*
900	* Return pre[n].
901	*/
902	uint32 t;
903	t = specifics[n];
904	specifics[n] = specifics[preptr-1];
905	specifics[preptr-1] = t; /* (not really needed) */
906	preptr--;
907	*ret = t;
908	} else {
909	/*
910	* Random case.
911	* Sign and exponent:
912	* - FIXME
913	* Significand:
914	* - with prob 1/5, a totally random bit pattern
915	* - with prob 1/5, all 1s down to some point and then random
916	* - with prob 1/5, all 1s up to some point and then random
917	* - with prob 1/5, all 0s down to some point and then random
918	* - with prob 1/5, all 0s up to some point and then random
919	*/
920	n = random32() % 5;
921	f = random32(); /* some random bits */
922	bits = random32() % 22 + 1; /* 1-22 */
923	switch (n) {
924	case 0:
925	break; /* leave f alone */
926	case 1:
927	f |= (1<<bits)-1;
928	break;
929	case 2:
930	f &= ~((1<<bits)-1);
931	break;
932	case 3:
933	f |= ~((1<<bits)-1);
934	break;
935	case 4:
936	f &= (1<<bits)-1;
937	break;
938	}
939	f &= 0x7FFFFF;
940	f |= (random32() & 0xFF800000);/* FIXME - do better */
941	*ret = f;
942	}
943	}
944	static void float64_case(uint32 *ret) {
945	int n, bits;
946	uint32 f, g;
947	static int premax, preptr;
948	static uint32 (*specifics)[2] = NULL;
949
950	if (!ret) {
951	if (specifics)
952	free(ptr: specifics);
953	specifics = NULL;
954	premax = preptr = 0;
955	return;
956	}
957
958	if (!specifics) {
959	int exps[] = {
960	-1023, -1022, -1021, -129, -128, -127, -126, -53, -4, -3, -2,
961	-1, 0, 1, 2, 3, 4, 29, 30, 31, 32, 53, 61, 62, 63, 64, 127,
962	128, 129, 1022, 1023, 1024
963	};
964	int sign, eptr;
965	uint32 se, j;
966	/*
967	* We want a cross product of:
968	* - each of two sign bits (2)
969	* - each of the above (unbiased) exponents (32)
970	* - the following list of fraction parts:
971	* * zero (1)
972	* * all bits (1)
973	* * one-bit-set (52)
974	* * one-bit-clear (52)
975	* * one-bit-and-above (49: 3 are duplicates)
976	* * one-bit-and-below (49: 3 are duplicates)
977	* (total 204)
978	* (total 13056)
979	*/
980	specifics = malloc(size: 13056 * sizeof(*specifics));
981	preptr = 0;
982	for (sign = 0; sign <= 1; sign++) {
983	for (eptr = 0; eptr < sizeof(exps)/sizeof(*exps); eptr++) {
984	se = (sign ? 0x80000000 : 0) | ((exps[eptr]+1023) << 20);
985	/*
986	* Zero.
987	*/
988	specifics[preptr][0] = 0;
989	specifics[preptr][1] = 0;
990	specifics[preptr++][0] |= se;
991	/*
992	* All bits.
993	*/
994	specifics[preptr][0] = 0xFFFFF;
995	specifics[preptr][1] = ~0;
996	specifics[preptr++][0] |= se;
997	/*
998	* One-bit-set.
999	*/
1000	for (j = 1; j && j <= 0x80000000; j <<= 1) {
1001	specifics[preptr][0] = 0;
1002	specifics[preptr][1] = j;
1003	specifics[preptr++][0] |= se;
1004	if (j & 0xFFFFF) {
1005	specifics[preptr][0] = j;
1006	specifics[preptr][1] = 0;
1007	specifics[preptr++][0] |= se;
1008	}
1009	}
1010	/*
1011	* One-bit-clear.
1012	*/
1013	for (j = 1; j && j <= 0x80000000; j <<= 1) {
1014	specifics[preptr][0] = 0xFFFFF;
1015	specifics[preptr][1] = ~j;
1016	specifics[preptr++][0] |= se;
1017	if (j & 0xFFFFF) {
1018	specifics[preptr][0] = 0xFFFFF ^ j;
1019	specifics[preptr][1] = ~0;
1020	specifics[preptr++][0] |= se;
1021	}
1022	}
1023	/*
1024	* One-bit-and-everything-below.
1025	*/
1026	for (j = 2; j && j <= 0x80000000; j <<= 1) {
1027	specifics[preptr][0] = 0;
1028	specifics[preptr][1] = 2*j-1;
1029	specifics[preptr++][0] |= se;
1030	}
1031	for (j = 1; j && j <= 0x20000; j <<= 1) {
1032	specifics[preptr][0] = 2*j-1;
1033	specifics[preptr][1] = ~0;
1034	specifics[preptr++][0] |= se;
1035	}
1036	/*
1037	* One-bit-and-everything-above.
1038	*/
1039	for (j = 4; j && j <= 0x80000000; j <<= 1) {
1040	specifics[preptr][0] = 0xFFFFF;
1041	specifics[preptr][1] = ~(j-1);
1042	specifics[preptr++][0] |= se;
1043	}
1044	for (j = 1; j && j <= 0x40000; j <<= 1) {
1045	specifics[preptr][0] = 0xFFFFF ^ (j-1);
1046	specifics[preptr][1] = 0;
1047	specifics[preptr++][0] |= se;
1048	}
1049	/*
1050	* Done.
1051	*/
1052	}
1053	}
1054	assert(preptr == 13056);
1055	premax = preptr;
1056	}
1057
1058	/*
1059	* Decide whether to return a pre or a random case.
1060	*/
1061	n = (uint32) random32() % (uint32) (premax+1);
1062	if (n < preptr) {
1063	/*
1064	* Return pre[n].
1065	*/
1066	uint32 t;
1067	t = specifics[n][0];
1068	specifics[n][0] = specifics[preptr-1][0];
1069	specifics[preptr-1][0] = t; /* (not really needed) */
1070	ret[0] = t;
1071	t = specifics[n][1];
1072	specifics[n][1] = specifics[preptr-1][1];
1073	specifics[preptr-1][1] = t; /* (not really needed) */
1074	ret[1] = t;
1075	preptr--;
1076	} else {
1077	/*
1078	* Random case.
1079	* Sign and exponent:
1080	* - FIXME
1081	* Significand:
1082	* - with prob 1/5, a totally random bit pattern
1083	* - with prob 1/5, all 1s down to some point and then random
1084	* - with prob 1/5, all 1s up to some point and then random
1085	* - with prob 1/5, all 0s down to some point and then random
1086	* - with prob 1/5, all 0s up to some point and then random
1087	*/
1088	n = random32() % 5;
1089	f = random32(); /* some random bits */
1090	g = random32(); /* some random bits */
1091	bits = random32() % 51 + 1; /* 1-51 */
1092	switch (n) {
1093	case 0:
1094	break; /* leave f alone */
1095	case 1:
1096	if (bits <= 32)
1097	f |= (1<<bits)-1;
1098	else {
1099	bits -= 32;
1100	g |= (1<<bits)-1;
1101	f = ~0;
1102	}
1103	break;
1104	case 2:
1105	if (bits <= 32)
1106	f &= ~((1<<bits)-1);
1107	else {
1108	bits -= 32;
1109	g &= ~((1<<bits)-1);
1110	f = 0;
1111	}
1112	break;
1113	case 3:
1114	if (bits <= 32)
1115	g &= (1<<bits)-1;
1116	else {
1117	bits -= 32;
1118	f &= (1<<bits)-1;
1119	g = 0;
1120	}
1121	break;
1122	case 4:
1123	if (bits <= 32)
1124	g |= ~((1<<bits)-1);
1125	else {
1126	bits -= 32;
1127	f |= ~((1<<bits)-1);
1128	g = ~0;
1129	}
1130	break;
1131	}
1132	g &= 0xFFFFF;
1133	g |= (random32() & 0xFFF00000);/* FIXME - do better */
1134	ret[0] = g;
1135	ret[1] = f;
1136	}
1137	}
1138
1139	static void cases_biased(uint32 *out, uint32 lowbound,
1140	uint32 highbound) {
1141	do {
1142	out[0] = highbound - random_upto_biased(limit: highbound-lowbound, bias: 8);
1143	out[1] = random_upto(limit: 0xFFFFFFFF);
1144	out[0] |= random_sign;
1145	} while (iszero(x: out)); /* rule out zero */
1146	}
1147
1148	static void cases_biased_positive(uint32 *out, uint32 lowbound,
1149	uint32 highbound) {
1150	do {
1151	out[0] = highbound - random_upto_biased(limit: highbound-lowbound, bias: 8);
1152	out[1] = random_upto(limit: 0xFFFFFFFF);
1153	} while (iszero(x: out)); /* rule out zero */
1154	}
1155
1156	static void cases_biased_float(uint32 *out, uint32 lowbound,
1157	uint32 highbound) {
1158	do {
1159	out[0] = highbound - random_upto_biased(limit: highbound-lowbound, bias: 8);
1160	out[1] = 0;
1161	out[0] |= random_sign;
1162	} while (iszero(x: out)); /* rule out zero */
1163	}
1164
1165	static void cases_semi1(uint32 *out, uint32 param1,
1166	uint32 param2) {
1167	float64_case(ret: out);
1168	}
1169
1170	static void cases_semi1_float(uint32 *out, uint32 param1,
1171	uint32 param2) {
1172	float32_case(ret: out);
1173	}
1174
1175	static void cases_semi2(uint32 *out, uint32 param1,
1176	uint32 param2) {
1177	float64_case(ret: out);
1178	float64_case(ret: out+2);
1179	}
1180
1181	static void cases_semi2_float(uint32 *out, uint32 param1,
1182	uint32 param2) {
1183	float32_case(ret: out);
1184	float32_case(ret: out+2);
1185	}
1186
1187	static void cases_ldexp(uint32 *out, uint32 param1,
1188	uint32 param2) {
1189	float64_case(ret: out);
1190	out[2] = random_upto(limit: 2048)-1024;
1191	}
1192
1193	static void cases_ldexp_float(uint32 *out, uint32 param1,
1194	uint32 param2) {
1195	float32_case(ret: out);
1196	out[2] = random_upto(limit: 256)-128;
1197	}
1198
1199	static void cases_uniform(uint32 *out, uint32 lowbound,
1200	uint32 highbound) {
1201	do {
1202	out[0] = highbound - random_upto(limit: highbound-lowbound);
1203	out[1] = random_upto(limit: 0xFFFFFFFF);
1204	out[0] |= random_sign;
1205	} while (iszero(x: out)); /* rule out zero */
1206	}
1207	static void cases_uniform_float(uint32 *out, uint32 lowbound,
1208	uint32 highbound) {
1209	do {
1210	out[0] = highbound - random_upto(limit: highbound-lowbound);
1211	out[1] = 0;
1212	out[0] |= random_sign;
1213	} while (iszero(x: out)); /* rule out zero */
1214	}
1215
1216	static void cases_uniform_positive(uint32 *out, uint32 lowbound,
1217	uint32 highbound) {
1218	do {
1219	out[0] = highbound - random_upto(limit: highbound-lowbound);
1220	out[1] = random_upto(limit: 0xFFFFFFFF);
1221	} while (iszero(x: out)); /* rule out zero */
1222	}
1223	static void cases_uniform_float_positive(uint32 *out, uint32 lowbound,
1224	uint32 highbound) {
1225	do {
1226	out[0] = highbound - random_upto(limit: highbound-lowbound);
1227	out[1] = 0;
1228	} while (iszero(x: out)); /* rule out zero */
1229	}
1230
1231
1232	static void log_cases(uint32 *out, uint32 param1,
1233	uint32 param2) {
1234	do {
1235	out[0] = random_upto(limit: 0x7FEFFFFF);
1236	out[1] = random_upto(limit: 0xFFFFFFFF);
1237	} while (iszero(x: out)); /* rule out zero */
1238	}
1239
1240	static void log_cases_float(uint32 *out, uint32 param1,
1241	uint32 param2) {
1242	do {
1243	out[0] = random_upto(limit: 0x7F7FFFFF);
1244	out[1] = 0;
1245	} while (iszero(x: out)); /* rule out zero */
1246	}
1247
1248	static void log1p_cases(uint32 *out, uint32 param1, uint32 param2)
1249	{
1250	uint32 sign = random_sign;
1251	if (sign == 0) {
1252	cases_uniform_positive(out, lowbound: 0x3c700000, highbound: 0x43400000);
1253	} else {
1254	cases_uniform_positive(out, lowbound: 0x3c000000, highbound: 0x3ff00000);
1255	}
1256	out[0] |= sign;
1257	}
1258
1259	static void log1p_cases_float(uint32 *out, uint32 param1, uint32 param2)
1260	{
1261	uint32 sign = random_sign;
1262	if (sign == 0) {
1263	cases_uniform_float_positive(out, lowbound: 0x32000000, highbound: 0x4c000000);
1264	} else {
1265	cases_uniform_float_positive(out, lowbound: 0x30000000, highbound: 0x3f800000);
1266	}
1267	out[0] |= sign;
1268	}
1269
1270	static void minmax_cases(uint32 *out, uint32 param1, uint32 param2)
1271	{
1272	do {
1273	out[0] = random_upto(limit: 0x7FEFFFFF);
1274	out[1] = random_upto(limit: 0xFFFFFFFF);
1275	out[0] |= random_sign;
1276	out[2] = random_upto(limit: 0x7FEFFFFF);
1277	out[3] = random_upto(limit: 0xFFFFFFFF);
1278	out[2] |= random_sign;
1279	} while (iszero(x: out)); /* rule out zero */
1280	}
1281
1282	static void minmax_cases_float(uint32 *out, uint32 param1, uint32 param2)
1283	{
1284	do {
1285	out[0] = random_upto(limit: 0x7F7FFFFF);
1286	out[1] = 0;
1287	out[0] |= random_sign;
1288	out[2] = random_upto(limit: 0x7F7FFFFF);
1289	out[3] = 0;
1290	out[2] |= random_sign;
1291	} while (iszero(x: out)); /* rule out zero */
1292	}
1293
1294	static void rred_cases(uint32 *out, uint32 param1,
1295	uint32 param2) {
1296	do {
1297	out[0] = ((0x3fc00000 + random_upto(limit: 0x036fffff)) |
1298	(random_upto(limit: 1) << 31));
1299	out[1] = random_upto(limit: 0xFFFFFFFF);
1300	} while (iszero(x: out)); /* rule out zero */
1301	}
1302
1303	static void rred_cases_float(uint32 *out, uint32 param1,
1304	uint32 param2) {
1305	do {
1306	out[0] = ((0x3e000000 + random_upto(limit: 0x0cffffff)) |
1307	(random_upto(limit: 1) << 31));
1308	out[1] = 0; /* for iszero */
1309	} while (iszero(x: out)); /* rule out zero */
1310	}
1311
1312	static void atan2_cases(uint32 *out, uint32 param1,
1313	uint32 param2) {
1314	do {
1315	int expdiff = random_upto(limit: 101)-51;
1316	int swap;
1317	if (expdiff < 0) {
1318	expdiff = -expdiff;
1319	swap = 2;
1320	} else
1321	swap = 0;
1322	out[swap ^ 0] = random_upto(limit: 0x7FEFFFFF-((expdiff+1)<<20));
1323	out[swap ^ 2] = random_upto(limit: ((expdiff+1)<<20)-1) + out[swap ^ 0];
1324	out[1] = random_upto(limit: 0xFFFFFFFF);
1325	out[3] = random_upto(limit: 0xFFFFFFFF);
1326	out[0] |= random_sign;
1327	out[2] |= random_sign;
1328	} while (iszero(x: out) || iszero(x: out+2));/* rule out zero */
1329	}
1330
1331	static void atan2_cases_float(uint32 *out, uint32 param1,
1332	uint32 param2) {
1333	do {
1334	int expdiff = random_upto(limit: 44)-22;
1335	int swap;
1336	if (expdiff < 0) {
1337	expdiff = -expdiff;
1338	swap = 2;
1339	} else
1340	swap = 0;
1341	out[swap ^ 0] = random_upto(limit: 0x7F7FFFFF-((expdiff+1)<<23));
1342	out[swap ^ 2] = random_upto(limit: ((expdiff+1)<<23)-1) + out[swap ^ 0];
1343	out[0] |= random_sign;
1344	out[2] |= random_sign;
1345	out[1] = out[3] = 0; /* for iszero */
1346	} while (iszero(x: out) || iszero(x: out+2));/* rule out zero */
1347	}
1348
1349	static void pow_cases(uint32 *out, uint32 param1,
1350	uint32 param2) {
1351	/*
1352	* Pick an exponent e (-0x33 to +0x7FE) for x, and here's the
1353	* range of numbers we can use as y:
1354	*
1355	* For e < 0x3FE, the range is [-0x400/(0x3FE-e),+0x432/(0x3FE-e)]
1356	* For e > 0x3FF, the range is [-0x432/(e-0x3FF),+0x400/(e-0x3FF)]
1357	*
1358	* For e == 0x3FE or e == 0x3FF, the range gets infinite at one
1359	* end or the other, so we have to be cleverer: pick a number n
1360	* of useful bits in the mantissa (1 thru 52, so 1 must imply
1361	* 0x3ff00000.00000001 whereas 52 is anything at least as big
1362	* as 0x3ff80000.00000000; for e == 0x3fe, 1 necessarily means
1363	* 0x3fefffff.ffffffff and 52 is anything at most as big as
1364	* 0x3fe80000.00000000). Then, as it happens, a sensible
1365	* maximum power is 2^(63-n) for e == 0x3fe, and 2^(62-n) for
1366	* e == 0x3ff.
1367	*
1368	* We inevitably get some overflows in approximating the log
1369	* curves by these nasty step functions, but that's all right -
1370	* we do want _some_ overflows to be tested.
1371	*
1372	* Having got that, then, it's just a matter of inventing a
1373	* probability distribution for all of this.
1374	*/
1375	int e, n;
1376	uint32 dmin, dmax;
1377	const uint32 pmin = 0x3e100000;
1378
1379	/*
1380	* Generate exponents in a slightly biased fashion.
1381	*/
1382	e = (random_upto(limit: 1) ? /* is exponent small or big? */
1383	0x3FE - random_upto_biased(limit: 0x431,bias: 2) : /* small */
1384	0x3FF + random_upto_biased(limit: 0x3FF,bias: 2)); /* big */
1385
1386	/*
1387	* Now split into cases.
1388	*/
1389	if (e < 0x3FE || e > 0x3FF) {
1390	uint32 imin, imax;
1391	if (e < 0x3FE)
1392	imin = 0x40000 / (0x3FE - e), imax = 0x43200 / (0x3FE - e);
1393	else
1394	imin = 0x43200 / (e - 0x3FF), imax = 0x40000 / (e - 0x3FF);
1395	/* Power range runs from -imin to imax. Now convert to doubles */
1396	dmin = doubletop(x: imin, scale: -8);
1397	dmax = doubletop(x: imax, scale: -8);
1398	/* Compute the number of mantissa bits. */
1399	n = (e > 0 ? 53 : 52+e);
1400	} else {
1401	/* Critical exponents. Generate a top bit index. */
1402	n = 52 - random_upto_biased(limit: 51, bias: 4);
1403	if (e == 0x3FE)
1404	dmax = 63 - n;
1405	else
1406	dmax = 62 - n;
1407	dmax = (dmax << 20) + 0x3FF00000;
1408	dmin = dmax;
1409	}
1410	/* Generate a mantissa. */
1411	if (n <= 32) {
1412	out[0] = 0;
1413	out[1] = random_upto(limit: (1 << (n-1)) - 1) + (1 << (n-1));
1414	} else if (n == 33) {
1415	out[0] = 1;
1416	out[1] = random_upto(limit: 0xFFFFFFFF);
1417	} else if (n > 33) {
1418	out[0] = random_upto(limit: (1 << (n-33)) - 1) + (1 << (n-33));
1419	out[1] = random_upto(limit: 0xFFFFFFFF);
1420	}
1421	/* Negate the mantissa if e == 0x3FE. */
1422	if (e == 0x3FE) {
1423	out[1] = -out[1];
1424	out[0] = -out[0];
1425	if (out[1]) out[0]--;
1426	}
1427	/* Put the exponent on. */
1428	out[0] &= 0xFFFFF;
1429	out[0] |= ((e > 0 ? e : 0) << 20);
1430	/* Generate a power. Powers don't go below 2^-30. */
1431	if (random_upto(limit: 1)) {
1432	/* Positive power */
1433	out[2] = dmax - random_upto_biased(limit: dmax-pmin, bias: 10);
1434	} else {
1435	/* Negative power */
1436	out[2] = (dmin - random_upto_biased(limit: dmin-pmin, bias: 10)) | 0x80000000;
1437	}
1438	out[3] = random_upto(limit: 0xFFFFFFFF);
1439	}
1440	static void pow_cases_float(uint32 *out, uint32 param1,
1441	uint32 param2) {
1442	/*
1443	* Pick an exponent e (-0x16 to +0xFE) for x, and here's the
1444	* range of numbers we can use as y:
1445	*
1446	* For e < 0x7E, the range is [-0x80/(0x7E-e),+0x95/(0x7E-e)]
1447	* For e > 0x7F, the range is [-0x95/(e-0x7F),+0x80/(e-0x7F)]
1448	*
1449	* For e == 0x7E or e == 0x7F, the range gets infinite at one
1450	* end or the other, so we have to be cleverer: pick a number n
1451	* of useful bits in the mantissa (1 thru 23, so 1 must imply
1452	* 0x3f800001 whereas 23 is anything at least as big as
1453	* 0x3fc00000; for e == 0x7e, 1 necessarily means 0x3f7fffff
1454	* and 23 is anything at most as big as 0x3f400000). Then, as
1455	* it happens, a sensible maximum power is 2^(31-n) for e ==
1456	* 0x7e, and 2^(30-n) for e == 0x7f.
1457	*
1458	* We inevitably get some overflows in approximating the log
1459	* curves by these nasty step functions, but that's all right -
1460	* we do want _some_ overflows to be tested.
1461	*
1462	* Having got that, then, it's just a matter of inventing a
1463	* probability distribution for all of this.
1464	*/
1465	int e, n;
1466	uint32 dmin, dmax;
1467	const uint32 pmin = 0x38000000;
1468
1469	/*
1470	* Generate exponents in a slightly biased fashion.
1471	*/
1472	e = (random_upto(limit: 1) ? /* is exponent small or big? */
1473	0x7E - random_upto_biased(limit: 0x94,bias: 2) : /* small */
1474	0x7F + random_upto_biased(limit: 0x7f,bias: 2)); /* big */
1475
1476	/*
1477	* Now split into cases.
1478	*/
1479	if (e < 0x7E || e > 0x7F) {
1480	uint32 imin, imax;
1481	if (e < 0x7E)
1482	imin = 0x8000 / (0x7e - e), imax = 0x9500 / (0x7e - e);
1483	else
1484	imin = 0x9500 / (e - 0x7f), imax = 0x8000 / (e - 0x7f);
1485	/* Power range runs from -imin to imax. Now convert to doubles */
1486	dmin = floatval(x: imin, scale: -8);
1487	dmax = floatval(x: imax, scale: -8);
1488	/* Compute the number of mantissa bits. */
1489	n = (e > 0 ? 24 : 23+e);
1490	} else {
1491	/* Critical exponents. Generate a top bit index. */
1492	n = 23 - random_upto_biased(limit: 22, bias: 4);
1493	if (e == 0x7E)
1494	dmax = 31 - n;
1495	else
1496	dmax = 30 - n;
1497	dmax = (dmax << 23) + 0x3F800000;
1498	dmin = dmax;
1499	}
1500	/* Generate a mantissa. */
1501	out[0] = random_upto(limit: (1 << (n-1)) - 1) + (1 << (n-1));
1502	out[1] = 0;
1503	/* Negate the mantissa if e == 0x7E. */
1504	if (e == 0x7E) {
1505	out[0] = -out[0];
1506	}
1507	/* Put the exponent on. */
1508	out[0] &= 0x7FFFFF;
1509	out[0] |= ((e > 0 ? e : 0) << 23);
1510	/* Generate a power. Powers don't go below 2^-15. */
1511	if (random_upto(limit: 1)) {
1512	/* Positive power */
1513	out[2] = dmax - random_upto_biased(limit: dmax-pmin, bias: 10);
1514	} else {
1515	/* Negative power */
1516	out[2] = (dmin - random_upto_biased(limit: dmin-pmin, bias: 10)) | 0x80000000;
1517	}
1518	out[3] = 0;
1519	}
1520
1521	void vet_for_decline(Testable *fn, uint32 *args, uint32 *result, int got_errno_in) {
1522	int declined = 0;
1523
1524	switch (fn->type) {
1525	case args1:
1526	case rred:
1527	case semi1:
1528	case t_frexp:
1529	case t_modf:
1530	case classify:
1531	case t_ldexp:
1532	declined |= lib_fo && is_dhard(args+0);
1533	break;
1534	case args1f:
1535	case rredf:
1536	case semi1f:
1537	case t_frexpf:
1538	case t_modff:
1539	case classifyf:
1540	declined |= lib_fo && is_shard(args+0);
1541	break;
1542	case args2:
1543	case semi2:
1544	case args1c:
1545	case args1cr:
1546	case compare:
1547	declined |= lib_fo && is_dhard(args+0);
1548	declined |= lib_fo && is_dhard(args+2);
1549	break;
1550	case args2f:
1551	case semi2f:
1552	case t_ldexpf:
1553	case comparef:
1554	case args1fc:
1555	case args1fcr:
1556	declined |= lib_fo && is_shard(args+0);
1557	declined |= lib_fo && is_shard(args+2);
1558	break;
1559	case args2c:
1560	declined |= lib_fo && is_dhard(args+0);
1561	declined |= lib_fo && is_dhard(args+2);
1562	declined |= lib_fo && is_dhard(args+4);
1563	declined |= lib_fo && is_dhard(args+6);
1564	break;
1565	case args2fc:
1566	declined |= lib_fo && is_shard(args+0);
1567	declined |= lib_fo && is_shard(args+2);
1568	declined |= lib_fo && is_shard(args+4);
1569	declined |= lib_fo && is_shard(args+6);
1570	break;
1571	}
1572
1573	switch (fn->type) {
1574	case args1: /* return an extra-precise result */
1575	case args2:
1576	case rred:
1577	case semi1: /* return a double result */
1578	case semi2:
1579	case t_ldexp:
1580	case t_frexp: /* return double * int */
1581	case args1cr:
1582	declined |= lib_fo && is_dhard(result);
1583	break;
1584	case args1f:
1585	case args2f:
1586	case rredf:
1587	case semi1f:
1588	case semi2f:
1589	case t_ldexpf:
1590	case args1fcr:
1591	declined |= lib_fo && is_shard(result);
1592	break;
1593	case t_modf: /* return double * double */
1594	declined |= lib_fo && is_dhard(result+0);
1595	declined |= lib_fo && is_dhard(result+2);
1596	break;
1597	case t_modff: /* return float * float */
1598	declined |= lib_fo && is_shard(result+2);
1599	/* fall through */
1600	case t_frexpf: /* return float * int */
1601	declined |= lib_fo && is_shard(result+0);
1602	break;
1603	case args1c:
1604	case args2c:
1605	declined |= lib_fo && is_dhard(result+0);
1606	declined |= lib_fo && is_dhard(result+4);
1607	break;
1608	case args1fc:
1609	case args2fc:
1610	declined |= lib_fo && is_shard(result+0);
1611	declined |= lib_fo && is_shard(result+4);
1612	break;
1613	}
1614
1615	/* Expect basic arithmetic tests to be declined if the command
1616	* line said that would happen */
1617	declined |= (lib_no_arith && (fn->func == (funcptr)mpc_add ||
1618	fn->func == (funcptr)mpc_sub ||
1619	fn->func == (funcptr)mpc_mul ||
1620	fn->func == (funcptr)mpc_div));
1621
1622	if (!declined) {
1623	if (got_errno_in)
1624	ntests++;
1625	else
1626	ntests += 3;
1627	}
1628	}
1629
1630	void docase(Testable *fn, uint32 *args) {
1631	uint32 result[8]; /* real part in first 4, imaginary part in last 4 */
1632	char *errstr = NULL;
1633	mpfr_t a, b, r;
1634	mpc_t ac, bc, rc;
1635	int rejected, printextra;
1636	wrapperctx ctx;
1637
1638	mpfr_init2(a, MPFR_PREC);
1639	mpfr_init2(b, MPFR_PREC);
1640	mpfr_init2(r, MPFR_PREC);
1641	mpc_init2(ac, MPFR_PREC);
1642	mpc_init2(bc, MPFR_PREC);
1643	mpc_init2(rc, MPFR_PREC);
1644
1645	printf(format: "func=%s", fn->name);
1646
1647	rejected = 0; /* FIXME */
1648
1649	switch (fn->type) {
1650	case args1:
1651	case rred:
1652	case semi1:
1653	case t_frexp:
1654	case t_modf:
1655	case classify:
1656	printf(format: " op1=%08x.%08x", args[0], args[1]);
1657	break;
1658	case args1f:
1659	case rredf:
1660	case semi1f:
1661	case t_frexpf:
1662	case t_modff:
1663	case classifyf:
1664	printf(format: " op1=%08x", args[0]);
1665	break;
1666	case args2:
1667	case semi2:
1668	case compare:
1669	printf(format: " op1=%08x.%08x", args[0], args[1]);
1670	printf(format: " op2=%08x.%08x", args[2], args[3]);
1671	break;
1672	case args2f:
1673	case semi2f:
1674	case t_ldexpf:
1675	case comparef:
1676	printf(format: " op1=%08x", args[0]);
1677	printf(format: " op2=%08x", args[2]);
1678	break;
1679	case t_ldexp:
1680	printf(format: " op1=%08x.%08x", args[0], args[1]);
1681	printf(format: " op2=%08x", args[2]);
1682	break;
1683	case args1c:
1684	case args1cr:
1685	printf(format: " op1r=%08x.%08x", args[0], args[1]);
1686	printf(format: " op1i=%08x.%08x", args[2], args[3]);
1687	break;
1688	case args2c:
1689	printf(format: " op1r=%08x.%08x", args[0], args[1]);
1690	printf(format: " op1i=%08x.%08x", args[2], args[3]);
1691	printf(format: " op2r=%08x.%08x", args[4], args[5]);
1692	printf(format: " op2i=%08x.%08x", args[6], args[7]);
1693	break;
1694	case args1fc:
1695	case args1fcr:
1696	printf(format: " op1r=%08x", args[0]);
1697	printf(format: " op1i=%08x", args[2]);
1698	break;
1699	case args2fc:
1700	printf(format: " op1r=%08x", args[0]);
1701	printf(format: " op1i=%08x", args[2]);
1702	printf(format: " op2r=%08x", args[4]);
1703	printf(format: " op2i=%08x", args[6]);
1704	break;
1705	default:
1706	fprintf(stderr, format: "internal inconsistency?!\n");
1707	abort();
1708	}
1709
1710	if (rejected == 2) {
1711	printf(format: " - test case rejected\n");
1712	goto cleanup;
1713	}
1714
1715	wrapper_init(ctx: &ctx);
1716
1717	if (rejected == 0) {
1718	switch (fn->type) {
1719	case args1:
1720	set_mpfr_d(x: a, h: args[0], l: args[1]);
1721	wrapper_op_real(ctx: &ctx, r: a, size: 2, ieee: args);
1722	((testfunc1)(fn->func))(r, a, GMP_RNDN);
1723	get_mpfr_d(x: r, h: &result[0], l: &result[1], extra: &result[2]);
1724	wrapper_result_real(ctx: &ctx, r, size: 2, ieee: result);
1725	if (wrapper_run(ctx: &ctx, wrappers: fn->wrappers))
1726	get_mpfr_d(x: r, h: &result[0], l: &result[1], extra: &result[2]);
1727	break;
1728	case args1cr:
1729	set_mpc_d(z: ac, rh: args[0], rl: args[1], ih: args[2], il: args[3]);
1730	wrapper_op_complex(ctx: &ctx, c: ac, size: 2, ieee: args);
1731	((testfunc1cr)(fn->func))(r, ac, GMP_RNDN);
1732	get_mpfr_d(x: r, h: &result[0], l: &result[1], extra: &result[2]);
1733	wrapper_result_real(ctx: &ctx, r, size: 2, ieee: result);
1734	if (wrapper_run(ctx: &ctx, wrappers: fn->wrappers))
1735	get_mpfr_d(x: r, h: &result[0], l: &result[1], extra: &result[2]);
1736	break;
1737	case args1f:
1738	set_mpfr_f(x: a, f: args[0]);
1739	wrapper_op_real(ctx: &ctx, r: a, size: 1, ieee: args);
1740	((testfunc1)(fn->func))(r, a, GMP_RNDN);
1741	get_mpfr_f(x: r, f: &result[0], extra: &result[1]);
1742	wrapper_result_real(ctx: &ctx, r, size: 1, ieee: result);
1743	if (wrapper_run(ctx: &ctx, wrappers: fn->wrappers))
1744	get_mpfr_f(x: r, f: &result[0], extra: &result[1]);
1745	break;
1746	case args1fcr:
1747	set_mpc_f(z: ac, r: args[0], i: args[2]);
1748	wrapper_op_complex(ctx: &ctx, c: ac, size: 1, ieee: args);
1749	((testfunc1cr)(fn->func))(r, ac, GMP_RNDN);
1750	get_mpfr_f(x: r, f: &result[0], extra: &result[1]);
1751	wrapper_result_real(ctx: &ctx, r, size: 1, ieee: result);
1752	if (wrapper_run(ctx: &ctx, wrappers: fn->wrappers))
1753	get_mpfr_f(x: r, f: &result[0], extra: &result[1]);
1754	break;
1755	case args2:
1756	set_mpfr_d(x: a, h: args[0], l: args[1]);
1757	wrapper_op_real(ctx: &ctx, r: a, size: 2, ieee: args);
1758	set_mpfr_d(x: b, h: args[2], l: args[3]);
1759	wrapper_op_real(ctx: &ctx, r: b, size: 2, ieee: args+2);
1760	((testfunc2)(fn->func))(r, a, b, GMP_RNDN);
1761	get_mpfr_d(x: r, h: &result[0], l: &result[1], extra: &result[2]);
1762	wrapper_result_real(ctx: &ctx, r, size: 2, ieee: result);
1763	if (wrapper_run(ctx: &ctx, wrappers: fn->wrappers))
1764	get_mpfr_d(x: r, h: &result[0], l: &result[1], extra: &result[2]);
1765	break;
1766	case args2f:
1767	set_mpfr_f(x: a, f: args[0]);
1768	wrapper_op_real(ctx: &ctx, r: a, size: 1, ieee: args);
1769	set_mpfr_f(x: b, f: args[2]);
1770	wrapper_op_real(ctx: &ctx, r: b, size: 1, ieee: args+2);
1771	((testfunc2)(fn->func))(r, a, b, GMP_RNDN);
1772	get_mpfr_f(x: r, f: &result[0], extra: &result[1]);
1773	wrapper_result_real(ctx: &ctx, r, size: 1, ieee: result);
1774	if (wrapper_run(ctx: &ctx, wrappers: fn->wrappers))
1775	get_mpfr_f(x: r, f: &result[0], extra: &result[1]);
1776	break;
1777	case rred:
1778	set_mpfr_d(x: a, h: args[0], l: args[1]);
1779	wrapper_op_real(ctx: &ctx, r: a, size: 2, ieee: args);
1780	((testrred)(fn->func))(r, a, (int *)&result[3]);
1781	get_mpfr_d(x: r, h: &result[0], l: &result[1], extra: &result[2]);
1782	wrapper_result_real(ctx: &ctx, r, size: 2, ieee: result);
1783	/* We never need to mess about with the integer auxiliary
1784	* output. */
1785	if (wrapper_run(ctx: &ctx, wrappers: fn->wrappers))
1786	get_mpfr_d(x: r, h: &result[0], l: &result[1], extra: &result[2]);
1787	break;
1788	case rredf:
1789	set_mpfr_f(x: a, f: args[0]);
1790	wrapper_op_real(ctx: &ctx, r: a, size: 1, ieee: args);
1791	((testrred)(fn->func))(r, a, (int *)&result[3]);
1792	get_mpfr_f(x: r, f: &result[0], extra: &result[1]);
1793	wrapper_result_real(ctx: &ctx, r, size: 1, ieee: result);
1794	/* We never need to mess about with the integer auxiliary
1795	* output. */
1796	if (wrapper_run(ctx: &ctx, wrappers: fn->wrappers))
1797	get_mpfr_f(x: r, f: &result[0], extra: &result[1]);
1798	break;
1799	case semi1:
1800	case semi1f:
1801	errstr = ((testsemi1)(fn->func))(args, result);
1802	break;
1803	case semi2:
1804	case compare:
1805	errstr = ((testsemi2)(fn->func))(args, args+2, result);
1806	break;
1807	case semi2f:
1808	case comparef:
1809	case t_ldexpf:
1810	errstr = ((testsemi2f)(fn->func))(args, args+2, result);
1811	break;
1812	case t_ldexp:
1813	errstr = ((testldexp)(fn->func))(args, args+2, result);
1814	break;
1815	case t_frexp:
1816	errstr = ((testfrexp)(fn->func))(args, result, result+2);
1817	break;
1818	case t_frexpf:
1819	errstr = ((testfrexp)(fn->func))(args, result, result+2);
1820	break;
1821	case t_modf:
1822	errstr = ((testmodf)(fn->func))(args, result, result+2);
1823	break;
1824	case t_modff:
1825	errstr = ((testmodf)(fn->func))(args, result, result+2);
1826	break;
1827	case classify:
1828	errstr = ((testclassify)(fn->func))(args, &result[0]);
1829	break;
1830	case classifyf:
1831	errstr = ((testclassifyf)(fn->func))(args, &result[0]);
1832	break;
1833	case args1c:
1834	set_mpc_d(z: ac, rh: args[0], rl: args[1], ih: args[2], il: args[3]);
1835	wrapper_op_complex(ctx: &ctx, c: ac, size: 2, ieee: args);
1836	((testfunc1c)(fn->func))(rc, ac, MPC_RNDNN);
1837	get_mpc_d(z: rc, rh: &result[0], rl: &result[1], rextra: &result[2], ih: &result[4], il: &result[5], iextra: &result[6]);
1838	wrapper_result_complex(ctx: &ctx, c: rc, size: 2, ieee: result);
1839	if (wrapper_run(ctx: &ctx, wrappers: fn->wrappers))
1840	get_mpc_d(z: rc, rh: &result[0], rl: &result[1], rextra: &result[2], ih: &result[4], il: &result[5], iextra: &result[6]);
1841	break;
1842	case args2c:
1843	set_mpc_d(z: ac, rh: args[0], rl: args[1], ih: args[2], il: args[3]);
1844	wrapper_op_complex(ctx: &ctx, c: ac, size: 2, ieee: args);
1845	set_mpc_d(z: bc, rh: args[4], rl: args[5], ih: args[6], il: args[7]);
1846	wrapper_op_complex(ctx: &ctx, c: bc, size: 2, ieee: args+4);
1847	((testfunc2c)(fn->func))(rc, ac, bc, MPC_RNDNN);
1848	get_mpc_d(z: rc, rh: &result[0], rl: &result[1], rextra: &result[2], ih: &result[4], il: &result[5], iextra: &result[6]);
1849	wrapper_result_complex(ctx: &ctx, c: rc, size: 2, ieee: result);
1850	if (wrapper_run(ctx: &ctx, wrappers: fn->wrappers))
1851	get_mpc_d(z: rc, rh: &result[0], rl: &result[1], rextra: &result[2], ih: &result[4], il: &result[5], iextra: &result[6]);
1852	break;
1853	case args1fc:
1854	set_mpc_f(z: ac, r: args[0], i: args[2]);
1855	wrapper_op_complex(ctx: &ctx, c: ac, size: 1, ieee: args);
1856	((testfunc1c)(fn->func))(rc, ac, MPC_RNDNN);
1857	get_mpc_f(z: rc, r: &result[0], rextra: &result[1], i: &result[4], iextra: &result[5]);
1858	wrapper_result_complex(ctx: &ctx, c: rc, size: 1, ieee: result);
1859	if (wrapper_run(ctx: &ctx, wrappers: fn->wrappers))
1860	get_mpc_f(z: rc, r: &result[0], rextra: &result[1], i: &result[4], iextra: &result[5]);
1861	break;
1862	case args2fc:
1863	set_mpc_f(z: ac, r: args[0], i: args[2]);
1864	wrapper_op_complex(ctx: &ctx, c: ac, size: 1, ieee: args);
1865	set_mpc_f(z: bc, r: args[4], i: args[6]);
1866	wrapper_op_complex(ctx: &ctx, c: bc, size: 1, ieee: args+4);
1867	((testfunc2c)(fn->func))(rc, ac, bc, MPC_RNDNN);
1868	get_mpc_f(z: rc, r: &result[0], rextra: &result[1], i: &result[4], iextra: &result[5]);
1869	wrapper_result_complex(ctx: &ctx, c: rc, size: 1, ieee: result);
1870	if (wrapper_run(ctx: &ctx, wrappers: fn->wrappers))
1871	get_mpc_f(z: rc, r: &result[0], rextra: &result[1], i: &result[4], iextra: &result[5]);
1872	break;
1873	default:
1874	fprintf(stderr, format: "internal inconsistency?!\n");
1875	abort();
1876	}
1877	}
1878
1879	switch (fn->type) {
1880	case args1: /* return an extra-precise result */
1881	case args2:
1882	case args1cr:
1883	case rred:
1884	printextra = 1;
1885	if (rejected == 0) {
1886	errstr = NULL;
1887	if (!mpfr_zero_p(a)) {
1888	if ((result[0] & 0x7FFFFFFF) == 0 && result[1] == 0) {
1889	/*
1890	* If the output is +0 or -0 apart from the extra
1891	* precision in result[2], then there's a tricky
1892	* judgment call about what we require in the
1893	* output. If we output the extra bits and set
1894	* errstr="?underflow" then mathtest will tolerate
1895	* the function under test rounding down to zero
1896	* _or_ up to the minimum denormal; whereas if we
1897	* suppress the extra bits and set
1898	* errstr="underflow", then mathtest will enforce
1899	* that the function really does underflow to zero.
1900	*
1901	* But where to draw the line? It seems clear to
1902	* me that numbers along the lines of
1903	* 00000000.00000000.7ff should be treated
1904	* similarly to 00000000.00000000.801, but on the
1905	* other hand, we must surely be prepared to
1906	* enforce a genuine underflow-to-zero in _some_
1907	* case where the true mathematical output is
1908	* nonzero but absurdly tiny.
1909	*
1910	* I think a reasonable place to draw the
1911	* distinction is at 00000000.00000000.400, i.e.
1912	* one quarter of the minimum positive denormal.
1913	* If a value less than that rounds up to the
1914	* minimum denormal, that must mean the function
1915	* under test has managed to make an error of an
1916	* entire factor of two, and that's something we
1917	* should fix. Above that, you can misround within
1918	* the limits of your accuracy bound if you have
1919	* to.
1920	*/
1921	if (result[2] < 0x40000000) {
1922	/* Total underflow (ERANGE + UFL) is required,
1923	* and we suppress the extra bits to make
1924	* mathtest enforce that the output is really
1925	* zero. */
1926	errstr = "underflow";
1927	printextra = 0;
1928	} else {
1929	/* Total underflow is not required, but if the
1930	* function rounds down to zero anyway, then
1931	* we should be prepared to tolerate it. */
1932	errstr = "?underflow";
1933	}
1934	} else if (!(result[0] & 0x7ff00000)) {
1935	/*
1936	* If the output is denormal, we usually expect a
1937	* UFL exception, warning the user of partial
1938	* underflow. The exception is if the denormal
1939	* being returned is just one of the input values,
1940	* unchanged even in principle. I bodgily handle
1941	* this by just special-casing the functions in
1942	* question below.
1943	*/
1944	if (!strcmp(s1: fn->name, s2: "fmax") ||
1945	!strcmp(s1: fn->name, s2: "fmin") ||
1946	!strcmp(s1: fn->name, s2: "creal") ||
1947	!strcmp(s1: fn->name, s2: "cimag")) {
1948	/* no error expected */
1949	} else {
1950	errstr = "u";
1951	}
1952	} else if ((result[0] & 0x7FFFFFFF) > 0x7FEFFFFF) {
1953	/*
1954	* Infinite results are usually due to overflow,
1955	* but one exception is lgamma of a negative
1956	* integer.
1957	*/
1958	if (!strcmp(s1: fn->name, s2: "lgamma") &&
1959	(args[0] & 0x80000000) != 0 && /* negative */
1960	is_dinteger(in: args)) {
1961	errstr = "ERANGE status=z";
1962	} else {
1963	errstr = "overflow";
1964	}
1965	printextra = 0;
1966	}
1967	} else {
1968	/* lgamma(0) is also a pole. */
1969	if (!strcmp(s1: fn->name, s2: "lgamma")) {
1970	errstr = "ERANGE status=z";
1971	printextra = 0;
1972	}
1973	}
1974	}
1975
1976	if (!printextra || (rejected && !(rejected==1 && result[2]!=0))) {
1977	printf(format: " result=%08x.%08x",
1978	result[0], result[1]);
1979	} else {
1980	printf(format: " result=%08x.%08x.%03x",
1981	result[0], result[1], (result[2] >> 20) & 0xFFF);
1982	}
1983	if (fn->type == rred) {
1984	printf(format: " res2=%08x", result[3]);
1985	}
1986	break;
1987	case args1f:
1988	case args2f:
1989	case args1fcr:
1990	case rredf:
1991	printextra = 1;
1992	if (rejected == 0) {
1993	errstr = NULL;
1994	if (!mpfr_zero_p(a)) {
1995	if ((result[0] & 0x7FFFFFFF) == 0) {
1996	/*
1997	* Decide whether to print the extra bits based on
1998	* just how close to zero the number is. See the
1999	* big comment in the double-precision case for
2000	* discussion.
2001	*/
2002	if (result[1] < 0x40000000) {
2003	errstr = "underflow";
2004	printextra = 0;
2005	} else {
2006	errstr = "?underflow";
2007	}
2008	} else if (!(result[0] & 0x7f800000)) {
2009	/*
2010	* Functions which do not report partial overflow
2011	* are listed here as special cases. (See the
2012	* corresponding double case above for a fuller
2013	* comment.)
2014	*/
2015	if (!strcmp(s1: fn->name, s2: "fmaxf") ||
2016	!strcmp(s1: fn->name, s2: "fminf") ||
2017	!strcmp(s1: fn->name, s2: "crealf") ||
2018	!strcmp(s1: fn->name, s2: "cimagf")) {
2019	/* no error expected */
2020	} else {
2021	errstr = "u";
2022	}
2023	} else if ((result[0] & 0x7FFFFFFF) > 0x7F7FFFFF) {
2024	/*
2025	* Infinite results are usually due to overflow,
2026	* but one exception is lgamma of a negative
2027	* integer.
2028	*/
2029	if (!strcmp(s1: fn->name, s2: "lgammaf") &&
2030	(args[0] & 0x80000000) != 0 && /* negative */
2031	is_sinteger(in: args)) {
2032	errstr = "ERANGE status=z";
2033	} else {
2034	errstr = "overflow";
2035	}
2036	printextra = 0;
2037	}
2038	} else {
2039	/* lgamma(0) is also a pole. */
2040	if (!strcmp(s1: fn->name, s2: "lgammaf")) {
2041	errstr = "ERANGE status=z";
2042	printextra = 0;
2043	}
2044	}
2045	}
2046
2047	if (!printextra || (rejected && !(rejected==1 && result[1]!=0))) {
2048	printf(format: " result=%08x",
2049	result[0]);
2050	} else {
2051	printf(format: " result=%08x.%03x",
2052	result[0], (result[1] >> 20) & 0xFFF);
2053	}
2054	if (fn->type == rredf) {
2055	printf(format: " res2=%08x", result[3]);
2056	}
2057	break;
2058	case semi1: /* return a double result */
2059	case semi2:
2060	case t_ldexp:
2061	printf(format: " result=%08x.%08x", result[0], result[1]);
2062	break;
2063	case semi1f:
2064	case semi2f:
2065	case t_ldexpf:
2066	printf(format: " result=%08x", result[0]);
2067	break;
2068	case t_frexp: /* return double * int */
2069	printf(format: " result=%08x.%08x res2=%08x", result[0], result[1],
2070	result[2]);
2071	break;
2072	case t_modf: /* return double * double */
2073	printf(format: " result=%08x.%08x res2=%08x.%08x",
2074	result[0], result[1], result[2], result[3]);
2075	break;
2076	case t_modff: /* return float * float */
2077	/* fall through */
2078	case t_frexpf: /* return float * int */
2079	printf(format: " result=%08x res2=%08x", result[0], result[2]);
2080	break;
2081	case classify:
2082	case classifyf:
2083	case compare:
2084	case comparef:
2085	printf(format: " result=%x", result[0]);
2086	break;
2087	case args1c:
2088	case args2c:
2089	if (0/* errstr */) {
2090	printf(format: " resultr=%08x.%08x", result[0], result[1]);
2091	printf(format: " resulti=%08x.%08x", result[4], result[5]);
2092	} else {
2093	printf(format: " resultr=%08x.%08x.%03x",
2094	result[0], result[1], (result[2] >> 20) & 0xFFF);
2095	printf(format: " resulti=%08x.%08x.%03x",
2096	result[4], result[5], (result[6] >> 20) & 0xFFF);
2097	}
2098	/* Underflow behaviour doesn't seem to be specified for complex arithmetic */
2099	errstr = "?underflow";
2100	break;
2101	case args1fc:
2102	case args2fc:
2103	if (0/* errstr */) {
2104	printf(format: " resultr=%08x", result[0]);
2105	printf(format: " resulti=%08x", result[4]);
2106	} else {
2107	printf(format: " resultr=%08x.%03x",
2108	result[0], (result[1] >> 20) & 0xFFF);
2109	printf(format: " resulti=%08x.%03x",
2110	result[4], (result[5] >> 20) & 0xFFF);
2111	}
2112	/* Underflow behaviour doesn't seem to be specified for complex arithmetic */
2113	errstr = "?underflow";
2114	break;
2115	}
2116
2117	if (errstr && *(errstr+1) == '\0') {
2118	printf(format: " errno=0 status=%c",*errstr);
2119	} else if (errstr && *errstr == '?') {
2120	printf(format: " maybeerror=%s", errstr+1);
2121	} else if (errstr && errstr[0] == 'E') {
2122	printf(format: " errno=%s", errstr);
2123	} else {
2124	printf(format: " error=%s", errstr && *errstr ? errstr : "0");
2125	}
2126
2127	printf(format: "\n");
2128
2129	vet_for_decline(fn, args, result, got_errno_in: 0);
2130
2131	cleanup:
2132	mpfr_clear(a);
2133	mpfr_clear(b);
2134	mpfr_clear(r);
2135	mpc_clear(ac);
2136	mpc_clear(bc);
2137	mpc_clear(rc);
2138	}
2139
2140	void gencases(Testable *fn, int number) {
2141	int i;
2142	uint32 args[8];
2143
2144	float32_case(NULL);
2145	float64_case(NULL);
2146
2147	printf(format: "random=on\n"); /* signal to runtests.pl that the following tests are randomly generated */
2148	for (i = 0; i < number; i++) {
2149	/* generate test point */
2150	fn->cases(args, fn->caseparam1, fn->caseparam2);
2151	docase(fn, args);
2152	}
2153	printf(format: "random=off\n");
2154	}
2155
2156	static uint32 doubletop(int x, int scale) {
2157	int e = 0x412 + scale;
2158	while (!(x & 0x100000))
2159	x <<= 1, e--;
2160	return (e << 20) + x;
2161	}
2162
2163	static uint32 floatval(int x, int scale) {
2164	int e = 0x95 + scale;
2165	while (!(x & 0x800000))
2166	x <<= 1, e--;
2167	return (e << 23) + x;
2168	}
2169