Hi,
I have just begun learning Julia and, as an exercise, I tried to
re-implement an algorithm that I had previously implemented in
Mathematica. The problem consists of calculating the following expression
Q(t)=-e^{-t-1}+sum _{n=0}^{floor(t)} e^{n-t} (t-n) L_n^1(t-n)/(n+1)
for large t (here L_n^1 are generalized Laguerre
polynomials<http://en.wikipedia.org/wiki/Laguerre_polynomials>
).
I was pleasantly surprised to find that my straight forward implementation
(included below, see the definition of Q2) in Julia runs twice as fast as
the Mathematica code for small t (t<1500). As a benchmark I tried to
compute Q(1500.0) which takes ~0.18s in Julia and ~0.32s in Mathematica.
For larger t, I have to use BigFloat (or I get NaNs otherwise) which
significantly decreases performance of Julia code (e.g. it takes ~4.2s to
compute Q(BigFloat(1500))). In my implementation the use of BigFloat seems
to be necessary as I would like to compute values of Q for t as large as
10^4. On the other hand, Mathematica still gives valid results for large t
using only 64bit floats (maybe by employing a more sophisticated algorithm
to compute Laguerre polynomials or by using higher precision internally and
only where it is actually needed?). So, effectively, as t increases
Mathematica, with its 64bit float, becomes much faster than my Julia code,
which has to use BigFloat. For example, Q(BigFloat(3000)) in Julia takes
~17s while in Mathematica it takes only ~3s to compute Q(3000.0).
I also tested the speed of generating the Laguerre polynomials themselves,
which gives
Mathematica:
Timing[LaguerreL[#, 1, #] & /@ Range[1000];]
{0.456666, Null}
Timing[LaguerreL[#, 1, #] & /@ Range[2000];]
{2.366666, Null}
Julia (implementation of laguerre_recursive included below)
@elapsed [laguerrel_recursive(n,1,n) for n=1:1000]
0.10840264
@elapsed [laguerrel_recursive(n,1,BigFloat(n)) for n=1:2000] #here I have
to use BigFloat or I get NaNs
8.334336062
So it seems that BigFloat is the major drawback. I really wander how
Mathematica handles the generation of Laguerre polynomials. I am aware of
the GSL bindings and the function sf_laguerre_n, but it also returns NaNs
for sf_laguerre_n(2000,1,2000). In fact, sf_laguerre_n generates Laguerre
polynomials in the same way as laguerre_recursive.
I would like to know if it is possible to improve the performance of my
code to a point where it is, at least, as fast as Mathematica (possibly by
removing the need for BigFloat). My implementation of Q(t) in Julia:
function Q(t::Number)
s=sum([exp(-(t-n))*(t-n)/(1+n)*laguerrel_recursive(n,1,t-n) for
n=0:int(floor(t))])-exp(-t-1)
return(s)
end
# implementation via recursive definition of laguerre polynomials
function laguerrel_recursive(n::Integer, alpha::Number, x::Number)
l0 = 1 # L_{0}
l1 = 1+alpha-x # L_{1}
if(n==0)
return(l0)
elseif(n==1)
return(l1)
else
ln = l1 # initial L_{n}
ln1= l0 # initial L_{n-1}
for i=2:n
ln2=ln1 # new L_{n-2}
ln1=ln # new L_{n-1}
ln=((2*i+alpha-1-x)*ln1-(i+alpha-1)*ln2)/i # recursive formula
for L_{n}
end
return(ln)
end
end
And the analogous implementation in Mathematica
Q[t_] := Sum[Exp[-(t - n)] (t - n)/(1 + n) LaguerreL[n, 1, t - n], {n, 0,
Floor[t]}] - Exp[-t - 1]
Paweł Biernat
