[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #18 from Jerry DeLisle --- Created attachment 64206 --> https://gcc.gnu.org/bugzilla/attachment.cgi?id=64206&action=edit A modified test case that is less sensitive to variations around 475.78 I played around a bit to see if the test could be de-senitized. This tests fine here with quite a bit of variation in epsilon and the value around 475. Seems to work OK here, but one can experiment with it.
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #17 from Steve Kargl --- On Tue, Apr 14, 2026 at 09:14:31PM +, jvdelisle at gcc dot gnu.org wrote: > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 > > --- Comment #16 from Jerry DeLisle --- > (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.
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #16 from Jerry DeLisle --- (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.
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #15 from Christopher Albert --- Created attachment 64204 --> https://gcc.gnu.org/bugzilla/attachment.cgi?id=64204&action=edit 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.
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 Jerry DeLisle changed: What|Removed |Added Last reconfirmed||2026-04-14 Ever confirmed|0 |1 Status|UNCONFIRMED |NEW --- Comment #14 from Jerry DeLisle --- (In reply to Christopher Albert from comment #13) --- snip --- > > So I tend to choose a different value to test after giving it some thought. Steve has the most experiance with these numerics so I concur. Lets change the test value to 475. I can do it if all agree.
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #13 from Christopher Albert --- (In reply to anlauf from comment #12) > bessel_6.f90 has already received multiple updates. > > Let's make a choice: increase the tolerance as done previously (multiple > times), > or choose a different value to test (-> comment#11). a) What speaks for increasing tolerance and keeping edge case: Subtle behaviour changes from gcc code will be picked up with high sensitivity. b) What speaks for choosing a different value to test: Subtle behavior changes from hardware floating point arithmetics will not break tests anymore. So what's realistically more probable? From the history of this problem I tend to believe that when we have a test that is affected by FMA vs non-FMA, subtle hardware behavior is more risky than breaking exactly this one test-case but nothing else by code changes in gcc. So I tend to choose a different value to test after giving it some thought.
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #12 from anlauf at gcc dot gnu.org --- bessel_6.f90 has already received multiple updates. Let's make a choice: increase the tolerance as done previously (multiple times), or choose a different value to test (-> comment#11).
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819
--- Comment #11 from Steve Kargl ---
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.00 2.601273e-02
475.00 2.601272e-02
ulp = 2.665015
475.790833 3.706664e-07
475.790833 3.726422e-07
ulp = 69519.074463
% cat j0.c
#include
#include
#include
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;
}
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #10 from anlauf at gcc dot gnu.org --- (In reply to Jerry DeLisle from comment #8) > (In reply to Christopher Albert from comment #7) > > Created attachment 64177 [details] > > testsuite-only patch: relax bessel_6 tolerance near x=475.78 > > > > The testcase point at x = 475.78 is worth keeping, but the final tolerance > > entry is too tight for the generic non-FMA path. > > > > This testsuite-only patch changes the last factor from 15 to 32 and keeps > > the same evaluation point. > > > > Tested with the existing build via direct runtest: > > - `cd /tmp/gcc-bessel-bisect-build/gcc/testsuite/gfortran && runtest --tool > > gfortran dg.exp=bessel_6.f90` -> 12 expected passes > > - `cd /tmp/gcc-bessel-bisect-build/gcc/testsuite/gfortran && runtest --tool > > gfortran dg.exp=bessel_7.f90` -> 12 expected passes > > If there is no failure in the testsuite now, there is no point to changing > the testcase. Platforms w/o FMA? https://gcc.gnu.org/pipermail/gcc-testresults/2026-April/873251.html
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #9 from Christopher Albert --- (In reply to Jerry DeLisle from comment #8) > (In reply to Christopher Albert from comment #7) > > Created attachment 64177 [details] > > testsuite-only patch: relax bessel_6 tolerance near x=475.78 > > > > The testcase point at x = 475.78 is worth keeping, but the final tolerance > > entry is too tight for the generic non-FMA path. > > > > This testsuite-only patch changes the last factor from 15 to 32 and keeps > > the same evaluation point. > > > > Tested with the existing build via direct runtest: > > - `cd /tmp/gcc-bessel-bisect-build/gcc/testsuite/gfortran && runtest --tool > > gfortran dg.exp=bessel_6.f90` -> 12 expected passes > > - `cd /tmp/gcc-bessel-bisect-build/gcc/testsuite/gfortran && runtest --tool > > gfortran dg.exp=bessel_7.f90` -> 12 expected passes > > If there is no failure in the testsuite now, there is no point to changing > the testcase. The problem is, I still frequently hit this problem when compiling and testing locally on my machine - no idea why this configuration is so sensitive to it. Should I just apply FMA flags then per default?
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #8 from Jerry DeLisle --- (In reply to Christopher Albert from comment #7) > Created attachment 64177 [details] > testsuite-only patch: relax bessel_6 tolerance near x=475.78 > > The testcase point at x = 475.78 is worth keeping, but the final tolerance > entry is too tight for the generic non-FMA path. > > This testsuite-only patch changes the last factor from 15 to 32 and keeps > the same evaluation point. > > Tested with the existing build via direct runtest: > - `cd /tmp/gcc-bessel-bisect-build/gcc/testsuite/gfortran && runtest --tool > gfortran dg.exp=bessel_6.f90` -> 12 expected passes > - `cd /tmp/gcc-bessel-bisect-build/gcc/testsuite/gfortran && runtest --tool > gfortran dg.exp=bessel_7.f90` -> 12 expected passes If there is no failure in the testsuite now, there is no point to changing the testcase.
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #1 from anlauf at gcc dot gnu.org --- I do not see the issue on Suse Leap 16.0 with % /lib64/libc.so.6 | head -7 GNU C Library (GNU libc) stable release version 2.40 (git ef321e23c2). Copyright (C) 2024 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. Configured for x86_64-suse-linux. Compiled by GNU CC version 13.4.0. Do you have another compiler on your system that uses the same libc?
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #6 from Christopher Albert --- I checked this further and do not see a GCC source-regression commit to bisect. Between releases/gcc-15.2.0 and current trunk, the only relevant change is 77a8ed6de20b adding `dg-do run` to bessel_6.f90. `libgfortran/generated/bessel_r4.c` has no functional change in that range. A local build from releases/gcc-15.2.0 already fails with the same standalone reproducer when linked against its locally built libgfortran. I also isolated this to code generation of bessel_r4.c on my system: - generic build: FAIL - `-march=x86-64-v3 -mno-fma`: FAIL - `-march=x86-64-v3 -ffp-contract=off`: FAIL - `-mavx2 -mfma`: PASS - `-march=x86-64-v3 -mfma`: PASS The passing variants use fused multiply-subtract in `_gfortran_bessel_jn_r4`, while the failing ones use separate multiply/subtract operations. So this looks like a numerical-sensitivity issue in the REAL(4) recurrence, not a trunk source regression. The testsuite change exposed it, but did not introduce it.
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #3 from anlauf at gcc dot gnu.org --- (In reply to Christopher Albert from comment #2) > I tested other compilers on the same system and libc. > > Results: > - system gfortran 15.2.1: passes > - nvfortran 26.3: passes > - local GCC 16 object linked against system libgfortran: passes > - system gfortran object linked against local GCC 16 libgfortran: fails > > So the failure follows the local libgfortran, not the frontend or glibc > alone. That points directly at the GCC 16 libgfortran transformational > REAL(4) BESSEL_JN runtime path. The only difference I see then is the CFLAGS. (I use "-g -march=native -O2" on an Intel i5-8250U and also compile with --disable-bootstrap.)
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #7 from Christopher Albert --- Created attachment 64177 --> https://gcc.gnu.org/bugzilla/attachment.cgi?id=64177&action=edit testsuite-only patch: relax bessel_6 tolerance near x=475.78 The testcase point at x = 475.78 is worth keeping, but the final tolerance entry is too tight for the generic non-FMA path. This testsuite-only patch changes the last factor from 15 to 32 and keeps the same evaluation point. Tested with the existing build via direct runtest: - `cd /tmp/gcc-bessel-bisect-build/gcc/testsuite/gfortran && runtest --tool gfortran dg.exp=bessel_6.f90` -> 12 expected passes - `cd /tmp/gcc-bessel-bisect-build/gcc/testsuite/gfortran && runtest --tool gfortran dg.exp=bessel_7.f90` -> 12 expected passes
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #4 from Drea Pinski --- (In reply to anlauf from comment #3) > (In reply to Christopher Albert from comment #2) > > I tested other compilers on the same system and libc. > > > > Results: > > - system gfortran 15.2.1: passes > > - nvfortran 26.3: passes > > - local GCC 16 object linked against system libgfortran: passes > > - system gfortran object linked against local GCC 16 libgfortran: fails > > > > So the failure follows the local libgfortran, not the frontend or glibc > > alone. That points directly at the GCC 16 libgfortran transformational > > REAL(4) BESSEL_JN runtime path. > > The only difference I see then is the CFLAGS. > > (I use "-g -march=native -O2" on an Intel i5-8250U and also compile with > --disable-bootstrap.) Hmm, that might mean using fma more here: ret->base_addr[i*stride] = x2rev * (i+1+n1) * last2 - last1; THough maybe it is not more accuracy but rather different rounding. Note I have not looked into the source enough to see the difference.
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 --- Comment #5 from Steve Kargl --- > > 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).
[Bug libfortran/124819] libgfortran: REAL(4) transformational BESSEL_JN loses accuracy against scalar path on x86_64-linux-gnu
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819 Christopher Albert changed: What|Removed |Added CC||albert at tugraz dot at, ||jvdelisle at gcc dot gnu.org --- Comment #2 from Christopher Albert --- I tested other compilers on the same system and libc. Results: - system gfortran 15.2.1: passes - nvfortran 26.3: passes - local GCC 16 object linked against system libgfortran: passes - system gfortran object linked against local GCC 16 libgfortran: fails So the failure follows the local libgfortran, not the frontend or glibc alone. That points directly at the GCC 16 libgfortran transformational REAL(4) BESSEL_JN runtime path.
