BesselK, BesselI, BesselJ, BesselY are unnecessarily slow
- To: mathgroup at smc.vnet.net
- Subject: [mg128413] BesselK, BesselI, BesselJ, BesselY are unnecessarily slow
- From: yenleeloh1 at gmail.com
- Date: Sun, 14 Oct 2012 00:14:42 -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
Dear Mathgroup readers:
I did timings on Mathematica's implementation of BesselK[n,x], BesselI, BesselJ, and BesselY[n,x] for n=0 and 0<x<10 using the following code:
xlist = Range[0, 10, 0.0001];
nxmax = Length@xlist;
First@Timing[ BesselK[0, xlist] ]/nxmax/10^-6 \[Mu]s // SetAccuracy[#, 2] &
First@Timing[ BesselI[0, xlist] ]/nxmax/10^-6 \[Mu]s // SetAccuracy[#, 2] &
First@Timing[ BesselJ[0, xlist] ]/nxmax/10^-6 \[Mu]s // SetAccuracy[#, 2] &
First@Timing[ BesselY[0, xlist] ]/nxmax/10^-6 \[Mu]s // SetAccuracy[#, 2] &
I am using Mathematica 8.0.1.0 on Mac OS X x86 (32-bit, 64-bit kernel). The timings were
BesselK[] 2.6 microseconds/call
BesselI[] 17.3 microseconds/call
BesselJ[] 19.0 microseconds/call
BesselY[] 144.1 microseconds/call
The implementation in R (R 2.15.0 GUI 1.51 Leopard build 64-bit) is much faster. In particular, R's besselY() function is 50 times faster than Mathematica's BesselY[] function.
besselK() 0.75 microsecs/call
besselI() 0.73 microsecs/call
besselJ() 0.75 microsecs/call
besselY() 3.04 microsecs/call
Is Mathematica's implementation really that lame, or am I missing something? Does anyone have info on timings with GSL, Numerical Recipes, etc? Please reply by email. Thank you.
My R code is below:
xlist <- seq(0.0001,10,0.0001)
nxmax <- length(xlist)
sprintf ("%1$g microsecs", system.time ( besselK(xlist, 0) )[1]/nxmax/10^-6);
sprintf ("%1$g microsecs", system.time ( besselI(xlist, 0) )[1]/nxmax/10^-6);
sprintf ("%1$g microsecs", system.time ( besselJ(xlist, 0) )[1]/nxmax/10^-6);
sprintf ("%1$g microsecs", system.time ( besselY(xlist, 0) )[1]/nxmax/10^-6);