数值计算阶乘和多项式的组合 [英] Numerically calculate combinations of factorials and polynomials

问题描述

我正在尝试编写一个简短的C ++例程，以针对给定的整数j> i(通常位于0到100之间)和复数z(以| z |为界)计算以下函数F(i，j，z).< 100)，其中L是

问题是我希望可以在CUDA内核中调用此函数(即具有 __ device __ 属性).因此，标准库/Boost/etc函数是不可能的，除非它们足够简单以至于我自己可以重新实现-这尤其与Boost和C ++ 17中存在的Laguerre多项式有关.无论我是否设法为Laguerre多项式包装任何标准函数，我都仍然具有类似的前置因子来计算形式(z ^ j/j！).

问题:如何在不引入明显的数值不稳定性的情况下实现这种功能的相对简单的实现?

到目前为止，

我的想法是独立计算L及其前置因子.我将通过首先从0到j-i循环计算并计算(z ^ 1 * z ^ 2/2 * ... * z ^(j-1)/(j-i！)！).然后，我将计算剩余因子exp(-| z | ^ 2/2)*(j-i)！* sqrt(i！/j！)(以类似方式或通过CUDA数学中实现的Gamma函数).然后的想法是找到一种最小的算法来计算相关的Laguerre多项式，除非我设法从一个例子中包装一个实现.Boost或GNU C ++.

编辑/旁注:对于某些i/j值，F的表达式实际上在数字上会爆炸.在我得到它的源中，它是错误推导的，而关联的Laguerre多项式的索引应该改为L_i ^(j-i).但这不会使答案/评论中建议的方法无效.

解决方案

我建议为Laguerre多项式的系数找到递归关系:

  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))//自己验证:) 

这使您在计算多项式时避免了大多数阶乘.

此后，我将遵循Severin的对数计算方法为了不使双浮点数范围超载:

  log(F)= log(sqrt(i！/j！))-| z | ^ 2 +(j-i)* log(-z)+ log(L(| z | ^ 2))log(L)= log((2 * j-i)！)+ log(sum)//其中使用上述递归关系计算总和 

并利用以下事实:

  log(a！)= sum(k = 1..a，log(k)) 

还有:

  log(z)= log(| z |)+ I * arg(z)对于复数zlog(-z)= log(| z |)+ I * arg(-z)log(-z)= log(| z |)-I * arg(z) 

对于 log(sqrt(i！/j！))部分，我会做(假设j> = i):

 日志(sqrt(i！/j！))= 0.5 *(log(i！)-log(j！))= -0.5 * sum(k == i + 1..j，log(k)) 

我还没有尝试过，所以在这里和那里肯定会有一些小错误.这个答案更多是关于技术，而不是复制粘贴准备好答案

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:

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.

解决方案

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

这篇关于数值计算阶乘和多项式的组合的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！