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 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. */ 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 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 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 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 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 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 } 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. */ 194 static inline int T(qnanpropagation) (struct T(args) a) 195 { 196 return T(reduce) (a, isnan, ||) && !T(reduce) (a, RT(issignaling), ||); 197 } 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. */ 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 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