[Rd] two bessel function bugs for nu<0

[Martin Maechler]

Thank you both, Hiroyuki and Robin,
 
>>>>> "Robin" == Robin Hankin <r.hankin at noc.soton.ac.uk>
>>>>>     on Tue, 19 Jun 2007 10:25:27 +0100 writes:

    Robin> I can reproduce both these bugs and confirm that the suggested fix
    Robin> agrees with Mathematica and Maple for a few trial values.

    Robin> I can confirm that Hiroyuki's algebra is indeed
    Robin> consistent with AMS-55 equation 9.1.2
    Robin> and the old source isn't.   I'd need more
    Robin> time to look at equation 9.6.2.

    Robin> I'm not sure why, in bessel_i.c, we are using a float ("expo")
    Robin> and a long ("ize") as a Boolean [flag to indicate whether or not
    Robin> to return scaled function values].

    Robin> PS1:
    Robin> My first thought was to check against the GSL library
    Robin> but  this doesn't allow non-integer orders for besselI()

    Robin> PS2: The source code apologizes for the method used,
    Robin> suggesting that it may be numerically and computationally
    Robin> "sub-optimal".

I had wanted to address the very succinct and fine bug report by
Hiroyuki, but have been diverted by other matters. 
Thank you for the confirmation.

I'm currently testing the changes and plan to commit them, also
for the 2.5.1 release candidate.

Thank you, once more.

With best regards,
Martin 


    Robin> 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
    >> 
    >> ______________________________________________

    Robin> --
    Robin> Robin Hankin
    Robin> Uncertainty Analyst
    Robin> National Oceanography Centre, Southampton
    Robin> European Way, Southampton SO14 3ZH, UK
    Robin> tel  023-8059-7743