853 lines · 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-exception7 */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_MPFR21# define USE_MPFR 022#endif23#if USE_MPFR24# include <mpfr.h>25#endif26 27#ifndef WANT_VMATH28/* Enable the build of vector math code. */29# define WANT_VMATH 130#endif31 32static inline uint64_t33asuint64 (double f)34{35 union36 {37 double f;38 uint64_t i;39 } u = {f};40 return u.i;41}42 43static inline double44asdouble (uint64_t i)45{46 union47 {48 uint64_t i;49 double f;50 } u = {i};51 return u.f;52}53 54static inline uint32_t55asuint (float f)56{57 union58 {59 float f;60 uint32_t i;61 } u = {f};62 return u.i;63}64 65static inline float66asfloat (uint32_t i)67{68 union69 {70 uint32_t i;71 float f;72 } u = {i};73 return u.f;74}75 76static uint64_t seed = 0x0123456789abcdef;77static uint64_t78rand64 (void)79{80 seed = 6364136223846793005ull * seed + 1;81 return seed ^ (seed >> 32);82}83 84/* Uniform random in [0,n]. */85static uint64_t86randn (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 104struct gen105{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 115struct args_f1116{117 float x;118};119 120struct args_f2121{122 float x;123 float x2;124};125 126struct args_d1127{128 double x;129};130 131struct args_d2132{133 double x;134 double x2;135};136 137/* result = y + tail*2^ulpexp. */138struct ret_f139{140 float y;141 double tail;142 int ulpexp;143 int ex;144 int ex_may;145};146 147struct ret_d148{149 double y;150 double tail;151 int ulpexp;152 int ex;153 int ex_may;154};155 156static inline uint64_t157next1 (struct gen *g)158{159 /* For single argument use randomized incremental steps,160 that produce dense sampling without collisions and allow161 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 169static inline uint64_t170next2 (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 178static struct args_f1179next_f1 (void *g)180{181 return (struct args_f1){asfloat (next1 (g))};182}183 184static struct args_f2185next_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 192static struct args_d1193next_d1 (void *g)194{195 return (struct args_d1){asdouble (next1 (g))};196}197 198static struct args_d2199next_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 206struct conf207{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. */219static float sincosf_sinf(float x) {(void)cosf(x); return sinf(x);}220static float sincosf_cosf(float x) {(void)sinf(x); return cosf(x);}221static double sincos_sin(double x) {(void)cos(x); return sin(x);}222static double sincos_cos(double x) {(void)sin(x); return cos(x);}223#if USE_MPFR224static 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); }225static 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#endif227 228/* A bit of a hack: call vector functions twice with the same229 input in lane 0 but a different value in other lanes: once230 with an in-range value and then with a special case value. */231static int secondcall;232 233/* Wrappers for vector functions. */234#if __aarch64__ && WANT_VMATH235typedef __f32x4_t v_float;236typedef __f64x2_t v_double;237static const float fv[2] = {1.0f, -INFINITY};238static const double dv[2] = {1.0, -INFINITY};239static inline v_float argf(float x) { return (v_float){x,x,x,fv[secondcall]}; }240static inline v_double argd(double x) { return (v_double){x,dv[secondcall]}; }241 242static float v_sinf(float x) { return __v_sinf(argf(x))[0]; }243static float v_cosf(float x) { return __v_cosf(argf(x))[0]; }244static float v_expf_1u(float x) { return __v_expf_1u(argf(x))[0]; }245static float v_expf(float x) { return __v_expf(argf(x))[0]; }246static float v_exp2f_1u(float x) { return __v_exp2f_1u(argf(x))[0]; }247static float v_exp2f(float x) { return __v_exp2f(argf(x))[0]; }248static float v_logf(float x) { return __v_logf(argf(x))[0]; }249static float v_powf(float x, float y) { return __v_powf(argf(x),argf(y))[0]; }250static double v_sin(double x) { return __v_sin(argd(x))[0]; }251static double v_cos(double x) { return __v_cos(argd(x))[0]; }252static double v_exp(double x) { return __v_exp(argd(x))[0]; }253static double v_log(double x) { return __v_log(argd(x))[0]; }254static double v_pow(double x, double y) { return __v_pow(argd(x),argd(y))[0]; }255#ifdef __vpcs256static float vn_sinf(float x) { return __vn_sinf(argf(x))[0]; }257static float vn_cosf(float x) { return __vn_cosf(argf(x))[0]; }258static float vn_expf_1u(float x) { return __vn_expf_1u(argf(x))[0]; }259static float vn_expf(float x) { return __vn_expf(argf(x))[0]; }260static float vn_exp2f_1u(float x) { return __vn_exp2f_1u(argf(x))[0]; }261static float vn_exp2f(float x) { return __vn_exp2f(argf(x))[0]; }262static float vn_logf(float x) { return __vn_logf(argf(x))[0]; }263static float vn_powf(float x, float y) { return __vn_powf(argf(x),argf(y))[0]; }264static double vn_sin(double x) { return __vn_sin(argd(x))[0]; }265static double vn_cos(double x) { return __vn_cos(argd(x))[0]; }266static double vn_exp(double x) { return __vn_exp(argd(x))[0]; }267static double vn_log(double x) { return __vn_log(argd(x))[0]; }268static double vn_pow(double x, double y) { return __vn_pow(argd(x),argd(y))[0]; }269static float Z_sinf(float x) { return _ZGVnN4v_sinf(argf(x))[0]; }270static float Z_cosf(float x) { return _ZGVnN4v_cosf(argf(x))[0]; }271static float Z_expf(float x) { return _ZGVnN4v_expf(argf(x))[0]; }272static float Z_exp2f(float x) { return _ZGVnN4v_exp2f(argf(x))[0]; }273static float Z_logf(float x) { return _ZGVnN4v_logf(argf(x))[0]; }274static float Z_powf(float x, float y) { return _ZGVnN4vv_powf(argf(x),argf(y))[0]; }275static double Z_sin(double x) { return _ZGVnN2v_sin(argd(x))[0]; }276static double Z_cos(double x) { return _ZGVnN2v_cos(argd(x))[0]; }277static double Z_exp(double x) { return _ZGVnN2v_exp(argd(x))[0]; }278static double Z_log(double x) { return _ZGVnN2v_log(argd(x))[0]; }279static double Z_pow(double x, double y) { return _ZGVnN2vv_pow(argd(x),argd(y))[0]; }280#endif281#endif282 283struct fun284{285 const char *name;286 int arity;287 int singleprec;288 int twice;289 union290 {291 float (*f1) (float);292 float (*f2) (float, float);293 double (*d1) (double);294 double (*d2) (double, double);295 } fun;296 union297 {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_MPFR304 union305 {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#endif312};313 314static const struct fun fun[] = {315#if USE_MPFR316# 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#else319# 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#endif322#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_VMATH341 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 __vpcs369 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#endif394#endif395#endif396#undef F397#undef F1398#undef F2399#undef D1400#undef D2401 {0}};402 403/* Boilerplate for generic calls. */404 405static inline int406ulpscale_f (float x)407{408 int e = asuint (x) >> 23 & 0xff;409 if (!e)410 e++;411 return e - 0x7f - 23;412}413static inline int414ulpscale_d (double x)415{416 int e = asuint64 (x) >> 52 & 0x7ff;417 if (!e)418 e++;419 return e - 0x3ff - 52;420}421static inline float422call_f1 (const struct fun *f, struct args_f1 a)423{424 return f->fun.f1 (a.x);425}426static inline float427call_f2 (const struct fun *f, struct args_f2 a)428{429 return f->fun.f2 (a.x, a.x2);430}431 432static inline double433call_d1 (const struct fun *f, struct args_d1 a)434{435 return f->fun.d1 (a.x);436}437static inline double438call_d2 (const struct fun *f, struct args_d2 a)439{440 return f->fun.d2 (a.x, a.x2);441}442static inline double443call_long_f1 (const struct fun *f, struct args_f1 a)444{445 return f->fun_long.f1 (a.x);446}447static inline double448call_long_f2 (const struct fun *f, struct args_f2 a)449{450 return f->fun_long.f2 (a.x, a.x2);451}452static inline long double453call_long_d1 (const struct fun *f, struct args_d1 a)454{455 return f->fun_long.d1 (a.x);456}457static inline long double458call_long_d2 (const struct fun *f, struct args_d2 a)459{460 return f->fun_long.d2 (a.x, a.x2);461}462static inline void463printcall_f1 (const struct fun *f, struct args_f1 a)464{465 printf ("%s(%a)", f->name, a.x);466}467static inline void468printcall_f2 (const struct fun *f, struct args_f2 a)469{470 printf ("%s(%a, %a)", f->name, a.x, a.x2);471}472static inline void473printcall_d1 (const struct fun *f, struct args_d1 a)474{475 printf ("%s(%a)", f->name, a.x);476}477static inline void478printcall_d2 (const struct fun *f, struct args_d2 a)479{480 printf ("%s(%a, %a)", f->name, a.x, a.x2);481}482static inline void483printgen_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}488static inline void489printgen_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}495static inline void496printgen_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}501static inline void502printgen_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_SNAN515# define IEEE_754_2008_SNAN 1516#endif517static inline int518issignaling_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}525static inline int526issignaling_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_MPFR535static mpfr_rnd_t536rmap (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 50553#define prec_mpfr_d 80554#define prec_f 24555#define prec_d 53556#define emin_f -148557#define emin_d -1073558#define emax_f 128559#define emax_d 1024560static inline int561call_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}567static inline int568call_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}576static inline int577call_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}583static inline int584call_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#endif593 594#define float_f float595#define double_f double596#define copysign_f copysignf597#define nextafter_f nextafterf598#define fabs_f fabsf599#define asuint_f asuint600#define asfloat_f asfloat601#define scalbn_f scalbnf602#define lscalbn_f scalbn603#define halfinf_f 0x1p127f604#define min_normal_f 0x1p-126f605 606#define float_d double607#define double_d long double608#define copysign_d copysign609#define nextafter_d nextafter610#define fabs_d fabs611#define asuint_d asuint64612#define asfloat_d asdouble613#define scalbn_d scalbn614#define lscalbn_d scalbnl615#define halfinf_d 0x1p1023616#define min_normal_d 0x1p-1022617 618#define NEW_RT619#define RT(x) x##_f620#define T(x) x##_f1621#include "ulp.h"622#undef T623#define T(x) x##_f2624#include "ulp.h"625#undef T626#undef RT627 628#define NEW_RT629#define RT(x) x##_d630#define T(x) x##_d1631#include "ulp.h"632#undef T633#define T(x) x##_d2634#include "ulp.h"635#undef T636#undef RT637 638static void639usage (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 653static int654cmp (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 else666 usage ();667 return r;668}669 670static uint64_t671getnum (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 sign699 | (singleprec ? asuint (strtof (s, 0)) : asuint64 (strtod (s, 0)));700}701 702static void703parsegen (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 else755 usage ();756}757 758int759main (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}853