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