使用 numpy 乘法后多项式值错误

Question

在下面的脚本中，我将许多多项式相乘并尝试通过设置 x=1855 来获得正确的值。

答案必须为零，因为乘法中的一个方程是 x-1855 但在乘以第七个 times.The 结果显示多项式值和多项式后答案远大于零。

我在想它是否在某处溢出但我不确定。我需要help.Thank你~

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

多项式（前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

乘法

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

结果

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

Answer 1

不是溢出，是fp分辨率不够

让我们来看看您首先看到错误结果的七个因素：

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

让我们关注最后一个系数。我们可以使用 np.nextafter 在任一方向上获得下一个可表示的数字。

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

正如我们所见，分辨率比 1 差很多，这意味着系数实际上是不准确的。

如果你只对整数系数和参数感兴趣，可以使用Python整数来避免这个问题（以牺牲性能为代价）：

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

如果您需要浮点数，请考虑使用类似 bigfloat 的东西（sympy 也可以）。

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

使用 numpy 乘法后多项式值错误

Wrong polynomial value after multiplication by using numpy

python

numpy

multiplication

polynomials