Bakuriu Bakuriu - 20 days ago 10
Python Question

Multiplying polynomials with numpy.convolve return wrong result

I'm trying to multiply two polynomials using the

numpy.convolve
function.
I thought that this would be really easy, but I found out that it does not always return the correct product.

The implementation of multiplication is quite simple:

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


An example that gives a wrong result:

>>> 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


The exponentiation is probably not a problem, since I use exactly the same algorithm with an other polynomial implementation that gives the correct result. And by "exactly the same algorithm" I mean that the polynomial class that uses
numpy.convolve
is a subclass of the other class, reimplementing only
__mul__
, and
_mod()
[in
_mod()
I simply call the other
_mod()
method and then delete the 0 coefficients for the terms greater than the polynomial degree].

An other example in which this fail:

>>> 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)


The wrong results only appeard with "big numbers", and also the exponent should be "big"[I couldn't find an example with exponent < 96].

I have no clue of why this is happening. I'm using
numpy.convolve
because it was suggested in another my question, and I'm just trying to compare the speed of my approach to the numpy approach. Maybe I'm using
numpy.convolve
in the wrong way?

Later I'll try to do some more debugging, trying to understand where and when exactly things get wrong.

Answer

NumPy is a library that implements fast operation on arrays of fixed-size numeric data types. It does not implement arbitrary precision arithmetic. So what you are seeing here is integer overflow: NumPy is representing your coefficients as 64-bit integers and when these are large enough, they overflow in numpy.convolve.

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

If you need arithmetic on arbitrary-precision integers, you need to implement your own convolution so that it can use Python's arbitrary-precision integers. For example,

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