我正在尝试编写一个简短的 C++ 例程来计算给定整数 j > i （通常它们位于 0 到 100 之间）和复数 z （以 |z| 关联拉盖尔多项式:
问题是我希望这个函数可以从 CUDA 内核中调用（即使用__device__属性）。因此，标准库/Boost/等函数是不可能的，除非它们足够简单，可以自己重新实现——这尤其与 Boost 和 C++17 中存在的拉盖尔多项式有关。无论我是否设法包装拉盖尔多项式的任何标准函数，我仍然有一个类似的预因子来计算形式 (z^j/j!)。
问题：如何在不引入显着数值不稳定的情况下相对简单地实现这样的函数？
My idea到目前为止就是独立计算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)。这并不会使答案/评论中建议的方法无效。
我建议找到拉盖尔多项式系数的递推关系：
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)
for the 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))
我还没有尝试过，所以肯定会有一些小错误。这个答案更多的是关于技术，而不是复制粘贴就绪的答案
本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)
