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.