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