Неверное значение полинома после умножения с использованием numpy

Question

В следующем сценарии я умножаю много полиномов вместе и пытаюсь получить его правильное значение, установив x = 1855.

Ответ должен быть равен нулю, потому что одно из уравнений в умножении - x-1855но ответ намного больше нуля после умножения в седьмой раз. Результат показывает полиномиальное значение и полином.

Я думаю, что это где-то переполнение, но я не уверен.Мне нужна помощь. Спасибо ~

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)