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

From: Robin Hankin <r.hankin_at_noc.soton.ac.uk>
Date: Tue, 19 Jun 2007 10:25:27 +0100

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

______________________________________________
R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Received on Tue 19 Jun 2007 - 09:29:35 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 Fri 22 Jun 2007 - 19:35:11 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.