一个用于计算“多项式系数”的numpy/scipy函数

Question

有没有Python函数（可能来自numpy或scipy）可以计算在扩展式(1+x+x**2+x**3+...+x**(k-1))**n中x**r的系数，其中k>=1，n>=0，且0<=r<=n(k-1)？

有时这被称为多项式系数(PC)（例如，参见此处）。

如果没有，您能想到一种有效的方法来计算它吗？（我不感兴趣的是朴素/贪婪的方法。）

- Bach

---

我认为这个问题可以用封闭形式解决，也许更适合在数学交流网站上讨论。 - wim

或者查看http://docs.scipy.org/doc/numpy/reference/generated/numpy.polynomial.polynomial.polyfit.html - Louis

@Louis - 你确定你提供的链接是相关的吗？ - Bach

在您看来，什么是天真/贪心的方式？ - Fred Foo

@larsmans - 一种天真的方法是循环遍历从0,...,k-1长度为n的所有组合，并计算每个组合的总和是否为r。 - Bach

哦，对了，我误读了问题。我以为你只是想评估那个公式。 - Fred Foo

Answer 1

你实际上正在对 [1, 1, 1, ..., 1, 1, 1] 进行 n 次卷积。
因此，你是否考虑在足够补零的数组上使用 FFT，将其元素提高到幂次为 n，然后使用反向 FFT 来恢复所有系数。

(1+x+x**2+x**3+...+x**(k-1))**n

然后只需阅读你感兴趣的那些内容吗？

更新:

由于FFT是循环的，您需要一个数组，其大小不小于要在其中执行FFT的项数。

(1+x+x**2+x**3+...+x**(k-1))**n

或者说，(k-1)*n+1，这样结果就不会在末尾环绕（或者至少是只在受影响的元素上添加零）。通常其长度也应该是2的幂次方，因为这是FFT算法所要求的（如果实现不需要，则会在输入中填充零直到满足此要求）。

在类C的伪代码中：

unsigned int m = 1;
while(m<(k-1)*n+1) m *= 2;

complex c[m];
for(unsigned int i=0;i!=k;++i) c[i] = complex(1.0, 0.0);
for(unsigned int i=k;i!=m;++i) c[i] = complex(0.0, 0.0);

c = fft(c);
for(unsigned int i=0;i!=m;++i) c[i] = pow(c[i], double(n));
c = inv_fft(c);

在此过程结束时，复数数组 c 的第 r 个元素具有实部等于 x ** r 的系数，虚部为零。
现在，由于所有操作都是浮点运算，因此您需要注意这些元素将累计舍入误差。您可以通过将它们四舍五入到最接近的整数来部分纠正这种误差，但请注意，对于足够大的 k 和 n ，这些误差将超过0.5，因此可能会产生一些小的相对误差。
在网上快速搜索后发现，numpy有FFT及其反转的实现，分别是 numpy.fft.rfft 和 numpy.fft.irfft ，当输入数据为实数时，您可以使用它们。

- thus spake a.k.

---

谢谢你的回复。我在傅里叶变换领域的知识很差，所以我不确定如何将你的想法转化为代码... - Bach

@Bach：我已经更正了一个错误，之前我说你应该将项乘以“n”，但实际上你应该将它们提高到“n”的幂。我还添加了一些伪代码，希望能让它更清晰明了。 - thus spake a.k.