1 /* 2 * ULP error checking tool for math functions. 3 * 4 * Copyright (c) 2019, 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 D1 (exp) 335 D1 (exp2) 336 D1 (log) 337 D1 (log2) 338 D2 (pow) 339 #if WANT_VMATH 340 F (__s_sinf, __s_sinf, sin, mpfr_sin, 1, 1, f1, 0) 341 F (__s_cosf, __s_cosf, cos, mpfr_cos, 1, 1, f1, 0) 342 F (__s_expf_1u, __s_expf_1u, exp, mpfr_exp, 1, 1, f1, 0) 343 F (__s_expf, __s_expf, exp, mpfr_exp, 1, 1, f1, 0) 344 F (__s_exp2f_1u, __s_exp2f_1u, exp2, mpfr_exp2, 1, 1, f1, 0) 345 F (__s_exp2f, __s_exp2f, exp2, mpfr_exp2, 1, 1, f1, 0) 346 F (__s_powf, __s_powf, pow, mpfr_pow, 2, 1, f2, 0) 347 F (__s_logf, __s_logf, log, mpfr_log, 1, 1, f1, 0) 348 F (__s_sin, __s_sin, sinl, mpfr_sin, 1, 0, d1, 0) 349 F (__s_cos, __s_cos, cosl, mpfr_cos, 1, 0, d1, 0) 350 F (__s_exp, __s_exp, expl, mpfr_exp, 1, 0, d1, 0) 351 F (__s_log, __s_log, logl, mpfr_log, 1, 0, d1, 0) 352 F (__s_pow, __s_pow, powl, mpfr_pow, 2, 0, d2, 0) 353 #if __aarch64__ 354 F (__v_sinf, v_sinf, sin, mpfr_sin, 1, 1, f1, 1) 355 F (__v_cosf, v_cosf, cos, mpfr_cos, 1, 1, f1, 1) 356 F (__v_expf_1u, v_expf_1u, exp, mpfr_exp, 1, 1, f1, 1) 357 F (__v_expf, v_expf, exp, mpfr_exp, 1, 1, f1, 1) 358 F (__v_exp2f_1u, v_exp2f_1u, exp2, mpfr_exp2, 1, 1, f1, 1) 359 F (__v_exp2f, v_exp2f, exp2, mpfr_exp2, 1, 1, f1, 1) 360 F (__v_logf, v_logf, log, mpfr_log, 1, 1, f1, 1) 361 F (__v_powf, v_powf, pow, mpfr_pow, 2, 1, f2, 1) 362 F (__v_sin, v_sin, sinl, mpfr_sin, 1, 0, d1, 1) 363 F (__v_cos, v_cos, cosl, mpfr_cos, 1, 0, d1, 1) 364 F (__v_exp, v_exp, expl, mpfr_exp, 1, 0, d1, 1) 365 F (__v_log, v_log, logl, mpfr_log, 1, 0, d1, 1) 366 F (__v_pow, v_pow, powl, mpfr_pow, 2, 0, d2, 1) 367 #ifdef __vpcs 368 F (__vn_sinf, vn_sinf, sin, mpfr_sin, 1, 1, f1, 1) 369 F (__vn_cosf, vn_cosf, cos, mpfr_cos, 1, 1, f1, 1) 370 F (__vn_expf_1u, vn_expf_1u, exp, mpfr_exp, 1, 1, f1, 1) 371 F (__vn_expf, vn_expf, exp, mpfr_exp, 1, 1, f1, 1) 372 F (__vn_exp2f_1u, vn_exp2f_1u, exp2, mpfr_exp2, 1, 1, f1, 1) 373 F (__vn_exp2f, vn_exp2f, exp2, mpfr_exp2, 1, 1, f1, 1) 374 F (__vn_logf, vn_logf, log, mpfr_log, 1, 1, f1, 1) 375 F (__vn_powf, vn_powf, pow, mpfr_pow, 2, 1, f2, 1) 376 F (__vn_sin, vn_sin, sinl, mpfr_sin, 1, 0, d1, 1) 377 F (__vn_cos, vn_cos, cosl, mpfr_cos, 1, 0, d1, 1) 378 F (__vn_exp, vn_exp, expl, mpfr_exp, 1, 0, d1, 1) 379 F (__vn_log, vn_log, logl, mpfr_log, 1, 0, d1, 1) 380 F (__vn_pow, vn_pow, powl, mpfr_pow, 2, 0, d2, 1) 381 F (_ZGVnN4v_sinf, Z_sinf, sin, mpfr_sin, 1, 1, f1, 1) 382 F (_ZGVnN4v_cosf, Z_cosf, cos, mpfr_cos, 1, 1, f1, 1) 383 F (_ZGVnN4v_expf, Z_expf, exp, mpfr_exp, 1, 1, f1, 1) 384 F (_ZGVnN4v_exp2f, Z_exp2f, exp2, mpfr_exp2, 1, 1, f1, 1) 385 F (_ZGVnN4v_logf, Z_logf, log, mpfr_log, 1, 1, f1, 1) 386 F (_ZGVnN4vv_powf, Z_powf, pow, mpfr_pow, 2, 1, f2, 1) 387 F (_ZGVnN2v_sin, Z_sin, sinl, mpfr_sin, 1, 0, d1, 1) 388 F (_ZGVnN2v_cos, Z_cos, cosl, mpfr_cos, 1, 0, d1, 1) 389 F (_ZGVnN2v_exp, Z_exp, expl, mpfr_exp, 1, 0, d1, 1) 390 F (_ZGVnN2v_log, Z_log, logl, mpfr_log, 1, 0, d1, 1) 391 F (_ZGVnN2vv_pow, Z_pow, powl, mpfr_pow, 2, 0, d2, 1) 392 #endif 393 #endif 394 #endif 395 #undef F 396 #undef F1 397 #undef F2 398 #undef D1 399 #undef D2 400 {0}}; 401 402 /* Boilerplate for generic calls. */ 403 404 static inline int 405 ulpscale_f (float x) 406 { 407 int e = asuint (x) >> 23 & 0xff; 408 if (!e) 409 e++; 410 return e - 0x7f - 23; 411 } 412 static inline int 413 ulpscale_d (double x) 414 { 415 int e = asuint64 (x) >> 52 & 0x7ff; 416 if (!e) 417 e++; 418 return e - 0x3ff - 52; 419 } 420 static inline float 421 call_f1 (const struct fun *f, struct args_f1 a) 422 { 423 return f->fun.f1 (a.x); 424 } 425 static inline float 426 call_f2 (const struct fun *f, struct args_f2 a) 427 { 428 return f->fun.f2 (a.x, a.x2); 429 } 430 431 static inline double 432 call_d1 (const struct fun *f, struct args_d1 a) 433 { 434 return f->fun.d1 (a.x); 435 } 436 static inline double 437 call_d2 (const struct fun *f, struct args_d2 a) 438 { 439 return f->fun.d2 (a.x, a.x2); 440 } 441 static inline double 442 call_long_f1 (const struct fun *f, struct args_f1 a) 443 { 444 return f->fun_long.f1 (a.x); 445 } 446 static inline double 447 call_long_f2 (const struct fun *f, struct args_f2 a) 448 { 449 return f->fun_long.f2 (a.x, a.x2); 450 } 451 static inline long double 452 call_long_d1 (const struct fun *f, struct args_d1 a) 453 { 454 return f->fun_long.d1 (a.x); 455 } 456 static inline long double 457 call_long_d2 (const struct fun *f, struct args_d2 a) 458 { 459 return f->fun_long.d2 (a.x, a.x2); 460 } 461 static inline void 462 printcall_f1 (const struct fun *f, struct args_f1 a) 463 { 464 printf ("%s(%a)", f->name, a.x); 465 } 466 static inline void 467 printcall_f2 (const struct fun *f, struct args_f2 a) 468 { 469 printf ("%s(%a, %a)", f->name, a.x, a.x2); 470 } 471 static inline void 472 printcall_d1 (const struct fun *f, struct args_d1 a) 473 { 474 printf ("%s(%a)", f->name, a.x); 475 } 476 static inline void 477 printcall_d2 (const struct fun *f, struct args_d2 a) 478 { 479 printf ("%s(%a, %a)", f->name, a.x, a.x2); 480 } 481 static inline void 482 printgen_f1 (const struct fun *f, struct gen *gen) 483 { 484 printf ("%s in [%a;%a]", f->name, asfloat (gen->start), 485 asfloat (gen->start + gen->len)); 486 } 487 static inline void 488 printgen_f2 (const struct fun *f, struct gen *gen) 489 { 490 printf ("%s in [%a;%a] x [%a;%a]", f->name, asfloat (gen->start), 491 asfloat (gen->start + gen->len), asfloat (gen->start2), 492 asfloat (gen->start2 + gen->len2)); 493 } 494 static inline void 495 printgen_d1 (const struct fun *f, struct gen *gen) 496 { 497 printf ("%s in [%a;%a]", f->name, asdouble (gen->start), 498 asdouble (gen->start + gen->len)); 499 } 500 static inline void 501 printgen_d2 (const struct fun *f, struct gen *gen) 502 { 503 printf ("%s in [%a;%a] x [%a;%a]", f->name, asdouble (gen->start), 504 asdouble (gen->start + gen->len), asdouble (gen->start2), 505 asdouble (gen->start2 + gen->len2)); 506 } 507 508 #define reduce_f1(a, f, op) (f (a.x)) 509 #define reduce_f2(a, f, op) (f (a.x) op f (a.x2)) 510 #define reduce_d1(a, f, op) (f (a.x)) 511 #define reduce_d2(a, f, op) (f (a.x) op f (a.x2)) 512 513 #ifndef IEEE_754_2008_SNAN 514 # define IEEE_754_2008_SNAN 1 515 #endif 516 static inline int 517 issignaling_f (float x) 518 { 519 uint32_t ix = asuint (x); 520 if (!IEEE_754_2008_SNAN) 521 return (ix & 0x7fc00000) == 0x7fc00000; 522 return 2 * (ix ^ 0x00400000) > 2u * 0x7fc00000; 523 } 524 static inline int 525 issignaling_d (double x) 526 { 527 uint64_t ix = asuint64 (x); 528 if (!IEEE_754_2008_SNAN) 529 return (ix & 0x7ff8000000000000) == 0x7ff8000000000000; 530 return 2 * (ix ^ 0x0008000000000000) > 2 * 0x7ff8000000000000ULL; 531 } 532 533 #if USE_MPFR 534 static mpfr_rnd_t 535 rmap (int r) 536 { 537 switch (r) 538 { 539 case FE_TONEAREST: 540 return MPFR_RNDN; 541 case FE_TOWARDZERO: 542 return MPFR_RNDZ; 543 case FE_UPWARD: 544 return MPFR_RNDU; 545 case FE_DOWNWARD: 546 return MPFR_RNDD; 547 } 548 return -1; 549 } 550 551 #define prec_mpfr_f 50 552 #define prec_mpfr_d 80 553 #define prec_f 24 554 #define prec_d 53 555 #define emin_f -148 556 #define emin_d -1073 557 #define emax_f 128 558 #define emax_d 1024 559 static inline int 560 call_mpfr_f1 (mpfr_t y, const struct fun *f, struct args_f1 a, mpfr_rnd_t r) 561 { 562 MPFR_DECL_INIT (x, prec_f); 563 mpfr_set_flt (x, a.x, MPFR_RNDN); 564 return f->fun_mpfr.f1 (y, x, r); 565 } 566 static inline int 567 call_mpfr_f2 (mpfr_t y, const struct fun *f, struct args_f2 a, mpfr_rnd_t r) 568 { 569 MPFR_DECL_INIT (x, prec_f); 570 MPFR_DECL_INIT (x2, prec_f); 571 mpfr_set_flt (x, a.x, MPFR_RNDN); 572 mpfr_set_flt (x2, a.x2, MPFR_RNDN); 573 return f->fun_mpfr.f2 (y, x, x2, r); 574 } 575 static inline int 576 call_mpfr_d1 (mpfr_t y, const struct fun *f, struct args_d1 a, mpfr_rnd_t r) 577 { 578 MPFR_DECL_INIT (x, prec_d); 579 mpfr_set_d (x, a.x, MPFR_RNDN); 580 return f->fun_mpfr.d1 (y, x, r); 581 } 582 static inline int 583 call_mpfr_d2 (mpfr_t y, const struct fun *f, struct args_d2 a, mpfr_rnd_t r) 584 { 585 MPFR_DECL_INIT (x, prec_d); 586 MPFR_DECL_INIT (x2, prec_d); 587 mpfr_set_d (x, a.x, MPFR_RNDN); 588 mpfr_set_d (x2, a.x2, MPFR_RNDN); 589 return f->fun_mpfr.d2 (y, x, x2, r); 590 } 591 #endif 592 593 #define float_f float 594 #define double_f double 595 #define copysign_f copysignf 596 #define nextafter_f nextafterf 597 #define fabs_f fabsf 598 #define asuint_f asuint 599 #define asfloat_f asfloat 600 #define scalbn_f scalbnf 601 #define lscalbn_f scalbn 602 #define halfinf_f 0x1p127f 603 #define min_normal_f 0x1p-126f 604 605 #define float_d double 606 #define double_d long double 607 #define copysign_d copysign 608 #define nextafter_d nextafter 609 #define fabs_d fabs 610 #define asuint_d asuint64 611 #define asfloat_d asdouble 612 #define scalbn_d scalbn 613 #define lscalbn_d scalbnl 614 #define halfinf_d 0x1p1023 615 #define min_normal_d 0x1p-1022 616 617 #define NEW_RT 618 #define RT(x) x##_f 619 #define T(x) x##_f1 620 #include "ulp.h" 621 #undef T 622 #define T(x) x##_f2 623 #include "ulp.h" 624 #undef T 625 #undef RT 626 627 #define NEW_RT 628 #define RT(x) x##_d 629 #define T(x) x##_d1 630 #include "ulp.h" 631 #undef T 632 #define T(x) x##_d2 633 #include "ulp.h" 634 #undef T 635 #undef RT 636 637 static void 638 usage (void) 639 { 640 puts ("./ulp [-q] [-m] [-f] [-r nudz] [-l soft-ulplimit] [-e ulplimit] func " 641 "lo [hi [x lo2 hi2] [count]]"); 642 puts ("Compares func against a higher precision implementation in [lo; hi]."); 643 puts ("-q: quiet."); 644 puts ("-m: use mpfr even if faster method is available."); 645 puts ("-f: disable fenv testing (rounding modes and exceptions)."); 646 puts ("Supported func:"); 647 for (const struct fun *f = fun; f->name; f++) 648 printf ("\t%s\n", f->name); 649 exit (1); 650 } 651 652 static int 653 cmp (const struct fun *f, struct gen *gen, const struct conf *conf) 654 { 655 int r = 1; 656 if (f->arity == 1 && f->singleprec) 657 r = cmp_f1 (f, gen, conf); 658 else if (f->arity == 2 && f->singleprec) 659 r = cmp_f2 (f, gen, conf); 660 else if (f->arity == 1 && !f->singleprec) 661 r = cmp_d1 (f, gen, conf); 662 else if (f->arity == 2 && !f->singleprec) 663 r = cmp_d2 (f, gen, conf); 664 else 665 usage (); 666 return r; 667 } 668 669 static uint64_t 670 getnum (const char *s, int singleprec) 671 { 672 // int i; 673 uint64_t sign = 0; 674 // char buf[12]; 675 676 if (s[0] == '+') 677 s++; 678 else if (s[0] == '-') 679 { 680 sign = singleprec ? 1ULL << 31 : 1ULL << 63; 681 s++; 682 } 683 /* 0xXXXX is treated as bit representation, '-' flips the sign bit. */ 684 if (s[0] == '0' && tolower (s[1]) == 'x' && strchr (s, 'p') == 0) 685 return sign ^ strtoull (s, 0, 0); 686 // /* SNaN, QNaN, NaN, Inf. */ 687 // for (i=0; s[i] && i < sizeof buf; i++) 688 // buf[i] = tolower(s[i]); 689 // buf[i] = 0; 690 // if (strcmp(buf, "snan") == 0) 691 // return sign | (singleprec ? 0x7fa00000 : 0x7ff4000000000000); 692 // if (strcmp(buf, "qnan") == 0 || strcmp(buf, "nan") == 0) 693 // return sign | (singleprec ? 0x7fc00000 : 0x7ff8000000000000); 694 // if (strcmp(buf, "inf") == 0 || strcmp(buf, "infinity") == 0) 695 // return sign | (singleprec ? 0x7f800000 : 0x7ff0000000000000); 696 /* Otherwise assume it's a floating-point literal. */ 697 return sign 698 | (singleprec ? asuint (strtof (s, 0)) : asuint64 (strtod (s, 0))); 699 } 700 701 static void 702 parsegen (struct gen *g, int argc, char *argv[], const struct fun *f) 703 { 704 int singleprec = f->singleprec; 705 int arity = f->arity; 706 uint64_t a, b, a2, b2, n; 707 if (argc < 1) 708 usage (); 709 b = a = getnum (argv[0], singleprec); 710 n = 0; 711 if (argc > 1 && strcmp (argv[1], "x") == 0) 712 { 713 argc -= 2; 714 argv += 2; 715 } 716 else if (argc > 1) 717 { 718 b = getnum (argv[1], singleprec); 719 if (argc > 2 && strcmp (argv[2], "x") == 0) 720 { 721 argc -= 3; 722 argv += 3; 723 } 724 } 725 b2 = a2 = getnum (argv[0], singleprec); 726 if (argc > 1) 727 b2 = getnum (argv[1], singleprec); 728 if (argc > 2) 729 n = strtoull (argv[2], 0, 0); 730 if (argc > 3) 731 usage (); 732 //printf("ab %lx %lx ab2 %lx %lx n %lu\n", a, b, a2, b2, n); 733 if (arity == 1) 734 { 735 g->start = a; 736 g->len = b - a; 737 if (n - 1 > b - a) 738 n = b - a + 1; 739 g->off = 0; 740 g->step = n ? (g->len + 1) / n : 1; 741 g->start2 = g->len2 = 0; 742 g->cnt = n; 743 } 744 else if (arity == 2) 745 { 746 g->start = a; 747 g->len = b - a; 748 g->off = g->step = 0; 749 g->start2 = a2; 750 g->len2 = b2 - a2; 751 g->cnt = n; 752 } 753 else 754 usage (); 755 } 756 757 int 758 main (int argc, char *argv[]) 759 { 760 const struct fun *f; 761 struct gen gen; 762 struct conf conf; 763 conf.rc = 'n'; 764 conf.quiet = 0; 765 conf.mpfr = 0; 766 conf.fenv = 1; 767 conf.softlim = 0; 768 conf.errlim = INFINITY; 769 for (;;) 770 { 771 argc--; 772 argv++; 773 if (argc < 1) 774 usage (); 775 if (argv[0][0] != '-') 776 break; 777 switch (argv[0][1]) 778 { 779 case 'e': 780 argc--; 781 argv++; 782 if (argc < 1) 783 usage (); 784 conf.errlim = strtod (argv[0], 0); 785 break; 786 case 'f': 787 conf.fenv = 0; 788 break; 789 case 'l': 790 argc--; 791 argv++; 792 if (argc < 1) 793 usage (); 794 conf.softlim = strtod (argv[0], 0); 795 break; 796 case 'm': 797 conf.mpfr = 1; 798 break; 799 case 'q': 800 conf.quiet = 1; 801 break; 802 case 'r': 803 conf.rc = argv[0][2]; 804 if (!conf.rc) 805 { 806 argc--; 807 argv++; 808 if (argc < 1) 809 usage (); 810 conf.rc = argv[0][0]; 811 } 812 break; 813 default: 814 usage (); 815 } 816 } 817 switch (conf.rc) 818 { 819 case 'n': 820 conf.r = FE_TONEAREST; 821 break; 822 case 'u': 823 conf.r = FE_UPWARD; 824 break; 825 case 'd': 826 conf.r = FE_DOWNWARD; 827 break; 828 case 'z': 829 conf.r = FE_TOWARDZERO; 830 break; 831 default: 832 usage (); 833 } 834 for (f = fun; f->name; f++) 835 if (strcmp (argv[0], f->name) == 0) 836 break; 837 if (!f->name) 838 usage (); 839 if (!f->singleprec && LDBL_MANT_DIG == DBL_MANT_DIG) 840 conf.mpfr = 1; /* Use mpfr if long double has no extra precision. */ 841 if (!USE_MPFR && conf.mpfr) 842 { 843 puts ("mpfr is not available."); 844 return 0; 845 } 846 argc--; 847 argv++; 848 parsegen (&gen, argc, argv, f); 849 conf.n = gen.cnt; 850 return cmp (f, &gen, &conf); 851 } 852