libc/AOR_v20.02/math/test/ulp.h

   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
  17 static 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.  */
  41 static 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 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 (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
  90 static 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
  97 static int RT(isok_nofenv) (RT(float) ygot, RT(float) ywant)
  98 {
  99   return RT(asuint) (ygot) == RT(asuint) (ywant);
 100 }
 101 #endif
 102
 103 static 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
 115 static 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
 122 static 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: 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, 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   else
 145     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'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 }
 155 static 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   else
 167     p->tail = RT(lscalbn) (yl - p->y, -p->ulpexp);
 168   return 0;
 169 }
 170
 171 /* There are nan input args and all quiet.  */
 172 static inline int T(qnanpropagation) (struct T(args) a)
 173 {
 174   return T(reduce) (a, isnan, ||) && !T(reduce) (a, RT(issignaling), ||);
 175 }
 176 static 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.  */
 182 static 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
 245 static 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       else
 268         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         else
 274           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_mpfr
 286                  ? 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 handling
 318               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   else
 352     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 }