1 /*
2 * Generic functions for ULP error estimation.
3 *
4 * Copyright (c) 2019-2024, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8 /* For each different math function type,
9 T(x) should add a different suffix to x.
10 RT(x) should add a return type specific suffix to x. */
11
12 #ifdef NEW_RT
13 #undef NEW_RT
14
15 # if USE_MPFR
RT(ulpscale_mpfr)16 static int RT(ulpscale_mpfr) (mpfr_t x, int t)
17 {
18 /* TODO: pow of 2 cases. */
19 if (mpfr_regular_p (x))
20 {
21 mpfr_exp_t e = mpfr_get_exp (x) - RT(prec);
22 if (e < RT(emin))
23 e = RT(emin) - 1;
24 if (e > RT(emax) - RT(prec))
25 e = RT(emax) - RT(prec);
26 return e;
27 }
28 if (mpfr_zero_p (x))
29 return RT(emin) - 1;
30 if (mpfr_inf_p (x))
31 return RT(emax) - RT(prec);
32 /* NaN. */
33 return 0;
34 }
35 # endif
36
37 /* Difference between exact result and closest real number that
38 gets rounded to got, i.e. error before rounding, for a correctly
39 rounded result the difference is 0. */
RT(ulperr)40 static double RT (ulperr) (RT (float) got, const struct RT (ret) * p, int r,
41 int ignore_zero_sign)
42 {
43 RT(float) want = p->y;
44 RT(float) d;
45 double e;
46
47 if (RT(asuint) (got) == RT(asuint) (want))
48 return 0.0;
49 if (isnan (got) && isnan (want))
50 /* Ignore sign of NaN, and signalling-ness for MPFR. */
51 # if USE_MPFR
52 return 0;
53 # else
54 return RT (issignaling) (got) == RT (issignaling) (want) ? 0 : INFINITY;
55 # endif
56 if (signbit (got) != signbit (want))
57 {
58 /* Fall through to ULP calculation if ignoring sign of zero and at
59 exactly one of want and got is non-zero. */
60 if (ignore_zero_sign && want == got)
61 return 0.0;
62 if (!ignore_zero_sign || (want != 0 && got != 0))
63 return INFINITY;
64 }
65 if (!isfinite (want) || !isfinite (got))
66 {
67 if (isnan (got) != isnan (want))
68 return INFINITY;
69 if (isnan (want))
70 return 0;
71 if (isinf (got))
72 {
73 got = RT(copysign) (RT(halfinf), got);
74 want *= 0.5f;
75 }
76 if (isinf (want))
77 {
78 want = RT(copysign) (RT(halfinf), want);
79 got *= 0.5f;
80 }
81 }
82 if (r == FE_TONEAREST)
83 {
84 // TODO: incorrect when got vs want cross a powof2 boundary
85 /* error = got > want
86 ? got - want - tail ulp - 0.5 ulp
87 : got - want - tail ulp + 0.5 ulp. */
88 d = got - want;
89 e = d > 0 ? -p->tail - 0.5 : -p->tail + 0.5;
90 }
91 else
92 {
93 if ((r == FE_DOWNWARD && got < want) || (r == FE_UPWARD && got > want)
94 || (r == FE_TOWARDZERO && fabs (got) < fabs (want)))
95 got = RT(nextafter) (got, want);
96 d = got - want;
97 e = -p->tail;
98 }
99 return RT(scalbn) (d, -p->ulpexp) + e;
100 }
101
RT(isok)102 static int RT(isok) (RT(float) ygot, int exgot, RT(float) ywant, int exwant,
103 int exmay)
104 {
105 return RT(asuint) (ygot) == RT(asuint) (ywant)
106 && ((exgot ^ exwant) & ~exmay) == 0;
107 }
108
RT(isok_nofenv)109 static int RT(isok_nofenv) (RT(float) ygot, RT(float) ywant)
110 {
111 return RT(asuint) (ygot) == RT(asuint) (ywant);
112 }
113 #endif
114
T(call_fenv)115 static inline void T (call_fenv) (const struct fun *f, struct T (args) a,
116 int r, RT (float) * y, int *ex,
117 const struct conf *conf)
118 {
119 if (r != FE_TONEAREST)
120 fesetround (r);
121 feclearexcept (FE_ALL_EXCEPT);
122 *y = T (call) (f, a, conf);
123 *ex = fetestexcept (FE_ALL_EXCEPT);
124 if (r != FE_TONEAREST)
125 fesetround (FE_TONEAREST);
126 }
127
T(call_nofenv)128 static inline void T (call_nofenv) (const struct fun *f, struct T (args) a,
129 int r, RT (float) * y, int *ex,
130 const struct conf *conf)
131 {
132 if (r != FE_TONEAREST)
133 fesetround (r);
134 *y = T (call) (f, a, conf);
135 *ex = 0;
136 if (r != FE_TONEAREST)
137 fesetround (FE_TONEAREST);
138 }
139
T(call_long_fenv)140 static inline int T (call_long_fenv) (const struct fun *f, struct T (args) a,
141 int r, struct RT (ret) * p,
142 RT (float) ygot, int exgot)
143 {
144 if (r != FE_TONEAREST)
145 fesetround (r);
146 feclearexcept (FE_ALL_EXCEPT);
147 volatile struct T(args) va = a; // TODO: barrier
148 a = va;
149 RT(double) yl = T(call_long) (f, a);
150 p->y = (RT(float)) yl;
151 volatile RT(float) vy = p->y; // TODO: barrier
152 (void) vy;
153 p->ex = fetestexcept (FE_ALL_EXCEPT);
154 if (r != FE_TONEAREST)
155 fesetround (FE_TONEAREST);
156 p->ex_may = FE_INEXACT;
157 if (RT(isok) (ygot, exgot, p->y, p->ex, p->ex_may))
158 return 1;
159 p->ulpexp = RT(ulpscale) (p->y);
160 if (isinf (p->y))
161 p->tail = RT(lscalbn) (yl - (RT(double)) 2 * RT(halfinf), -p->ulpexp);
162 else
163 p->tail = RT(lscalbn) (yl - p->y, -p->ulpexp);
164 if (RT(fabs) (p->y) < RT(min_normal))
165 {
166 /* TODO: subnormal result is treated as undeflow even if it's
167 exact since call_long may not raise inexact correctly. */
168 if (p->y != 0 || (p->ex & FE_INEXACT))
169 p->ex |= FE_UNDERFLOW | FE_INEXACT;
170 }
171 return 0;
172 }
T(call_long_nofenv)173 static inline int T(call_long_nofenv) (const struct fun *f, struct T(args) a,
174 int r, struct RT(ret) * p,
175 RT(float) ygot, int exgot)
176 {
177 if (r != FE_TONEAREST)
178 fesetround (r);
179 RT(double) yl = T(call_long) (f, a);
180 p->y = (RT(float)) yl;
181 if (r != FE_TONEAREST)
182 fesetround (FE_TONEAREST);
183 if (RT(isok_nofenv) (ygot, p->y))
184 return 1;
185 p->ulpexp = RT(ulpscale) (p->y);
186 if (isinf (p->y))
187 p->tail = RT(lscalbn) (yl - (RT(double)) 2 * RT(halfinf), -p->ulpexp);
188 else
189 p->tail = RT(lscalbn) (yl - p->y, -p->ulpexp);
190 return 0;
191 }
192
193 /* There are nan input args and all quiet. */
T(qnanpropagation)194 static inline int T(qnanpropagation) (struct T(args) a)
195 {
196 return T(reduce) (a, isnan, ||) && !T(reduce) (a, RT(issignaling), ||);
197 }
T(sum)198 static inline RT(float) T(sum) (struct T(args) a)
199 {
200 return T(reduce) (a, , +);
201 }
202
203 /* returns 1 if the got result is ok. */
T(call_mpfr_fix)204 static inline int T(call_mpfr_fix) (const struct fun *f, struct T(args) a,
205 int r_fenv, struct RT(ret) * p,
206 RT(float) ygot, int exgot)
207 {
208 #if USE_MPFR
209 int t, t2;
210 mpfr_rnd_t r = rmap (r_fenv);
211 MPFR_DECL_INIT(my, RT(prec_mpfr));
212 MPFR_DECL_INIT(mr, RT(prec));
213 MPFR_DECL_INIT(me, RT(prec_mpfr));
214 mpfr_clear_flags ();
215 t = T(call_mpfr) (my, f, a, r);
216 /* Double rounding. */
217 t2 = mpfr_set (mr, my, r);
218 if (t2)
219 t = t2;
220 mpfr_set_emin (RT(emin));
221 mpfr_set_emax (RT(emax));
222 t = mpfr_check_range (mr, t, r);
223 t = mpfr_subnormalize (mr, t, r);
224 mpfr_set_emax (MPFR_EMAX_DEFAULT);
225 mpfr_set_emin (MPFR_EMIN_DEFAULT);
226 p->y = mpfr_get_d (mr, r);
227 p->ex = t ? FE_INEXACT : 0;
228 p->ex_may = FE_INEXACT;
229 if (mpfr_underflow_p () && (p->ex & FE_INEXACT))
230 /* TODO: handle before and after rounding uflow cases. */
231 p->ex |= FE_UNDERFLOW;
232 if (mpfr_overflow_p ())
233 p->ex |= FE_OVERFLOW | FE_INEXACT;
234 if (mpfr_divby0_p ())
235 p->ex |= FE_DIVBYZERO;
236 //if (mpfr_erangeflag_p ())
237 // p->ex |= FE_INVALID;
238 if (!mpfr_nanflag_p () && RT(isok) (ygot, exgot, p->y, p->ex, p->ex_may))
239 return 1;
240 if (mpfr_nanflag_p () && !T(qnanpropagation) (a))
241 p->ex |= FE_INVALID;
242 p->ulpexp = RT(ulpscale_mpfr) (my, t);
243 if (!isfinite (p->y))
244 {
245 p->tail = 0;
246 if (isnan (p->y))
247 {
248 /* If an input was nan keep its sign. */
249 p->y = T(sum) (a);
250 if (!isnan (p->y))
251 p->y = (p->y - p->y) / (p->y - p->y);
252 return RT(isok) (ygot, exgot, p->y, p->ex, p->ex_may);
253 }
254 mpfr_set_si_2exp (mr, signbit (p->y) ? -1 : 1, 1024, MPFR_RNDN);
255 if (mpfr_cmpabs (my, mr) >= 0)
256 return RT(isok) (ygot, exgot, p->y, p->ex, p->ex_may);
257 }
258 mpfr_sub (me, my, mr, MPFR_RNDN);
259 mpfr_mul_2si (me, me, -p->ulpexp, MPFR_RNDN);
260 p->tail = mpfr_get_d (me, MPFR_RNDN);
261 return 0;
262 #else
263 abort ();
264 #endif
265 }
266
T(cmp)267 static int T(cmp) (const struct fun *f, struct gen *gen,
268 const struct conf *conf)
269 {
270 double maxerr = 0;
271 uint64_t cnt = 0;
272 uint64_t cnt1 = 0;
273 uint64_t cnt2 = 0;
274 uint64_t cntfail = 0;
275 int r = conf->r;
276 int use_mpfr = conf->mpfr;
277 int fenv = conf->fenv;
278
279 for (;;)
280 {
281 struct RT(ret) want;
282 struct T(args) a = T(next) (gen);
283 int exgot;
284 int exgot2;
285 RT(float) ygot;
286 RT(float) ygot2;
287 int fail = 0;
288 if (fenv)
289 T (call_fenv) (f, a, r, &ygot, &exgot, conf);
290 else
291 T (call_nofenv) (f, a, r, &ygot, &exgot, conf);
292 if (f->twice) {
293 secondcall = 1;
294 if (fenv)
295 T (call_fenv) (f, a, r, &ygot2, &exgot2, conf);
296 else
297 T (call_nofenv) (f, a, r, &ygot2, &exgot2, conf);
298 secondcall = 0;
299 if (RT(asuint) (ygot) != RT(asuint) (ygot2))
300 {
301 fail = 1;
302 cntfail++;
303 T(printcall) (f, a);
304 printf (" got %a then %a for same input\n", ygot, ygot2);
305 }
306 }
307 cnt++;
308 int ok = use_mpfr
309 ? T(call_mpfr_fix) (f, a, r, &want, ygot, exgot)
310 : (fenv ? T(call_long_fenv) (f, a, r, &want, ygot, exgot)
311 : T(call_long_nofenv) (f, a, r, &want, ygot, exgot));
312 if (!ok)
313 {
314 int print = 0;
315 double err = RT (ulperr) (ygot, &want, r, conf->ignore_zero_sign);
316 double abserr = fabs (err);
317 // TODO: count errors below accuracy limit.
318 if (abserr > 0)
319 cnt1++;
320 if (abserr > 1)
321 cnt2++;
322 if (abserr > conf->errlim)
323 {
324 print = 1;
325 if (!fail)
326 {
327 fail = 1;
328 cntfail++;
329 }
330 }
331 if (abserr > maxerr)
332 {
333 maxerr = abserr;
334 if (!conf->quiet && abserr > conf->softlim)
335 print = 1;
336 }
337 if (print)
338 {
339 T(printcall) (f, a);
340 // TODO: inf ulp handling
341 printf (" got %a want %a %+g ulp err %g\n", ygot, want.y,
342 want.tail, err);
343 }
344 int diff = fenv ? exgot ^ want.ex : 0;
345 if (fenv && (diff & ~want.ex_may))
346 {
347 if (!fail)
348 {
349 fail = 1;
350 cntfail++;
351 }
352 T(printcall) (f, a);
353 printf (" is %a %+g ulp, got except 0x%0x", want.y, want.tail,
354 exgot);
355 if (diff & exgot)
356 printf (" wrongly set: 0x%x", diff & exgot);
357 if (diff & ~exgot)
358 printf (" wrongly clear: 0x%x", diff & ~exgot);
359 putchar ('\n');
360 }
361 }
362 if (cnt >= conf->n)
363 break;
364 if (!conf->quiet && cnt % 0x100000 == 0)
365 printf ("progress: %6.3f%% cnt %llu cnt1 %llu cnt2 %llu cntfail %llu "
366 "maxerr %g\n",
367 100.0 * cnt / conf->n, (unsigned long long) cnt,
368 (unsigned long long) cnt1, (unsigned long long) cnt2,
369 (unsigned long long) cntfail, maxerr);
370 }
371 double cc = cnt;
372 if (cntfail)
373 printf ("FAIL ");
374 else
375 printf ("PASS ");
376 T(printgen) (f, gen);
377 printf (" round %c errlim %g maxerr %g %s cnt %llu cnt1 %llu %g%% cnt2 %llu "
378 "%g%% cntfail %llu %g%%\n",
379 conf->rc, conf->errlim,
380 maxerr, conf->r == FE_TONEAREST ? "+0.5" : "+1.0",
381 (unsigned long long) cnt,
382 (unsigned long long) cnt1, 100.0 * cnt1 / cc,
383 (unsigned long long) cnt2, 100.0 * cnt2 / cc,
384 (unsigned long long) cntfail, 100.0 * cntfail / cc);
385 return !!cntfail;
386 }
387