https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819
--- Comment #5 from Steve Kargl <kargl at gcc dot gnu.org> --- > > The scalar side lowers to repeated jnf calls, while the transformational > REAL(4) side goes through libgfortran/generated/bessel_r4.c:bessel_jn_r4 and > uses a REAL(4) backward recurrence seeded from jnf(n2, x) and jnf(n2-1, x). > That recurrence seems to accumulate enough rounding error on this system to > miss the test tolerance. > This is not an accumulation of rounding errors issue. Downward recursion is stable for bessel functions. Bessel function implementations near the zeros of the function are notoriously inaccurate. jnf(30,477.77999) = 3.96262854e-04, and a zero occurs near x = 475.791. This is simply a cancellation issue where there are not enough bits. fma might improve the accuracy a bit, but won't eliminate the cancellation. Using a wider precision type as an intermediate will work for REAL(4) and REAL(8), but becomes cumbersome for REAL(10), and REAL(16).
