On Thursday, March 20, 2014 8:26:56 AM UTC-4, Paweł Biernat wrote:
>
> 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.
>
The basic problem here is that laguerre_n(2000,1,2000) is about 6.55e431,
which is larger than the maximum representable Float64 (about 1e308), so it
is literally impossible to compute this function in double precision.
However, you don't actually need to compute laguerre_n. Your summand is:
laguerre_n(n,1,t-n) * exp(-(t-n))*(t-n)/(1+n)
which includes an exponential scaling that cancels the exponential growth
of the Laguerre polynomial.
What you should do is to rewrite your Laguerre recurrence to compute the
exponentially scaled Laguerre polynomial, i.e. including the exp(-(t-n))
factor, rather than the Laguerre polynomial by itself. This will require
some thought; it's not just a matter of multiplying by exp(-(t-n)), because
that will underflow.
(This is the kind of thought process that leads to all sorts of rescaled
special functions, such as erfcx, exp1p, and lngamma.)
