Bill Allombert on Mon, 13 Feb 2012 15:00:24 +0100


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: bessel function algorithm


On Mon, Feb 13, 2012 at 10:24:45AM -0200, Marcel Bezerra wrote:
> 
> Henri, I would like to thank you again for the help.
> 
> Do you know where I could find the algorithm used in PARI for Bessel
> function of the first kind?
> 
> I want to study that algorithm, but I couldn't find anything about it in
> the PARI's documentation pages.

I think we just use the standard formula:

J_alpha(x) = sum_{m=0}^oo (-1)^m/(m!*Gamma(m+alpha+1)) (x/2)^(2*m+alpha)

You can check the source code if you like (in src/basemath/trans3.c).

> I am aware there are some MATLAB programs that intend to calculate it. I
> made some comparisons with Hongxue Cai's algorithm and PARI, and I noticed
> the results are different after the fourth decimal place.
> 
> The values I tested were: besselj(1+I,0); besselj(1+I,1+I); besselj(1+I,I).
> 
> The results I got from MATLAB using the Hongxue's cbessel(nu,z) algorithm
> are, respectively:
> 
> NaN+NaN*i; 3.841170714997851e-001 -9.914709109073297e-002i;
> 1.464168620011011e-001 +5.274551029538346e-002i

This suggests that these values are computed using C double complex. PARI use
multiprecision and is vastly more accurate. You can increase the precision and
check that the results matches.

Cheers,
Bill.