Bakuriu - 1 year ago 200
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.

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
``````
Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download