4

I am trying to write a short C++ routine to calculate the following function F(i,j,z) for given integers j > i (typically they lie between 0 and 100) and complex number z (bounded by |z| < 100), where L are the associated Laguerre Polynomials:

enter image description here

The issue is that I want this function to be callable from within a CUDA kernel (i.e. with a __device__ attribute). Standard library/Boost/etc functions are therefore out of the questions, unless they are simple enough to re-implement on my own - this especially relates to the Laguerre polynomials which exist in Boost and C++17. Regardless if I manage to wrap any standard function for Laguerre polynomials, I still have a similar pre-factor to calculate of the form (z^j/j!).

Question: How can I do a relatively simple implementation of such a function, without introducing significant numerical instability?

My idea so far is to calculate L and its pre-factor independently. The pre-factor I will calculate by first looping from 0 to j-i and calculate (z^1 * z^2/2 * ... * z^(j-1)/(j-i)!). I will then calculate the remaining factor exp(-|z|^2/2) *(j-i)! * sqrt(i!/j!) (either in a similar way, or through the Gamma-function, which is implemented in CUDA math). The idea is then to find a minimal algorithm to calculate the associated Laguerre polynomial, unless I manage to wrap an implementation from e.g. Boost or GNU C++.

Edit/side note: The expression for F actually blows up numerically for some values of i/j. It was derived wrong in the source where I got it, and the indices of the associated Laguerre polynomials should instead be L_i^(j-i). That does not invalidate the approaches suggested in the answers/comments.

fromGiants
  • 474
  • 1
  • 7
  • 16

2 Answers2

2

I recommend finding a recurrence relation for the coefficients of the Laguerre Polynomial:

C(k+1) = g(k)C(k)
g(k) = C(k+1) / C(k)
g(k) = -z * (j - k) / ((j - i + k + 1) * (k + 1)) //Verify this yourself :)

This allows you to avoid most of factorials in computing the polynomial.

After that I would follow Severin's idea of doing the calculations in logarithms so as to not overload the double floating point range:

log(F) = log(sqrt(i!/j!)) - |z|^2 + (j-i) * log(-z) + log(L(|z|^2))
log(L) = log((2*j - i)!) + log(sum) // where the summation is computed using the recurrence relation above

and using the fact that:

log(a!) = sum(k=1..a, log(k))

and also:

log(z) = log(|z|) + I * arg(z) for complex z
log(-z) = log(|z|) + I * arg(-z)
log(-z) = log(|z|) - I * arg(z)

for the log(sqrt(i!/j!)) part I would do (assuming that j >= i):

  log(sqrt(i!/j!))
= 0.5 * (log(i!) - log(j!))
= -0.5 * sum(k==i+1..j, log(k))

I haven't tried this out so there could definitely be little mistakes here and there. This answer is more about the technique rather than a copy-paste-ready answer

SirGuy
  • 10,660
  • 2
  • 36
  • 66
1

Well, what you should do is to logarithm it

Assuming natural logarithm,

q = log(z^j/j!) = log(z^j) - log(j!) = j*log(z) - log(Gamma(j+1))

First term is simple, second term is standard C++ function lgamma(x) (or you could use GSL).

compute value of q and return cexp(q)

You could fold exponent in this method as well

Severin Pappadeux
  • 18,636
  • 3
  • 38
  • 64