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