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