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