https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819
--- Comment #17 from Steve Kargl <sgk at troutmask dot apl.washington.edu> --- On Tue, Apr 14, 2026 at 09:14:31PM +0000, jvdelisle at gcc dot gnu.org wrote: > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 > > --- Comment #16 from Jerry DeLisle <jvdelisle at gcc dot gnu.org> --- > (In reply to Christopher Albert from comment #15) > > Created attachment 64204 [details] > > alternative patch: shift bessel_6 evaluation point away from the zero > > > > Alternative patch: instead of increasing the tolerance, this shifts the > > bessel_6 evaluation point away from the near-zero cancellation region while > > keeping the existing tolerance. > > I was experimenting with this myself. The value 475.0 is too far and fails. > Your 475.5 must be a sweet spot. It works here. > After a quick peek at the testcase, I'm surprised there aren't other issues. It appears that the testcase sets x and chooses some maximum order mymax and then does a comparison of results for that x for all orders 0 <= n <= mymax. The zeros of Bessel functions can be quite fascinating (https://dlmf.nist.gov/10.21). The transformational BESSEL_JN(0, mymax, X) is mapped to gfortran's _gfortran_bessel_jn_r4() and BESSEL_JN(i, X) is mapped to __builtin_jnf(). The builtin will ultimately lead to jnf() in libm, and likely has a more robust algorithm in choice of starting values in the downward recursion and the normalization of value. I am aware of a paper that tries to deal with the zeros of Bessels (https://www.cl.cam.ac.uk/~jrh13/papers/bessel.pdf) I've never quite worked out how to implement his ideas for the other floating-point precisions. PS: IMO, J3 should never have added the Bessel functions to the Fortran standard.
