Addressing polynomial multiplication and division "overflow" issue

Question

I have a list of the coefficient to degree 1 polynomials, with a[i][0]*x^1 + a[i][1]

a = np.array([[ 1.        , 77.48514702],
          [ 1.        ,  0.        ],
             [ 1.        ,  2.4239275 ],
           [ 1.        ,  1.21848739],
            [ 1.        ,  0.        ],
            [ 1.        ,  1.18181818],
           [ 1.        ,  1.375     ],
           [ 1.        ,  2.        ],
          [ 1.        ,  2.        ],
          [ 1.        ,  2.        ]])

And running into issues with the following operation,

np.polydiv(reduce(np.polymul, a), a[0])[0] != reduce(np.polymul, a[1:])

where

In [185]: reduce(np.polymul, a[1:])
Out[185]:
array([  1.        ,  12.19923307,  63.08691612, 179.21045388,
       301.91486027, 301.5756213 , 165.35814595,  38.39582615,
         0.        ,   0.        ])

and

In [186]: np.polydiv(reduce(np.polymul, a), a[0])[0]
Out[186]:
array([ 1.00000000e+00,  1.21992331e+01,  6.30869161e+01,  1.79210454e+02,
        3.01914860e+02,  3.01575621e+02,  1.65358169e+02,  3.83940472e+01,
        1.37845155e-01, -1.06809521e+01])

First of all the remainder of np.polydiv(reduce(np.polymul, a), a[0]) is way bigger than 0, 827.61514239 to be exact, and secondly, the last two terms to quotient should be 0, but way larger from 0. 1.37845155e-01, -1.06809521e+01.

I'm wondering what are my options to improve the accuracy?

Answer 1

There is a slightly complicated way to keep the product first and then divide structure.

By first employ n points and evaluate on a.

xs = np.linspace(0, 1., 10) 
ys = np.array([np.prod(list(map(lambda r: np.polyval(r, x), a))) for x in xs])

then do the division on ys instead of coefficients.

ys = ys/np.array([np.polyval(a[0], x) for x in xs])

finally recover the coefficient using polynomial interpolation with xs and ys

from scipy.interpolate import lagrange
lagrange(xs, ys)

Found a second solution

b = a[:,::-1]
np.polydiv(reduce(np.polymul, b), b[0])[0][::-1]