Re: [Rd] two bessel function bugs for nu<0

From: Martin Maechler <maechler_at_stat.math.ethz.ch>
Date: Fri, 22 Jun 2007 18:38:08 +0200

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

______________________________________________
R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Fri 22 Jun 2007 - 16:52:16 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sat 23 Jun 2007 - 01:35:20 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.