解决多项式乘法和除法“溢出”问题问题

Question

我有一个 1 次多项式的系数列表，其中 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.        ]])

并遇到以下操作的问题，

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

其中

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

和

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

首先，np.polydiv(reduce(np.polymul, a), a[0]) 的余数远大于 0，准确地说是 827.61514239，其次，商的最后两项应该为 0，但比 0 更大。1.37845155e-01, -1.06809521e+01。

我想知道我有哪些选择可以提高准确性？

Answer 1

有一个稍微复杂的方法，先保留产品，然后划分结构。

首先使用n个点并对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])

然后对 ys 而不是系数进行除法。

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

最后使用 xs 和 ys 使用多项式插值恢复系数

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

找到了第二个解决方案

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