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

source code of libc/AOR_v20.02/math/test/ulp.h