Wrong polynomial value after multiplication by using numpy

Question

In the following script, I multiply many polynomials together and try to get its correct value by setting x=1855.

The answer has to be zero because one of the equation in the multiplication is x-1855 but the answer is way much more than zero after multiply the seventh times.The result shows the polynomial value and the polynomial.

I'm thinking that whether it is overflow somewhere but I am not sure. I need help.Thank you~

from matplotlib import pyplot as plt
import csv
import math
import numpy as np


x=np.array([1860, 1855, 1844, 1828, 1550, 1507, 1496, 1486, 1485, 1480, 1475, 1442,
1032,  950,  593,  543,  499,  243,  211,  200,  189,  173,  168,  152, 141,  140, 
124,   97,  87,   76,   70,   65,   59,   49,   43,   38,27,   22,   11,],'d')


p1=1
for i in range(1,39):
         p2 = np.poly1d([1,-x[i]])
         p1 = np.polymul(p1,p2)
         print(p2)
         print(p1(1855))

polynomials(first 8)

1 x - 1855, 1 x - 1844 , 1 x - 1828 , 1 x - 1550 ,1 x - 1507 ,1 x - 1496, 1 x - 1486,1 x - 1485 

multiplication

1 x - 1855                                                           -->first
   2
1 x - 3699 x + 3.421e+06                                             -->second
   3        2
1 x - 5527 x + 1.018e+07 x - 6.253e+09                               -->third
   4        3             2
1 x - 7077 x + 1.875e+07 x - 2.204e+10 x + 9.692e+12                 -->forth
   5        4             3             2
1 x - 8584 x + 2.941e+07 x - 5.029e+10 x + 4.29e+13 x - 1.461e+16    -->fifth
   6             5             4             3             2
1 x - 1.008e+04 x + 4.226e+07 x - 9.429e+10 x + 1.181e+14 x - 7.878e+16 x + 2.185e+19 
-->sixth
   7             6             5             4             3
1 x - 1.157e+04 x + 5.723e+07 x - 1.571e+11 x + 2.583e+14 x
              2
 - 2.543e+17 x + 1.389e+20 x - 3.247e+22                             --seventh
   8             7             6             5             4
1 x - 1.305e+04 x + 7.441e+07 x - 2.421e+11 x + 4.915e+14 x
              3             2
 - 6.378e+17 x + 5.166e+20 x - 2.388e+23 x + 4.822e+25               -->eighth
 ....

Result

when x=1855
first one=first two multiplication=...=first six multiplication=0
first seven multiplication=37748736.0
first eight multiplication=558345748480.0

Answer 1

It is not overflow, it is insufficient fp resolution.

Let's have a look at seven factors where you first see the wrong result:

>>> import numpy as np
>>> import functools as ft
>>> 
>>> x=np.array([1860, 1855, 1844, 1828, 1550, 1507, 1496, 1486, 1485, 1480, 1475, 1442,
... 1032,  950,  593,  543,  499,  243,  211,  200,  189,  173,  168,  152, 141,  140, 
... 124,   97,  87,   76,   70,   65,   59,   49,   43,   38,27,   22,   11,],'d')
>>> 
>>> L = np.c_[np.ones_like(x), -x]
>>> first7 = ft.reduce(np.polymul, L[:7])
>>> first7
array([ 1.00000000e+00, -1.19400000e+04,  6.10047450e+07, -1.72890531e+11,
        2.93522255e+14, -2.98513911e+17,  1.68387944e+20, -4.06415732e+22])

Let us focus on the last coefficient. We can use np.nextafter to obtain the next representable number in either direction.

>>> first7[-1] - np.nextafter(first7[-1], 0)
-8388608.0    
>>> first7[-1] - np.nextafter(first7[-1], 2 * first7[-1])
8388608.0

As we can see the resolution is much worse than 1, meaning that the coefficients are in practical terms guaranteed to be inaccurate.

If you are only interested in integer coefficients and arguments, you can use Python integers to avoid this problem (at the expense of performance):

>>> Li = L.astype(int).astype(object)
>>> fullprod = ft.reduce(np.polymul, Li)
>>> np.polyval(fullprod, np.array(1855, dtype=object))
0

If you need floating point, then consider using something like bigfloat (sympy might also work).

>>> from bigfloat import BigFloat, setcontext, Context
# set precision (1000 is probably way too much)
>>> setcontext(Context(precision=1000))
# convert array elements
>>> Lb = np.frompyfunc(BigFloat, 1, 1)(L)
# compute product
>>> fullprod = ft.reduce(np.polymul, Lb)
# evaluate
>>> np.polyval(fullprod, BigFloat(1855))
BigFloat.exact('0', precision=1000)