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