https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819
--- Comment #11 from Steve Kargl <kargl at gcc dot gnu.org> ---
We're fighting a losing battle with the evaluation of jn(30,x)
near one of its zeros. One can either adjust the expected
precision as is done in Christopher's patch or change the
value tested from 475.78 to 475.0, which has a ulp of 2.66.
% cc -o z -O3 j0.c -lm && ./z
475.000000 2.601273e-02
475.000000 2.601272e-02
ulp = 2.665015
475.790833 3.706664e-07
475.790833 3.726422e-07
ulp = 69519.074463
% cat j0.c
#include <float.h>
#include <math.h>
#include <stdio.h>
double
ulp(float app, double acc)
{
int n;
double f;
f = frexp(acc, &n);
f = fabs(acc - app);
f = ldexp(f, FLT_MANT_DIG - n);
return (f);
}
int
main(void)
{
float f, x;
double d, u, um, xm;
x = 475.;
f = jnf(30, x);
d = jn(30,(double)x);
u = ulp(f, d);
printf("% f %e\n", x, jnf(30, x));
printf("% f %le\n", x, jn(30, (double)x));
printf("ulp = %f\n", u);
um = -1.;
do {
f = jnf(30, x);
d = jn(30,(double)x);
u = ulp(f, d);
if (u > um) {
um = u;
xm = x;
}
x = nextafterf(x,1.e30f);
} while (x < 476);
printf("% f %e\n", xm, jnf(30, xm));
printf("% f %le\n", xm, jn(30, (double)xm));
printf("ulp = %f\n", um);
return 0;
}
