MathGroup Archive 2012

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

Search the Archive

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);






  • Prev by Date: BesselJ[] with machine-precision input gives sub-machine-precision output
  • Next by Date: Re: plot function of two variables
  • Previous by thread: BesselJ[] with machine-precision input gives sub-machine-precision output
  • Next by thread: It's worse than you think. Re: BesselJ[] with machine-precision input