numpyfloating-accuracypolynomialspolynomial-math

# Addressing polynomial multiplication and division "overflow" issue

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?

Solution

• 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]
``````