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