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
asuint64(double f)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
asdouble(uint64_t i)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
asuint(float f)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
asfloat(uint32_t i)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
rand64(void)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
randn(uint64_t n)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
next1(struct gen * g)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
next2(uint64_t * x2,struct gen * g)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
next_f1(void * g)179 next_f1 (void *g)
180 {
181 return (struct args_f1){asfloat (next1 (g))};
182 }
183
184 static struct args_f2
next_f2(void * g)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
next_d1(void * g)193 next_d1 (void *g)
194 {
195 return (struct args_d1){asdouble (next1 (g))};
196 }
197
198 static struct args_d2
next_d2(void * g)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. */
sincosf_sinf(float x)219 static float sincosf_sinf(float x) {(void)cosf(x); return sinf(x);}
sincosf_cosf(float x)220 static float sincosf_cosf(float x) {(void)sinf(x); return cosf(x);}
sincos_sin(double x)221 static double sincos_sin(double x) {(void)cos(x); return sin(x);}
sincos_cos(double x)222 static double sincos_cos(double x) {(void)sin(x); return cos(x);}
223 #if USE_MPFR
sincos_mpfr_sin(mpfr_t y,const mpfr_t x,mpfr_rnd_t r)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); }
sincos_mpfr_cos(mpfr_t y,const mpfr_t x,mpfr_rnd_t 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};
argf(float x)239 static inline v_float argf(float x) { return (v_float){x,x,x,fv[secondcall]}; }
argd(double x)240 static inline v_double argd(double x) { return (v_double){x,dv[secondcall]}; }
241
v_sinf(float x)242 static float v_sinf(float x) { return __v_sinf(argf(x))[0]; }
v_cosf(float x)243 static float v_cosf(float x) { return __v_cosf(argf(x))[0]; }
v_expf_1u(float x)244 static float v_expf_1u(float x) { return __v_expf_1u(argf(x))[0]; }
v_expf(float x)245 static float v_expf(float x) { return __v_expf(argf(x))[0]; }
v_exp2f_1u(float x)246 static float v_exp2f_1u(float x) { return __v_exp2f_1u(argf(x))[0]; }
v_exp2f(float x)247 static float v_exp2f(float x) { return __v_exp2f(argf(x))[0]; }
v_logf(float x)248 static float v_logf(float x) { return __v_logf(argf(x))[0]; }
v_powf(float x,float y)249 static float v_powf(float x, float y) { return __v_powf(argf(x),argf(y))[0]; }
v_sin(double x)250 static double v_sin(double x) { return __v_sin(argd(x))[0]; }
v_cos(double x)251 static double v_cos(double x) { return __v_cos(argd(x))[0]; }
v_exp(double x)252 static double v_exp(double x) { return __v_exp(argd(x))[0]; }
v_log(double x)253 static double v_log(double x) { return __v_log(argd(x))[0]; }
v_pow(double x,double y)254 static double v_pow(double x, double y) { return __v_pow(argd(x),argd(y))[0]; }
255 #ifdef __vpcs
vn_sinf(float x)256 static float vn_sinf(float x) { return __vn_sinf(argf(x))[0]; }
vn_cosf(float x)257 static float vn_cosf(float x) { return __vn_cosf(argf(x))[0]; }
vn_expf_1u(float x)258 static float vn_expf_1u(float x) { return __vn_expf_1u(argf(x))[0]; }
vn_expf(float x)259 static float vn_expf(float x) { return __vn_expf(argf(x))[0]; }
vn_exp2f_1u(float x)260 static float vn_exp2f_1u(float x) { return __vn_exp2f_1u(argf(x))[0]; }
vn_exp2f(float x)261 static float vn_exp2f(float x) { return __vn_exp2f(argf(x))[0]; }
vn_logf(float x)262 static float vn_logf(float x) { return __vn_logf(argf(x))[0]; }
vn_powf(float x,float y)263 static float vn_powf(float x, float y) { return __vn_powf(argf(x),argf(y))[0]; }
vn_sin(double x)264 static double vn_sin(double x) { return __vn_sin(argd(x))[0]; }
vn_cos(double x)265 static double vn_cos(double x) { return __vn_cos(argd(x))[0]; }
vn_exp(double x)266 static double vn_exp(double x) { return __vn_exp(argd(x))[0]; }
vn_log(double x)267 static double vn_log(double x) { return __vn_log(argd(x))[0]; }
vn_pow(double x,double y)268 static double vn_pow(double x, double y) { return __vn_pow(argd(x),argd(y))[0]; }
Z_sinf(float x)269 static float Z_sinf(float x) { return _ZGVnN4v_sinf(argf(x))[0]; }
Z_cosf(float x)270 static float Z_cosf(float x) { return _ZGVnN4v_cosf(argf(x))[0]; }
Z_expf(float x)271 static float Z_expf(float x) { return _ZGVnN4v_expf(argf(x))[0]; }
Z_exp2f(float x)272 static float Z_exp2f(float x) { return _ZGVnN4v_exp2f(argf(x))[0]; }
Z_logf(float x)273 static float Z_logf(float x) { return _ZGVnN4v_logf(argf(x))[0]; }
Z_powf(float x,float y)274 static float Z_powf(float x, float y) { return _ZGVnN4vv_powf(argf(x),argf(y))[0]; }
Z_sin(double x)275 static double Z_sin(double x) { return _ZGVnN2v_sin(argd(x))[0]; }
Z_cos(double x)276 static double Z_cos(double x) { return _ZGVnN2v_cos(argd(x))[0]; }
Z_exp(double x)277 static double Z_exp(double x) { return _ZGVnN2v_exp(argd(x))[0]; }
Z_log(double x)278 static double Z_log(double x) { return _ZGVnN2v_log(argd(x))[0]; }
Z_pow(double x,double y)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
ulpscale_f(float x)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
ulpscale_d(double x)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
call_f1(const struct fun * f,struct args_f1 a)422 call_f1 (const struct fun *f, struct args_f1 a)
423 {
424 return f->fun.f1 (a.x);
425 }
426 static inline float
call_f2(const struct fun * f,struct args_f2 a)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
call_d1(const struct fun * f,struct args_d1 a)433 call_d1 (const struct fun *f, struct args_d1 a)
434 {
435 return f->fun.d1 (a.x);
436 }
437 static inline double
call_d2(const struct fun * f,struct args_d2 a)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
call_long_f1(const struct fun * f,struct args_f1 a)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
call_long_f2(const struct fun * f,struct args_f2 a)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
call_long_d1(const struct fun * f,struct args_d1 a)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
call_long_d2(const struct fun * f,struct args_d2 a)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
printcall_f1(const struct fun * f,struct args_f1 a)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
printcall_f2(const struct fun * f,struct args_f2 a)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
printcall_d1(const struct fun * f,struct args_d1 a)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
printcall_d2(const struct fun * f,struct args_d2 a)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
printgen_f1(const struct fun * f,struct gen * gen)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
printgen_f2(const struct fun * f,struct gen * gen)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
printgen_d1(const struct fun * f,struct gen * gen)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
printgen_d2(const struct fun * f,struct gen * gen)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
issignaling_f(float x)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
issignaling_d(double x)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
rmap(int r)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
call_mpfr_f1(mpfr_t y,const struct fun * f,struct args_f1 a,mpfr_rnd_t r)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
call_mpfr_f2(mpfr_t y,const struct fun * f,struct args_f2 a,mpfr_rnd_t r)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
call_mpfr_d1(mpfr_t y,const struct fun * f,struct args_d1 a,mpfr_rnd_t r)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
call_mpfr_d2(mpfr_t y,const struct fun * f,struct args_d2 a,mpfr_rnd_t r)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
usage(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
cmp(const struct fun * f,struct gen * gen,const struct conf * conf)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
getnum(const char * s,int singleprec)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
parsegen(struct gen * g,int argc,char * argv[],const struct fun * f)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
main(int argc,char * argv[])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 }
853