two bessel function bugs for nu<0
I can reproduce both these bugs and confirm that the suggested fix
agrees with Mathematica and Maple for a few trial values.
I can confirm that Hiroyuki's algebra is indeed
consistent with AMS-55 equation 9.1.2
and the old source isn't. I'd need more
time to look at equation 9.6.2.
I'm not sure why, in bessel_i.c, we are using a float ("expo")
and a long ("ize") as a Boolean [flag to indicate whether or not
to return scaled function values].
PS1:
My first thought was to check against the GSL library
but this doesn't allow non-integer orders for besselI()
PS2: The source code apologizes for the method used,
suggesting that it may be numerically and computationally
"sub-optimal".
Best wishes
rksh
On 18 Jun 2007, at 23:33, Hiroyuki Kawakatsu wrote:
#bug 1: besselI() for nu<0 and expon.scaled=TRUE #tested with R-devel (2007-06-17 r41981) x <- 2.3 nu <- -0.4 print(paste(besselI(x, nu, TRUE), "=", exp(-x)*besselI(x, nu, FALSE))) #fix: #$ diff bessel_i_old.c bessel_i_new.c #57c57 #< bessel_k(x, -alpha, expo) * ((ize == 1)? 2. : 2.*exp(-x))/M_PI #--- #> bessel_k(x, -alpha, expo) * ((ize == 1)? 2. : 2.*exp(-2.0*x))/ M_PI #bug 2: besselY() for nu<0 #don't know how to check in R; a few random checks against mathematica 5.2 #fix: #$ diff bessel_y_old.c bessel_y_new.c #55c55 #< return(bessel_y(x, -alpha) + bessel_j(x, -alpha) * sin(-M_PI * alpha)); #--- #> return(bessel_y(x, -alpha) * cos(M_PI * alpha) - bessel_j(x, -alpha) * sin(M_PI * alpha)); h. -- ---------------------------------- Hiroyuki Kawakatsu Business School Dublin City University Dublin 9, Ireland Tel +353 (0)1 700 7496
______________________________________________ R-devel at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
-- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743