libc/AOR_v20.02/math/test/ulp.c

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