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