2

我正在尝试使用该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)

错误的结果只出现在“大数字”中，而且指数应该是“大”[我找不到指数 < 96 的例子]。

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

稍后我会尝试做更多的调试，试图了解究竟在哪里以及什么时候出了问题。

4

2 回答 2

4

NumPy 是一个在固定大小的数值数据类型数组上实现快速操作的库。它不实现任意精度算术。所以你在这里看到的是整数溢出：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
于 2012-08-26T10:48:23.470 回答
0

如果在执行重复卷积时速度是一个问题dtype=object（如您的情况），您应该考虑使用karatsuba模块而不是numpy.convolve. 它旨在使用高级对象（例如 with dtype=object），但依赖于计划来预先计算卷积。它对于单个卷积是无用的，但如果需要许多类似的卷积则非常有效。

于 2016-11-20T17:33:32.543 回答