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