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