数值计算阶乘和多项式的组合

Question

我正在尝试编写一个简短的 C++ 例程来为给定的整数 j > i（通常它们位于 0 和 100 之间）和复数 z（以 | 为界）计算以下函数 F(i,j,z) z| < 100)，其中 L 是 associated Laguerre Polynomials:

问题是我希望这个函数可以从 CUDA 内核中调用（即具有 __device__ 属性）。因此，标准 library/Boost/etc 函数是不可能的，除非它们足够简单到我自己 re-implement - 这尤其与存在于 Boost 和 C++17 中的拉盖尔多项式有关。不管我是否设法为拉盖尔多项式包装任何标准函数，我仍然有一个类似的 pre-factor 来计算形式 (z^j/j!).

问题：我怎样才能相对简单地实现这样一个函数，而不引入明显的数值不稳定性？

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

Edit/side 注意： 对于 i/j 的某些值，F 的表达式实际上在数值上爆炸了。它在我得到它的来源中是错误的，相关的拉盖尔多项式的索引应该是 L_i^(j-i)。这并不会使 answers/comments.

中建议的方法无效

Answer 1

嗯，你应该做的是对它取对数

假设自然对数，

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

第一项很简单，第二项是标准 C++ 函数 lgamma(x)（或者您可以使用 GSL）。

计算 q 和 return cexp(q)

的值

你也可以用这个方法折叠指数

Answer 2

我建议找到拉盖尔多项式系数的递归关系：

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 :)

这可以让您在计算多项式时避免大部分阶乘。

之后我会按照 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) // where the summation is computed using the recurrence relation above

并利用以下事实：

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

还有：

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

对于 log(sqrt(i!/j!)) 我会做的部分（假设 j >= i）：

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

我还没有试过这个，所以肯定会有一些小错误。这个答案更多的是关于技术而不是复制粘贴就绪的答案

数值计算阶乘和多项式的组合

Numerically calculate combinations of factorials and polynomials

c++

algorithm

numerical-methods

factorial

polynomials