1 /*
2 * ULP error checking tool for math functions.
3 *
4 * Copyright (c) 2019-2023, 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 struct conf
202 {
203 int r;
204 int rc;
205 int quiet;
206 int mpfr;
207 int fenv;
208 unsigned long long n;
209 double softlim;
210 double errlim;
211 int ignore_zero_sign;
212 };
213
214 /* A bit of a hack: call vector functions twice with the same
215 input in lane 0 but a different value in other lanes: once
216 with an in-range value and then with a special case value. */
217 static int secondcall;
218
219 /* Wrappers for vector functions. */
220 #ifdef __vpcs
221 typedef __f32x4_t v_float;
222 typedef __f64x2_t v_double;
223 /* First element of fv and dv may be changed by -c argument. */
224 static float fv[2] = {1.0f, -INFINITY};
225 static double dv[2] = {1.0, -INFINITY};
argf(float x)226 static inline v_float argf(float x) { return (v_float){x,x,x,fv[secondcall]}; }
argd(double x)227 static inline v_double argd(double x) { return (v_double){x,dv[secondcall]}; }
228 #if WANT_SVE_MATH
229 #include <arm_sve.h>
230 typedef __SVFloat32_t sv_float;
231 typedef __SVFloat64_t sv_double;
232
svargf(float x)233 static inline sv_float svargf(float x) {
234 int n = svcntw();
235 float base[n];
236 for (int i=0; i<n; i++)
237 base[i] = (float)x;
238 base[n-1] = (float) fv[secondcall];
239 return svld1(svptrue_b32(), base);
240 }
svargd(double x)241 static inline sv_double svargd(double x) {
242 int n = svcntd();
243 double base[n];
244 for (int i=0; i<n; i++)
245 base[i] = x;
246 base[n-1] = dv[secondcall];
247 return svld1(svptrue_b64(), base);
248 }
svretf(sv_float vec)249 static inline float svretf(sv_float vec) {
250 int n = svcntw();
251 float res[n];
252 svst1(svptrue_b32(), res, vec);
253 return res[0];
254 }
svretd(sv_double vec)255 static inline double svretd(sv_double vec) {
256 int n = svcntd();
257 double res[n];
258 svst1(svptrue_b64(), res, vec);
259 return res[0];
260 }
261 #endif
262 #endif
263
264 #include "test/ulp_wrappers.h"
265
266 struct fun
267 {
268 const char *name;
269 int arity;
270 int singleprec;
271 int twice;
272 union
273 {
274 float (*f1) (float);
275 float (*f2) (float, float);
276 double (*d1) (double);
277 double (*d2) (double, double);
278 } fun;
279 union
280 {
281 double (*f1) (double);
282 double (*f2) (double, double);
283 long double (*d1) (long double);
284 long double (*d2) (long double, long double);
285 } fun_long;
286 #if USE_MPFR
287 union
288 {
289 int (*f1) (mpfr_t, const mpfr_t, mpfr_rnd_t);
290 int (*f2) (mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t);
291 int (*d1) (mpfr_t, const mpfr_t, mpfr_rnd_t);
292 int (*d2) (mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t);
293 } fun_mpfr;
294 #endif
295 };
296
297 static const struct fun fun[] = {
298 #if USE_MPFR
299 # define F(x, x_wrap, x_long, x_mpfr, a, s, t, twice) \
300 {#x, a, s, twice, {.t = x_wrap}, {.t = x_long}, {.t = x_mpfr}},
301 #else
302 # define F(x, x_wrap, x_long, x_mpfr, a, s, t, twice) \
303 {#x, a, s, twice, {.t = x_wrap}, {.t = x_long}},
304 #endif
305 #define F1(x) F (x##f, x##f, x, mpfr_##x, 1, 1, f1, 0)
306 #define F2(x) F (x##f, x##f, x, mpfr_##x, 2, 1, f2, 0)
307 #define D1(x) F (x, x, x##l, mpfr_##x, 1, 0, d1, 0)
308 #define D2(x) F (x, x, x##l, mpfr_##x, 2, 0, d2, 0)
309 /* Neon routines. */
310 #define VF1(x) F (__v_##x##f, v_##x##f, x, mpfr_##x, 1, 1, f1, 0)
311 #define VF2(x) F (__v_##x##f, v_##x##f, x, mpfr_##x, 2, 1, f2, 0)
312 #define VD1(x) F (__v_##x, v_##x, x##l, mpfr_##x, 1, 0, d1, 0)
313 #define VD2(x) F (__v_##x, v_##x, x##l, mpfr_##x, 2, 0, d2, 0)
314 #define VNF1(x) F (__vn_##x##f, vn_##x##f, x, mpfr_##x, 1, 1, f1, 0)
315 #define VNF2(x) F (__vn_##x##f, vn_##x##f, x, mpfr_##x, 2, 1, f2, 0)
316 #define VND1(x) F (__vn_##x, vn_##x, x##l, mpfr_##x, 1, 0, d1, 0)
317 #define VND2(x) F (__vn_##x, vn_##x, x##l, mpfr_##x, 2, 0, d2, 0)
318 #define ZVF1(x) F (_ZGVnN4v_##x##f, Z_##x##f, x, mpfr_##x, 1, 1, f1, 0)
319 #define ZVF2(x) F (_ZGVnN4vv_##x##f, Z_##x##f, x, mpfr_##x, 2, 1, f2, 0)
320 #define ZVD1(x) F (_ZGVnN2v_##x, Z_##x, x##l, mpfr_##x, 1, 0, d1, 0)
321 #define ZVD2(x) F (_ZGVnN2vv_##x, Z_##x, x##l, mpfr_##x, 2, 0, d2, 0)
322 #define ZVNF1(x) VNF1 (x) ZVF1 (x)
323 #define ZVNF2(x) VNF2 (x) ZVF2 (x)
324 #define ZVND1(x) VND1 (x) ZVD1 (x)
325 #define ZVND2(x) VND2 (x) ZVD2 (x)
326 /* SVE routines. */
327 #define SVF1(x) F (__sv_##x##f, sv_##x##f, x, mpfr_##x, 1, 1, f1, 0)
328 #define SVF2(x) F (__sv_##x##f, sv_##x##f, x, mpfr_##x, 2, 1, f2, 0)
329 #define SVD1(x) F (__sv_##x, sv_##x, x##l, mpfr_##x, 1, 0, d1, 0)
330 #define SVD2(x) F (__sv_##x, sv_##x, x##l, mpfr_##x, 2, 0, d2, 0)
331 #define ZSVF1(x) F (_ZGVsMxv_##x##f, Z_sv_##x##f, x, mpfr_##x, 1, 1, f1, 0)
332 #define ZSVF2(x) F (_ZGVsMxvv_##x##f, Z_sv_##x##f, x, mpfr_##x, 2, 1, f2, 0)
333 #define ZSVD1(x) F (_ZGVsMxv_##x, Z_sv_##x, x##l, mpfr_##x, 1, 0, d1, 0)
334 #define ZSVD2(x) F (_ZGVsMxvv_##x, Z_sv_##x, x##l, mpfr_##x, 2, 0, d2, 0)
335
336 #include "test/ulp_funcs.h"
337
338 #undef F
339 #undef F1
340 #undef F2
341 #undef D1
342 #undef D2
343 #undef SVF1
344 #undef SVF2
345 #undef SVD1
346 #undef SVD2
347 {0}};
348
349 /* Boilerplate for generic calls. */
350
351 static inline int
ulpscale_f(float x)352 ulpscale_f (float x)
353 {
354 int e = asuint (x) >> 23 & 0xff;
355 if (!e)
356 e++;
357 return e - 0x7f - 23;
358 }
359 static inline int
ulpscale_d(double x)360 ulpscale_d (double x)
361 {
362 int e = asuint64 (x) >> 52 & 0x7ff;
363 if (!e)
364 e++;
365 return e - 0x3ff - 52;
366 }
367 static inline float
call_f1(const struct fun * f,struct args_f1 a)368 call_f1 (const struct fun *f, struct args_f1 a)
369 {
370 return f->fun.f1 (a.x);
371 }
372 static inline float
call_f2(const struct fun * f,struct args_f2 a)373 call_f2 (const struct fun *f, struct args_f2 a)
374 {
375 return f->fun.f2 (a.x, a.x2);
376 }
377
378 static inline double
call_d1(const struct fun * f,struct args_d1 a)379 call_d1 (const struct fun *f, struct args_d1 a)
380 {
381 return f->fun.d1 (a.x);
382 }
383 static inline double
call_d2(const struct fun * f,struct args_d2 a)384 call_d2 (const struct fun *f, struct args_d2 a)
385 {
386 return f->fun.d2 (a.x, a.x2);
387 }
388 static inline double
call_long_f1(const struct fun * f,struct args_f1 a)389 call_long_f1 (const struct fun *f, struct args_f1 a)
390 {
391 return f->fun_long.f1 (a.x);
392 }
393 static inline double
call_long_f2(const struct fun * f,struct args_f2 a)394 call_long_f2 (const struct fun *f, struct args_f2 a)
395 {
396 return f->fun_long.f2 (a.x, a.x2);
397 }
398 static inline long double
call_long_d1(const struct fun * f,struct args_d1 a)399 call_long_d1 (const struct fun *f, struct args_d1 a)
400 {
401 return f->fun_long.d1 (a.x);
402 }
403 static inline long double
call_long_d2(const struct fun * f,struct args_d2 a)404 call_long_d2 (const struct fun *f, struct args_d2 a)
405 {
406 return f->fun_long.d2 (a.x, a.x2);
407 }
408 static inline void
printcall_f1(const struct fun * f,struct args_f1 a)409 printcall_f1 (const struct fun *f, struct args_f1 a)
410 {
411 printf ("%s(%a)", f->name, a.x);
412 }
413 static inline void
printcall_f2(const struct fun * f,struct args_f2 a)414 printcall_f2 (const struct fun *f, struct args_f2 a)
415 {
416 printf ("%s(%a, %a)", f->name, a.x, a.x2);
417 }
418 static inline void
printcall_d1(const struct fun * f,struct args_d1 a)419 printcall_d1 (const struct fun *f, struct args_d1 a)
420 {
421 printf ("%s(%a)", f->name, a.x);
422 }
423 static inline void
printcall_d2(const struct fun * f,struct args_d2 a)424 printcall_d2 (const struct fun *f, struct args_d2 a)
425 {
426 printf ("%s(%a, %a)", f->name, a.x, a.x2);
427 }
428 static inline void
printgen_f1(const struct fun * f,struct gen * gen)429 printgen_f1 (const struct fun *f, struct gen *gen)
430 {
431 printf ("%s in [%a;%a]", f->name, asfloat (gen->start),
432 asfloat (gen->start + gen->len));
433 }
434 static inline void
printgen_f2(const struct fun * f,struct gen * gen)435 printgen_f2 (const struct fun *f, struct gen *gen)
436 {
437 printf ("%s in [%a;%a] x [%a;%a]", f->name, asfloat (gen->start),
438 asfloat (gen->start + gen->len), asfloat (gen->start2),
439 asfloat (gen->start2 + gen->len2));
440 }
441 static inline void
printgen_d1(const struct fun * f,struct gen * gen)442 printgen_d1 (const struct fun *f, struct gen *gen)
443 {
444 printf ("%s in [%a;%a]", f->name, asdouble (gen->start),
445 asdouble (gen->start + gen->len));
446 }
447 static inline void
printgen_d2(const struct fun * f,struct gen * gen)448 printgen_d2 (const struct fun *f, struct gen *gen)
449 {
450 printf ("%s in [%a;%a] x [%a;%a]", f->name, asdouble (gen->start),
451 asdouble (gen->start + gen->len), asdouble (gen->start2),
452 asdouble (gen->start2 + gen->len2));
453 }
454
455 #define reduce_f1(a, f, op) (f (a.x))
456 #define reduce_f2(a, f, op) (f (a.x) op f (a.x2))
457 #define reduce_d1(a, f, op) (f (a.x))
458 #define reduce_d2(a, f, op) (f (a.x) op f (a.x2))
459
460 #ifndef IEEE_754_2008_SNAN
461 # define IEEE_754_2008_SNAN 1
462 #endif
463 static inline int
issignaling_f(float x)464 issignaling_f (float x)
465 {
466 uint32_t ix = asuint (x);
467 if (!IEEE_754_2008_SNAN)
468 return (ix & 0x7fc00000) == 0x7fc00000;
469 return 2 * (ix ^ 0x00400000) > 2u * 0x7fc00000;
470 }
471 static inline int
issignaling_d(double x)472 issignaling_d (double x)
473 {
474 uint64_t ix = asuint64 (x);
475 if (!IEEE_754_2008_SNAN)
476 return (ix & 0x7ff8000000000000) == 0x7ff8000000000000;
477 return 2 * (ix ^ 0x0008000000000000) > 2 * 0x7ff8000000000000ULL;
478 }
479
480 #if USE_MPFR
481 static mpfr_rnd_t
rmap(int r)482 rmap (int r)
483 {
484 switch (r)
485 {
486 case FE_TONEAREST:
487 return MPFR_RNDN;
488 case FE_TOWARDZERO:
489 return MPFR_RNDZ;
490 case FE_UPWARD:
491 return MPFR_RNDU;
492 case FE_DOWNWARD:
493 return MPFR_RNDD;
494 }
495 return -1;
496 }
497
498 #define prec_mpfr_f 50
499 #define prec_mpfr_d 80
500 #define prec_f 24
501 #define prec_d 53
502 #define emin_f -148
503 #define emin_d -1073
504 #define emax_f 128
505 #define emax_d 1024
506 static inline int
call_mpfr_f1(mpfr_t y,const struct fun * f,struct args_f1 a,mpfr_rnd_t r)507 call_mpfr_f1 (mpfr_t y, const struct fun *f, struct args_f1 a, mpfr_rnd_t r)
508 {
509 MPFR_DECL_INIT (x, prec_f);
510 mpfr_set_flt (x, a.x, MPFR_RNDN);
511 return f->fun_mpfr.f1 (y, x, r);
512 }
513 static inline int
call_mpfr_f2(mpfr_t y,const struct fun * f,struct args_f2 a,mpfr_rnd_t r)514 call_mpfr_f2 (mpfr_t y, const struct fun *f, struct args_f2 a, mpfr_rnd_t r)
515 {
516 MPFR_DECL_INIT (x, prec_f);
517 MPFR_DECL_INIT (x2, prec_f);
518 mpfr_set_flt (x, a.x, MPFR_RNDN);
519 mpfr_set_flt (x2, a.x2, MPFR_RNDN);
520 return f->fun_mpfr.f2 (y, x, x2, r);
521 }
522 static inline int
call_mpfr_d1(mpfr_t y,const struct fun * f,struct args_d1 a,mpfr_rnd_t r)523 call_mpfr_d1 (mpfr_t y, const struct fun *f, struct args_d1 a, mpfr_rnd_t r)
524 {
525 MPFR_DECL_INIT (x, prec_d);
526 mpfr_set_d (x, a.x, MPFR_RNDN);
527 return f->fun_mpfr.d1 (y, x, r);
528 }
529 static inline int
call_mpfr_d2(mpfr_t y,const struct fun * f,struct args_d2 a,mpfr_rnd_t r)530 call_mpfr_d2 (mpfr_t y, const struct fun *f, struct args_d2 a, mpfr_rnd_t r)
531 {
532 MPFR_DECL_INIT (x, prec_d);
533 MPFR_DECL_INIT (x2, prec_d);
534 mpfr_set_d (x, a.x, MPFR_RNDN);
535 mpfr_set_d (x2, a.x2, MPFR_RNDN);
536 return f->fun_mpfr.d2 (y, x, x2, r);
537 }
538 #endif
539
540 #define float_f float
541 #define double_f double
542 #define copysign_f copysignf
543 #define nextafter_f nextafterf
544 #define fabs_f fabsf
545 #define asuint_f asuint
546 #define asfloat_f asfloat
547 #define scalbn_f scalbnf
548 #define lscalbn_f scalbn
549 #define halfinf_f 0x1p127f
550 #define min_normal_f 0x1p-126f
551
552 #define float_d double
553 #define double_d long double
554 #define copysign_d copysign
555 #define nextafter_d nextafter
556 #define fabs_d fabs
557 #define asuint_d asuint64
558 #define asfloat_d asdouble
559 #define scalbn_d scalbn
560 #define lscalbn_d scalbnl
561 #define halfinf_d 0x1p1023
562 #define min_normal_d 0x1p-1022
563
564 #define NEW_RT
565 #define RT(x) x##_f
566 #define T(x) x##_f1
567 #include "ulp.h"
568 #undef T
569 #define T(x) x##_f2
570 #include "ulp.h"
571 #undef T
572 #undef RT
573
574 #define NEW_RT
575 #define RT(x) x##_d
576 #define T(x) x##_d1
577 #include "ulp.h"
578 #undef T
579 #define T(x) x##_d2
580 #include "ulp.h"
581 #undef T
582 #undef RT
583
584 static void
usage(void)585 usage (void)
586 {
587 puts ("./ulp [-q] [-m] [-f] [-r {n|u|d|z}] [-l soft-ulplimit] [-e ulplimit] func "
588 "lo [hi [x lo2 hi2] [count]]");
589 puts ("Compares func against a higher precision implementation in [lo; hi].");
590 puts ("-q: quiet.");
591 puts ("-m: use mpfr even if faster method is available.");
592 puts ("-f: disable fenv exceptions testing.");
593 #ifdef ___vpcs
594 puts ("-c: neutral 'control value' to test behaviour when one lane can affect another. \n"
595 " This should be different from tested input in other lanes, and non-special \n"
596 " (i.e. should not trigger fenv exceptions). Default is 1.");
597 #endif
598 puts ("-z: ignore sign of 0.");
599 puts ("Supported func:");
600 for (const struct fun *f = fun; f->name; f++)
601 printf ("\t%s\n", f->name);
602 exit (1);
603 }
604
605 static int
cmp(const struct fun * f,struct gen * gen,const struct conf * conf)606 cmp (const struct fun *f, struct gen *gen, const struct conf *conf)
607 {
608 int r = 1;
609 if (f->arity == 1 && f->singleprec)
610 r = cmp_f1 (f, gen, conf);
611 else if (f->arity == 2 && f->singleprec)
612 r = cmp_f2 (f, gen, conf);
613 else if (f->arity == 1 && !f->singleprec)
614 r = cmp_d1 (f, gen, conf);
615 else if (f->arity == 2 && !f->singleprec)
616 r = cmp_d2 (f, gen, conf);
617 else
618 usage ();
619 return r;
620 }
621
622 static uint64_t
getnum(const char * s,int singleprec)623 getnum (const char *s, int singleprec)
624 {
625 // int i;
626 uint64_t sign = 0;
627 // char buf[12];
628
629 if (s[0] == '+')
630 s++;
631 else if (s[0] == '-')
632 {
633 sign = singleprec ? 1ULL << 31 : 1ULL << 63;
634 s++;
635 }
636 /* 0xXXXX is treated as bit representation, '-' flips the sign bit. */
637 if (s[0] == '0' && tolower (s[1]) == 'x' && strchr (s, 'p') == 0)
638 return sign ^ strtoull (s, 0, 0);
639 // /* SNaN, QNaN, NaN, Inf. */
640 // for (i=0; s[i] && i < sizeof buf; i++)
641 // buf[i] = tolower(s[i]);
642 // buf[i] = 0;
643 // if (strcmp(buf, "snan") == 0)
644 // return sign | (singleprec ? 0x7fa00000 : 0x7ff4000000000000);
645 // if (strcmp(buf, "qnan") == 0 || strcmp(buf, "nan") == 0)
646 // return sign | (singleprec ? 0x7fc00000 : 0x7ff8000000000000);
647 // if (strcmp(buf, "inf") == 0 || strcmp(buf, "infinity") == 0)
648 // return sign | (singleprec ? 0x7f800000 : 0x7ff0000000000000);
649 /* Otherwise assume it's a floating-point literal. */
650 return sign
651 | (singleprec ? asuint (strtof (s, 0)) : asuint64 (strtod (s, 0)));
652 }
653
654 static void
parsegen(struct gen * g,int argc,char * argv[],const struct fun * f)655 parsegen (struct gen *g, int argc, char *argv[], const struct fun *f)
656 {
657 int singleprec = f->singleprec;
658 int arity = f->arity;
659 uint64_t a, b, a2, b2, n;
660 if (argc < 1)
661 usage ();
662 b = a = getnum (argv[0], singleprec);
663 n = 0;
664 if (argc > 1 && strcmp (argv[1], "x") == 0)
665 {
666 argc -= 2;
667 argv += 2;
668 }
669 else if (argc > 1)
670 {
671 b = getnum (argv[1], singleprec);
672 if (argc > 2 && strcmp (argv[2], "x") == 0)
673 {
674 argc -= 3;
675 argv += 3;
676 }
677 }
678 b2 = a2 = getnum (argv[0], singleprec);
679 if (argc > 1)
680 b2 = getnum (argv[1], singleprec);
681 if (argc > 2)
682 n = strtoull (argv[2], 0, 0);
683 if (argc > 3)
684 usage ();
685 //printf("ab %lx %lx ab2 %lx %lx n %lu\n", a, b, a2, b2, n);
686 if (arity == 1)
687 {
688 g->start = a;
689 g->len = b - a;
690 if (n - 1 > b - a)
691 n = b - a + 1;
692 g->off = 0;
693 g->step = n ? (g->len + 1) / n : 1;
694 g->start2 = g->len2 = 0;
695 g->cnt = n;
696 }
697 else if (arity == 2)
698 {
699 g->start = a;
700 g->len = b - a;
701 g->off = g->step = 0;
702 g->start2 = a2;
703 g->len2 = b2 - a2;
704 g->cnt = n;
705 }
706 else
707 usage ();
708 }
709
710 int
main(int argc,char * argv[])711 main (int argc, char *argv[])
712 {
713 const struct fun *f;
714 struct gen gen;
715 struct conf conf;
716 conf.rc = 'n';
717 conf.quiet = 0;
718 conf.mpfr = 0;
719 conf.fenv = 1;
720 conf.softlim = 0;
721 conf.errlim = INFINITY;
722 conf.ignore_zero_sign = 0;
723 for (;;)
724 {
725 argc--;
726 argv++;
727 if (argc < 1)
728 usage ();
729 if (argv[0][0] != '-')
730 break;
731 switch (argv[0][1])
732 {
733 case 'e':
734 argc--;
735 argv++;
736 if (argc < 1)
737 usage ();
738 conf.errlim = strtod (argv[0], 0);
739 break;
740 case 'f':
741 conf.fenv = 0;
742 break;
743 case 'l':
744 argc--;
745 argv++;
746 if (argc < 1)
747 usage ();
748 conf.softlim = strtod (argv[0], 0);
749 break;
750 case 'm':
751 conf.mpfr = 1;
752 break;
753 case 'q':
754 conf.quiet = 1;
755 break;
756 case 'r':
757 conf.rc = argv[0][2];
758 if (!conf.rc)
759 {
760 argc--;
761 argv++;
762 if (argc < 1 || argv[0][1] != '\0')
763 usage ();
764 conf.rc = argv[0][0];
765 }
766 break;
767 case 'z':
768 conf.ignore_zero_sign = 1;
769 break;
770 #ifdef __vpcs
771 case 'c':
772 argc--;
773 argv++;
774 fv[0] = strtof(argv[0], 0);
775 dv[0] = strtod(argv[0], 0);
776 break;
777 #endif
778 default:
779 usage ();
780 }
781 }
782 switch (conf.rc)
783 {
784 case 'n':
785 conf.r = FE_TONEAREST;
786 break;
787 case 'u':
788 conf.r = FE_UPWARD;
789 break;
790 case 'd':
791 conf.r = FE_DOWNWARD;
792 break;
793 case 'z':
794 conf.r = FE_TOWARDZERO;
795 break;
796 default:
797 usage ();
798 }
799 for (f = fun; f->name; f++)
800 if (strcmp (argv[0], f->name) == 0)
801 break;
802 if (!f->name)
803 {
804 #ifndef __vpcs
805 /* Ignore vector math functions if vector math is not supported. */
806 if (strncmp (argv[0], "_ZGVnN", 6) == 0)
807 exit (0);
808 #endif
809 #if !WANT_SVE_MATH
810 if (strncmp (argv[0], "_ZGVsMxv", 8) == 0)
811 exit (0);
812 #endif
813 printf ("math function %s not supported\n", argv[0]);
814 exit (1);
815 }
816 if (!f->singleprec && LDBL_MANT_DIG == DBL_MANT_DIG)
817 conf.mpfr = 1; /* Use mpfr if long double has no extra precision. */
818 if (!USE_MPFR && conf.mpfr)
819 {
820 puts ("mpfr is not available.");
821 return 0;
822 }
823 argc--;
824 argv++;
825 parsegen (&gen, argc, argv, f);
826 conf.n = gen.cnt;
827 return cmp (f, &gen, &conf);
828 }
829