https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63488
kargl at gcc dot gnu.org changed:
What |Removed |Added
----------------------------------------------------------------------------
CC| |kargl at gcc dot gnu.org
--- Comment #5 from kargl at gcc dot gnu.org ---
(In reply to Paul Zimmermann from comment #0)
> For the input aa below (with 36 digits, we can recover the exact 113-bit
> binary value by rounding to nearest) mpfr_y0 gets the result cc (more
> precisely, the corresponding 113-bit binary value) which should be correctly
> rounded, and y0q gets the result dd which differs from 180688 ulps.
>
> aa=8.935761195587725798762818805462843676e-01
> cc=-7.446238819993339201075302662649733975e-07 [MPFR]
> dd=-7.446238819993339201075302662815669696e-07
> ulp difference = 180688.000000
>
> For other functions and one million random inputs per function, I got much
> smaller differences (at most 42 ulps).
>
I suppose I don't understand the intent of this PR. The
approximation of a special function with either a polynomial
or a rational approximation is know to have issues near the
zeros of the special function. Given that the zeros of
y0(x) lies close to x=(n+1/4)*pi with n=0,1,2..., I'm
curious to see what you would find if you test in small
intervals about x.
Testing 10000 values in a small interval about the lowest
10 zeros, the double precision y0() on FreeBSD (which comes
from fdlibm) has
ULP: 5918, x = 8.9358398197930655e-01
ULP: 15390, x = 3.9576788857941221e+00
ULP: 17833, x = 7.0860502172517021e+00
ULP: 2510, x = 1.0222342340788490e+01
ULP: 4396, x = 1.3361094710349882e+01
ULP: 1203, x = 1.6500927187922073e+01
ULP: 14157, x = 1.9641309720499763e+01
ULP: 2141, x = 2.2782032287080856e+01
ULP: 1048, x = 2.5922964874664050e+01
ULP: 2219, x = 2.9064027475248540e+01
Note, I did not try to find the worse case.
