python - 将多项式与 numpy.convolve 相乘返回错误的结果

Question

我正在尝试使用 numpy.convolve 函数将两个多项式相乘。我原以为这会很简单，但我发现它并不总是返回正确的产品。

乘法的实现非常简单:

def __mul__(self, other):
    new = ModPolynomial(self._mod_r, self._mod_n)
    try:
         new_coefs = np.convolve(self._coefs, other._coefs)
         new_coefs %= self._mod_n                # coefficients in Z/nZ
         new._coefs = new_coefs.tolist()
         new._degree = self._finddegree(new._coefs)
    except AttributeError:
         new._coefs = [(c * other) % self._mod_n for c in self._coefs]
         new._degree = self._finddegree(new._coefs)
    new._mod()       #the polynomial mod x^r-1
    return new

一个给出错误结果的例子:

>>> N = 5915587277
>>> r = 1091
>>> pol = polynomials.ModPolynomial(r,N)
>>> pol += 1
>>> r_minus_1 = (pol**(r-1)).coefficients
>>> # (x+1)^r = (x+1)^(r-1) * (x+1) = (x+1)^(r-1) * x + (x+1)^(r-1) * 1
... shifted = [r_minus_1[-1]] + list(r_minus_1[:-1])   #(x+1)^(r-1) * x mod x^r-1
>>> res = [(a+b)%N for a,b in zip(shifted, r_minus_1)]
>>> tuple(res) == tuple((pol**r).coefficients)
False

求幂可能不是问题，因为我使用完全相同的算法和另一个给出正确结果的多项式实现。 “完全相同的算法”是指使用 numpy.convolve 的多项式类是另一个类的子类，仅重新实现 __mul__ 和 _mod ()[在_mod()中我简单地调用了另一个_mod()方法，然后删除大于多项式次数的项的0系数]。

另一个失败的例子:

>>> pol = polynomials.ModPolynomial(r, N)   #r,N same as before
>>> pol += 1
>>> (pol**96).coefficients
(1, 96, 4560, 142880, 3321960, 61124064, 927048304, 88017926, 2458096246, 1029657217, 4817106694, 4856395617, 384842111, 2941717277, 5186464497, 5873440931, 526082533, 39852453, 1160839201, 1963430115, 3122515485, 3694777161, 1571327669, 5827174319, 2249616721, 501768628, 5713942687, 1107927924, 3954439741, 1346794697, 4091850467, 2576431255, 94278252, 5838836826, 3146740571, 1898930862, 4578131646, 1668290724, 2073150016, 2424971880, 1386950302, 1658296694, 5652662386, 1467437114, 2496056685, 1276577534, 4774523858, 5138784090, 4607975862, 5138784090, 4774523858, 1276577534, 2496056685, 1467437114, 5652662386, 1658296694, 1386950302, 2424971880, 2073150016, 1668290724, 4578131646, 1898930862, 3146740571, 5838836826, 94278252, 2576431255, 4091850467, 1346794697, 3954439741, 1107927924, 5713942687, 501768628, 2249616721, 5827174319, 1571327669, 3694777161, 3122515485, 1963430115, 1160839201, 39852453, 526082533, 5873440931, 5186464497, 2941717277, 384842111, 4856395617, 4817106694, 1029657217, 2458096246, 88017926, 927048304, 61124064, 3321960, 142880, 4560, 96, 1)
#the correct result[taken from wolframalpha]:
(1, 96, 4560, 142880, 3321960, 61124064, 927048304, 88017926, 2458096246, 1029657217, 4817106694, 4856395617, 384842111, 2941717277, 5186464497, 5873440931, 526082533, 39852453, 1160839201L, 1963430115L, 3122515485L, 3694777161L, 1571327669L, 5827174319L, 1209974072L, 5377713256L, 4674300038L, 68285275L, 2914797092L, 307152048L, 3052207818L, 1536788606L, 4970222880L, 4799194177L, 2107097922L, 859288213L, 4578131646L, 1668290724L, 1033507367L, 1385329231L, 347307653L, 618654045L, 4613019737L, 427794465L, 1456414036L, 236934885L, 3734881209L, 4099141441L, 3568333213L, 4099141441L, 3734881209L, 236934885L, 1456414036L, 427794465L, 4613019737L, 618654045L, 347307653L, 1385329231L, 1033507367L, 1668290724L, 4578131646L, 859288213L, 2107097922L, 4799194177L, 4970222880L, 1536788606L, 3052207818L, 307152048L, 2914797092L, 68285275L, 4674300038L, 5377713256L, 1209974072L, 5827174319L, 1571327669L, 3694777161L, 3122515485L, 1963430115L, 1160839201L, 39852453, 526082533, 5873440931, 5186464497, 2941717277, 384842111, 4856395617, 4817106694, 1029657217, 2458096246, 88017926, 927048304, 61124064, 3321960, 142880, 4560, 96, 1)

错误的结果只出现在“big numbers”的情况下，指数也应该是“big”[我找不到指数<96的例子]。

我不知道为什么会这样。我正在使用 numpy.convolve 因为它在我的另一个问题中被建议，我只是想比较我的方法与 numpy 方法的速度。也许我以错误的方式使用了 numpy.convolve？

稍后我将尝试进行更多调试，尝试了解出错的确切位置和时间。

Answer 1

NumPy 是一个实现对 fixed-size numeric data types 数组快速运算的库.它不实现任意精度算术。所以你在这里看到的是整数溢出:NumPy 将你的系数表示为 64 位整数，当这些系数足够大时，它们会在 numpy.convolve 中溢出。

>>> import numpy
>>> a = numpy.convolve([1,1], [1,1])
>>> type(a[0])
<type 'numpy.int64'>

如果您需要对任意精度整数进行算术运算，则需要实现自己的卷积，以便它可以使用 Python 的任意精度整数。例如，

def convolve(a, b):
    """
    Generate the discrete, linear convolution of two one-dimensional sequences.
    """
    return [sum(a[j] * b[i - j] for j in range(i + 1)
                if j < len(a) and i - j < len(b))
            for i in range(len(a) + len(b) - 1)]

>>> a = [1,1]
>>> for i in range(95): a = convolve(a, [1,1])
... 
>>> from math import factorial as f
>>> all(a[i] == f(96) / f(i) / f(96 - i) for i in range(97))
True