>>>>> Mikael Jagan
>>>>> on Mon, 29 Dec 2025 16:10:24 -0500 writes:
> These lines (561-567) of src/nmath/bessel_j.c seem relevant. :-)
> L250:
> /* ---------------------------------------------------
> Normalize. Divide all b[N] by sum.
> ---------------------------------------------------*/
> /* if (nu + 1. != 1.) poor test */
> if(fabs(nu) > 1e-15)
> sum *= (Rf_gamma_cody(nu) * pow(.5* *x, -nu));
> Mikael
Indeed; thank you, Mikael, and also Leo Mada, Ben Bolker, Richard O'Keefe.
And yes, clearly a buglet.
Yes, R's bessel{IJKH}() algorithms are not perfect, even when
based on published code/algortithms and have already been
tweaked to be more "robust" than the original code.
This "imperfectness" where the main reason I had created the
{Bessel} CRAN package ( https://cran.R-project.org/package=Bessel )
a long time ago in order to provide alternative algorithms -
notably to work with *complex* 'x', but also to explore pure R
code implementations of asymptotic approximations, (almost?)
always related to |x| -> oo or |nu| --> oo ("oo" = `Inf`).
The case here, |nu| --> 0 , is a new not-yet reported
bug which should be addressed and fixed.
Note that the Bessel's package versions do not seem to suffer
from the problem here.
(nuSml <- 2^-c(seq(30, 53, by=1/2), 75, 100, 300, 800, 1022, 1074, Inf))
## bug in R :
Jsml <- sapply(nuSml, \(nu) besselJ(pi/2, nu))
options(digits = 14) # show more digits
tail(cbind(nuSml, Jsml), 20)
## nuSml Jsml
## [35,] 7.1054273576010e-15 4.7200121576824e-01
## [36,] 5.0242958677881e-15 4.7200121576824e-01
## [37,] 3.5527136788005e-15 4.7200121576824e-01
## [38,] 2.5121479338940e-15 4.7200121576824e-01
## [39,] 1.7763568394003e-15 4.7200121576823e-01
## [40,] 1.2560739669470e-15 4.7200121576823e-01
## [41,] 8.8817841970013e-16 5.3142612486306e+14
## [42,] 6.2803698347351e-16 7.5155003318072e+14
## [43,] 4.4408920985006e-16 1.0628522497261e+15
## [44,] 3.1401849173676e-16 1.5031000663614e+15
## [45,] 2.2204460492503e-16 2.1257044994522e+15
## [46,] 1.5700924586838e-16 3.0062001327229e+15
## [47,] 1.1102230246252e-16 1.0000000000000e+00
## [48,] 2.6469779601697e-23 1.0000000000000e+00
## [49,] 7.8886090522101e-31 1.0000000000000e+00
## [50,] 4.9090934652977e-91 1.0000000000000e+00
## [51,] 1.4996968138956e-241 1.0000000000000e+00
## [52,] 2.2250738585072e-308 1.0000000000000e+00
## [53,] 4.9406564584125e-324 1.0000000000000e+00
## [54,] 0.0000000000000e+00 4.7200121576823e-01
if(!require("Bessel")) { install.packages("Bessel"); library(Bessel) }
bJsml <- sapply(nuSml, \(nu) BesselJ(pi/2, nu))
tail(cbind(nuSml, bJsml), 20) # shows that all is fine here
-------------------------
I'm looking into fixing this,
a first good step is probably to just set nu=0 when |nu| is too
small, but I have to first find, i
but you (readers) could help well
by exploring the nu -> 0 case over "all" x instead of just
x = pi/2. (notably large x).
Also, I've just experienced my first really positive interaction
with an LLM, starting from mathoverflow / stackoverflow, ending
here:
https://stackoverflow.com/ai-assist/shared/46dd1d0b-7da7-46ad-b2c5-e63eca35c721
In any case, thank you once more for the report!
Martin
--
Martin Maechler
ETH Zurich and R Core team
______________________________________________
[email protected] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide https://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
