Thanks Andrew,
On linux, glibc must be doing something defensive at zeros
of J0(x). On FreeBSD, with -O0, -O1, and -O2 I get
% gfcx -o z -O0 bessel.f90 && ./z
900691. 182038. 6333. 48537. 203435. 364229. 14664. 5813. 18993. 1306.
% gfcx -o z -O1 bessel.f90 && ./z
900691. 182038. 6333. 48537. 203435. 364229. 14664. 5813. 18993. 1306.
% gfcx -o z -O2 bessel.f90 && ./z
900691. 182038. 6333. 48537. 203435. 364229. 14664. 5813. 18993. 1306.
which is what I expected. The above values are estimates of ULP.
900691 corresponds to a single correct digit!
With -O3, I get
% gfcx -o z -O3 bessel.f90 -fdump-tree-optimized && ./z
0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
which is not at all what I expected! In fact, if I look for symbols in
z there is no call to j0f in the math library.
% gfcx -o z -O bessel.f90
% nm z | grep -i j0
U j0f@FBSD_1.0
% gfcx -o z -O3 bessel.f90
% nm z | grep -i j0
So, the Bessel function is optimized out at -O3.
Scratching head???
--
steve
On 9/11/26 22:17, Andrew Benson wrote:
Hi Steve
Using gfortran 16.0.1 on Linux with -O2 or lower:
8. 0. 0. 0. 0. 0. 0. 0. 0. 0.
Using -O3:
0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
-Andrew
--
* Andrew Benson: http://users.obs.carnegiescience.edu/abenson <http://
users.obs.carnegiescience.edu/abenson>
* Galacticus: https://github.com/galacticusorg/galacticus <https://
github.com/galacticusorg/galacticus>
On Fri, Sep 11, 2026, 9:56 PM Steve Kargl <[email protected]
<mailto:[email protected]>> wrote:
All,
I can test FreeBSD, but have no easy access to other
operating systems. Can others compile the following
program, run the executable, and either post your
results here or send them to directly.
Thanks.
program besseltest
implicit none
integer i
integer n(10)
real(4) x(10), f(10)
real(4), parameter :: e(10) = [ &
& 5.64343985e-8, 2.44766589e-8, 1.03553056e-7, -3.53017149e-9, &
& -6.48150511e-9, 5.15323162e-9, -1.50233475e-7, 1.25245705e-7, &
& 5.43311032e-8, 1.22055454e-7]
x = [ 2.40482556, 5.52007811, 8.65372791, 11.7915344, 14.9309177, &
& 18.0710640, 21.2116366, 24.3524715, 27.4934791, 30.6346065]
f = bessel_j0(x)
n = exponent(e);
f = abs(e - f);
f = scale(f, 24 - n);
print '(*(F0.0,1x))', f;
end program besseltest
--
steve
