1	/*
2	* mathtest.c - test rig for mathlib
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 <assert.h>
10	#include <stdio.h>
11	#include <stdlib.h>
12	#include <string.h>
13	#include <setjmp.h>
14	#include <ctype.h>
15	#include <math.h>
16	#include <errno.h>
17	#include <limits.h>
18	#include <fenv.h>
19	#include "mathlib.h"
20
21	#ifndef math_errhandling
22	# define math_errhandling 0
23	#endif
24
25	#ifdef __cplusplus
26	#define EXTERN_C extern "C"
27	#else
28	#define EXTERN_C extern
29	#endif
30
31	#ifndef TRUE
32	#define TRUE 1
33	#endif
34	#ifndef FALSE
35	#define FALSE 0
36	#endif
37
38	#ifdef IMPORT_SYMBOL
39	#define STR2(x) #x
40	#define STR(x) STR2(x)
41	_Pragma(STR(import IMPORT_SYMBOL))
42	#endif
43
44	int dmsd, dlsd;
45	int quiet = 0;
46	int doround = 0;
47	unsigned statusmask = FE_ALL_EXCEPT;
48
49	#define EXTRABITS (12)
50	#define ULPUNIT (1<<EXTRABITS)
51
52	typedef int (*test) (void);
53
54	/*
55	struct to hold info about a function (which could actually be a macro)
56	*/
57	typedef struct {
58	enum {
59	t_func, t_macro
60	} type;
61	enum {
62	at_d, at_s, /* double or single precision float */
63	at_d2, at_s2, /* same, but taking two args */
64	at_di, at_si, /* double/single and an int */
65	at_dip, at_sip, /* double/single and an int ptr */
66	at_ddp, at_ssp, /* d/s and a d/s ptr */
67	at_dc, at_sc, /* double or single precision complex */
68	at_dc2, at_sc2 /* same, but taking two args */
69	} argtype;
70	enum {
71	rt_d, rt_s, rt_i, /* double, single, int */
72	rt_dc, rt_sc, /* double, single precision complex */
73	rt_d2, rt_s2 /* also use res2 */
74	} rettype;
75	union {
76	void* ptr;
77	double (*d_d_ptr)(double);
78	float (*s_s_ptr)(float);
79	int (*d_i_ptr)(double);
80	int (*s_i_ptr)(float);
81	double (*d2_d_ptr)(double, double);
82	float (*s2_s_ptr)(float, float);
83	double (*di_d_ptr)(double,int);
84	float (*si_s_ptr)(float,int);
85	double (*dip_d_ptr)(double,int*);
86	float (*sip_s_ptr)(float,int*);
87	double (*ddp_d_ptr)(double,double*);
88	float (*ssp_s_ptr)(float,float*);
89	} func;
90	enum {
91	m_none,
92	m_isfinite, m_isfinitef,
93	m_isgreater, m_isgreaterequal,
94	m_isgreaterequalf, m_isgreaterf,
95	m_isinf, m_isinff,
96	m_isless, m_islessequal,
97	m_islessequalf, m_islessf,
98	m_islessgreater, m_islessgreaterf,
99	m_isnan, m_isnanf,
100	m_isnormal, m_isnormalf,
101	m_isunordered, m_isunorderedf,
102	m_fpclassify, m_fpclassifyf,
103	m_signbit, m_signbitf,
104	/* not actually a macro, but makes things easier */
105	m_rred, m_rredf,
106	m_cadd, m_csub, m_cmul, m_cdiv,
107	m_caddf, m_csubf, m_cmulf, m_cdivf
108	} macro_name; /* only used if a macro/something that can't be done using func */
109	long long tolerance;
110	const char* name;
111	} test_func;
112
113	/* used in qsort */
114	int compare_tfuncs(const void* a, const void* b) {
115	return strcmp(s1: ((test_func*)a)->name, s2: ((test_func*)b)->name);
116	}
117
118	int is_double_argtype(int argtype) {
119	switch(argtype) {
120	case at_d:
121	case at_d2:
122	case at_dc:
123	case at_dc2:
124	return 1;
125	default:
126	return 0;
127	}
128	}
129
130	int is_single_argtype(int argtype) {
131	switch(argtype) {
132	case at_s:
133	case at_s2:
134	case at_sc:
135	case at_sc2:
136	return 1;
137	default:
138	return 0;
139	}
140	}
141
142	int is_double_rettype(int rettype) {
143	switch(rettype) {
144	case rt_d:
145	case rt_dc:
146	case rt_d2:
147	return 1;
148	default:
149	return 0;
150	}
151	}
152
153	int is_single_rettype(int rettype) {
154	switch(rettype) {
155	case rt_s:
156	case rt_sc:
157	case rt_s2:
158	return 1;
159	default:
160	return 0;
161	}
162	}
163
164	int is_complex_argtype(int argtype) {
165	switch(argtype) {
166	case at_dc:
167	case at_sc:
168	case at_dc2:
169	case at_sc2:
170	return 1;
171	default:
172	return 0;
173	}
174	}
175
176	int is_complex_rettype(int rettype) {
177	switch(rettype) {
178	case rt_dc:
179	case rt_sc:
180	return 1;
181	default:
182	return 0;
183	}
184	}
185
186	/*
187	* Special-case flags indicating that some functions' error
188	* tolerance handling is more complicated than a fixed relative
189	* error bound.
190	*/
191	#define ABSLOWERBOUND 0x4000000000000000LL
192	#define PLUSMINUSPIO2 0x1000000000000000LL
193
194	#define ARM_PREFIX(x) x
195
196	#define TFUNC(arg,ret,name,tolerance) { t_func, arg, ret, (void*)&name, m_none, tolerance, #name }
197	#define TFUNCARM(arg,ret,name,tolerance) { t_func, arg, ret, (void*)& ARM_PREFIX(name), m_none, tolerance, #name }
198	#define MFUNC(arg,ret,name,tolerance) { t_macro, arg, ret, NULL, m_##name, tolerance, #name }
199
200	/* sincosf wrappers for easier testing. */
201	static float sincosf_sinf(float x) { float s,c; sincosf(x: x, sinx: &s, cosx: &c); return s; }
202	static float sincosf_cosf(float x) { float s,c; sincosf(x: x, sinx: &s, cosx: &c); return c; }
203
204	test_func tfuncs[] = {
205	/* trigonometric */
206	TFUNC(at_d,rt_d, acos, 4*ULPUNIT),
207	TFUNC(at_d,rt_d, asin, 4*ULPUNIT),
208	TFUNC(at_d,rt_d, atan, 4*ULPUNIT),
209	TFUNC(at_d2,rt_d, atan2, 4*ULPUNIT),
210
211	TFUNC(at_d,rt_d, tan, 2*ULPUNIT),
212	TFUNC(at_d,rt_d, sin, 2*ULPUNIT),
213	TFUNC(at_d,rt_d, cos, 2*ULPUNIT),
214
215	TFUNC(at_s,rt_s, acosf, 4*ULPUNIT),
216	TFUNC(at_s,rt_s, asinf, 4*ULPUNIT),
217	TFUNC(at_s,rt_s, atanf, 4*ULPUNIT),
218	TFUNC(at_s2,rt_s, atan2f, 4*ULPUNIT),
219	TFUNCARM(at_s,rt_s, tanf, 4*ULPUNIT),
220	TFUNCARM(at_s,rt_s, sinf, 3*ULPUNIT/4),
221	TFUNCARM(at_s,rt_s, cosf, 3*ULPUNIT/4),
222	TFUNCARM(at_s,rt_s, sincosf_sinf, 3*ULPUNIT/4),
223	TFUNCARM(at_s,rt_s, sincosf_cosf, 3*ULPUNIT/4),
224
225	/* hyperbolic */
226	TFUNC(at_d, rt_d, atanh, 4*ULPUNIT),
227	TFUNC(at_d, rt_d, asinh, 4*ULPUNIT),
228	TFUNC(at_d, rt_d, acosh, 4*ULPUNIT),
229	TFUNC(at_d,rt_d, tanh, 4*ULPUNIT),
230	TFUNC(at_d,rt_d, sinh, 4*ULPUNIT),
231	TFUNC(at_d,rt_d, cosh, 4*ULPUNIT),
232
233	TFUNC(at_s, rt_s, atanhf, 4*ULPUNIT),
234	TFUNC(at_s, rt_s, asinhf, 4*ULPUNIT),
235	TFUNC(at_s, rt_s, acoshf, 4*ULPUNIT),
236	TFUNC(at_s,rt_s, tanhf, 4*ULPUNIT),
237	TFUNC(at_s,rt_s, sinhf, 4*ULPUNIT),
238	TFUNC(at_s,rt_s, coshf, 4*ULPUNIT),
239
240	/* exponential and logarithmic */
241	TFUNC(at_d,rt_d, log, 3*ULPUNIT/4),
242	TFUNC(at_d,rt_d, log10, 3*ULPUNIT),
243	TFUNC(at_d,rt_d, log2, 3*ULPUNIT/4),
244	TFUNC(at_d,rt_d, log1p, 2*ULPUNIT),
245	TFUNC(at_d,rt_d, exp, 3*ULPUNIT/4),
246	TFUNC(at_d,rt_d, exp2, 3*ULPUNIT/4),
247	TFUNC(at_d,rt_d, expm1, ULPUNIT),
248	TFUNCARM(at_s,rt_s, logf, ULPUNIT),
249	TFUNC(at_s,rt_s, log10f, 3*ULPUNIT),
250	TFUNCARM(at_s,rt_s, log2f, ULPUNIT),
251	TFUNC(at_s,rt_s, log1pf, 2*ULPUNIT),
252	TFUNCARM(at_s,rt_s, expf, 3*ULPUNIT/4),
253	TFUNCARM(at_s,rt_s, exp2f, 3*ULPUNIT/4),
254	TFUNC(at_s,rt_s, expm1f, ULPUNIT),
255
256	/* power */
257	TFUNC(at_d2,rt_d, pow, 3*ULPUNIT/4),
258	TFUNC(at_d,rt_d, sqrt, ULPUNIT/2),
259	TFUNC(at_d,rt_d, cbrt, 2*ULPUNIT),
260	TFUNC(at_d2, rt_d, hypot, 4*ULPUNIT),
261
262	TFUNCARM(at_s2,rt_s, powf, ULPUNIT),
263	TFUNC(at_s,rt_s, sqrtf, ULPUNIT/2),
264	TFUNC(at_s,rt_s, cbrtf, 2*ULPUNIT),
265	TFUNC(at_s2, rt_s, hypotf, 4*ULPUNIT),
266
267	/* error function */
268	TFUNC(at_d,rt_d, erf, 16*ULPUNIT),
269	TFUNC(at_s,rt_s, erff, 16*ULPUNIT),
270	TFUNC(at_d,rt_d, erfc, 16*ULPUNIT),
271	TFUNC(at_s,rt_s, erfcf, 16*ULPUNIT),
272
273	/* gamma functions */
274	TFUNC(at_d,rt_d, tgamma, 16*ULPUNIT),
275	TFUNC(at_s,rt_s, tgammaf, 16*ULPUNIT),
276	TFUNC(at_d,rt_d, lgamma, 16*ULPUNIT | ABSLOWERBOUND),
277	TFUNC(at_s,rt_s, lgammaf, 16*ULPUNIT | ABSLOWERBOUND),
278
279	TFUNC(at_d,rt_d, ceil, 0),
280	TFUNC(at_s,rt_s, ceilf, 0),
281	TFUNC(at_d2,rt_d, copysign, 0),
282	TFUNC(at_s2,rt_s, copysignf, 0),
283	TFUNC(at_d,rt_d, floor, 0),
284	TFUNC(at_s,rt_s, floorf, 0),
285	TFUNC(at_d2,rt_d, fmax, 0),
286	TFUNC(at_s2,rt_s, fmaxf, 0),
287	TFUNC(at_d2,rt_d, fmin, 0),
288	TFUNC(at_s2,rt_s, fminf, 0),
289	TFUNC(at_d2,rt_d, fmod, 0),
290	TFUNC(at_s2,rt_s, fmodf, 0),
291	MFUNC(at_d, rt_i, fpclassify, 0),
292	MFUNC(at_s, rt_i, fpclassifyf, 0),
293	TFUNC(at_dip,rt_d, frexp, 0),
294	TFUNC(at_sip,rt_s, frexpf, 0),
295	MFUNC(at_d, rt_i, isfinite, 0),
296	MFUNC(at_s, rt_i, isfinitef, 0),
297	MFUNC(at_d, rt_i, isgreater, 0),
298	MFUNC(at_d, rt_i, isgreaterequal, 0),
299	MFUNC(at_s, rt_i, isgreaterequalf, 0),
300	MFUNC(at_s, rt_i, isgreaterf, 0),
301	MFUNC(at_d, rt_i, isinf, 0),
302	MFUNC(at_s, rt_i, isinff, 0),
303	MFUNC(at_d, rt_i, isless, 0),
304	MFUNC(at_d, rt_i, islessequal, 0),
305	MFUNC(at_s, rt_i, islessequalf, 0),
306	MFUNC(at_s, rt_i, islessf, 0),
307	MFUNC(at_d, rt_i, islessgreater, 0),
308	MFUNC(at_s, rt_i, islessgreaterf, 0),
309	MFUNC(at_d, rt_i, isnan, 0),
310	MFUNC(at_s, rt_i, isnanf, 0),
311	MFUNC(at_d, rt_i, isnormal, 0),
312	MFUNC(at_s, rt_i, isnormalf, 0),
313	MFUNC(at_d, rt_i, isunordered, 0),
314	MFUNC(at_s, rt_i, isunorderedf, 0),
315	TFUNC(at_di,rt_d, ldexp, 0),
316	TFUNC(at_si,rt_s, ldexpf, 0),
317	TFUNC(at_ddp,rt_d2, modf, 0),
318	TFUNC(at_ssp,rt_s2, modff, 0),
319	#ifndef BIGRANGERED
320	MFUNC(at_d, rt_d, rred, 2*ULPUNIT),
321	#else
322	MFUNC(at_d, rt_d, m_rred, ULPUNIT),
323	#endif
324	MFUNC(at_d, rt_i, signbit, 0),
325	MFUNC(at_s, rt_i, signbitf, 0),
326	};
327
328	/*
329	* keywords are: func size op1 op2 result res2 errno op1r op1i op2r op2i resultr resulti
330	* also we ignore: wrongresult wrongres2 wrongerrno
331	* op1 equivalent to op1r, same with op2 and result
332	*/
333
334	typedef struct {
335	test_func *func;
336	unsigned op1r[2]; /* real part, also used for non-complex numbers */
337	unsigned op1i[2]; /* imaginary part */
338	unsigned op2r[2];
339	unsigned op2i[2];
340	unsigned resultr[3];
341	unsigned resulti[3];
342	enum {
343	rc_none, rc_zero, rc_infinity, rc_nan, rc_finite
344	} resultc; /* special complex results, rc_none means use resultr and resulti as normal */
345	unsigned res2[2];
346	unsigned status; /* IEEE status return, if any */
347	unsigned maybestatus; /* for optional status, or allowance for spurious */
348	int nresult; /* number of result words */
349	int in_err, in_err_limit;
350	int err;
351	int maybeerr;
352	int valid;
353	int comment;
354	int random;
355	} testdetail;
356
357	enum { /* keywords */
358	k_errno, k_errno_in, k_error, k_func, k_maybeerror, k_maybestatus, k_op1, k_op1i, k_op1r, k_op2, k_op2i, k_op2r,
359	k_random, k_res2, k_result, k_resultc, k_resulti, k_resultr, k_status,
360	k_wrongres2, k_wrongresult, k_wrongstatus, k_wrongerrno
361	};
362	char *keywords[] = {
363	"errno", "errno_in", "error", "func", "maybeerror", "maybestatus", "op1", "op1i", "op1r", "op2", "op2i", "op2r",
364	"random", "res2", "result", "resultc", "resulti", "resultr", "status",
365	"wrongres2", "wrongresult", "wrongstatus", "wrongerrno"
366	};
367
368	enum {
369	e_0, e_EDOM, e_ERANGE,
370
371	/*
372	* This enum makes sure that we have the right number of errnos in the
373	* errno[] array
374	*/
375	e_number_of_errnos
376	};
377	char *errnos[] = {
378	"0", "EDOM", "ERANGE"
379	};
380
381	enum {
382	e_none, e_divbyzero, e_domain, e_overflow, e_underflow
383	};
384	char *errors[] = {
385	"0", "divbyzero", "domain", "overflow", "underflow"
386	};
387
388	static int verbose, fo, strict;
389
390	/* state toggled by random=on / random=off */
391	static int randomstate;
392
393	/* Canonify a double NaN: SNaNs all become 7FF00000.00000001 and QNaNs
394	* all become 7FF80000.00000001 */
395	void canon_dNaN(unsigned a[2]) {
396	if ((a[0] & 0x7FF00000) != 0x7FF00000)
397	return; /* not Inf or NaN */
398	if (!(a[0] & 0xFFFFF) && !a[1])
399	return; /* Inf */
400	a[0] &= 0x7FF80000; /* canonify top word */
401	a[1] = 0x00000001; /* canonify bottom word */
402	}
403
404	/* Canonify a single NaN: SNaNs all become 7F800001 and QNaNs
405	* all become 7FC00001. Returns classification of the NaN. */
406	void canon_sNaN(unsigned a[1]) {
407	if ((a[0] & 0x7F800000) != 0x7F800000)
408	return; /* not Inf or NaN */
409	if (!(a[0] & 0x7FFFFF))
410	return; /* Inf */
411	a[0] &= 0x7FC00000; /* canonify most bits */
412	a[0] |= 0x00000001; /* canonify bottom bit */
413	}
414
415	/*
416	* Detect difficult operands for FO mode.
417	*/
418	int is_dhard(unsigned a[2])
419	{
420	if ((a[0] & 0x7FF00000) == 0x7FF00000)
421	return TRUE; /* inf or NaN */
422	if ((a[0] & 0x7FF00000) == 0 &&
423	((a[0] & 0x7FFFFFFF) | a[1]) != 0)
424	return TRUE; /* denormal */
425	return FALSE;
426	}
427	int is_shard(unsigned a[1])
428	{
429	if ((a[0] & 0x7F800000) == 0x7F800000)
430	return TRUE; /* inf or NaN */
431	if ((a[0] & 0x7F800000) == 0 &&
432	(a[0] & 0x7FFFFFFF) != 0)
433	return TRUE; /* denormal */
434	return FALSE;
435	}
436
437	/*
438	* Normalise all zeroes into +0, for FO mode.
439	*/
440	void dnormzero(unsigned a[2])
441	{
442	if (a[0] == 0x80000000 && a[1] == 0)
443	a[0] = 0;
444	}
445	void snormzero(unsigned a[1])
446	{
447	if (a[0] == 0x80000000)
448	a[0] = 0;
449	}
450
451	static int find(char *word, char **array, int asize) {
452	int i, j;
453
454	asize /= sizeof(char *);
455
456	i = -1; j = asize; /* strictly between i and j */
457	while (j-i > 1) {
458	int k = (i+j) / 2;
459	int c = strcmp(s1: word, s2: array[k]);
460	if (c > 0)
461	i = k;
462	else if (c < 0)
463	j = k;
464	else /* found it! */
465	return k;
466	}
467	return -1; /* not found */
468	}
469
470	static test_func* find_testfunc(char *word) {
471	int i, j, asize;
472
473	asize = sizeof(tfuncs)/sizeof(test_func);
474
475	i = -1; j = asize; /* strictly between i and j */
476	while (j-i > 1) {
477	int k = (i+j) / 2;
478	int c = strcmp(s1: word, s2: tfuncs[k].name);
479	if (c > 0)
480	i = k;
481	else if (c < 0)
482	j = k;
483	else /* found it! */
484	return tfuncs + k;
485	}
486	return NULL; /* not found */
487	}
488
489	static long long calc_error(unsigned a[2], unsigned b[3], int shift, int rettype) {
490	unsigned r0, r1, r2;
491	int sign, carry;
492	long long result;
493
494	/*
495	* If either number is infinite, require exact equality. If
496	* either number is NaN, require that both are NaN. If either
497	* of these requirements is broken, return INT_MAX.
498	*/
499	if (is_double_rettype(rettype)) {
500	if ((a[0] & 0x7FF00000) == 0x7FF00000 ||
501	(b[0] & 0x7FF00000) == 0x7FF00000) {
502	if (((a[0] & 0x800FFFFF) || a[1]) &&
503	((b[0] & 0x800FFFFF) || b[1]) &&
504	(a[0] & 0x7FF00000) == 0x7FF00000 &&
505	(b[0] & 0x7FF00000) == 0x7FF00000)
506	return 0; /* both NaN - OK */
507	if (!((a[0] & 0xFFFFF) || a[1]) &&
508	!((b[0] & 0xFFFFF) || b[1]) &&
509	a[0] == b[0])
510	return 0; /* both same sign of Inf - OK */
511	return LLONG_MAX;
512	}
513	} else {
514	if ((a[0] & 0x7F800000) == 0x7F800000 ||
515	(b[0] & 0x7F800000) == 0x7F800000) {
516	if ((a[0] & 0x807FFFFF) &&
517	(b[0] & 0x807FFFFF) &&
518	(a[0] & 0x7F800000) == 0x7F800000 &&
519	(b[0] & 0x7F800000) == 0x7F800000)
520	return 0; /* both NaN - OK */
521	if (!(a[0] & 0x7FFFFF) &&
522	!(b[0] & 0x7FFFFF) &&
523	a[0] == b[0])
524	return 0; /* both same sign of Inf - OK */
525	return LLONG_MAX;
526	}
527	}
528
529	/*
530	* Both finite. Return INT_MAX if the signs differ.
531	*/
532	if ((a[0] ^ b[0]) & 0x80000000)
533	return LLONG_MAX;
534
535	/*
536	* Now it's just straight multiple-word subtraction.
537	*/
538	if (is_double_rettype(rettype)) {
539	r2 = -b[2]; carry = (r2 == 0);
540	r1 = a[1] + ~b[1] + carry; carry = (r1 < a[1] || (carry && r1 == a[1]));
541	r0 = a[0] + ~b[0] + carry;
542	} else {
543	r2 = -b[1]; carry = (r2 == 0);
544	r1 = a[0] + ~b[0] + carry; carry = (r1 < a[0] || (carry && r1 == a[0]));
545	r0 = ~0 + carry;
546	}
547
548	/*
549	* Forgive larger errors in specialised cases.
550	*/
551	if (shift > 0) {
552	if (shift > 32*3)
553	return 0; /* all errors are forgiven! */
554	while (shift >= 32) {
555	r2 = r1;
556	r1 = r0;
557	r0 = -(r0 >> 31);
558	shift -= 32;
559	}
560
561	if (shift > 0) {
562	r2 = (r2 >> shift) | (r1 << (32-shift));
563	r1 = (r1 >> shift) | (r0 << (32-shift));
564	r0 = (r0 >> shift) | ((-(r0 >> 31)) << (32-shift));
565	}
566	}
567
568	if (r0 & 0x80000000) {
569	sign = 1;
570	r2 = ~r2; carry = (r2 == 0);
571	r1 = 0 + ~r1 + carry; carry = (carry && (r2 == 0));
572	r0 = 0 + ~r0 + carry;
573	} else {
574	sign = 0;
575	}
576
577	if (r0 >= (1LL<<(31-EXTRABITS)))
578	return LLONG_MAX; /* many ulps out */
579
580	result = (r2 >> (32-EXTRABITS)) & (ULPUNIT-1);
581	result |= r1 << EXTRABITS;
582	result |= (long long)r0 << (32+EXTRABITS);
583	if (sign)
584	result = -result;
585	return result;
586	}
587
588	/* special named operands */
589
590	typedef struct {
591	unsigned op1, op2;
592	char* name;
593	} special_op;
594
595	static special_op special_ops_double[] = {
596	{0x00000000,0x00000000,"0"},
597	{0x3FF00000,0x00000000,"1"},
598	{0x7FF00000,0x00000000,"inf"},
599	{0x7FF80000,0x00000001,"qnan"},
600	{0x7FF00000,0x00000001,"snan"},
601	{0x3ff921fb,0x54442d18,"pi2"},
602	{0x400921fb,0x54442d18,"pi"},
603	{0x3fe921fb,0x54442d18,"pi4"},
604	{0x4002d97c,0x7f3321d2,"3pi4"},
605	};
606
607	static special_op special_ops_float[] = {
608	{0x00000000,0,"0"},
609	{0x3f800000,0,"1"},
610	{0x7f800000,0,"inf"},
611	{0x7fc00000,0,"qnan"},
612	{0x7f800001,0,"snan"},
613	{0x3fc90fdb,0,"pi2"},
614	{0x40490fdb,0,"pi"},
615	{0x3f490fdb,0,"pi4"},
616	{0x4016cbe4,0,"3pi4"},
617	};
618
619	/*
620	This is what is returned by the below functions.
621	We need it to handle the sign of the number
622	*/
623	static special_op tmp_op = {0,0,0};
624
625	special_op* find_special_op_from_op(unsigned op1, unsigned op2, int is_double) {
626	int i;
627	special_op* sop;
628	if(is_double) {
629	sop = special_ops_double;
630	} else {
631	sop = special_ops_float;
632	}
633	for(i = 0; i < sizeof(special_ops_double)/sizeof(special_op); i++) {
634	if(sop->op1 == (op1&0x7fffffff) && sop->op2 == op2) {
635	if(tmp_op.name) free(ptr: tmp_op.name);
636	tmp_op.name = malloc(size: strlen(s: sop->name)+2);
637	if(op1>>31) {
638	sprintf(s: tmp_op.name,format: "-%s",sop->name);
639	} else {
640	strcpy(dest: tmp_op.name,src: sop->name);
641	}
642	return &tmp_op;
643	}
644	sop++;
645	}
646	return NULL;
647	}
648
649	special_op* find_special_op_from_name(const char* name, int is_double) {
650	int i, neg=0;
651	special_op* sop;
652	if(is_double) {
653	sop = special_ops_double;
654	} else {
655	sop = special_ops_float;
656	}
657	if(*name=='-') {
658	neg=1;
659	name++;
660	} else if(*name=='+') {
661	name++;
662	}
663	for(i = 0; i < sizeof(special_ops_double)/sizeof(special_op); i++) {
664	if(0 == strcmp(s1: name,s2: sop->name)) {
665	tmp_op.op1 = sop->op1;
666	if(neg) {
667	tmp_op.op1 |= 0x80000000;
668	}
669	tmp_op.op2 = sop->op2;
670	return &tmp_op;
671	}
672	sop++;
673	}
674	return NULL;
675	}
676
677	/*
678	helper function for the below
679	type=0 for single, 1 for double, 2 for no sop
680	*/
681	int do_op(char* q, unsigned* op, const char* name, int num, int sop_type) {
682	int i;
683	int n=num;
684	special_op* sop = NULL;
685	for(i = 0; i < num; i++) {
686	op[i] = 0;
687	}
688	if(sop_type<2) {
689	sop = find_special_op_from_name(name: q,is_double: sop_type);
690	}
691	if(sop != NULL) {
692	op[0] = sop->op1;
693	op[1] = sop->op2;
694	} else {
695	switch(num) {
696	case 1: n = sscanf(s: q, format: "%x", &op[0]); break;
697	case 2: n = sscanf(s: q, format: "%x.%x", &op[0], &op[1]); break;
698	case 3: n = sscanf(s: q, format: "%x.%x.%x", &op[0], &op[1], &op[2]); break;
699	default: return -1;
700	}
701	}
702	if (verbose) {
703	printf(format: "%s=",name);
704	for (i = 0; (i < n); ++i) printf(format: "%x.", op[i]);
705	printf(format: " (n=%d)\n", n);
706	}
707	return n;
708	}
709
710	testdetail parsetest(char *testbuf, testdetail oldtest) {
711	char *p; /* Current part of line: Option name */
712	char *q; /* Current part of line: Option value */
713	testdetail ret; /* What we return */
714	int k; /* Function enum from k_* */
715	int n; /* Used as returns for scanfs */
716	int argtype=2, rettype=2; /* for do_op */
717
718	/* clear ret */
719	memset(s: &ret, c: 0, n: sizeof(ret));
720
721	if (verbose) printf(format: "Parsing line: %s\n", testbuf);
722	while (*testbuf && isspace(*testbuf)) testbuf++;
723	if (testbuf[0] == ';' || testbuf[0] == '#' || testbuf[0] == '!' ||
724	testbuf[0] == '>' || testbuf[0] == '\0') {
725	ret.comment = 1;
726	if (verbose) printf(format: "Line is a comment\n");
727	return ret;
728	}
729	ret.comment = 0;
730
731	if (*testbuf == '+') {
732	if (oldtest.valid) {
733	ret = oldtest; /* structure copy */
734	} else {
735	fprintf(stderr, format: "copy from invalid: ignored\n");
736	}
737	testbuf++;
738	}
739
740	ret.random = randomstate;
741
742	ret.in_err = 0;
743	ret.in_err_limit = e_number_of_errnos;
744
745	p = strtok(s: testbuf, delim: " \t");
746	while (p != NULL) {
747	q = strchr(s: p, c: '=');
748	if (!q)
749	goto balderdash;
750	*q++ = '\0';
751	k = find(word: p, array: keywords, asize: sizeof(keywords));
752	switch (k) {
753	case k_random:
754	randomstate = (!strcmp(s1: q, s2: "on"));
755	ret.comment = 1;
756	return ret; /* otherwise ignore this line */
757	case k_func:
758	if (verbose) printf(format: "func=%s ", q);
759	//ret.func = find(q, funcs, sizeof(funcs));
760	ret.func = find_testfunc(word: q);
761	if (ret.func == NULL)
762	{
763	if (verbose) printf(format: "(id=unknown)\n");
764	goto balderdash;
765	}
766	if(is_single_argtype(argtype: ret.func->argtype))
767	argtype = 0;
768	else if(is_double_argtype(argtype: ret.func->argtype))
769	argtype = 1;
770	if(is_single_rettype(rettype: ret.func->rettype))
771	rettype = 0;
772	else if(is_double_rettype(rettype: ret.func->rettype))
773	rettype = 1;
774	//ret.size = sizes[ret.func];
775	if (verbose) printf(format: "(name=%s) (size=%d)\n", ret.func->name, ret.func->argtype);
776	break;
777	case k_op1:
778	case k_op1r:
779	n = do_op(q,op: ret.op1r,name: "op1r",num: 2,sop_type: argtype);
780	if (n < 1)
781	goto balderdash;
782	break;
783	case k_op1i:
784	n = do_op(q,op: ret.op1i,name: "op1i",num: 2,sop_type: argtype);
785	if (n < 1)
786	goto balderdash;
787	break;
788	case k_op2:
789	case k_op2r:
790	n = do_op(q,op: ret.op2r,name: "op2r",num: 2,sop_type: argtype);
791	if (n < 1)
792	goto balderdash;
793	break;
794	case k_op2i:
795	n = do_op(q,op: ret.op2i,name: "op2i",num: 2,sop_type: argtype);
796	if (n < 1)
797	goto balderdash;
798	break;
799	case k_resultc:
800	puts(s: q);
801	if(strncmp(s1: q,s2: "inf",n: 3)==0) {
802	ret.resultc = rc_infinity;
803	} else if(strcmp(s1: q,s2: "zero")==0) {
804	ret.resultc = rc_zero;
805	} else if(strcmp(s1: q,s2: "nan")==0) {
806	ret.resultc = rc_nan;
807	} else if(strcmp(s1: q,s2: "finite")==0) {
808	ret.resultc = rc_finite;
809	} else {
810	goto balderdash;
811	}
812	break;
813	case k_result:
814	case k_resultr:
815	n = (do_op)(q,op: ret.resultr,name: "resultr",num: 3,sop_type: rettype);
816	if (n < 1)
817	goto balderdash;
818	ret.nresult = n; /* assume real and imaginary have same no. words */
819	break;
820	case k_resulti:
821	n = do_op(q,op: ret.resulti,name: "resulti",num: 3,sop_type: rettype);
822	if (n < 1)
823	goto balderdash;
824	break;
825	case k_res2:
826	n = do_op(q,op: ret.res2,name: "res2",num: 2,sop_type: rettype);
827	if (n < 1)
828	goto balderdash;
829	break;
830	case k_status:
831	while (*q) {
832	if (*q == 'i') ret.status |= FE_INVALID;
833	if (*q == 'z') ret.status |= FE_DIVBYZERO;
834	if (*q == 'o') ret.status |= FE_OVERFLOW;
835	if (*q == 'u') ret.status |= FE_UNDERFLOW;
836	q++;
837	}
838	break;
839	case k_maybeerror:
840	n = find(word: q, array: errors, asize: sizeof(errors));
841	if (n < 0)
842	goto balderdash;
843	if(math_errhandling&MATH_ERREXCEPT) {
844	switch(n) {
845	case e_domain: ret.maybestatus |= FE_INVALID; break;
846	case e_divbyzero: ret.maybestatus |= FE_DIVBYZERO; break;
847	case e_overflow: ret.maybestatus |= FE_OVERFLOW; break;
848	case e_underflow: ret.maybestatus |= FE_UNDERFLOW; break;
849	}
850	}
851	{
852	switch(n) {
853	case e_domain:
854	ret.maybeerr = e_EDOM; break;
855	case e_divbyzero:
856	case e_overflow:
857	case e_underflow:
858	ret.maybeerr = e_ERANGE; break;
859	}
860	}
861	case k_maybestatus:
862	while (*q) {
863	if (*q == 'i') ret.maybestatus |= FE_INVALID;
864	if (*q == 'z') ret.maybestatus |= FE_DIVBYZERO;
865	if (*q == 'o') ret.maybestatus |= FE_OVERFLOW;
866	if (*q == 'u') ret.maybestatus |= FE_UNDERFLOW;
867	q++;
868	}
869	break;
870	case k_error:
871	n = find(word: q, array: errors, asize: sizeof(errors));
872	if (n < 0)
873	goto balderdash;
874	if(math_errhandling&MATH_ERREXCEPT) {
875	switch(n) {
876	case e_domain: ret.status |= FE_INVALID; break;
877	case e_divbyzero: ret.status |= FE_DIVBYZERO; break;
878	case e_overflow: ret.status |= FE_OVERFLOW; break;
879	case e_underflow: ret.status |= FE_UNDERFLOW; break;
880	}
881	}
882	if(math_errhandling&MATH_ERRNO) {
883	switch(n) {
884	case e_domain:
885	ret.err = e_EDOM; break;
886	case e_divbyzero:
887	case e_overflow:
888	case e_underflow:
889	ret.err = e_ERANGE; break;
890	}
891	}
892	if(!(math_errhandling&MATH_ERRNO)) {
893	switch(n) {
894	case e_domain:
895	ret.maybeerr = e_EDOM; break;
896	case e_divbyzero:
897	case e_overflow:
898	case e_underflow:
899	ret.maybeerr = e_ERANGE; break;
900	}
901	}
902	break;
903	case k_errno:
904	ret.err = find(word: q, array: errnos, asize: sizeof(errnos));
905	if (ret.err < 0)
906	goto balderdash;
907	break;
908	case k_errno_in:
909	ret.in_err = find(word: q, array: errnos, asize: sizeof(errnos));
910	if (ret.err < 0)
911	goto balderdash;
912	ret.in_err_limit = ret.in_err + 1;
913	break;
914	case k_wrongresult:
915	case k_wrongstatus:
916	case k_wrongres2:
917	case k_wrongerrno:
918	/* quietly ignore these keys */
919	break;
920	default:
921	goto balderdash;
922	}
923	p = strtok(NULL, delim: " \t");
924	}
925	ret.valid = 1;
926	return ret;
927
928	/* come here from almost any error */
929	balderdash:
930	ret.valid = 0;
931	return ret;
932	}
933
934	typedef enum {
935	test_comment, /* deliberately not a test */
936	test_invalid, /* accidentally not a test */
937	test_decline, /* was a test, and wasn't run */
938	test_fail, /* was a test, and failed */
939	test_pass /* was a test, and passed */
940	} testresult;
941
942	char failtext[512];
943
944	typedef union {
945	unsigned i[2];
946	double f;
947	double da[2];
948	} dbl;
949
950	typedef union {
951	unsigned i;
952	float f;
953	float da[2];
954	} sgl;
955
956	/* helper function for runtest */
957	void print_error(int rettype, unsigned *result, char* text, char** failp) {
958	special_op *sop;
959	char *str;
960
961	if(result) {
962	*failp += sprintf(s: *failp,format: " %s=",text);
963	sop = find_special_op_from_op(op1: result[0],op2: result[1],is_double: is_double_rettype(rettype));
964	if(sop) {
965	*failp += sprintf(s: *failp,format: "%s",sop->name);
966	} else {
967	if(is_double_rettype(rettype)) {
968	str="%08x.%08x";
969	} else {
970	str="%08x";
971	}
972	*failp += sprintf(s: *failp,format: str,result[0],result[1]);
973	}
974	}
975	}
976
977
978	void print_ulps_helper(const char *name, long long ulps, char** failp) {
979	if(ulps == LLONG_MAX) {
980	*failp += sprintf(s: *failp, format: " %s=HUGE", name);
981	} else {
982	*failp += sprintf(s: *failp, format: " %s=%.3f", name, (double)ulps / ULPUNIT);
983	}
984	}
985
986	/* for complex args make ulpsr or ulpsri = 0 to not print */
987	void print_ulps(int rettype, long long ulpsr, long long ulpsi, char** failp) {
988	if(is_complex_rettype(rettype)) {
989	if (ulpsr) print_ulps_helper(name: "ulpsr",ulps: ulpsr,failp);
990	if (ulpsi) print_ulps_helper(name: "ulpsi",ulps: ulpsi,failp);
991	} else {
992	if (ulpsr) print_ulps_helper(name: "ulps",ulps: ulpsr,failp);
993	}
994	}
995
996	int runtest(testdetail t) {
997	int err, status;
998
999	dbl d_arg1, d_arg2, d_res, d_res2;
1000	sgl s_arg1, s_arg2, s_res, s_res2;
1001
1002	int deferred_decline = FALSE;
1003	char *failp = failtext;
1004
1005	unsigned int intres=0;
1006
1007	int res2_adjust = 0;
1008
1009	if (t.comment)
1010	return test_comment;
1011	if (!t.valid)
1012	return test_invalid;
1013
1014	/* Set IEEE status to mathlib-normal */
1015	feclearexcept(FE_ALL_EXCEPT);
1016
1017	/* Deal with operands */
1018	#define DO_DOP(arg,op) arg.i[dmsd] = t.op[0]; arg.i[dlsd] = t.op[1]
1019	DO_DOP(d_arg1,op1r);
1020	DO_DOP(d_arg2,op2r);
1021	s_arg1.i = t.op1r[0]; s_arg2.i = t.op2r[0];
1022
1023	/*
1024	* Detect NaNs, infinities and denormals on input, and set a
1025	* deferred decline flag if we're in FO mode.
1026	*
1027	* (We defer the decline rather than doing it immediately
1028	* because even in FO mode the operation is not permitted to
1029	* crash or tight-loop; so we _run_ the test, and then ignore
1030	* all the results.)
1031	*/
1032	if (fo) {
1033	if (is_double_argtype(argtype: t.func->argtype) && is_dhard(a: t.op1r))
1034	deferred_decline = TRUE;
1035	if (t.func->argtype==at_d2 && is_dhard(a: t.op2r))
1036	deferred_decline = TRUE;
1037	if (is_single_argtype(argtype: t.func->argtype) && is_shard(a: t.op1r))
1038	deferred_decline = TRUE;
1039	if (t.func->argtype==at_s2 && is_shard(a: t.op2r))
1040	deferred_decline = TRUE;
1041	if (is_double_rettype(rettype: t.func->rettype) && is_dhard(a: t.resultr))
1042	deferred_decline = TRUE;
1043	if (t.func->rettype==rt_d2 && is_dhard(a: t.res2))
1044	deferred_decline = TRUE;
1045	if (is_single_argtype(argtype: t.func->rettype) && is_shard(a: t.resultr))
1046	deferred_decline = TRUE;
1047	if (t.func->rettype==rt_s2 && is_shard(a: t.res2))
1048	deferred_decline = TRUE;
1049	if (t.err == e_ERANGE)
1050	deferred_decline = TRUE;
1051	}
1052
1053	/*
1054	* Perform the operation
1055	*/
1056
1057	errno = t.in_err == e_EDOM ? EDOM : t.in_err == e_ERANGE ? ERANGE : 0;
1058	if (t.err == e_0)
1059	t.err = t.in_err;
1060	if (t.maybeerr == e_0)
1061	t.maybeerr = t.in_err;
1062
1063	if(t.func->type == t_func) {
1064	switch(t.func->argtype) {
1065	case at_d: d_res.f = t.func->func.d_d_ptr(d_arg1.f); break;
1066	case at_s: s_res.f = t.func->func.s_s_ptr(s_arg1.f); break;
1067	case at_d2: d_res.f = t.func->func.d2_d_ptr(d_arg1.f, d_arg2.f); break;
1068	case at_s2: s_res.f = t.func->func.s2_s_ptr(s_arg1.f, s_arg2.f); break;
1069	case at_di: d_res.f = t.func->func.di_d_ptr(d_arg1.f, d_arg2.i[dmsd]); break;
1070	case at_si: s_res.f = t.func->func.si_s_ptr(s_arg1.f, s_arg2.i); break;
1071	case at_dip: d_res.f = t.func->func.dip_d_ptr(d_arg1.f, (int*)&intres); break;
1072	case at_sip: s_res.f = t.func->func.sip_s_ptr(s_arg1.f, (int*)&intres); break;
1073	case at_ddp: d_res.f = t.func->func.ddp_d_ptr(d_arg1.f, &d_res2.f); break;
1074	case at_ssp: s_res.f = t.func->func.ssp_s_ptr(s_arg1.f, &s_res2.f); break;
1075	default:
1076	printf(format: "unhandled function: %s\n",t.func->name);
1077	return test_fail;
1078	}
1079	} else {
1080	/* printf("macro: name=%s, num=%i, s1.i=0x%08x s1.f=%f\n",t.func->name, t.func->macro_name, s_arg1.i, (double)s_arg1.f); */
1081	switch(t.func->macro_name) {
1082	case m_isfinite: intres = isfinite(d_arg1.f); break;
1083	case m_isinf: intres = isinf(d_arg1.f); break;
1084	case m_isnan: intres = isnan(d_arg1.f); break;
1085	case m_isnormal: intres = isnormal(d_arg1.f); break;
1086	case m_signbit: intres = signbit(d_arg1.f); break;
1087	case m_fpclassify: intres = fpclassify(d_arg1.f); break;
1088	case m_isgreater: intres = isgreater(d_arg1.f, d_arg2.f); break;
1089	case m_isgreaterequal: intres = isgreaterequal(d_arg1.f, d_arg2.f); break;
1090	case m_isless: intres = isless(d_arg1.f, d_arg2.f); break;
1091	case m_islessequal: intres = islessequal(d_arg1.f, d_arg2.f); break;
1092	case m_islessgreater: intres = islessgreater(d_arg1.f, d_arg2.f); break;
1093	case m_isunordered: intres = isunordered(d_arg1.f, d_arg2.f); break;
1094
1095	case m_isfinitef: intres = isfinite(s_arg1.f); break;
1096	case m_isinff: intres = isinf(s_arg1.f); break;
1097	case m_isnanf: intres = isnan(s_arg1.f); break;
1098	case m_isnormalf: intres = isnormal(s_arg1.f); break;
1099	case m_signbitf: intres = signbit(s_arg1.f); break;
1100	case m_fpclassifyf: intres = fpclassify(s_arg1.f); break;
1101	case m_isgreaterf: intres = isgreater(s_arg1.f, s_arg2.f); break;
1102	case m_isgreaterequalf: intres = isgreaterequal(s_arg1.f, s_arg2.f); break;
1103	case m_islessf: intres = isless(s_arg1.f, s_arg2.f); break;
1104	case m_islessequalf: intres = islessequal(s_arg1.f, s_arg2.f); break;
1105	case m_islessgreaterf: intres = islessgreater(s_arg1.f, s_arg2.f); break;
1106	case m_isunorderedf: intres = isunordered(s_arg1.f, s_arg2.f); break;
1107
1108	default:
1109	printf(format: "unhandled macro: %s\n",t.func->name);
1110	return test_fail;
1111	}
1112	}
1113
1114	/*
1115	* Decline the test if the deferred decline flag was set above.
1116	*/
1117	if (deferred_decline)
1118	return test_decline;
1119
1120	/* printf("intres=%i\n",intres); */
1121
1122	/* Clear the fail text (indicating a pass unless we change it) */
1123	failp[0] = '\0';
1124
1125	/* Check the IEEE status bits (except INX, which we disregard).
1126	* We don't bother with this for complex numbers, because the
1127	* complex functions are hard to get exactly right and we don't
1128	* have to anyway (C99 annex G is only informative). */
1129	if (!(is_complex_argtype(argtype: t.func->argtype) || is_complex_rettype(rettype: t.func->rettype))) {
1130	status = fetestexcept(FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW);
1131	if ((status|t.maybestatus|~statusmask) != (t.status|t.maybestatus|~statusmask)) {
1132	if (quiet) failtext[0]='x';
1133	else {
1134	failp += sprintf(s: failp,
1135	format: " wrongstatus=%s%s%s%s%s",
1136	(status & FE_INVALID ? "i" : ""),
1137	(status & FE_DIVBYZERO ? "z" : ""),
1138	(status & FE_OVERFLOW ? "o" : ""),
1139	(status & FE_UNDERFLOW ? "u" : ""),
1140	(status ? "" : "OK"));
1141	}
1142	}
1143	}
1144
1145	/* Check the result */
1146	{
1147	unsigned resultr[2], resulti[2];
1148	unsigned tresultr[3], tresulti[3], wres;
1149
1150	switch(t.func->rettype) {
1151	case rt_d:
1152	case rt_d2:
1153	tresultr[0] = t.resultr[0];
1154	tresultr[1] = t.resultr[1];
1155	resultr[0] = d_res.i[dmsd]; resultr[1] = d_res.i[dlsd];
1156	wres = 2;
1157	break;
1158	case rt_i:
1159	tresultr[0] = t.resultr[0];
1160	resultr[0] = intres;
1161	wres = 1;
1162	break;
1163	case rt_s:
1164	case rt_s2:
1165	tresultr[0] = t.resultr[0];
1166	resultr[0] = s_res.i;
1167	wres = 1;
1168	break;
1169	default:
1170	puts(s: "unhandled rettype in runtest");
1171	wres = 0;
1172	}
1173	if(t.resultc != rc_none) {
1174	int err = 0;
1175	switch(t.resultc) {
1176	case rc_zero:
1177	if(resultr[0] != 0 || resulti[0] != 0 ||
1178	(wres==2 && (resultr[1] != 0 || resulti[1] != 0))) {
1179	err = 1;
1180	}
1181	break;
1182	case rc_infinity:
1183	if(wres==1) {
1184	if(!((resultr[0]&0x7fffffff)==0x7f800000 ||
1185	(resulti[0]&0x7fffffff)==0x7f800000)) {
1186	err = 1;
1187	}
1188	} else {
1189	if(!(((resultr[0]&0x7fffffff)==0x7ff00000 && resultr[1]==0) ||
1190	((resulti[0]&0x7fffffff)==0x7ff00000 && resulti[1]==0))) {
1191	err = 1;
1192	}
1193	}
1194	break;
1195	case rc_nan:
1196	if(wres==1) {
1197	if(!((resultr[0]&0x7fffffff)>0x7f800000 ||
1198	(resulti[0]&0x7fffffff)>0x7f800000)) {
1199	err = 1;
1200	}
1201	} else {
1202	canon_dNaN(a: resultr);
1203	canon_dNaN(a: resulti);
1204	if(!(((resultr[0]&0x7fffffff)>0x7ff00000 && resultr[1]==1) ||
1205	((resulti[0]&0x7fffffff)>0x7ff00000 && resulti[1]==1))) {
1206	err = 1;
1207	}
1208	}
1209	break;
1210	case rc_finite:
1211	if(wres==1) {
1212	if(!((resultr[0]&0x7fffffff)<0x7f800000 ||
1213	(resulti[0]&0x7fffffff)<0x7f800000)) {
1214	err = 1;
1215	}
1216	} else {
1217	if(!((resultr[0]&0x7fffffff)<0x7ff00000 ||
1218	(resulti[0]&0x7fffffff)<0x7ff00000)) {
1219	err = 1;
1220	}
1221	}
1222	break;
1223	default:
1224	break;
1225	}
1226	if(err) {
1227	print_error(rettype: t.func->rettype,result: resultr,text: "wrongresultr",failp: &failp);
1228	print_error(rettype: t.func->rettype,result: resulti,text: "wrongresulti",failp: &failp);
1229	}
1230	} else if (t.nresult > wres) {
1231	/*
1232	* The test case data has provided the result to more
1233	* than double precision. Instead of testing exact
1234	* equality, we test against our maximum error
1235	* tolerance.
1236	*/
1237	int rshift, ishift;
1238	long long ulpsr, ulpsi, ulptolerance;
1239
1240	tresultr[wres] = t.resultr[wres] << (32-EXTRABITS);
1241	tresulti[wres] = t.resulti[wres] << (32-EXTRABITS);
1242	if(strict) {
1243	ulptolerance = 4096; /* one ulp */
1244	} else {
1245	ulptolerance = t.func->tolerance;
1246	}
1247	rshift = ishift = 0;
1248	if (ulptolerance & ABSLOWERBOUND) {
1249	/*
1250	* Hack for the lgamma functions, which have an
1251	* error behaviour that can't conveniently be
1252	* characterised in pure ULPs. Really, we want to
1253	* say that the error in lgamma is "at most N ULPs,
1254	* or at most an absolute error of X, whichever is
1255	* larger", for appropriately chosen N,X. But since
1256	* these two functions are the only cases where it
1257	* arises, I haven't bothered to do it in a nice way
1258	* in the function table above.
1259	*
1260	* (The difficult cases arise with negative input
1261	* values such that |gamma(x)| is very near to 1; in
1262	* this situation implementations tend to separately
1263	* compute lgamma(|x|) and the log of the correction
1264	* term from the Euler reflection formula, and
1265	* subtract - which catastrophically loses
1266	* significance.)
1267	*
1268	* As far as I can tell, nobody cares about this:
1269	* GNU libm doesn't get those cases right either,
1270	* and OpenCL explicitly doesn't state a ULP error
1271	* limit for lgamma. So my guess is that this is
1272	* simply considered acceptable error behaviour for
1273	* this particular function, and hence I feel free
1274	* to allow for it here.
1275	*/
1276	ulptolerance &= ~ABSLOWERBOUND;
1277	if (t.op1r[0] & 0x80000000) {
1278	if (t.func->rettype == rt_d)
1279	rshift = 0x400 - ((tresultr[0] >> 20) & 0x7ff);
1280	else if (t.func->rettype == rt_s)
1281	rshift = 0x80 - ((tresultr[0] >> 23) & 0xff);
1282	if (rshift < 0)
1283	rshift = 0;
1284	}
1285	}
1286	if (ulptolerance & PLUSMINUSPIO2) {
1287	ulptolerance &= ~PLUSMINUSPIO2;
1288	/*
1289	* Hack for range reduction, which can reduce
1290	* borderline cases in the wrong direction, i.e.
1291	* return a value just outside one end of the interval
1292	* [-pi/4,+pi/4] when it could have returned a value
1293	* just inside the other end by subtracting an
1294	* adjacent multiple of pi/2.
1295	*
1296	* We tolerate this, up to a point, because the
1297	* trigonometric functions making use of the output of
1298	* rred can cope and because making the range reducer
1299	* do the exactly right thing in every case would be
1300	* more expensive.
1301	*/
1302	if (wres == 1) {
1303	/* Upper bound of overshoot derived in rredf.h */
1304	if ((resultr[0]&0x7FFFFFFF) <= 0x3f494b02 &&
1305	(resultr[0]&0x7FFFFFFF) > 0x3f490fda &&
1306	(resultr[0]&0x80000000) != (tresultr[0]&0x80000000)) {
1307	unsigned long long val;
1308	val = tresultr[0];
1309	val = (val << 32) | tresultr[1];
1310	/*
1311	* Compute the alternative permitted result by
1312	* subtracting from the sum of the extended
1313	* single-precision bit patterns of +pi/4 and
1314	* -pi/4. This is a horrible hack which only
1315	* works because we can be confident that
1316	* numbers in this range all have the same
1317	* exponent!
1318	*/
1319	val = 0xfe921fb54442d184ULL - val;
1320	tresultr[0] = val >> 32;
1321	tresultr[1] = (val >> (32-EXTRABITS)) << (32-EXTRABITS);
1322	/*
1323	* Also, expect a correspondingly different
1324	* value of res2 as a result of this change.
1325	* The adjustment depends on whether we just
1326	* flipped the result from + to - or vice
1327	* versa.
1328	*/
1329	if (resultr[0] & 0x80000000) {
1330	res2_adjust = +1;
1331	} else {
1332	res2_adjust = -1;
1333	}
1334	}
1335	}
1336	}
1337	ulpsr = calc_error(a: resultr, b: tresultr, shift: rshift, rettype: t.func->rettype);
1338	if(is_complex_rettype(rettype: t.func->rettype)) {
1339	ulpsi = calc_error(a: resulti, b: tresulti, shift: ishift, rettype: t.func->rettype);
1340	} else {
1341	ulpsi = 0;
1342	}
1343	unsigned *rr = (ulpsr > ulptolerance || ulpsr < -ulptolerance) ? resultr : NULL;
1344	unsigned *ri = (ulpsi > ulptolerance || ulpsi < -ulptolerance) ? resulti : NULL;
1345	/* printf("tolerance=%i, ulpsr=%i, ulpsi=%i, rr=%p, ri=%p\n",ulptolerance,ulpsr,ulpsi,rr,ri); */
1346	if (rr || ri) {
1347	if (quiet) failtext[0]='x';
1348	else {
1349	print_error(rettype: t.func->rettype,result: rr,text: "wrongresultr",failp: &failp);
1350	print_error(rettype: t.func->rettype,result: ri,text: "wrongresulti",failp: &failp);
1351	print_ulps(rettype: t.func->rettype,ulpsr: rr ? ulpsr : 0, ulpsi: ri ? ulpsi : 0,failp: &failp);
1352	}
1353	}
1354	} else {
1355	if(is_complex_rettype(rettype: t.func->rettype))
1356	/*
1357	* Complex functions are not fully supported,
1358	* this is unreachable, but prevents warnings.
1359	*/
1360	abort();
1361	/*
1362	* The test case data has provided the result in
1363	* exactly the output precision. Therefore we must
1364	* complain about _any_ violation.
1365	*/
1366	switch(t.func->rettype) {
1367	case rt_dc:
1368	canon_dNaN(a: tresulti);
1369	canon_dNaN(a: resulti);
1370	if (fo) {
1371	dnormzero(a: tresulti);
1372	dnormzero(a: resulti);
1373	}
1374	/* deliberate fall-through */
1375	case rt_d:
1376	canon_dNaN(a: tresultr);
1377	canon_dNaN(a: resultr);
1378	if (fo) {
1379	dnormzero(a: tresultr);
1380	dnormzero(a: resultr);
1381	}
1382	break;
1383	case rt_sc:
1384	canon_sNaN(a: tresulti);
1385	canon_sNaN(a: resulti);
1386	if (fo) {
1387	snormzero(a: tresulti);
1388	snormzero(a: resulti);
1389	}
1390	/* deliberate fall-through */
1391	case rt_s:
1392	canon_sNaN(a: tresultr);
1393	canon_sNaN(a: resultr);
1394	if (fo) {
1395	snormzero(a: tresultr);
1396	snormzero(a: resultr);
1397	}
1398	break;
1399	default:
1400	break;
1401	}
1402	if(is_complex_rettype(rettype: t.func->rettype)) {
1403	unsigned *rr, *ri;
1404	if(resultr[0] != tresultr[0] ||
1405	(wres > 1 && resultr[1] != tresultr[1])) {
1406	rr = resultr;
1407	} else {
1408	rr = NULL;
1409	}
1410	if(resulti[0] != tresulti[0] ||
1411	(wres > 1 && resulti[1] != tresulti[1])) {
1412	ri = resulti;
1413	} else {
1414	ri = NULL;
1415	}
1416	if(rr || ri) {
1417	if (quiet) failtext[0]='x';
1418	print_error(rettype: t.func->rettype,result: rr,text: "wrongresultr",failp: &failp);
1419	print_error(rettype: t.func->rettype,result: ri,text: "wrongresulti",failp: &failp);
1420	}
1421	} else if (resultr[0] != tresultr[0] ||
1422	(wres > 1 && resultr[1] != tresultr[1])) {
1423	if (quiet) failtext[0]='x';
1424	print_error(rettype: t.func->rettype,result: resultr,text: "wrongresult",failp: &failp);
1425	}
1426	}
1427	/*
1428	* Now test res2, for those functions (frexp, modf, rred)
1429	* which use it.
1430	*/
1431	if (t.func->func.ptr == &frexp || t.func->func.ptr == &frexpf ||
1432	t.func->macro_name == m_rred || t.func->macro_name == m_rredf) {
1433	unsigned tres2 = t.res2[0];
1434	if (res2_adjust) {
1435	/* Fix for range reduction, propagated from further up */
1436	tres2 = (tres2 + res2_adjust) & 3;
1437	}
1438	if (tres2 != intres) {
1439	if (quiet) failtext[0]='x';
1440	else {
1441	failp += sprintf(s: failp,
1442	format: " wrongres2=%08x", intres);
1443	}
1444	}
1445	} else if (t.func->func.ptr == &modf || t.func->func.ptr == &modff) {
1446	tresultr[0] = t.res2[0];
1447	tresultr[1] = t.res2[1];
1448	if (is_double_rettype(rettype: t.func->rettype)) {
1449	canon_dNaN(a: tresultr);
1450	resultr[0] = d_res2.i[dmsd];
1451	resultr[1] = d_res2.i[dlsd];
1452	canon_dNaN(a: resultr);
1453	if (fo) {
1454	dnormzero(a: tresultr);
1455	dnormzero(a: resultr);
1456	}
1457	} else {
1458	canon_sNaN(a: tresultr);
1459	resultr[0] = s_res2.i;
1460	resultr[1] = s_res2.i;
1461	canon_sNaN(a: resultr);
1462	if (fo) {
1463	snormzero(a: tresultr);
1464	snormzero(a: resultr);
1465	}
1466	}
1467	if (resultr[0] != tresultr[0] ||
1468	(wres > 1 && resultr[1] != tresultr[1])) {
1469	if (quiet) failtext[0]='x';
1470	else {
1471	if (is_double_rettype(rettype: t.func->rettype))
1472	failp += sprintf(s: failp, format: " wrongres2=%08x.%08x",
1473	resultr[0], resultr[1]);
1474	else
1475	failp += sprintf(s: failp, format: " wrongres2=%08x",
1476	resultr[0]);
1477	}
1478	}
1479	}
1480	}
1481
1482	/* Check errno */
1483	err = (errno == EDOM ? e_EDOM : errno == ERANGE ? e_ERANGE : e_0);
1484	if (err != t.err && err != t.maybeerr) {
1485	if (quiet) failtext[0]='x';
1486	else {
1487	failp += sprintf(s: failp, format: " wrongerrno=%s expecterrno=%s ", errnos[err], errnos[t.err]);
1488	}
1489	}
1490
1491	return *failtext ? test_fail : test_pass;
1492	}
1493
1494	int passed, failed, declined;
1495
1496	void runtests(char *name, FILE *fp) {
1497	char testbuf[512], linebuf[512];
1498	int lineno = 1;
1499	testdetail test;
1500
1501	test.valid = 0;
1502
1503	if (verbose) printf(format: "runtests: %s\n", name);
1504	while (fgets(s: testbuf, n: sizeof(testbuf), stream: fp)) {
1505	int res, print_errno;
1506	testbuf[strcspn(s: testbuf, reject: "\r\n")] = '\0';
1507	strcpy(dest: linebuf, src: testbuf);
1508	test = parsetest(testbuf, oldtest: test);
1509	print_errno = 0;
1510	while (test.in_err < test.in_err_limit) {
1511	res = runtest(t: test);
1512	if (res == test_pass) {
1513	if (verbose)
1514	printf(format: "%s:%d: pass\n", name, lineno);
1515	++passed;
1516	} else if (res == test_decline) {
1517	if (verbose)
1518	printf(format: "%s:%d: declined\n", name, lineno);
1519	++declined;
1520	} else if (res == test_fail) {
1521	if (!quiet)
1522	printf(format: "%s:%d: FAIL%s: %s%s%s%s\n", name, lineno,
1523	test.random ? " (random)" : "",
1524	linebuf,
1525	print_errno ? " errno_in=" : "",
1526	print_errno ? errnos[test.in_err] : "",
1527	failtext);
1528	++failed;
1529	} else if (res == test_invalid) {
1530	printf(format: "%s:%d: malformed: %s\n", name, lineno, linebuf);
1531	++failed;
1532	}
1533	test.in_err++;
1534	print_errno = 1;
1535	}
1536	lineno++;
1537	}
1538	}
1539
1540	int main(int ac, char **av) {
1541	char **files;
1542	int i, nfiles = 0;
1543	dbl d;
1544
1545	#ifdef MICROLIB
1546	/*
1547	* Invent argc and argv ourselves.
1548	*/
1549	char *argv[256];
1550	char args[256];
1551	{
1552	int sargs[2];
1553	char *p;
1554
1555	ac = 0;
1556
1557	sargs[0]=(int)args;
1558	sargs[1]=(int)sizeof(args);
1559	if (!__semihost(0x15, sargs)) {
1560	args[sizeof(args)-1] = '\0'; /* just in case */
1561	p = args;
1562	while (1) {
1563	while (*p == ' ' || *p == '\t') p++;
1564	if (!*p) break;
1565	argv[ac++] = p;
1566	while (*p && *p != ' ' && *p != '\t') p++;
1567	if (*p) *p++ = '\0';
1568	}
1569	}
1570
1571	av = argv;
1572	}
1573	#endif
1574
1575	/* Sort tfuncs */
1576	qsort(base: tfuncs, nmemb: sizeof(tfuncs)/sizeof(test_func), size: sizeof(test_func), compar: &compare_tfuncs);
1577
1578	/*
1579	* Autodetect the `double' endianness.
1580	*/
1581	dmsd = 0;
1582	d.f = 1.0; /* 0x3ff00000 / 0x00000000 */
1583	if (d.i[dmsd] == 0) {
1584	dmsd = 1;
1585	}
1586	/*
1587	* Now dmsd denotes what the compiler thinks we're at. Let's
1588	* check that it agrees with what the runtime thinks.
1589	*/
1590	d.i[0] = d.i[1] = 0x11111111;/* a random +ve number */
1591	d.f /= d.f; /* must now be one */
1592	if (d.i[dmsd] == 0) {
1593	fprintf(stderr, format: "YIKES! Compiler and runtime disagree on endianness"
1594	" of `double'. Bailing out\n");
1595	return 1;
1596	}
1597	dlsd = !dmsd;
1598
1599	/* default is terse */
1600	verbose = 0;
1601	fo = 0;
1602	strict = 0;
1603
1604	files = (char **)malloc(size: (ac+1) * sizeof(char *));
1605	if (!files) {
1606	fprintf(stderr, format: "initial malloc failed!\n");
1607	return 1;
1608	}
1609	#ifdef NOCMDLINE
1610	files[nfiles++] = "testfile";
1611	#endif
1612
1613	while (--ac) {
1614	char *p = *++av;
1615	if (*p == '-') {
1616	static char *options[] = {
1617	"-fo",
1618	#if 0
1619	"-noinexact",
1620	"-noround",
1621	#endif
1622	"-nostatus",
1623	"-quiet",
1624	"-strict",
1625	"-v",
1626	"-verbose",
1627	};
1628	enum {
1629	op_fo,
1630	#if 0
1631	op_noinexact,
1632	op_noround,
1633	#endif
1634	op_nostatus,
1635	op_quiet,
1636	op_strict,
1637	op_v,
1638	op_verbose,
1639	};
1640	switch (find(word: p, array: options, asize: sizeof(options))) {
1641	case op_quiet:
1642	quiet = 1;
1643	break;
1644	#if 0
1645	case op_noinexact:
1646	statusmask &= 0x0F; /* remove bit 4 */
1647	break;
1648	case op_noround:
1649	doround = 0;
1650	break;
1651	#endif
1652	case op_nostatus: /* no status word => noinx,noround */
1653	statusmask = 0;
1654	doround = 0;
1655	break;
1656	case op_v:
1657	case op_verbose:
1658	verbose = 1;
1659	break;
1660	case op_fo:
1661	fo = 1;
1662	break;
1663	case op_strict: /* tolerance is 1 ulp */
1664	strict = 1;
1665	break;
1666	default:
1667	fprintf(stderr, format: "unrecognised option: %s\n", p);
1668	break;
1669	}
1670	} else {
1671	files[nfiles++] = p;
1672	}
1673	}
1674
1675	passed = failed = declined = 0;
1676
1677	if (nfiles) {
1678	for (i = 0; i < nfiles; i++) {
1679	FILE *fp = fopen(filename: files[i], modes: "r");
1680	if (!fp) {
1681	fprintf(stderr, format: "Couldn't open %s\n", files[i]);
1682	} else
1683	runtests(name: files[i], fp);
1684	}
1685	} else
1686	runtests(name: "(stdin)", stdin);
1687
1688	printf(format: "Completed. Passed %d, failed %d (total %d",
1689	passed, failed, passed+failed);
1690	if (declined)
1691	printf(format: " plus %d declined", declined);
1692	printf(format: ")\n");
1693	if (failed || passed == 0)
1694	return 1;
1695	printf(format: "** TEST PASSED OK **\n");
1696	return 0;
1697	}
1698
1699	void undef_func() {
1700	failed++;
1701	puts(s: "ERROR: undefined function called");
1702	}
1703