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