2

関数を使用して 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()[in _mod()I では単に他のメソッドを呼び出して_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 に答える