Тип номера с точной информацией?

Question

Недавно кто-то хотел, чтобы функция ha sh без коллизий имела sh миллион значений до 32-битного значения ha sh. Если вы знаете парадокс дня рождения , вы знаете, что это вряд ли будет свободным от столкновений. Но, желая узнать вероятность, я вычислил ее следующим образом (начните с вероятности 1, затем для каждого из миллиона значений умножьте на вероятность, что это не одно из предыдущих):

>>> p = 1
>>> for i in range(10**6):
        p *= (2**32 - i) / 2**32

>>> p
2.7390147476139603e-51

Но я ' m умножая там миллион чисел с плавающей запятой, поэтому я беспокоюсь о том, что потеряю все больше и больше точности.

Есть ли числовой тип, который, в отличие от простых чисел с плавающей запятой, дает мне не только неточное число, но и говорит мне, насколько оно неточно является? Как диапазон вроде [2.73e-51, 2.74e-51] или с ошибкой 2.7390147476139603e-51 +/- 1e-54?

Или есть другой способ проверить точность результата?

Answer 1

Вот наихудший случай: при каждой операции (умножения или деления) явно умножайте результат на 1 + 2 ^ -52 или 1-2 ^ -52 и проверяйте (используя assert), что это действительно имело значение . Это должно оценить верхнюю границу неопределенности, и она все еще довольно мала - она ​​достигает конца без сбоев в утверждениях, а разница составляет одну часть из 10 ^ 9.

import sys

m_upper = (1 + 2**(1 - sys.float_info.mant_dig))
m_lower = (1 - 2**(1 - sys.float_info.mant_dig))

p_upper = p_lower = 1

for i in range(10**6):

    factor = (2**32 - i) / 2**32
    f_upper = factor * m_upper
    f_lower = factor * m_lower

    assert(f_upper > factor)
    assert(f_lower < factor)

    p_upper *= f_upper

    p_upper1 = p_upper * m_upper
    assert(p_upper1 > p_upper)
    p_upper = p_upper1
    
    p_lower *= f_lower

    p_lower1 = p_lower * m_lower
    assert(p_lower1 < p_lower)
    p_lower = p_lower1

print(p_upper, p_lower, p_upper - p_lower)

дает

2.739014748809663e-51 2.7390147464186476e-51 2.3910154124504752e-60

Обратите внимание, что если (1 - sys.float_info.mant_dig) заменить на -sys.float_info.mant_dig (т. Е. Использовать 2 ^ -53 вместо 2 ^ -52), утверждения начинают ошибаться.

Answer 2

(на основе ответа alaniwi )

Коэффициенты (2**32 - i) / 2**32 точны, т.е. они представлены как float в точности. Кроме того, стандарт с плавающей запятой гарантирует, что умножение даст наиболее точное значение float. Оно может быть ниже или выше реального продукта, но это ближайшее возможное значение float. Таким образом, если мы намеренно всегда отклоняемся к следующему большему значению float, оно никогда не будет меньше реального значения, то есть дает нам верхнюю границу. И мы получаем нижний связанный путем отклонения к следующему меньшему float значению.

Python 3.9 вводит math.nextafter , давайте использовать это:

>>> import math
>>> lower = upper = 1
>>> for i in range(10**6):
        factor = (2**32 - i) / 2**32
        lower = math.nextafter(lower * factor, -math.inf)
        upper = math.nextafter(upper * factor, math.inf)

>>> lower, upper
(2.739014747179961e-51, 2.739014748048138e-51)
>>> upper - lower
8.681767916298978e-61

Answer 3

Как прокомментировал Эри c Постпишил , это « интервальная арифметика c и связанные понятия».

Гугл python интервал arithmeti c находит PyInterval . Давайте попробуем это:

from interval import interval

p = interval[1]
for i in range(10**6):
    p *= (2**32 - i) / 2**32
print(p)

Вывод (запустить на repl.it ):

interval([2.7390147473969355e-51, 2.739014747831127e-51])

Давайте сравним это с границами из целочисленного вычисления :

interval upper 2.739014747831127e-51
integer upper   27390147476140722271150280539996691121583143640960
integer lower   27390147476140722271150280539996691121583143636646
interval lower 2.7390147473969355e-51

Таким образом, решение interval менее точное (это больший интервал, совпадают только первые десять цифр нижней и верхней границы), но оно правильное (реальное значение действительно внутри интервала). Думаю, в этом смысле он всегда будет правильным, хотя я не разбирался в том, как это работает.

Answer 4

Один из способов получить диапазон - использовать целые числа, масштабируя вероятность, скажем, на 10 ¹⁰⁰ . Для нижней границы всегда округлять в меньшую сторону, а для верхней границы всегда в большую:

>>> lower = 10**100
>>> for i in range(10**6):
        lower = lower * (2**32 - i) // 2**32

>>> lower
27390147476140722271150280539996691121583143636646

>>> upper = 10**100
>>> for i in range(10**6):
        upper = -(-upper * (2**32 - i) // 2**32)

>>> upper
27390147476140722271150280539996691121583143640960

Выравнивание:

upper  27390147476140722271150280539996691121583143640960
p     2.7390147476139603e-51
lower  27390147476140722271150280539996691121583143636646

Мы видим, что p ( float) на самом деле выходит за пределы реального диапазона, слишком мало. Но его первые двенадцать цифр верны, так что это кажется довольно хорошим.

Сравнивая lower и upper, мы также получаем довольно много совпадающих и, следовательно, правильных цифр: 2,73901474761407222711502805399966911215831436e-51. А с большим коэффициентом масштабирования мы можем получить еще больше.