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