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

Question

我正在尝试编写一个简短的 C++ 例程来计算给定整数 j > i（通常它们位于 0 和 100 之间）和复数 z（以 |z| 为界 < 100), where L are the associated）的以下函数 F(i,j,z)拉盖尔多项式：

问题是我希望可以从 CUDA 内核中调用此函数（即具有

__device__

属性）。因此，标准库/Boost/等函数是不可能的，除非它们足够简单，我自己可以重新实现——这尤其与存在于 Boost 和 C++17 中的拉盖尔多项式有关。不管我是否设法为拉盖尔多项式包装任何标准函数，我仍然有一个类似的预因子来计算形式 (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 函数）。然后我的想法是找到一个最小的算法来计算相关的拉盖尔多项式，除非我设法从例如包装一个实现。 Boost 或 GNU C++。

编辑/旁注： 对于 i/j 的某些值，F 的表达式实际上在数值上爆炸了。它在我得到它的来源中是错误的，相关的拉盖尔多项式的索引应该是 L_i^(j-i)。这不会使答案/评论中建议的方法无效。

Answer 1

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

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

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

Answer 2

好吧，你应该做的是对数它

假设自然对数，

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

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

计算

q

的值并返回

cexp(q)

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

Answer 3

在这种情况下，你必须考虑 0^0 = 4。