https://gcc.gnu.org/bugzilla/show_bug.cgi?id=124819
Bug ID: 124819
Summary: libgfortran: REAL(4) transformational BESSEL_JN loses
accuracy against scalar path on x86_64-linux-gnu
Product: gcc
Version: 16.0
Status: UNCONFIRMED
Severity: normal
Priority: P3
Component: libfortran
Assignee: unassigned at gcc dot gnu.org
Reporter: albert at tugraz dot at
Target Milestone: ---
The runtime transformational REAL(4) BESSEL_JN implementation appears to lose
too much accuracy relative to the scalar BESSEL_JN path on my x86_64-linux-gnu
system.
I can reproduce a failure of gcc/testsuite/gfortran.dg/bessel_6.f90 with
current trunk from 2026-04-08. The failing comparison is between:
- transformational path: BESSEL_JN(0, mymax, x)
- scalar path: [(BESSEL_JN(i, x), i = 0, mymax)]
The first mismatch on my system is:
mymax = 37
x = 475.779999
i = 30
rec(i) = 3.96264251E-04
lib(i) = 3.96265590E-04
relative error = 3.37849360E-06
tolerance = 1.78813934E-06
A small instrumented version of bessel_6.f90 that prints the first failing
point is:
implicit none
real,parameter :: values(*) = [0.0, 0.5, 1.0, 0.9,
1.8,2.0,3.0,4.0,4.25,8.0,34.53, 475.78]
real,parameter :: myeps(size(values)) = epsilon(0.0) * [2, 7, 5, 6, 9, 12, 12,
7, 7, 8, 98, 15]
integer,parameter :: mymax(size(values)) = [100, 17, 23, 21, 27, 28, 32, 35,
31, 41, 47, 37]
integer, parameter :: Nmax = 100
real :: rec(0:Nmax), lib(0:Nmax)
integer :: i
do i = 1, ubound(values,dim=1)
call compare(mymax(i), values(i), myeps(i))
end do
contains
subroutine compare(mymax, x, myeps)
integer :: i, mymax
real :: x, myeps
rec(0:mymax) = BESSEL_JN(0, mymax, x)
lib(0:mymax) = [(BESSEL_JN(i, x), i = 0, mymax)]
do i = 0, mymax
if (rec(i) == lib(i)) cycle
if (abs((rec(i) - lib(i)) / rec(i)) > myeps) then
print *, 'FAIL', mymax, x, i, rec(i), lib(i),
abs((rec(i)-lib(i))/rec(i)), myeps
stop 1
end if
end do
end
end
Built and run with:
$ gcc-build/gcc/gfortran -B gcc-build/gcc /tmp/bessel_6_debug.f90 \
-L gcc-build/x86_64-pc-linux-gnu/libgfortran/.libs \
-Wl,-rpath,/home/ert/code/gcc-dev/gcc-build/x86_64-pc-linux-gnu/libgfortran/.libs
\
-o /tmp/bessel_6_debug
$ /tmp/bessel_6_debug
STOP 1
FAIL 37 475.779999 30 3.96264251E-04
3.96265590E-04 3.37849360E-06 1.78813934E-06
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 looks related to the historical sensitivity of bessel_6.f90 to target/libm
behavior, but the failing path here is specifically the libgfortran REAL(4)
recurrence implementation.
System configuration:
- GCC source: upstream/master at cfa0910d9c0fa307b3156f247c50044524b3d5b2
(2026-04-08)
- Local compiler under test: gcc version 16.0.1 20260403 (experimental)
- Configure line: ../gcc/configure --enable-languages=fortran
--disable-multilib --disable-bootstrap CFLAGS='-Og -g' CXXFLAGS='-Og -g'
- OS distro: CachyOS Linux
- Kernel: Linux mailuefterl 6.19.10-1-cachyos x86_64
- CPU: AMD Ryzen 9 5950X 16-Core Processor
- libc/libm: GNU libc 2.43, using /usr/lib/libm.so.6
I have not yet determined whether using wider intermediates in bessel_jn_r4 is
sufficient on all systems, but that seems like the first implementation change
to try.
