Taking into account your suggestions I came up with an ad-hoc
implementation tailored for my needs. Instead of computing
Laguerre(n,1,t-n) I compute the whole expression
Laguerre(n,1,t-n)*exp(-(t-n))*(t-n) via a function laguerrel_xexpx.
Because I incorporated the exponential factor into the recursive formula
for laguerre polynomials the NaNs are no longer so significant. Even if
they occur they are due to underflows, which can be simply converted to
exact zeros. The Julia implementation is now roughly order of magnitude
faster than Mathematica and still works for arbitrarily large t. Many
thanks! Also, if you have any suggestions for further performance
improvement please let me know.
function Q2(t::Number)
s=sum([one(t)/(1+n)*laguerrel_xexpx(n,1,t-n) for
n=0:int(floor(t))])-exp(-t-1)
return(s)
end
# computes the LaguerreL(n,alpha,x)*exp(-x)*x
function laguerrel_xexpx(n::Integer, alpha::Number, x::Number)
l1 = float(1+alpha-x) # L_{1}
l0 = one(11) # L_{0}
if(n==0)
ln=l0*x*exp(-x)
elseif(n==1)
ln=l1*x*exp(-x)
else
ln = l1*exp(-x/n) # initial L_{n}
ln1= l0*exp(-x/n) # initial L_{n-1}
ln2= l0
for i=2:n
ln2=ln1*exp(-x/n) # new L_{n-2}
ln1=ln*exp(-x/n) # new L_{n-1}
ln=((2*i+alpha-1-x)*ln1-(i+alpha-1)*ln2)/i # recursive formula
for L_{n}
end
ln=ln*x
end
if(isnan(ln)) # probably underflow, so set ln to zero
return(zero(ln))
else
return(ln)
end
end
W dniu czwartek, 20 marca 2014 16:28:28 UTC+1 użytkownik Steven G. Johnson
napisał:
>
> For example, you may be able to use the analytical asymptotic form of the
> Laguerre polynomial for large n. Computing special functions efficiently
> is all about switching between different asymptotic forms, recurrences,
> etcetera, for different regions of the parameter space.
>
