2  * Generic functions for ULP error estimation.
4  * Copyright (c) 2019-2024, Arm Limited.
5  * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
21       mpfr_exp_t e = mpfr_get_exp (x) - RT(prec);  in RT()
23 	e = RT(emin) - 1;  in RT()
24       if (e > RT(emax) - RT(prec))  in RT()
25 	e = RT(emax) - RT(prec);  in RT()
29     return RT(emin) - 1;  in RT()
31     return RT(emax) - RT(prec);  in RT()
33   return 0;  in RT()
39    rounded result the difference is 0.  */
43   RT(float) want = p->y;  in RT()
50   /* Ignore sign of NaN, and signalling-ness for MPFR.  */  in RT()
52     return 0;  in RT()
54     return RT (issignaling) (got) == RT (issignaling) (want) ? 0 : INFINITY;  in RT()
58       /* Fall through to ULP calculation if ignoring sign of zero and at  in RT()
59 	 exactly one of want and got is non-zero.  */  in RT()
62       if (!ignore_zero_sign || (want != 0 && got != 0))  in RT()
70 	return 0;  in RT()
86 	      ? got - want - tail ulp - 0.5 ulp  in RT()
87 	      : got - want - tail ulp + 0.5 ulp.  */  in RT()
88       d = got - want;  in RT()
89       e = d > 0 ? -p->tail - 0.5 : -p->tail + 0.5;  in RT()
96       d = got - want;  in RT()
97       e = -p->tail;  in RT()
99   return RT(scalbn) (d, -p->ulpexp) + e;  in RT()
106 	 && ((exgot ^ exwant) & ~exmay) == 0;  in RT()
135   *ex = 0;  in T()
150   p->y = (RT(float)) yl;  in T()
151   volatile RT(float) vy = p->y; // TODO: barrier  in T()
153   p->ex = fetestexcept (FE_ALL_EXCEPT);  in T()
156   p->ex_may = FE_INEXACT;  in T()
157   if (RT(isok) (ygot, exgot, p->y, p->ex, p->ex_may))  in T()
159   p->ulpexp = RT(ulpscale) (p->y);  in T()
160   if (isinf (p->y))  in T()
161     p->tail = RT(lscalbn) (yl - (RT(double)) 2 * RT(halfinf), -p->ulpexp);  in T()
163     p->tail = RT(lscalbn) (yl - p->y, -p->ulpexp);  in T()
164   if (RT(fabs) (p->y) < RT(min_normal))  in T()
168       if (p->y != 0 || (p->ex & FE_INEXACT))  in T()
169 	p->ex |= FE_UNDERFLOW | FE_INEXACT;  in T()
171   return 0;  in T()
180   p->y = (RT(float)) yl;  in T()
183   if (RT(isok_nofenv) (ygot, p->y))  in T()
185   p->ulpexp = RT(ulpscale) (p->y);  in T()
186   if (isinf (p->y))  in T()
187     p->tail = RT(lscalbn) (yl - (RT(double)) 2 * RT(halfinf), -p->ulpexp);  in T()
189     p->tail = RT(lscalbn) (yl - p->y, -p->ulpexp);  in T()
190   return 0;  in T()
226   p->y = mpfr_get_d (mr, r);  in T()
227   p->ex = t ? FE_INEXACT : 0;  in T()
228   p->ex_may = FE_INEXACT;  in T()
229   if (mpfr_underflow_p () && (p->ex & FE_INEXACT))  in T()
231     p->ex |= FE_UNDERFLOW;  in T()
233     p->ex |= FE_OVERFLOW | FE_INEXACT;  in T()
235     p->ex |= FE_DIVBYZERO;  in T()
237   //  p->ex |= FE_INVALID;  in T()
238   if (!mpfr_nanflag_p () && RT(isok) (ygot, exgot, p->y, p->ex, p->ex_may))  in T()
241     p->ex |= FE_INVALID;  in T()
242   p->ulpexp = RT(ulpscale_mpfr) (my, t);  in T()
243   if (!isfinite (p->y))  in T()
245       p->tail = 0;  in T()
246       if (isnan (p->y))  in T()
249 	  p->y = T(sum) (a);  in T()
250 	  if (!isnan (p->y))  in T()
251 	    p->y = (p->y - p->y) / (p->y - p->y);  in T()
252 	  return RT(isok) (ygot, exgot, p->y, p->ex, p->ex_may);  in T()
254       mpfr_set_si_2exp (mr, signbit (p->y) ? -1 : 1, 1024, MPFR_RNDN);  in T()
255       if (mpfr_cmpabs (my, mr) >= 0)  in T()
256 	return RT(isok) (ygot, exgot, p->y, p->ex, p->ex_may);  in T()
259   mpfr_mul_2si (me, me, -p->ulpexp, MPFR_RNDN);  in T()
260   p->tail = mpfr_get_d (me, MPFR_RNDN);  in T()
261   return 0;  in T()
270   double maxerr = 0;  in T()
271   uint64_t cnt = 0;  in T()
272   uint64_t cnt1 = 0;  in T()
273   uint64_t cnt2 = 0;  in T()
274   uint64_t cntfail = 0;  in T()
275   int r = conf->r;  in T()
276   int use_mpfr = conf->mpfr;  in T()
277   int fenv = conf->fenv;  in T()
287       int fail = 0;  in T()
292       if (f->twice) {  in T()
298 	secondcall = 0;  in T()
314 	  int print = 0;  in T()
315 	  double err = RT (ulperr) (ygot, &want, r, conf->ignore_zero_sign);  in T()
318 	  if (abserr > 0)  in T()
322 	  if (abserr > conf->errlim)  in T()
334 	      if (!conf->quiet && abserr > conf->softlim)  in T()
340 	      // TODO: inf ulp handling  in T()
341 	      printf (" got %a want %a %+g ulp err %g\n", ygot, want.y,  in T()
344 	  int diff = fenv ? exgot ^ want.ex : 0;  in T()
353 	      printf (" is %a %+g ulp, got except 0x%0x", want.y, want.tail,  in T()
356 		printf (" wrongly set: 0x%x", diff & exgot);  in T()
358 		printf (" wrongly clear: 0x%x", diff & ~exgot);  in T()
362       if (cnt >= conf->n)  in T()
364       if (!conf->quiet && cnt % 0x100000 == 0)  in T()
367 		100.0 * cnt / conf->n, (unsigned long long) cnt,  in T()
379 	  conf->rc, conf->errlim,  in T()
380 	  maxerr, conf->r == FE_TONEAREST ? "+0.5" : "+1.0",  in T()