[Date Index] [Thread Index] [Author Index]

It's worse than you think. Re: BesselJ[] with machine-precision input

• To: mathgroup at smc.vnet.net
• Subject: [mg128417] It's worse than you think. Re: BesselJ[] with machine-precision input
• From: Richard Fateman <fateman at cs.berkeley.edu>
• Date: Sun, 14 Oct 2012 23:40:00 -0400 (EDT)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com
• Delivered-to: l-mathgroup@wolfram.com
• Delivered-to: mathgroup-newout@smc.vnet.net
• Delivered-to: mathgroup-newsend@smc.vnet.net
• References: <k5dem6$sfn$1@smc.vnet.net>

First of all, there's a lot more wrong with Bessel than you noticed.

try this on Mathematica 8:

q = 429515858585961022071539/6922263581864661506963;

a = N[BesselJ[0, q], 45];
b = N[BesselJ[0, q], 100];
c = N[BesselJ[0, q], 45];

You might expect a to be the same as c, but it isn't.

You might think that a or c would look like the first 45 digits of
b, but that's not true either.

That reveals another peculiarity of Mathematica, and not
directly something about Bessel:

  Note that q was chosen so that
Rationalize[N[q,45],0] == q.

Peculiarly, then
SetPrecision[N[q,45],100]-N[q,100] is not zero but about -3.9E-59


Secondly, if you are interested in fast evaluation of bessel function 
J_0(x) to double-float precision and no more, for 0<x<20, there are
much simpler computations.  Even for 30 digits correct, you can
do better. For more of an education, see for example ACM algorithm 236 
by Gautschi. For less of an education, see Numerical Recipes.

Now if you are comparing Mathematica to R, in arbitrary precision,
then perhaps you should look at Rmpfr:

http://cran.r-project.org/web/packages/Rmpfr/vignettes/Rmpfr-pkg.pdf

and then you would be more accurately comparing software with
similar capabilities: MPFR Bessel algorithms to Mathematica.

I think those MPFR algorithms can be considerably sped up. I suppose
Mathematica could use those.  If you want to make your own,
remember that they must work for any size or precision argument, so your 
interpolation formula will not work.

Mathematica's methods are not public, and according to your tests are 
apparently not fast or accurate; by my test the results are not
necessarily reproducible. I think we can summarize by saying these
programs could be improved.

RJF