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
32static inline uint64_t
33asuint64 (double f)
34{
35 union
36 {
37 double f;
38 uint64_t i;
39 } u = {f};
40 return u.i;
41}
42
43static inline double
44asdouble (uint64_t i)
45{
46 union
47 {
48 uint64_t i;
49 double f;
50 } u = {i};
51 return u.f;
52}
53
54static inline uint32_t
55asuint (float f)
56{
57 union
58 {
59 float f;
60 uint32_t i;
61 } u = {f};
62 return u.i;
63}
64
65static inline float
66asfloat (uint32_t i)
67{
68 union
69 {
70 uint32_t i;
71 float f;
72 } u = {i};
73 return u.f;
74}
75
76static uint64_t seed = 0x0123456789abcdef;
77static uint64_t
78rand64 (void)
79{
80 seed = 6364136223846793005ull * seed + 1;
81 return seed ^ (seed >> 32);
82}
83
84/* Uniform random in [0,n]. */
85static uint64_t
86randn (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 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
115struct args_f1
116{
117 float x;
118};
119
120struct args_f2
121{
122 float x;
123 float x2;
124};
125
126struct args_d1
127{
128 double x;
129};
130
131struct args_d2
132{
133 double x;
134 double x2;
135};
136
137/* result = y + tail*2^ulpexp. */
138struct ret_f
139{
140 float y;
141 double tail;
142 int ulpexp;
143 int ex;
144 int ex_may;
145};
146
147struct ret_d
148{
149 double y;
150 double tail;
151 int ulpexp;
152 int ex;
153 int ex_may;
154};
155
156static inline uint64_t
157next1 (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 (n: g->step / 2);
164 if (g->off > g->len)
165 g->off -= g->len; /* hack. */
166 return r;
167}
168
169static inline uint64_t
170next2 (uint64_t *x2, struct gen *g)
171{
172 /* For two arguments use uniform random sampling. */
173 uint64_t r = g->start + randn (n: g->len);
174 *x2 = g->start2 + randn (n: g->len2);
175 return r;
176}
177
178static struct args_f1
179next_f1 (void *g)
180{
181 return (struct args_f1){asfloat (i: next1 (g))};
182}
183
184static struct args_f2
185next_f2 (void *g)
186{
187 uint64_t x2;
188 uint64_t x = next2 (x2: &x2, g);
189 return (struct args_f2){asfloat (i: x), asfloat (i: x2)};
190}
191
192static struct args_d1
193next_d1 (void *g)
194{
195 return (struct args_d1){asdouble (i: next1 (g))};
196}
197
198static struct args_d2
199next_d2 (void *g)
200{
201 uint64_t x2;
202 uint64_t x = next2 (x2: &x2, g);
203 return (struct args_d2){asdouble (i: x), asdouble (i: x2)};
204}
205
206struct 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. */
219static float sincosf_sinf(float x) {(void)cosf(x: x); return sinf(x: x);}
220static float sincosf_cosf(float x) {(void)sinf(x: x); return cosf(x: x);}
221static double sincos_sin(double x) {(void)cos(x: x); return sin(x: x);}
222static double sincos_cos(double x) {(void)sin(x: x); return cos(x: x);}
223#if USE_MPFR
224static 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#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. */
231static int secondcall;
232
233/* Wrappers for vector functions. */
234#if __aarch64__ && WANT_VMATH
235typedef __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 __vpcs
256static 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#endif
281#endif
282
283struct 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
314static 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
405static inline int
406ulpscale_f (float x)
407{
408 int e = asuint (f: x) >> 23 & 0xff;
409 if (!e)
410 e++;
411 return e - 0x7f - 23;
412}
413static inline int
414ulpscale_d (double x)
415{
416 int e = asuint64 (f: x) >> 52 & 0x7ff;
417 if (!e)
418 e++;
419 return e - 0x3ff - 52;
420}
421static inline float
422call_f1 (const struct fun *f, struct args_f1 a)
423{
424 return f->fun.f1 (a.x);
425}
426static inline float
427call_f2 (const struct fun *f, struct args_f2 a)
428{
429 return f->fun.f2 (a.x, a.x2);
430}
431
432static inline double
433call_d1 (const struct fun *f, struct args_d1 a)
434{
435 return f->fun.d1 (a.x);
436}
437static inline double
438call_d2 (const struct fun *f, struct args_d2 a)
439{
440 return f->fun.d2 (a.x, a.x2);
441}
442static inline double
443call_long_f1 (const struct fun *f, struct args_f1 a)
444{
445 return f->fun_long.f1 (a.x);
446}
447static inline double
448call_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 double
453call_long_d1 (const struct fun *f, struct args_d1 a)
454{
455 return f->fun_long.d1 (a.x);
456}
457static inline long double
458call_long_d2 (const struct fun *f, struct args_d2 a)
459{
460 return f->fun_long.d2 (a.x, a.x2);
461}
462static inline void
463printcall_f1 (const struct fun *f, struct args_f1 a)
464{
465 printf (format: "%s(%a)", f->name, a.x);
466}
467static inline void
468printcall_f2 (const struct fun *f, struct args_f2 a)
469{
470 printf (format: "%s(%a, %a)", f->name, a.x, a.x2);
471}
472static inline void
473printcall_d1 (const struct fun *f, struct args_d1 a)
474{
475 printf (format: "%s(%a)", f->name, a.x);
476}
477static inline void
478printcall_d2 (const struct fun *f, struct args_d2 a)
479{
480 printf (format: "%s(%a, %a)", f->name, a.x, a.x2);
481}
482static inline void
483printgen_f1 (const struct fun *f, struct gen *gen)
484{
485 printf (format: "%s in [%a;%a]", f->name, asfloat (i: gen->start),
486 asfloat (i: gen->start + gen->len));
487}
488static inline void
489printgen_f2 (const struct fun *f, struct gen *gen)
490{
491 printf (format: "%s in [%a;%a] x [%a;%a]", f->name, asfloat (i: gen->start),
492 asfloat (i: gen->start + gen->len), asfloat (i: gen->start2),
493 asfloat (i: gen->start2 + gen->len2));
494}
495static inline void
496printgen_d1 (const struct fun *f, struct gen *gen)
497{
498 printf (format: "%s in [%a;%a]", f->name, asdouble (i: gen->start),
499 asdouble (i: gen->start + gen->len));
500}
501static inline void
502printgen_d2 (const struct fun *f, struct gen *gen)
503{
504 printf (format: "%s in [%a;%a] x [%a;%a]", f->name, asdouble (i: gen->start),
505 asdouble (i: gen->start + gen->len), asdouble (i: gen->start2),
506 asdouble (i: 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
517static inline int
518issignaling_f (float x)
519{
520 uint32_t ix = asuint (f: x);
521 if (!IEEE_754_2008_SNAN)
522 return (ix & 0x7fc00000) == 0x7fc00000;
523 return 2 * (ix ^ 0x00400000) > 2u * 0x7fc00000;
524}
525static inline int
526issignaling_d (double x)
527{
528 uint64_t ix = asuint64 (f: 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
535static mpfr_rnd_t
536rmap (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
560static inline int
561call_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 int
568call_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 int
577call_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 int
584call_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
638static void
639usage (void)
640{
641 puts (s: "./ulp [-q] [-m] [-f] [-r nudz] [-l soft-ulplimit] [-e ulplimit] func "
642 "lo [hi [x lo2 hi2] [count]]");
643 puts (s: "Compares func against a higher precision implementation in [lo; hi].");
644 puts (s: "-q: quiet.");
645 puts (s: "-m: use mpfr even if faster method is available.");
646 puts (s: "-f: disable fenv testing (rounding modes and exceptions).");
647 puts (s: "Supported func:");
648 for (const struct fun *f = fun; f->name; f++)
649 printf (format: "\t%s\n", f->name);
650 exit (status: 1);
651}
652
653static int
654cmp (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
670static uint64_t
671getnum (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 (c: s[1]) == 'x' && strchr (s: s, c: 'p') == 0)
686 return sign ^ strtoull (nptr: s, endptr: 0, base: 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 (f: strtof (nptr: s, endptr: 0)) : asuint64 (f: strtod (nptr: s, endptr: 0)));
700}
701
702static void
703parsegen (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 (s: argv[0], singleprec);
711 n = 0;
712 if (argc > 1 && strcmp (s1: argv[1], s2: "x") == 0)
713 {
714 argc -= 2;
715 argv += 2;
716 }
717 else if (argc > 1)
718 {
719 b = getnum (s: argv[1], singleprec);
720 if (argc > 2 && strcmp (s1: argv[2], s2: "x") == 0)
721 {
722 argc -= 3;
723 argv += 3;
724 }
725 }
726 b2 = a2 = getnum (s: argv[0], singleprec);
727 if (argc > 1)
728 b2 = getnum (s: argv[1], singleprec);
729 if (argc > 2)
730 n = strtoull (nptr: argv[2], endptr: 0, base: 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
758int
759main (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 (nptr: argv[0], endptr: 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 (nptr: argv[0], endptr: 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 (s1: argv[0], s2: 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 (s: "mpfr is not available.");
845 return 0;
846 }
847 argc--;
848 argv++;
849 parsegen (g: &gen, argc, argv, f);
850 conf.n = gen.cnt;
851 return cmp (f, gen: &gen, conf: &conf);
852}
853

source code of libc/AOR_v20.02/math/test/ulp.c