364 lines · c
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-exception7 */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_RT14#undef NEW_RT15 16# if USE_MPFR17static 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# endif37 38/* Difference between exact result and closest real number that39 gets rounded to got, i.e. error before rounding, for a correctly40 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) (got) == RT(asuint) (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), got);62 want *= 0.5f;63 }64 if (isinf (want))65 {66 want = RT(copysign) (RT(halfinf), want);67 got *= 0.5f;68 }69 }70 if (r == FE_TONEAREST)71 {72 // TODO: incorrect when got vs want cross a powof2 boundary73 /* error = got > want74 ? got - want - tail ulp - 0.5 ulp75 : 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 else80 {81 if ((r == FE_DOWNWARD && got < want) || (r == FE_UPWARD && got > want)82 || (r == FE_TOWARDZERO && fabs (got) < fabs (want)))83 got = RT(nextafter) (got, want);84 d = got - want;85 e = -p->tail;86 }87 return RT(scalbn) (d, -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) (ygot) == RT(asuint) (ywant)94 && ((exgot ^ exwant) & ~exmay) == 0;95}96 97static int RT(isok_nofenv) (RT(float) ygot, RT(float) ywant)98{99 return RT(asuint) (ygot) == RT(asuint) (ywant);100}101#endif102 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 (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 (r);128 feclearexcept (FE_ALL_EXCEPT);129 volatile struct T(args) va = a; // TODO: barrier130 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: barrier134 (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, p->y, p->ex, p->ex_may))140 return 1;141 p->ulpexp = RT(ulpscale) (p->y);142 if (isinf (p->y))143 p->tail = RT(lscalbn) (yl - (RT(double)) 2 * RT(halfinf), -p->ulpexp);144 else145 p->tail = RT(lscalbn) (yl - p->y, -p->ulpexp);146 if (RT(fabs) (p->y) < RT(min_normal))147 {148 /* TODO: subnormal result is treated as undeflow even if it's149 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, p->y))162 return 1;163 p->ulpexp = RT(ulpscale) (p->y);164 if (isinf (p->y))165 p->tail = RT(lscalbn) (yl - (RT(double)) 2 * RT(halfinf), -p->ulpexp);166 else167 p->tail = RT(lscalbn) (yl - p->y, -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_MPFR187 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#else241 abort ();242#endif243}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) (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, &ygot, &exgot);267 else268 T(call_nofenv) (f, a, r, &ygot, &exgot);269 if (f->twice) {270 secondcall = 1;271 if (fenv)272 T(call_fenv) (f, a, r, &ygot2, &exgot2);273 else274 T(call_nofenv) (f, a, r, &ygot2, &exgot2);275 secondcall = 0;276 if (RT(asuint) (ygot) != RT(asuint) (ygot2))277 {278 fail = 1;279 cntfail++;280 T(printcall) (f, a);281 printf (" got %a then %a for same input\n", ygot, ygot2);282 }283 }284 cnt++;285 int ok = use_mpfr286 ? T(call_mpfr_fix) (f, a, r, &want, ygot, exgot)287 : (fenv ? T(call_long_fenv) (f, a, r, &want, ygot, exgot)288 : T(call_long_nofenv) (f, a, r, &want, ygot, exgot));289 if (!ok)290 {291 int print = 0;292 double err = RT(ulperr) (ygot, &want, r);293 double abserr = fabs (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 handling318 printf (" 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 (" is %a %+g ulp, got except 0x%0x", want.y, want.tail,331 exgot);332 if (diff & exgot)333 printf (" wrongly set: 0x%x", diff & exgot);334 if (diff & ~exgot)335 printf (" wrongly clear: 0x%x", diff & ~exgot);336 putchar ('\n');337 }338 }339 if (cnt >= conf->n)340 break;341 if (!conf->quiet && cnt % 0x100000 == 0)342 printf ("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 ("FAIL ");351 else352 printf ("PASS ");353 T(printgen) (f, gen);354 printf (" 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