1 /* 2 * ULP error checking tool for math functions. 3 * 4 * Copyright (c) 2019-2020, Arm Limited. 5 * SPDX-License-Identifier: MIT 6 */ 7 8 #include <ctype.h> 9 #include <fenv.h> 10 #include <float.h> 11 #include <math.h> 12 #include <stdint.h> 13 #include <stdio.h> 14 #include <stdlib.h> 15 #include <string.h> 16 #include "mathlib.h" 17 18 /* Don't depend on mpfr by default. */ 19 #ifndef USE_MPFR 20 # define USE_MPFR 0 21 #endif 22 #if USE_MPFR 23 # include <mpfr.h> 24 #endif 25 26 #ifndef WANT_VMATH 27 /* Enable the build of vector math code. */ 28 # define WANT_VMATH 1 29 #endif 30 31 static inline uint64_t 32 asuint64 (double f) 33 { 34 union 35 { 36 double f; 37 uint64_t i; 38 } u = {f}; 39 return u.i; 40 } 41 42 static inline double 43 asdouble (uint64_t i) 44 { 45 union 46 { 47 uint64_t i; 48 double f; 49 } u = {i}; 50 return u.f; 51 } 52 53 static inline uint32_t 54 asuint (float f) 55 { 56 union 57 { 58 float f; 59 uint32_t i; 60 } u = {f}; 61 return u.i; 62 } 63 64 static inline float 65 asfloat (uint32_t i) 66 { 67 union 68 { 69 uint32_t i; 70 float f; 71 } u = {i}; 72 return u.f; 73 } 74 75 static uint64_t seed = 0x0123456789abcdef; 76 static uint64_t 77 rand64 (void) 78 { 79 seed = 6364136223846793005ull * seed + 1; 80 return seed ^ (seed >> 32); 81 } 82 83 /* Uniform random in [0,n]. */ 84 static uint64_t 85 randn (uint64_t n) 86 { 87 uint64_t r, m; 88 89 if (n == 0) 90 return 0; 91 n++; 92 if (n == 0) 93 return rand64 (); 94 for (;;) 95 { 96 r = rand64 (); 97 m = r % n; 98 if (r - m <= -n) 99 return m; 100 } 101 } 102 103 struct gen 104 { 105 uint64_t start; 106 uint64_t len; 107 uint64_t start2; 108 uint64_t len2; 109 uint64_t off; 110 uint64_t step; 111 uint64_t cnt; 112 }; 113 114 struct args_f1 115 { 116 float x; 117 }; 118 119 struct args_f2 120 { 121 float x; 122 float x2; 123 }; 124 125 struct args_d1 126 { 127 double x; 128 }; 129 130 struct args_d2 131 { 132 double x; 133 double x2; 134 }; 135 136 /* result = y + tail*2^ulpexp. */ 137 struct ret_f 138 { 139 float y; 140 double tail; 141 int ulpexp; 142 int ex; 143 int ex_may; 144 }; 145 146 struct ret_d 147 { 148 double y; 149 double tail; 150 int ulpexp; 151 int ex; 152 int ex_may; 153 }; 154 155 static inline uint64_t 156 next1 (struct gen *g) 157 { 158 /* For single argument use randomized incremental steps, 159 that produce dense sampling without collisions and allow 160 testing all inputs in a range. */ 161 uint64_t r = g->start + g->off; 162 g->off += g->step + randn (g->step / 2); 163 if (g->off > g->len) 164 g->off -= g->len; /* hack. */ 165 return r; 166 } 167 168 static inline uint64_t 169 next2 (uint64_t *x2, struct gen *g) 170 { 171 /* For two arguments use uniform random sampling. */ 172 uint64_t r = g->start + randn (g->len); 173 *x2 = g->start2 + randn (g->len2); 174 return r; 175 } 176 177 static struct args_f1 178 next_f1 (void *g) 179 { 180 return (struct args_f1){asfloat (next1 (g))}; 181 } 182 183 static struct args_f2 184 next_f2 (void *g) 185 { 186 uint64_t x2; 187 uint64_t x = next2 (&x2, g); 188 return (struct args_f2){asfloat (x), asfloat (x2)}; 189 } 190 191 static struct args_d1 192 next_d1 (void *g) 193 { 194 return (struct args_d1){asdouble (next1 (g))}; 195 } 196 197 static struct args_d2 198 next_d2 (void *g) 199 { 200 uint64_t x2; 201 uint64_t x = next2 (&x2, g); 202 return (struct args_d2){asdouble (x), asdouble (x2)}; 203 } 204 205 struct conf 206 { 207 int r; 208 int rc; 209 int quiet; 210 int mpfr; 211 int fenv; 212 unsigned long long n; 213 double softlim; 214 double errlim; 215 }; 216 217 /* Wrappers for sincos. */ 218 static float sincosf_sinf(float x) {(void)cosf(x); return sinf(x);} 219 static float sincosf_cosf(float x) {(void)sinf(x); return cosf(x);} 220 static double sincos_sin(double x) {(void)cos(x); return sin(x);} 221 static double sincos_cos(double x) {(void)sin(x); return cos(x);} 222 #if USE_MPFR 223 static int sincos_mpfr_sin(mpfr_t y, const mpfr_t x, mpfr_rnd_t r) { mpfr_cos(y,x,r); return mpfr_sin(y,x,r); } 224 static int sincos_mpfr_cos(mpfr_t y, const mpfr_t x, mpfr_rnd_t r) { mpfr_sin(y,x,r); return mpfr_cos(y,x,r); } 225 #endif 226 227 /* A bit of a hack: call vector functions twice with the same 228 input in lane 0 but a different value in other lanes: once 229 with an in-range value and then with a special case value. */ 230 static int secondcall; 231 232 /* Wrappers for vector functions. */ 233 #if __aarch64__ && WANT_VMATH 234 typedef __f32x4_t v_float; 235 typedef __f64x2_t v_double; 236 static const float fv[2] = {1.0f, -INFINITY}; 237 static const double dv[2] = {1.0, -INFINITY}; 238 static inline v_float argf(float x) { return (v_float){x,x,x,fv[secondcall]}; } 239 static inline v_double argd(double x) { return (v_double){x,dv[secondcall]}; } 240 241 static float v_sinf(float x) { return __v_sinf(argf(x))[0]; } 242 static float v_cosf(float x) { return __v_cosf(argf(x))[0]; } 243 static float v_expf_1u(float x) { return __v_expf_1u(argf(x))[0]; } 244 static float v_expf(float x) { return __v_expf(argf(x))[0]; } 245 static float v_exp2f_1u(float x) { return __v_exp2f_1u(argf(x))[0]; } 246 static float v_exp2f(float x) { return __v_exp2f(argf(x))[0]; } 247 static float v_logf(float x) { return __v_logf(argf(x))[0]; } 248 static float v_powf(float x, float y) { return __v_powf(argf(x),argf(y))[0]; } 249 static double v_sin(double x) { return __v_sin(argd(x))[0]; } 250 static double v_cos(double x) { return __v_cos(argd(x))[0]; } 251 static double v_exp(double x) { return __v_exp(argd(x))[0]; } 252 static double v_log(double x) { return __v_log(argd(x))[0]; } 253 static double v_pow(double x, double y) { return __v_pow(argd(x),argd(y))[0]; } 254 #ifdef __vpcs 255 static float vn_sinf(float x) { return __vn_sinf(argf(x))[0]; } 256 static float vn_cosf(float x) { return __vn_cosf(argf(x))[0]; } 257 static float vn_expf_1u(float x) { return __vn_expf_1u(argf(x))[0]; } 258 static float vn_expf(float x) { return __vn_expf(argf(x))[0]; } 259 static float vn_exp2f_1u(float x) { return __vn_exp2f_1u(argf(x))[0]; } 260 static float vn_exp2f(float x) { return __vn_exp2f(argf(x))[0]; } 261 static float vn_logf(float x) { return __vn_logf(argf(x))[0]; } 262 static float vn_powf(float x, float y) { return __vn_powf(argf(x),argf(y))[0]; } 263 static double vn_sin(double x) { return __vn_sin(argd(x))[0]; } 264 static double vn_cos(double x) { return __vn_cos(argd(x))[0]; } 265 static double vn_exp(double x) { return __vn_exp(argd(x))[0]; } 266 static double vn_log(double x) { return __vn_log(argd(x))[0]; } 267 static double vn_pow(double x, double y) { return __vn_pow(argd(x),argd(y))[0]; } 268 static float Z_sinf(float x) { return _ZGVnN4v_sinf(argf(x))[0]; } 269 static float Z_cosf(float x) { return _ZGVnN4v_cosf(argf(x))[0]; } 270 static float Z_expf(float x) { return _ZGVnN4v_expf(argf(x))[0]; } 271 static float Z_exp2f(float x) { return _ZGVnN4v_exp2f(argf(x))[0]; } 272 static float Z_logf(float x) { return _ZGVnN4v_logf(argf(x))[0]; } 273 static float Z_powf(float x, float y) { return _ZGVnN4vv_powf(argf(x),argf(y))[0]; } 274 static double Z_sin(double x) { return _ZGVnN2v_sin(argd(x))[0]; } 275 static double Z_cos(double x) { return _ZGVnN2v_cos(argd(x))[0]; } 276 static double Z_exp(double x) { return _ZGVnN2v_exp(argd(x))[0]; } 277 static double Z_log(double x) { return _ZGVnN2v_log(argd(x))[0]; } 278 static double Z_pow(double x, double y) { return _ZGVnN2vv_pow(argd(x),argd(y))[0]; } 279 #endif 280 #endif 281 282 struct fun 283 { 284 const char *name; 285 int arity; 286 int singleprec; 287 int twice; 288 union 289 { 290 float (*f1) (float); 291 float (*f2) (float, float); 292 double (*d1) (double); 293 double (*d2) (double, double); 294 } fun; 295 union 296 { 297 double (*f1) (double); 298 double (*f2) (double, double); 299 long double (*d1) (long double); 300 long double (*d2) (long double, long double); 301 } fun_long; 302 #if USE_MPFR 303 union 304 { 305 int (*f1) (mpfr_t, const mpfr_t, mpfr_rnd_t); 306 int (*f2) (mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t); 307 int (*d1) (mpfr_t, const mpfr_t, mpfr_rnd_t); 308 int (*d2) (mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t); 309 } fun_mpfr; 310 #endif 311 }; 312 313 static const struct fun fun[] = { 314 #if USE_MPFR 315 # define F(x, x_wrap, x_long, x_mpfr, a, s, t, twice) \ 316 {#x, a, s, twice, {.t = x_wrap}, {.t = x_long}, {.t = x_mpfr}}, 317 #else 318 # define F(x, x_wrap, x_long, x_mpfr, a, s, t, twice) \ 319 {#x, a, s, twice, {.t = x_wrap}, {.t = x_long}}, 320 #endif 321 #define F1(x) F (x##f, x##f, x, mpfr_##x, 1, 1, f1, 0) 322 #define F2(x) F (x##f, x##f, x, mpfr_##x, 2, 1, f2, 0) 323 #define D1(x) F (x, x, x##l, mpfr_##x, 1, 0, d1, 0) 324 #define D2(x) F (x, x, x##l, mpfr_##x, 2, 0, d2, 0) 325 F1 (sin) 326 F1 (cos) 327 F (sincosf_sinf, sincosf_sinf, sincos_sin, sincos_mpfr_sin, 1, 1, f1, 0) 328 F (sincosf_cosf, sincosf_cosf, sincos_cos, sincos_mpfr_cos, 1, 1, f1, 0) 329 F1 (exp) 330 F1 (exp2) 331 F1 (log) 332 F1 (log2) 333 F2 (pow) 334 F1 (erf) 335 D1 (exp) 336 D1 (exp2) 337 D1 (log) 338 D1 (log2) 339 D2 (pow) 340 D1 (erf) 341 #if WANT_VMATH 342 F (__s_sinf, __s_sinf, sin, mpfr_sin, 1, 1, f1, 0) 343 F (__s_cosf, __s_cosf, cos, mpfr_cos, 1, 1, f1, 0) 344 F (__s_expf_1u, __s_expf_1u, exp, mpfr_exp, 1, 1, f1, 0) 345 F (__s_expf, __s_expf, exp, mpfr_exp, 1, 1, f1, 0) 346 F (__s_exp2f_1u, __s_exp2f_1u, exp2, mpfr_exp2, 1, 1, f1, 0) 347 F (__s_exp2f, __s_exp2f, exp2, mpfr_exp2, 1, 1, f1, 0) 348 F (__s_powf, __s_powf, pow, mpfr_pow, 2, 1, f2, 0) 349 F (__s_logf, __s_logf, log, mpfr_log, 1, 1, f1, 0) 350 F (__s_sin, __s_sin, sinl, mpfr_sin, 1, 0, d1, 0) 351 F (__s_cos, __s_cos, cosl, mpfr_cos, 1, 0, d1, 0) 352 F (__s_exp, __s_exp, expl, mpfr_exp, 1, 0, d1, 0) 353 F (__s_log, __s_log, logl, mpfr_log, 1, 0, d1, 0) 354 F (__s_pow, __s_pow, powl, mpfr_pow, 2, 0, d2, 0) 355 #if __aarch64__ 356 F (__v_sinf, v_sinf, sin, mpfr_sin, 1, 1, f1, 1) 357 F (__v_cosf, v_cosf, cos, mpfr_cos, 1, 1, f1, 1) 358 F (__v_expf_1u, v_expf_1u, exp, mpfr_exp, 1, 1, f1, 1) 359 F (__v_expf, v_expf, exp, mpfr_exp, 1, 1, f1, 1) 360 F (__v_exp2f_1u, v_exp2f_1u, exp2, mpfr_exp2, 1, 1, f1, 1) 361 F (__v_exp2f, v_exp2f, exp2, mpfr_exp2, 1, 1, f1, 1) 362 F (__v_logf, v_logf, log, mpfr_log, 1, 1, f1, 1) 363 F (__v_powf, v_powf, pow, mpfr_pow, 2, 1, f2, 1) 364 F (__v_sin, v_sin, sinl, mpfr_sin, 1, 0, d1, 1) 365 F (__v_cos, v_cos, cosl, mpfr_cos, 1, 0, d1, 1) 366 F (__v_exp, v_exp, expl, mpfr_exp, 1, 0, d1, 1) 367 F (__v_log, v_log, logl, mpfr_log, 1, 0, d1, 1) 368 F (__v_pow, v_pow, powl, mpfr_pow, 2, 0, d2, 1) 369 #ifdef __vpcs 370 F (__vn_sinf, vn_sinf, sin, mpfr_sin, 1, 1, f1, 1) 371 F (__vn_cosf, vn_cosf, cos, mpfr_cos, 1, 1, f1, 1) 372 F (__vn_expf_1u, vn_expf_1u, exp, mpfr_exp, 1, 1, f1, 1) 373 F (__vn_expf, vn_expf, exp, mpfr_exp, 1, 1, f1, 1) 374 F (__vn_exp2f_1u, vn_exp2f_1u, exp2, mpfr_exp2, 1, 1, f1, 1) 375 F (__vn_exp2f, vn_exp2f, exp2, mpfr_exp2, 1, 1, f1, 1) 376 F (__vn_logf, vn_logf, log, mpfr_log, 1, 1, f1, 1) 377 F (__vn_powf, vn_powf, pow, mpfr_pow, 2, 1, f2, 1) 378 F (__vn_sin, vn_sin, sinl, mpfr_sin, 1, 0, d1, 1) 379 F (__vn_cos, vn_cos, cosl, mpfr_cos, 1, 0, d1, 1) 380 F (__vn_exp, vn_exp, expl, mpfr_exp, 1, 0, d1, 1) 381 F (__vn_log, vn_log, logl, mpfr_log, 1, 0, d1, 1) 382 F (__vn_pow, vn_pow, powl, mpfr_pow, 2, 0, d2, 1) 383 F (_ZGVnN4v_sinf, Z_sinf, sin, mpfr_sin, 1, 1, f1, 1) 384 F (_ZGVnN4v_cosf, Z_cosf, cos, mpfr_cos, 1, 1, f1, 1) 385 F (_ZGVnN4v_expf, Z_expf, exp, mpfr_exp, 1, 1, f1, 1) 386 F (_ZGVnN4v_exp2f, Z_exp2f, exp2, mpfr_exp2, 1, 1, f1, 1) 387 F (_ZGVnN4v_logf, Z_logf, log, mpfr_log, 1, 1, f1, 1) 388 F (_ZGVnN4vv_powf, Z_powf, pow, mpfr_pow, 2, 1, f2, 1) 389 F (_ZGVnN2v_sin, Z_sin, sinl, mpfr_sin, 1, 0, d1, 1) 390 F (_ZGVnN2v_cos, Z_cos, cosl, mpfr_cos, 1, 0, d1, 1) 391 F (_ZGVnN2v_exp, Z_exp, expl, mpfr_exp, 1, 0, d1, 1) 392 F (_ZGVnN2v_log, Z_log, logl, mpfr_log, 1, 0, d1, 1) 393 F (_ZGVnN2vv_pow, Z_pow, powl, mpfr_pow, 2, 0, d2, 1) 394 #endif 395 #endif 396 #endif 397 #undef F 398 #undef F1 399 #undef F2 400 #undef D1 401 #undef D2 402 {0}}; 403 404 /* Boilerplate for generic calls. */ 405 406 static inline int 407 ulpscale_f (float x) 408 { 409 int e = asuint (x) >> 23 & 0xff; 410 if (!e) 411 e++; 412 return e - 0x7f - 23; 413 } 414 static inline int 415 ulpscale_d (double x) 416 { 417 int e = asuint64 (x) >> 52 & 0x7ff; 418 if (!e) 419 e++; 420 return e - 0x3ff - 52; 421 } 422 static inline float 423 call_f1 (const struct fun *f, struct args_f1 a) 424 { 425 return f->fun.f1 (a.x); 426 } 427 static inline float 428 call_f2 (const struct fun *f, struct args_f2 a) 429 { 430 return f->fun.f2 (a.x, a.x2); 431 } 432 433 static inline double 434 call_d1 (const struct fun *f, struct args_d1 a) 435 { 436 return f->fun.d1 (a.x); 437 } 438 static inline double 439 call_d2 (const struct fun *f, struct args_d2 a) 440 { 441 return f->fun.d2 (a.x, a.x2); 442 } 443 static inline double 444 call_long_f1 (const struct fun *f, struct args_f1 a) 445 { 446 return f->fun_long.f1 (a.x); 447 } 448 static inline double 449 call_long_f2 (const struct fun *f, struct args_f2 a) 450 { 451 return f->fun_long.f2 (a.x, a.x2); 452 } 453 static inline long double 454 call_long_d1 (const struct fun *f, struct args_d1 a) 455 { 456 return f->fun_long.d1 (a.x); 457 } 458 static inline long double 459 call_long_d2 (const struct fun *f, struct args_d2 a) 460 { 461 return f->fun_long.d2 (a.x, a.x2); 462 } 463 static inline void 464 printcall_f1 (const struct fun *f, struct args_f1 a) 465 { 466 printf ("%s(%a)", f->name, a.x); 467 } 468 static inline void 469 printcall_f2 (const struct fun *f, struct args_f2 a) 470 { 471 printf ("%s(%a, %a)", f->name, a.x, a.x2); 472 } 473 static inline void 474 printcall_d1 (const struct fun *f, struct args_d1 a) 475 { 476 printf ("%s(%a)", f->name, a.x); 477 } 478 static inline void 479 printcall_d2 (const struct fun *f, struct args_d2 a) 480 { 481 printf ("%s(%a, %a)", f->name, a.x, a.x2); 482 } 483 static inline void 484 printgen_f1 (const struct fun *f, struct gen *gen) 485 { 486 printf ("%s in [%a;%a]", f->name, asfloat (gen->start), 487 asfloat (gen->start + gen->len)); 488 } 489 static inline void 490 printgen_f2 (const struct fun *f, struct gen *gen) 491 { 492 printf ("%s in [%a;%a] x [%a;%a]", f->name, asfloat (gen->start), 493 asfloat (gen->start + gen->len), asfloat (gen->start2), 494 asfloat (gen->start2 + gen->len2)); 495 } 496 static inline void 497 printgen_d1 (const struct fun *f, struct gen *gen) 498 { 499 printf ("%s in [%a;%a]", f->name, asdouble (gen->start), 500 asdouble (gen->start + gen->len)); 501 } 502 static inline void 503 printgen_d2 (const struct fun *f, struct gen *gen) 504 { 505 printf ("%s in [%a;%a] x [%a;%a]", f->name, asdouble (gen->start), 506 asdouble (gen->start + gen->len), asdouble (gen->start2), 507 asdouble (gen->start2 + gen->len2)); 508 } 509 510 #define reduce_f1(a, f, op) (f (a.x)) 511 #define reduce_f2(a, f, op) (f (a.x) op f (a.x2)) 512 #define reduce_d1(a, f, op) (f (a.x)) 513 #define reduce_d2(a, f, op) (f (a.x) op f (a.x2)) 514 515 #ifndef IEEE_754_2008_SNAN 516 # define IEEE_754_2008_SNAN 1 517 #endif 518 static inline int 519 issignaling_f (float x) 520 { 521 uint32_t ix = asuint (x); 522 if (!IEEE_754_2008_SNAN) 523 return (ix & 0x7fc00000) == 0x7fc00000; 524 return 2 * (ix ^ 0x00400000) > 2u * 0x7fc00000; 525 } 526 static inline int 527 issignaling_d (double x) 528 { 529 uint64_t ix = asuint64 (x); 530 if (!IEEE_754_2008_SNAN) 531 return (ix & 0x7ff8000000000000) == 0x7ff8000000000000; 532 return 2 * (ix ^ 0x0008000000000000) > 2 * 0x7ff8000000000000ULL; 533 } 534 535 #if USE_MPFR 536 static mpfr_rnd_t 537 rmap (int r) 538 { 539 switch (r) 540 { 541 case FE_TONEAREST: 542 return MPFR_RNDN; 543 case FE_TOWARDZERO: 544 return MPFR_RNDZ; 545 case FE_UPWARD: 546 return MPFR_RNDU; 547 case FE_DOWNWARD: 548 return MPFR_RNDD; 549 } 550 return -1; 551 } 552 553 #define prec_mpfr_f 50 554 #define prec_mpfr_d 80 555 #define prec_f 24 556 #define prec_d 53 557 #define emin_f -148 558 #define emin_d -1073 559 #define emax_f 128 560 #define emax_d 1024 561 static inline int 562 call_mpfr_f1 (mpfr_t y, const struct fun *f, struct args_f1 a, mpfr_rnd_t r) 563 { 564 MPFR_DECL_INIT (x, prec_f); 565 mpfr_set_flt (x, a.x, MPFR_RNDN); 566 return f->fun_mpfr.f1 (y, x, r); 567 } 568 static inline int 569 call_mpfr_f2 (mpfr_t y, const struct fun *f, struct args_f2 a, mpfr_rnd_t r) 570 { 571 MPFR_DECL_INIT (x, prec_f); 572 MPFR_DECL_INIT (x2, prec_f); 573 mpfr_set_flt (x, a.x, MPFR_RNDN); 574 mpfr_set_flt (x2, a.x2, MPFR_RNDN); 575 return f->fun_mpfr.f2 (y, x, x2, r); 576 } 577 static inline int 578 call_mpfr_d1 (mpfr_t y, const struct fun *f, struct args_d1 a, mpfr_rnd_t r) 579 { 580 MPFR_DECL_INIT (x, prec_d); 581 mpfr_set_d (x, a.x, MPFR_RNDN); 582 return f->fun_mpfr.d1 (y, x, r); 583 } 584 static inline int 585 call_mpfr_d2 (mpfr_t y, const struct fun *f, struct args_d2 a, mpfr_rnd_t r) 586 { 587 MPFR_DECL_INIT (x, prec_d); 588 MPFR_DECL_INIT (x2, prec_d); 589 mpfr_set_d (x, a.x, MPFR_RNDN); 590 mpfr_set_d (x2, a.x2, MPFR_RNDN); 591 return f->fun_mpfr.d2 (y, x, x2, r); 592 } 593 #endif 594 595 #define float_f float 596 #define double_f double 597 #define copysign_f copysignf 598 #define nextafter_f nextafterf 599 #define fabs_f fabsf 600 #define asuint_f asuint 601 #define asfloat_f asfloat 602 #define scalbn_f scalbnf 603 #define lscalbn_f scalbn 604 #define halfinf_f 0x1p127f 605 #define min_normal_f 0x1p-126f 606 607 #define float_d double 608 #define double_d long double 609 #define copysign_d copysign 610 #define nextafter_d nextafter 611 #define fabs_d fabs 612 #define asuint_d asuint64 613 #define asfloat_d asdouble 614 #define scalbn_d scalbn 615 #define lscalbn_d scalbnl 616 #define halfinf_d 0x1p1023 617 #define min_normal_d 0x1p-1022 618 619 #define NEW_RT 620 #define RT(x) x##_f 621 #define T(x) x##_f1 622 #include "ulp.h" 623 #undef T 624 #define T(x) x##_f2 625 #include "ulp.h" 626 #undef T 627 #undef RT 628 629 #define NEW_RT 630 #define RT(x) x##_d 631 #define T(x) x##_d1 632 #include "ulp.h" 633 #undef T 634 #define T(x) x##_d2 635 #include "ulp.h" 636 #undef T 637 #undef RT 638 639 static void 640 usage (void) 641 { 642 puts ("./ulp [-q] [-m] [-f] [-r nudz] [-l soft-ulplimit] [-e ulplimit] func " 643 "lo [hi [x lo2 hi2] [count]]"); 644 puts ("Compares func against a higher precision implementation in [lo; hi]."); 645 puts ("-q: quiet."); 646 puts ("-m: use mpfr even if faster method is available."); 647 puts ("-f: disable fenv testing (rounding modes and exceptions)."); 648 puts ("Supported func:"); 649 for (const struct fun *f = fun; f->name; f++) 650 printf ("\t%s\n", f->name); 651 exit (1); 652 } 653 654 static int 655 cmp (const struct fun *f, struct gen *gen, const struct conf *conf) 656 { 657 int r = 1; 658 if (f->arity == 1 && f->singleprec) 659 r = cmp_f1 (f, gen, conf); 660 else if (f->arity == 2 && f->singleprec) 661 r = cmp_f2 (f, gen, conf); 662 else if (f->arity == 1 && !f->singleprec) 663 r = cmp_d1 (f, gen, conf); 664 else if (f->arity == 2 && !f->singleprec) 665 r = cmp_d2 (f, gen, conf); 666 else 667 usage (); 668 return r; 669 } 670 671 static uint64_t 672 getnum (const char *s, int singleprec) 673 { 674 // int i; 675 uint64_t sign = 0; 676 // char buf[12]; 677 678 if (s[0] == '+') 679 s++; 680 else if (s[0] == '-') 681 { 682 sign = singleprec ? 1ULL << 31 : 1ULL << 63; 683 s++; 684 } 685 /* 0xXXXX is treated as bit representation, '-' flips the sign bit. */ 686 if (s[0] == '0' && tolower (s[1]) == 'x' && strchr (s, 'p') == 0) 687 return sign ^ strtoull (s, 0, 0); 688 // /* SNaN, QNaN, NaN, Inf. */ 689 // for (i=0; s[i] && i < sizeof buf; i++) 690 // buf[i] = tolower(s[i]); 691 // buf[i] = 0; 692 // if (strcmp(buf, "snan") == 0) 693 // return sign | (singleprec ? 0x7fa00000 : 0x7ff4000000000000); 694 // if (strcmp(buf, "qnan") == 0 || strcmp(buf, "nan") == 0) 695 // return sign | (singleprec ? 0x7fc00000 : 0x7ff8000000000000); 696 // if (strcmp(buf, "inf") == 0 || strcmp(buf, "infinity") == 0) 697 // return sign | (singleprec ? 0x7f800000 : 0x7ff0000000000000); 698 /* Otherwise assume it's a floating-point literal. */ 699 return sign 700 | (singleprec ? asuint (strtof (s, 0)) : asuint64 (strtod (s, 0))); 701 } 702 703 static void 704 parsegen (struct gen *g, int argc, char *argv[], const struct fun *f) 705 { 706 int singleprec = f->singleprec; 707 int arity = f->arity; 708 uint64_t a, b, a2, b2, n; 709 if (argc < 1) 710 usage (); 711 b = a = getnum (argv[0], singleprec); 712 n = 0; 713 if (argc > 1 && strcmp (argv[1], "x") == 0) 714 { 715 argc -= 2; 716 argv += 2; 717 } 718 else if (argc > 1) 719 { 720 b = getnum (argv[1], singleprec); 721 if (argc > 2 && strcmp (argv[2], "x") == 0) 722 { 723 argc -= 3; 724 argv += 3; 725 } 726 } 727 b2 = a2 = getnum (argv[0], singleprec); 728 if (argc > 1) 729 b2 = getnum (argv[1], singleprec); 730 if (argc > 2) 731 n = strtoull (argv[2], 0, 0); 732 if (argc > 3) 733 usage (); 734 //printf("ab %lx %lx ab2 %lx %lx n %lu\n", a, b, a2, b2, n); 735 if (arity == 1) 736 { 737 g->start = a; 738 g->len = b - a; 739 if (n - 1 > b - a) 740 n = b - a + 1; 741 g->off = 0; 742 g->step = n ? (g->len + 1) / n : 1; 743 g->start2 = g->len2 = 0; 744 g->cnt = n; 745 } 746 else if (arity == 2) 747 { 748 g->start = a; 749 g->len = b - a; 750 g->off = g->step = 0; 751 g->start2 = a2; 752 g->len2 = b2 - a2; 753 g->cnt = n; 754 } 755 else 756 usage (); 757 } 758 759 int 760 main (int argc, char *argv[]) 761 { 762 const struct fun *f; 763 struct gen gen; 764 struct conf conf; 765 conf.rc = 'n'; 766 conf.quiet = 0; 767 conf.mpfr = 0; 768 conf.fenv = 1; 769 conf.softlim = 0; 770 conf.errlim = INFINITY; 771 for (;;) 772 { 773 argc--; 774 argv++; 775 if (argc < 1) 776 usage (); 777 if (argv[0][0] != '-') 778 break; 779 switch (argv[0][1]) 780 { 781 case 'e': 782 argc--; 783 argv++; 784 if (argc < 1) 785 usage (); 786 conf.errlim = strtod (argv[0], 0); 787 break; 788 case 'f': 789 conf.fenv = 0; 790 break; 791 case 'l': 792 argc--; 793 argv++; 794 if (argc < 1) 795 usage (); 796 conf.softlim = strtod (argv[0], 0); 797 break; 798 case 'm': 799 conf.mpfr = 1; 800 break; 801 case 'q': 802 conf.quiet = 1; 803 break; 804 case 'r': 805 conf.rc = argv[0][2]; 806 if (!conf.rc) 807 { 808 argc--; 809 argv++; 810 if (argc < 1) 811 usage (); 812 conf.rc = argv[0][0]; 813 } 814 break; 815 default: 816 usage (); 817 } 818 } 819 switch (conf.rc) 820 { 821 case 'n': 822 conf.r = FE_TONEAREST; 823 break; 824 case 'u': 825 conf.r = FE_UPWARD; 826 break; 827 case 'd': 828 conf.r = FE_DOWNWARD; 829 break; 830 case 'z': 831 conf.r = FE_TOWARDZERO; 832 break; 833 default: 834 usage (); 835 } 836 for (f = fun; f->name; f++) 837 if (strcmp (argv[0], f->name) == 0) 838 break; 839 if (!f->name) 840 usage (); 841 if (!f->singleprec && LDBL_MANT_DIG == DBL_MANT_DIG) 842 conf.mpfr = 1; /* Use mpfr if long double has no extra precision. */ 843 if (!USE_MPFR && conf.mpfr) 844 { 845 puts ("mpfr is not available."); 846 return 0; 847 } 848 argc--; 849 argv++; 850 parsegen (&gen, argc, argv, f); 851 conf.n = gen.cnt; 852 return cmp (f, &gen, &conf); 853 } 854