1 /* 2 * ULP error checking tool for math functions. 3 * 4 * Copyright (c) 2019-2022, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 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 /* A bit of a hack: call vector functions twice with the same 218 input in lane 0 but a different value in other lanes: once 219 with an in-range value and then with a special case value. */ 220 static int secondcall; 221 222 /* Wrappers for vector functions. */ 223 #if __aarch64__ && WANT_VMATH 224 typedef __f32x4_t v_float; 225 typedef __f64x2_t v_double; 226 /* First element of fv and dv may be changed by -c argument. */ 227 static float fv[2] = {1.0f, -INFINITY}; 228 static double dv[2] = {1.0, -INFINITY}; 229 static inline v_float argf(float x) { return (v_float){x,x,x,fv[secondcall]}; } 230 static inline v_double argd(double x) { return (v_double){x,dv[secondcall]}; } 231 #if WANT_SVE_MATH 232 #include <arm_sve.h> 233 typedef __SVFloat32_t sv_float; 234 typedef __SVFloat64_t sv_double; 235 236 static inline sv_float svargf(float x) { 237 int n = svcntw(); 238 float base[n]; 239 for (int i=0; i<n; i++) 240 base[i] = (float)x; 241 base[n-1] = (float) fv[secondcall]; 242 return svld1(svptrue_b32(), base); 243 } 244 static inline sv_double svargd(double x) { 245 int n = svcntd(); 246 double base[n]; 247 for (int i=0; i<n; i++) 248 base[i] = x; 249 base[n-1] = dv[secondcall]; 250 return svld1(svptrue_b64(), base); 251 } 252 static inline float svretf(sv_float vec) { 253 int n = svcntw(); 254 float res[n]; 255 svst1(svptrue_b32(), res, vec); 256 return res[0]; 257 } 258 static inline double svretd(sv_double vec) { 259 int n = svcntd(); 260 double res[n]; 261 svst1(svptrue_b64(), res, vec); 262 return res[0]; 263 } 264 #endif 265 #endif 266 267 #if WANT_SVE_MATH 268 long double 269 dummyl (long double x) 270 { 271 return x; 272 } 273 274 double 275 dummy (double x) 276 { 277 return x; 278 } 279 280 static sv_double 281 __sv_dummy (sv_double x) 282 { 283 return x; 284 } 285 286 static sv_float 287 __sv_dummyf (sv_float x) 288 { 289 return x; 290 } 291 #endif 292 293 #include "test/ulp_wrappers.h" 294 295 /* Wrappers for SVE functions. */ 296 #if WANT_SVE_MATH 297 static double sv_dummy (double x) { return svretd (__sv_dummy (svargd (x))); } 298 static float sv_dummyf (float x) { return svretf (__sv_dummyf (svargf (x))); } 299 #endif 300 301 struct fun 302 { 303 const char *name; 304 int arity; 305 int singleprec; 306 int twice; 307 union 308 { 309 float (*f1) (float); 310 float (*f2) (float, float); 311 double (*d1) (double); 312 double (*d2) (double, double); 313 } fun; 314 union 315 { 316 double (*f1) (double); 317 double (*f2) (double, double); 318 long double (*d1) (long double); 319 long double (*d2) (long double, long double); 320 } fun_long; 321 #if USE_MPFR 322 union 323 { 324 int (*f1) (mpfr_t, const mpfr_t, mpfr_rnd_t); 325 int (*f2) (mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t); 326 int (*d1) (mpfr_t, const mpfr_t, mpfr_rnd_t); 327 int (*d2) (mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t); 328 } fun_mpfr; 329 #endif 330 }; 331 332 static const struct fun fun[] = { 333 #if USE_MPFR 334 # define F(x, x_wrap, x_long, x_mpfr, a, s, t, twice) \ 335 {#x, a, s, twice, {.t = x_wrap}, {.t = x_long}, {.t = x_mpfr}}, 336 #else 337 # define F(x, x_wrap, x_long, x_mpfr, a, s, t, twice) \ 338 {#x, a, s, twice, {.t = x_wrap}, {.t = x_long}}, 339 #endif 340 #define F1(x) F (x##f, x##f, x, mpfr_##x, 1, 1, f1, 0) 341 #define F2(x) F (x##f, x##f, x, mpfr_##x, 2, 1, f2, 0) 342 #define D1(x) F (x, x, x##l, mpfr_##x, 1, 0, d1, 0) 343 #define D2(x) F (x, x, x##l, mpfr_##x, 2, 0, d2, 0) 344 /* Neon routines. */ 345 #define VF1(x) F (__v_##x##f, v_##x##f, x, mpfr_##x, 1, 1, f1, 0) 346 #define VF2(x) F (__v_##x##f, v_##x##f, x, mpfr_##x, 2, 1, f2, 0) 347 #define VD1(x) F (__v_##x, v_##x, x##l, mpfr_##x, 1, 0, d1, 0) 348 #define VD2(x) F (__v_##x, v_##x, x##l, mpfr_##x, 2, 0, d2, 0) 349 #define VNF1(x) F (__vn_##x##f, vn_##x##f, x, mpfr_##x, 1, 1, f1, 0) 350 #define VNF2(x) F (__vn_##x##f, vn_##x##f, x, mpfr_##x, 2, 1, f2, 0) 351 #define VND1(x) F (__vn_##x, vn_##x, x##l, mpfr_##x, 1, 0, d1, 0) 352 #define VND2(x) F (__vn_##x, vn_##x, x##l, mpfr_##x, 2, 0, d2, 0) 353 #define ZVF1(x) F (_ZGVnN4v_##x##f, Z_##x##f, x, mpfr_##x, 1, 1, f1, 0) 354 #define ZVF2(x) F (_ZGVnN4vv_##x##f, Z_##x##f, x, mpfr_##x, 2, 1, f2, 0) 355 #define ZVD1(x) F (_ZGVnN2v_##x, Z_##x, x##l, mpfr_##x, 1, 0, d1, 0) 356 #define ZVD2(x) F (_ZGVnN2vv_##x, Z_##x, x##l, mpfr_##x, 2, 0, d2, 0) 357 #define ZVNF1(x) VNF1 (x) ZVF1 (x) 358 #define ZVNF2(x) VNF2 (x) ZVF2 (x) 359 #define ZVND1(x) VND1 (x) ZVD1 (x) 360 #define ZVND2(x) VND2 (x) ZVD2 (x) 361 #define SF1(x) F (__s_##x##f, __s_##x##f, x, mpfr_##x, 1, 1, f1, 0) 362 #define SF2(x) F (__s_##x##f, __s_##x##f, x, mpfr_##x, 2, 1, f2, 0) 363 #define SD1(x) F (__s_##x, __s_##x, x##l, mpfr_##x, 1, 0, d1, 0) 364 #define SD2(x) F (__s_##x, __s_##x, x##l, mpfr_##x, 2, 0, d2, 0) 365 /* SVE routines. */ 366 #define SVF1(x) F (__sv_##x##f, sv_##x##f, x, mpfr_##x, 1, 1, f1, 0) 367 #define SVF2(x) F (__sv_##x##f, sv_##x##f, x, mpfr_##x, 2, 1, f2, 0) 368 #define SVD1(x) F (__sv_##x, sv_##x, x##l, mpfr_##x, 1, 0, d1, 0) 369 #define SVD2(x) F (__sv_##x, sv_##x, x##l, mpfr_##x, 2, 0, d2, 0) 370 #define ZSVF1(x) F (_ZGVsMxv_##x##f, Z_sv_##x##f, x, mpfr_##x, 1, 1, f1, 0) 371 #define ZSVF2(x) F (_ZGVsMxvv_##x##f, Z_sv_##x##f, x, mpfr_##x, 2, 1, f2, 0) 372 #define ZSVD1(x) F (_ZGVsMxv_##x, Z_sv_##x, x##l, mpfr_##x, 1, 0, d1, 0) 373 #define ZSVD2(x) F (_ZGVsMxvv_##x, Z_sv_##x, x##l, mpfr_##x, 2, 0, d2, 0) 374 375 #include "test/ulp_funcs.h" 376 377 #if WANT_SVE_MATH 378 SVD1 (dummy) 379 SVF1 (dummy) 380 #endif 381 382 #undef F 383 #undef F1 384 #undef F2 385 #undef D1 386 #undef D2 387 #undef SVF1 388 #undef SVF2 389 #undef SVD1 390 #undef SVD2 391 {0}}; 392 393 /* Boilerplate for generic calls. */ 394 395 static inline int 396 ulpscale_f (float x) 397 { 398 int e = asuint (x) >> 23 & 0xff; 399 if (!e) 400 e++; 401 return e - 0x7f - 23; 402 } 403 static inline int 404 ulpscale_d (double x) 405 { 406 int e = asuint64 (x) >> 52 & 0x7ff; 407 if (!e) 408 e++; 409 return e - 0x3ff - 52; 410 } 411 static inline float 412 call_f1 (const struct fun *f, struct args_f1 a) 413 { 414 return f->fun.f1 (a.x); 415 } 416 static inline float 417 call_f2 (const struct fun *f, struct args_f2 a) 418 { 419 return f->fun.f2 (a.x, a.x2); 420 } 421 422 static inline double 423 call_d1 (const struct fun *f, struct args_d1 a) 424 { 425 return f->fun.d1 (a.x); 426 } 427 static inline double 428 call_d2 (const struct fun *f, struct args_d2 a) 429 { 430 return f->fun.d2 (a.x, a.x2); 431 } 432 static inline double 433 call_long_f1 (const struct fun *f, struct args_f1 a) 434 { 435 return f->fun_long.f1 (a.x); 436 } 437 static inline double 438 call_long_f2 (const struct fun *f, struct args_f2 a) 439 { 440 return f->fun_long.f2 (a.x, a.x2); 441 } 442 static inline long double 443 call_long_d1 (const struct fun *f, struct args_d1 a) 444 { 445 return f->fun_long.d1 (a.x); 446 } 447 static inline long double 448 call_long_d2 (const struct fun *f, struct args_d2 a) 449 { 450 return f->fun_long.d2 (a.x, a.x2); 451 } 452 static inline void 453 printcall_f1 (const struct fun *f, struct args_f1 a) 454 { 455 printf ("%s(%a)", f->name, a.x); 456 } 457 static inline void 458 printcall_f2 (const struct fun *f, struct args_f2 a) 459 { 460 printf ("%s(%a, %a)", f->name, a.x, a.x2); 461 } 462 static inline void 463 printcall_d1 (const struct fun *f, struct args_d1 a) 464 { 465 printf ("%s(%a)", f->name, a.x); 466 } 467 static inline void 468 printcall_d2 (const struct fun *f, struct args_d2 a) 469 { 470 printf ("%s(%a, %a)", f->name, a.x, a.x2); 471 } 472 static inline void 473 printgen_f1 (const struct fun *f, struct gen *gen) 474 { 475 printf ("%s in [%a;%a]", f->name, asfloat (gen->start), 476 asfloat (gen->start + gen->len)); 477 } 478 static inline void 479 printgen_f2 (const struct fun *f, struct gen *gen) 480 { 481 printf ("%s in [%a;%a] x [%a;%a]", f->name, asfloat (gen->start), 482 asfloat (gen->start + gen->len), asfloat (gen->start2), 483 asfloat (gen->start2 + gen->len2)); 484 } 485 static inline void 486 printgen_d1 (const struct fun *f, struct gen *gen) 487 { 488 printf ("%s in [%a;%a]", f->name, asdouble (gen->start), 489 asdouble (gen->start + gen->len)); 490 } 491 static inline void 492 printgen_d2 (const struct fun *f, struct gen *gen) 493 { 494 printf ("%s in [%a;%a] x [%a;%a]", f->name, asdouble (gen->start), 495 asdouble (gen->start + gen->len), asdouble (gen->start2), 496 asdouble (gen->start2 + gen->len2)); 497 } 498 499 #define reduce_f1(a, f, op) (f (a.x)) 500 #define reduce_f2(a, f, op) (f (a.x) op f (a.x2)) 501 #define reduce_d1(a, f, op) (f (a.x)) 502 #define reduce_d2(a, f, op) (f (a.x) op f (a.x2)) 503 504 #ifndef IEEE_754_2008_SNAN 505 # define IEEE_754_2008_SNAN 1 506 #endif 507 static inline int 508 issignaling_f (float x) 509 { 510 uint32_t ix = asuint (x); 511 if (!IEEE_754_2008_SNAN) 512 return (ix & 0x7fc00000) == 0x7fc00000; 513 return 2 * (ix ^ 0x00400000) > 2u * 0x7fc00000; 514 } 515 static inline int 516 issignaling_d (double x) 517 { 518 uint64_t ix = asuint64 (x); 519 if (!IEEE_754_2008_SNAN) 520 return (ix & 0x7ff8000000000000) == 0x7ff8000000000000; 521 return 2 * (ix ^ 0x0008000000000000) > 2 * 0x7ff8000000000000ULL; 522 } 523 524 #if USE_MPFR 525 static mpfr_rnd_t 526 rmap (int r) 527 { 528 switch (r) 529 { 530 case FE_TONEAREST: 531 return MPFR_RNDN; 532 case FE_TOWARDZERO: 533 return MPFR_RNDZ; 534 case FE_UPWARD: 535 return MPFR_RNDU; 536 case FE_DOWNWARD: 537 return MPFR_RNDD; 538 } 539 return -1; 540 } 541 542 #define prec_mpfr_f 50 543 #define prec_mpfr_d 80 544 #define prec_f 24 545 #define prec_d 53 546 #define emin_f -148 547 #define emin_d -1073 548 #define emax_f 128 549 #define emax_d 1024 550 static inline int 551 call_mpfr_f1 (mpfr_t y, const struct fun *f, struct args_f1 a, mpfr_rnd_t r) 552 { 553 MPFR_DECL_INIT (x, prec_f); 554 mpfr_set_flt (x, a.x, MPFR_RNDN); 555 return f->fun_mpfr.f1 (y, x, r); 556 } 557 static inline int 558 call_mpfr_f2 (mpfr_t y, const struct fun *f, struct args_f2 a, mpfr_rnd_t r) 559 { 560 MPFR_DECL_INIT (x, prec_f); 561 MPFR_DECL_INIT (x2, prec_f); 562 mpfr_set_flt (x, a.x, MPFR_RNDN); 563 mpfr_set_flt (x2, a.x2, MPFR_RNDN); 564 return f->fun_mpfr.f2 (y, x, x2, r); 565 } 566 static inline int 567 call_mpfr_d1 (mpfr_t y, const struct fun *f, struct args_d1 a, mpfr_rnd_t r) 568 { 569 MPFR_DECL_INIT (x, prec_d); 570 mpfr_set_d (x, a.x, MPFR_RNDN); 571 return f->fun_mpfr.d1 (y, x, r); 572 } 573 static inline int 574 call_mpfr_d2 (mpfr_t y, const struct fun *f, struct args_d2 a, mpfr_rnd_t r) 575 { 576 MPFR_DECL_INIT (x, prec_d); 577 MPFR_DECL_INIT (x2, prec_d); 578 mpfr_set_d (x, a.x, MPFR_RNDN); 579 mpfr_set_d (x2, a.x2, MPFR_RNDN); 580 return f->fun_mpfr.d2 (y, x, x2, r); 581 } 582 #endif 583 584 #define float_f float 585 #define double_f double 586 #define copysign_f copysignf 587 #define nextafter_f nextafterf 588 #define fabs_f fabsf 589 #define asuint_f asuint 590 #define asfloat_f asfloat 591 #define scalbn_f scalbnf 592 #define lscalbn_f scalbn 593 #define halfinf_f 0x1p127f 594 #define min_normal_f 0x1p-126f 595 596 #define float_d double 597 #define double_d long double 598 #define copysign_d copysign 599 #define nextafter_d nextafter 600 #define fabs_d fabs 601 #define asuint_d asuint64 602 #define asfloat_d asdouble 603 #define scalbn_d scalbn 604 #define lscalbn_d scalbnl 605 #define halfinf_d 0x1p1023 606 #define min_normal_d 0x1p-1022 607 608 #define NEW_RT 609 #define RT(x) x##_f 610 #define T(x) x##_f1 611 #include "ulp.h" 612 #undef T 613 #define T(x) x##_f2 614 #include "ulp.h" 615 #undef T 616 #undef RT 617 618 #define NEW_RT 619 #define RT(x) x##_d 620 #define T(x) x##_d1 621 #include "ulp.h" 622 #undef T 623 #define T(x) x##_d2 624 #include "ulp.h" 625 #undef T 626 #undef RT 627 628 static void 629 usage (void) 630 { 631 puts ("./ulp [-q] [-m] [-f] [-r nudz] [-l soft-ulplimit] [-e ulplimit] func " 632 "lo [hi [x lo2 hi2] [count]]"); 633 puts ("Compares func against a higher precision implementation in [lo; hi]."); 634 puts ("-q: quiet."); 635 puts ("-m: use mpfr even if faster method is available."); 636 puts ("-f: disable fenv testing (rounding modes and exceptions)."); 637 #if __aarch64__ && WANT_VMATH 638 puts ("-c: neutral 'control value' to test behaviour when one lane can affect another. \n" 639 " This should be different from tested input in other lanes, and non-special \n" 640 " (i.e. should not trigger fenv exceptions). Default is 1."); 641 #endif 642 puts ("Supported func:"); 643 for (const struct fun *f = fun; f->name; f++) 644 printf ("\t%s\n", f->name); 645 exit (1); 646 } 647 648 static int 649 cmp (const struct fun *f, struct gen *gen, const struct conf *conf) 650 { 651 int r = 1; 652 if (f->arity == 1 && f->singleprec) 653 r = cmp_f1 (f, gen, conf); 654 else if (f->arity == 2 && f->singleprec) 655 r = cmp_f2 (f, gen, conf); 656 else if (f->arity == 1 && !f->singleprec) 657 r = cmp_d1 (f, gen, conf); 658 else if (f->arity == 2 && !f->singleprec) 659 r = cmp_d2 (f, gen, conf); 660 else 661 usage (); 662 return r; 663 } 664 665 static uint64_t 666 getnum (const char *s, int singleprec) 667 { 668 // int i; 669 uint64_t sign = 0; 670 // char buf[12]; 671 672 if (s[0] == '+') 673 s++; 674 else if (s[0] == '-') 675 { 676 sign = singleprec ? 1ULL << 31 : 1ULL << 63; 677 s++; 678 } 679 /* 0xXXXX is treated as bit representation, '-' flips the sign bit. */ 680 if (s[0] == '0' && tolower (s[1]) == 'x' && strchr (s, 'p') == 0) 681 return sign ^ strtoull (s, 0, 0); 682 // /* SNaN, QNaN, NaN, Inf. */ 683 // for (i=0; s[i] && i < sizeof buf; i++) 684 // buf[i] = tolower(s[i]); 685 // buf[i] = 0; 686 // if (strcmp(buf, "snan") == 0) 687 // return sign | (singleprec ? 0x7fa00000 : 0x7ff4000000000000); 688 // if (strcmp(buf, "qnan") == 0 || strcmp(buf, "nan") == 0) 689 // return sign | (singleprec ? 0x7fc00000 : 0x7ff8000000000000); 690 // if (strcmp(buf, "inf") == 0 || strcmp(buf, "infinity") == 0) 691 // return sign | (singleprec ? 0x7f800000 : 0x7ff0000000000000); 692 /* Otherwise assume it's a floating-point literal. */ 693 return sign 694 | (singleprec ? asuint (strtof (s, 0)) : asuint64 (strtod (s, 0))); 695 } 696 697 static void 698 parsegen (struct gen *g, int argc, char *argv[], const struct fun *f) 699 { 700 int singleprec = f->singleprec; 701 int arity = f->arity; 702 uint64_t a, b, a2, b2, n; 703 if (argc < 1) 704 usage (); 705 b = a = getnum (argv[0], singleprec); 706 n = 0; 707 if (argc > 1 && strcmp (argv[1], "x") == 0) 708 { 709 argc -= 2; 710 argv += 2; 711 } 712 else if (argc > 1) 713 { 714 b = getnum (argv[1], singleprec); 715 if (argc > 2 && strcmp (argv[2], "x") == 0) 716 { 717 argc -= 3; 718 argv += 3; 719 } 720 } 721 b2 = a2 = getnum (argv[0], singleprec); 722 if (argc > 1) 723 b2 = getnum (argv[1], singleprec); 724 if (argc > 2) 725 n = strtoull (argv[2], 0, 0); 726 if (argc > 3) 727 usage (); 728 //printf("ab %lx %lx ab2 %lx %lx n %lu\n", a, b, a2, b2, n); 729 if (arity == 1) 730 { 731 g->start = a; 732 g->len = b - a; 733 if (n - 1 > b - a) 734 n = b - a + 1; 735 g->off = 0; 736 g->step = n ? (g->len + 1) / n : 1; 737 g->start2 = g->len2 = 0; 738 g->cnt = n; 739 } 740 else if (arity == 2) 741 { 742 g->start = a; 743 g->len = b - a; 744 g->off = g->step = 0; 745 g->start2 = a2; 746 g->len2 = b2 - a2; 747 g->cnt = n; 748 } 749 else 750 usage (); 751 } 752 753 int 754 main (int argc, char *argv[]) 755 { 756 const struct fun *f; 757 struct gen gen; 758 struct conf conf; 759 conf.rc = 'n'; 760 conf.quiet = 0; 761 conf.mpfr = 0; 762 conf.fenv = 1; 763 conf.softlim = 0; 764 conf.errlim = INFINITY; 765 for (;;) 766 { 767 argc--; 768 argv++; 769 if (argc < 1) 770 usage (); 771 if (argv[0][0] != '-') 772 break; 773 switch (argv[0][1]) 774 { 775 case 'e': 776 argc--; 777 argv++; 778 if (argc < 1) 779 usage (); 780 conf.errlim = strtod (argv[0], 0); 781 break; 782 case 'f': 783 conf.fenv = 0; 784 break; 785 case 'l': 786 argc--; 787 argv++; 788 if (argc < 1) 789 usage (); 790 conf.softlim = strtod (argv[0], 0); 791 break; 792 case 'm': 793 conf.mpfr = 1; 794 break; 795 case 'q': 796 conf.quiet = 1; 797 break; 798 case 'r': 799 conf.rc = argv[0][2]; 800 if (!conf.rc) 801 { 802 argc--; 803 argv++; 804 if (argc < 1) 805 usage (); 806 conf.rc = argv[0][0]; 807 } 808 break; 809 #if __aarch64__ && WANT_VMATH 810 case 'c': 811 argc--; 812 argv++; 813 fv[0] = strtof(argv[0], 0); 814 dv[0] = strtod(argv[0], 0); 815 break; 816 #endif 817 default: 818 usage (); 819 } 820 } 821 switch (conf.rc) 822 { 823 case 'n': 824 conf.r = FE_TONEAREST; 825 break; 826 case 'u': 827 conf.r = FE_UPWARD; 828 break; 829 case 'd': 830 conf.r = FE_DOWNWARD; 831 break; 832 case 'z': 833 conf.r = FE_TOWARDZERO; 834 break; 835 default: 836 usage (); 837 } 838 for (f = fun; f->name; f++) 839 if (strcmp (argv[0], f->name) == 0) 840 break; 841 if (!f->name) 842 usage (); 843 if (!f->singleprec && LDBL_MANT_DIG == DBL_MANT_DIG) 844 conf.mpfr = 1; /* Use mpfr if long double has no extra precision. */ 845 if (!USE_MPFR && conf.mpfr) 846 { 847 puts ("mpfr is not available."); 848 return 0; 849 } 850 argc--; 851 argv++; 852 parsegen (&gen, argc, argv, f); 853 conf.n = gen.cnt; 854 return cmp (f, &gen, &conf); 855 } 856