Ликбез о плавающей точке: сложение, катастрофическое сокращение и бабушка Кэхена / forpes.ru

Главная
Ликбез о плавающей точке: сложение, катастрофическое сокращение и бабушка Кэхена

Ликбез о плавающей точке: сложение, катастрофическое сокращение и бабушка Кэхена +13

18.11.2025 15:08

haqreu 0 650 Источник

Давайте продолжим обсуждение самой неоптимизированной в мире 32-битной библиотеки для работы с плавающей запятой TinyFloat. Библиотека написана на C++ и намеренно избегает встроенных типов плавающей запятой, полагаясь исключительно на 32-битные целые числа. Цель состоит в том, чтобы сделать код максимально читабельным — без бит-хаков и хитроумных уловок.

Кроме того, я хочу иметь подробную документацию о том, что происходит «под капотом». Оказалось, что лучший способ документировать код C++ — это полностью переписать его на Python :-)

Тема на сегодня - как работает сложение чисел с плавающей точкой и соответствующие подводные камни. Это третья статья из цикла, первые две читать тут и тут.

В этой статье я использую 8-битную плавающую запятую и игнорирую специальные значения (NaN, Inf). Я полагаюсь на нативный тип float в Python только для тестирования реализации.

Float8

С этого момента мы будем работать с числами с плавающей запятой со знаком, поэтому я добавляю последний бит, необходимый для получения класса Float8. Если вы читали две предыдущие главы, то знаете, что этот класс имеет целочисленные члены e и m (обозначающие экспоненту и мантиссу). Здесь я добавляю член s, который может принимать значение 1 или -1.

Эта тройка представляет вещественное число, которое можно записать в научной нотации $s \cdot m \cdot 2^{e-4}$ , проверьте строку 17 следующего листинга:

class Float8:
    def __init__(self, s, e, m):
        self.s, self.e, self.m = s, e, m

    @classmethod
    def from_uint8(cls, uint8):
        s = 1 if uint8<128 else -1
        e = (uint8 % 128) // 16 - 4
        m = (uint8 % 128) %  16
        if e == -4:
            e += 1  # if subnormal, increment the exponent, the hidden bit = 0
        else:
            m += 16 # if normal, recover the hidden bit = 1
        return cls(s, e, m)

    def __float__(self):
        return self.s * self.m * 2**(self.e-4)

    @classmethod
    def from_float(cls, f):
        dist = [ abs(float(Float8.from_uint8(i)) - abs(f)) for i in range(128) ] # distance from |f| to all non-negative Float8
        i = min(range(128), key=lambda j : dist[j])                              # take the (first) closest
        if i < 127 and dist[i] == dist[i+1]:                                     # if tie
            i += i % 2                                                           # round to even
        return Float8.from_uint8(i + (0 if f>=0 else 128))                       # do not forget to recover the sign

    def copy(self):
        return Float8(**self.__dict__)

Обратите внимание, что Float8 хранит нормализованные экспоненту и мантиссу: скрытый бит восстанавливается, а денормализованные числа преобразуются в нормализованную форму (строки 10-13).

Таким образом,

$s \in \{-1, 1\}$ ,
$e \in [-3\dots 3]$ ,
$m\in[0\dots 31]$ ,

и тройка кодирует число $s\cdot m \cdot 2^{e-4}$ без какой-либо остаточной логики обработки особых случаев.

Округление

Стандарт IEEE754 определяет ряд режимов округления, но режимом округления по умолчанию является режим округления-до-ближайшего- числа,-до-чётного-числа-при-равенстве (RNE). Ниже приведена иллюстрация режима округления RNE с использованием нашего 8-разрядного представления с плавающей запятой:

Вещественное число 4,65 ближе всего к значению Float8 4,75,
но 4,875 лежит ровно посередине между 4,75 и 5,0; в этом случае RNE требует округления до четного из двух представляемых значений.
Поскольку битовая маска 5,0 (0b01100100) является чётной, 4,875 округляется до 5,0. Аналогично, вещественное число 5,125 округляется до 5,0, а 5,375 — до 5,5.

В приведенном выше фрагменте функция from_float() принимает python float и округляет его до нашего Float8, используя режим RNE. Реализация проста: вычислить расстояние от вещественного числа до всех значений Float8 (строка 21), определить ближайшего(-их) кандидата(ов) и применить RNE в случае ничьей.

Наивное сложение

Давайте рассмотрим общую идею сложения чисел с плавающей запятой. Я хочу сложить два положительных числа и , представленных их целыми мантиссами и экспонентами :

$\begin{align}a &= a_m \cdot 2^{a_e - 4}\\b &= b_m \cdot 2^{b_e - 4}\end{align}$

Вот наивная реализация (строки 30-52):

class Float8:
    def __init__(self, s, e, m):
        self.s, self.e, self.m = s, e, m

    @classmethod
    def from_uint8(cls, uint8):
        s = 1 if uint8<128 else -1
        e = (uint8 % 128) // 16 - 4
        m = (uint8 % 128) %  16
        if e == -4:
            e += 1  # if subnormal, increment the exponent, the hidden bit = 0
        else:
            m += 16 # if normal, recover the hidden bit = 1
        return cls(s, e, m)

    def __float__(self):
        return self.s * self.m * 2**(self.e-4)

    @classmethod
    def from_float(cls, f):
        dist = [ abs(float(Float8.from_uint8(i)) - abs(f)) for i in range(128) ] # distance from |f| to all non-negative Float8
        i = min(range(128), key=lambda j : dist[j])                              # take the (first) closest
        if i < 127 and dist[i] == dist[i+1]:                                     # if tie
            i += i % 2                                                           # round to even
        return Float8.from_uint8(i + (0 if f>=0 else 128))                       # do not forget to recover the sign

    def copy(self):
        return Float8(**self.__dict__)

    def __add__(self, other):
        a, b = self.copy(), other.copy()
        if a.e < b.e:
            a, b = b, a 

        while a.e > b.e:                         # align exponents
            b.m //= 2
            b.e += 1

        sum = Float8(a.s if a.m >= b.m else b.s, # sign
                     a.e,                        # exponent
                     abs(a.s*a.m + b.s*b.m))     # mantissa
    
        while sum.m < 16 and sum.e > -3:         # normalize the result
            sum.m *= 2
            sum.e -= 1 

        while sum.m >= 32 and sum.e < 3:         # can't be more than one iteration
            sum.m //= 2
            sum.e += 1

        sum.m = min(sum.m, 31)                   # handle overflows (we have no infinities here)
        return sum

    def __sub__(self, other):
        return self + Float8(-other.s, other.e, other.m)

Под «наивным» я имею в виду, что на некоторое время я игнорирую правильное округление результата. Сложить два числа в научной нотации не так уж и сложно, давайте разделим это на два случая:

Равные показатели степени:

Если показатели степени равны, то мы можем сложить мантиссы напрямую:

$a + b = a_m \cdot 2^{a_e - 4} + b_m \cdot 2^{b_e - 4} = (a_m + b_m) \cdot 2^{a_e - 4}.$

Рассмотрим пример с и . В этом случае , и :

$\begin{align}a &= 4,25 = 17 \cdot 2^{2-4}\\b &= 5,25 = 21 \cdot 2^{2-4}\end{align}$

Тогда сумма может быть представлена экспонентой и мантиссой :

$a+b = (17+21) \cdot 2^{2-4} = 9,5.$

Хотя это действительно просто, нам нужно быть внимательными. Действительно, сумма может быть представлена показателем и мантиссой , но есть одна загвоздка.

Чтобы сохранить результат в Float8, мантисса должна укладываться в диапазон $[0\dots 31]$ , а 38 превышает 31.

Это не проблема: мы можем разделить мантиссу на 2 и увеличить экспоненту. Этот процесс называется нормализацией, см. строки 43-49 предыдущего листинга.

$9,5 = 38\cdot 2^{2-4} = 19 \cdot 2^{3-4}.$

Здесь сумма представлена мантиссой и экспонентой . Давайте протестируем нашу реализацию:

from float8 import *

r = Float8.from_float(4.5) + \
    Float8.from_float(5.25)
t = Float8.from_float(4.5 + 5.25)
print(f'The sum is {r.m}*2**({r.e}-4) = {float(r)}, expected {float(t)}')

Обратите внимание, что здесь я оцениваю правильность результата, вычисляя сумму с помощью нативного float и преобразуя ее в Float8.
В этом случае мы получаем правильный результат:

The sum is 19*2**(3-4) = 9.5, expected 9.5

Спойлер

Не забывайте, что экспонента также должна входить в диапазон $[-3\dots 3]$. Всякий раз, когда она превышает максимальное значение, стандарт IEEE754 гласит, что результатом является специальное значение `Inf`. [TinyFloat](https://github.com/ssloy/tinyfloat) корректно обрабатывает бесконечности, но в этой реализации на Python я округляю до максимального представляемого числа с плавающей запятой, т. е. $15.5$ (строка 51 в листинге).

Есть еще одна сложность, которую нужно решить. Попробуем сложить и :

from float8 import *

r = Float8.from_float(4.5) + \
    Float8.from_float(5.25)
t = Float8.from_float(4.5 + 5.25)
print(f'The sum is {r.m}*2**({r.e}-4) = {float(r)}, expected {float(t)}')

Как видите, наше наивное сложение отличается от истинного значения:

The sum is 19*2**(3-4) = 9.5, expected 10.0

Давайте внимательно рассмотрим, что происходит. Прежде всего, обратите внимание, что истинный результат не может быть точно представлен в нашем Float8, поэтому нам нужно округлить сумму.

До результирующей нормализации мантисса суммы была , а экспонента была равна :

$\begin{align}a = 4,5 & = 18 \cdot 2^{2-4}\\b = 5,25 & = 21 \cdot 2^{2-4}\\a + b & = 39 \cdot 2^{2-4}\end{align}$

Нормализация делит мантиссу пополам и увеличивает показатель степени:

$a + b = 39 \cdot 2^{2-4} = \frac{39}{2} \cdot 2^{3-4}$

Поскольку мантисса должна быть целым числом, наша наивная реализация усекает $\frac{39}2 = 19,5$ до , фактически теряя один бит точности. Согласно правилам RNE, должно быть округлено до .

Я отложу обсуждение округления на некоторое время, но помните, что нормализация может привести к потере одного бита информации.

Различные показатели:

Если показатели и различаются, мы можем предположить, что

, в противном случае мы просто меняем местами числа, см. строки 32-33 в листинге. С учетом этого мы можем выровнять показатели степени:

$b = b_m \cdot 2^{b_e - 4} = \frac{b_m}{2^{a_e - b_e}} \cdot 2^ {b_e - 4 + (a_e - b_e)} = \frac{b_m}{2^{a_e - b_e}} \cdot 2^{a_e - 4}.$

После выравнивания мы можем вычислить сумму так же, как и раньше, просто сложив (целые) мантиссы!

$a + b = \left(a_m + \frac{b_m}{2^{a_e - b_e}}\right) \cdot 2^{a_e - 4}.$

Это выравнивание выполняется в строках 35-37 приведенного выше листинга. Обратите внимание, что эта процедура очень похожа на нормализацию (строки 47-49), поэтому она имеет точно такую же проблему с округлением. Если $\frac{b_m}{2^{a_e - b_e}}$ является целым числом, то выравнивание не приводит к потере данных. Однако, если оно не является целым числом, нам нужно дополнительно потрудиться, чтобы правильно округлить результат.

Вот вам небольшой вопрос: почему мы сдвигаем вправо, а не сдвигаем влево? Теоретически, оба варианта приведут к выравниванию экспонент.

Спойлер

Сдвиг влево большой мантиссы может легко привести к переполнению мантиссы перед сложением. Сдвиг мантиссы вправо никогда не увеличивает величину; он только безопасно снижает точность. Следовательно, при сдвиге $b_m$ вправо мы теряем крошечные, незначительные биты (которые мы вскоре восстановим), тогда как сдвиг влево привел бы к потере *важных* старших битов.

GRS биты: правильное округление

В нашей простой реализации есть два момента, которые могут привести к неправильному окончательному округлению суммы:

деление мантиссы пополам при выравнивании экспонент
и деление мантиссы пополам при нормализации.

При сложении чисел с плавающей запятой во время вычислений используются дополнительные биты (справа от мантиссы). Когда мантисса сдвигается вправо, младшие её биты свдигатся в них. Оказывается, для правильного округления достаточно трех битов, но мы вернемся к этому позже.

Давайте вычислим вручную , уделяя пристальное внимание проблеме округления:

$\begin{alignat}{2}a & = 7 & = 28 \cdot 2^{2-4} \\b & =0,875 & = 14 \cdot 2^{0-4}\end{alignat}$

Экспоненты различаются,

, поэтому нам нужно их выровнять:

$b = 14 \cdot 2^{0-4} = 3,5 \cdot 2^{2-4}.$

В наивной реализации мы усекли 3,5 до 3, тем самым введя ошибку. На этот раз мы сохраняем дробную часть мантиссы. Выровняв показатели степени, мы можем просто сложить мантиссы:

$a + b = (28 + 3,5) \cdot 2^{2-4} = 31,5 \cdot 2^{2-4}$

Обратите внимание, что здесь мы не попадаем в замкнутый круг: для сложения 28 + 3,5 не требуется математика с плавающей запятой, оно выполняется с целыми числами (т. е. в представлении с фиксированной запятой, поскольку для округления достаточно трех дополнительных битов).

Теперь нам нужно выполнить нормализацию, поскольку превышает :

$a + b = 15,75 \cdot 2^{3-4}.$

Теперь нам нужно правильно округлить мантиссу, следуя правилу RNE. Решение принимается просто:

если дробная часть > 0,5 ИЛИ
   дробная часть = 0,5 И целая часть нечетная:
        округлить мантиссу в большую сторону
в противном случае:
        усечь мантиссу в меньшую сторону

Поскольку

, мы можем с уверенностью заключить, что

$a + b = 15,75 \cdot 2^{3-4} \approx 16 \cdot 2^{3-4}.$

Поздравляем, ближайшим к истинной сумме значением Float8 является !

Вот правильная реализация сложения:

    def __add__(self, other):
        a, b = self.copy(), other.copy()
        if a.e < b.e:
            a, b = b, a

        a.m *= 8                                 # reserve place for GRS bits
        b.m *= 8
        while a.e > b.e:                         # align exponents
            b.m = (b.m // 2) | (b.m % 2)         # LSB is sticky
            b.e += 1

        sum = Float8(a.s if a.m >= b.m else b.s, # sign
                     a.e,                        # exponent
                     abs(a.s*a.m + b.s*b.m))     # mantissa

        while sum.m < 16*8 and sum.e > -3:       # normalize the result
            sum.m *= 2
            sum.e -= 1

        while sum.m >= 32*8 and sum.e < 3:       # can't be more than one iteration
            sum.m = (sum.m // 2) | (sum.m % 2)   # do not forget the sticky bit
            sum.e += 1

        g = (sum.m // 4) % 2                     # guard bit
        r = (sum.m // 2) % 2                     # round bit
        s =  sum.m % 2                           # sticky bit
        sum.m //= 8

        if g and (r or s or sum.m % 2):          # round-to-nearest, even-on-ties
            sum.m += 1
            if sum.m == 32 and sum.e < 3:        # renormalize if necessary
                sum.m //= 2
                sum.e += 1

        sum.m = min(sum.m, 31)                   # handle overflows (we have no infinities here)
        return sum

Я зарезервировал три дополнительных бита точности в строках 35-36. В строках 53-56 я восстанавливаю целую часть и дробную часть мантиссы. Наконец, правило RNE применяется в строках 58-62.

Давайте попробуем код на приведенном выше примере:

from float8 import *

r = Float8.from_float(7) + \
    Float8.from_float(0.875)
t = Float8.from_float(7 + 0.875)
print(f'The sum is {r.m}*2**({r.e}-4) = {float(r)}, expected {float(t)}')

Результат действительно равен 8,0, как и ожидалось:

The sum is 16*2**(3-4) = 8.0, expected 8.0

Дополнительные три бита, которые я выделил, называются битами GRS (guard, round и sticky). Трех битов достаточно, чтобы гарантировать правильное округление RNE для всех размеров плавающих точек, включая 32-битные.

Домашнее задание № 1

Найдите пример, в котором нам понадобятся все три бита для определения округления.
Обратите внимание, что в приведенном выше примере мы не использовали LSB (sticky).

Спойлер

Во фремя фазы нормализации вам нужно сдвинуть мантиссу *влево*, а не вправо.

Домашнее задание № 2

Покажите, что трех битов достаточно для округления при одном условии: sticky-бит действительно липкий.

Когда мантисса сдвигается с конца, она сдвигается в GRS-биты.
В целом это работает как обычный сдвиг вправо, за исключением того, что в момент, когда любой взведённый бит сдвигается в "липкий" бит, он остается 1 с этого момента и до конца (это и делает его липким).

Проверьте строку 38 (а также строку 50) листинга:

            b.m = (b.m // 2) | (b.m % 2)         # LSB является липким

Мантисса сдвигается вправо, но как только младший бит взведён, он остается взведённым навсегда.

Спойлер

Выравнивание экспонент может потребовать более 3 сдвигов, поэтому мы не можем сохранить все биты, выпадающие из мантиссы, в битах GRS. Тем не менее, "липкое" поведение достаточно для окончательного раунда округления.

Численные ошибки

Как мы только что видели, результат арифметических операций с плавающей запятой подвержен ошибке округления, когда он не может быть точно представлен. Хотя современное оборудование и библиотеки дают правильно округленные результаты арифметических операций, эта ошибка округления может усиливаться при выполнении ряда операций, поскольку промежуточный результат должен быть округлен. При работе с конечной точностью следует тщательно учитывать числовые ошибки.

Катастрофическое сокращение

Представьте, что вы измеряете, сколько бензина вы используете в двух одинаковых поездках.

В первой поездке автомобиль использует 30,2 литра (ваша бензоколонка показывает с точностью до одного десятичного места, ±0,1 л).
Во время второй поездки, после небольшой настройки, вы используете 30,0 литров (опять же, точность $\pm 0,1\ L$ ).

Вы хотите выяснить, насколько меньше бензина вы использовали после тюнинга:

$\text{Экономия топлива} = 30,2\ L - 30,0\ L = 0,2\ L$

Но из-за ограниченной точности заправочной колонки каждое измерение может иметь погрешность $\pm 0,1\ L$ , поэтому фактическое количество использованного топлива могло быть в пределах:

Первая поездка: $30,1 - 30,3\ L$
Вторая поездка: $29,9 - 30,1\ L$

Таким образом, разница (экономия) может составлять от 0,0 л (30,1 – 30,1 долл.) до 0,4 л (30,3 – 29,9 долл.)! Ваша относительная погрешность огромна по сравнению с измеренной вами «экономией». Несмотря на то, что ваши исходные измерения казались точными до 0,1 л, вычитание почти равных чисел (для получения экономии) сделало ваш результат в основном просто шумом.

Это называется катастрофическим сокращением: когда вы вычитаете почти равные значения, любая погрешность или неопределенность в каждом из них увеличивается в вашем небольшом результате, независимо от использования компьютеров или физических инструментов.

Проблема заключается НЕ в плавающей запятой, а в том, что вы потеряли точность, вычитая два почти равных измерения. Катастрофическое сокращение может произойти в любой системе с ограниченной точностью (физические измерения, финансовые вычисления с ограниченным количеством десятичных знаков и т. д.).

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

Рассмотрим и .
Легко вычислить истинную разность квадратов .

Обратите внимание, что как , так и точно представляются в Float8. Однако, если мы вычислим в нашем представлении Float8, то получим в результате . Наиболее близким к значением Float8 является , а округляется до .

$\begin{alignat}{2}2.875^2 &= 8.265625 & \approx 8.5 \\2.75^2 &= 7.5625 & \approx 7.5\end{alignat}$

Вычитая одно из другого, мы получаем . Относительная погрешность нашего результата составляет $42\%$ !

$\frac{|1 - 0.703125|}{|0.703125|} \approx 0.42.$

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

В нашем примере разница может быть точно вычислена в нашем Float8, а сумма округляется до :

$\begin{alignat}{2}2,875 + 2,75 &= 5,625 & \approx 5,5 \\2,875 - 2,75 &= 0,125. &\end{alignat}$

Тогда произведение $5,5 \cdot 0,125 = 0,6875$ , которое снова точно вычисляется в нашем представлении Float8, имеет относительную погрешность всего $2\%$ :

$\frac{|0,6875 - 0,703125|}{|0,703125|} \approx 0,02.$

Ошибка в $2\%$ гораздо лучше, чем $42\%$ , которые были раньше!

Суммирование

Ладно, вычитание — дело сложное, но если мы складываем положительные числа, ничего не может пойти не так. Скажем, если я сложу раз $\frac{1}{128}$ , то получу , верно? Ну верно же?

Обратите внимание, что можно точно представить в Float8, поэтому когда я пишу Float8.from_float(1/128), нет потери информации.
Давайте проверим суммирование. Я начну с переменной-аккумулятора total, инициализированной как , и добавлю раз .

from float8 import *

total = Float8.from_float(0)
for _ in range(128):
    total += Float8.from_float(1/128)
print(float(total))

Барабанная дробь...

0.25

Важно понять, почему это происходит, поэтому давайте вручную сложим 0,25 и 0,0078125. Прежде всего, найдем их мантиссы и экспоненты:

$\begin{align}0,25 &= 16 \cdot 2^{-2-4}\\0,0078125 &= 1\cdot 2^{-3-4}\end{align}$

Нам нужно выровнять мантиссы:

$0,0078125 = 0,5 \cdot 2^{-2-4}$

Наконец, мы суммируем мантиссы и округляем их, используя правило RNE:

$0,25 + 0,0078125 = (16 + 0,5) \cdot 2^{-2-4} \approx 16 \cdot 2^{-2-4}.$

Другими словами, в представлении Float8, 0,25 + 0,0078125 = 0,25.

При сложении двух чисел с очень разными величинами меньшее число может быть проигнорировано или «потеряно» из-за ограниченной точности, что приводит к потере точности. Это происходит потому, что значимые цифры меньшего числа могут не поместиться в пределы допустимой точности при сложении с гораздо большим числом.

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

Проиллюстрируем эту идею рабочим кодом. Здесь я определяю две функции, возвращающие сумму массива. Первая — это простая последовательная реализация, а вторая реализует идею бинарного дерева:

def naive_sum(arr):
    s = Float8.from_float(0.)
    for v in arr:
        s += v
    return s

def pairwise_sum(arr):
    match len(arr):
        case 0:
            return Float8.from_float(0)
        case 1:
            return arr[0]
    mid = len(arr)//2
    return pairwise_sum(arr[:mid]) + \
           pairwise_sum(arr[mid:])

А теперь давайте вызовем обе функции на нашем тестовом массиве:

arr = [ Float8.from_float(1/128) ] * 128

print('naive sum:\t',    float(   naive_sum(arr)))
print('pairwise sum:\t', float(pairwise_sum(arr)))

Попарное сложение сильно улучшает результат:

naive sum:       0.25
pairwise sum:    1.0

Хотя попарное суммирование уменьшает ошибки округления при сложении больших последовательностей чисел, существуют еще более совершенные методы, такие как суммирование Кэхана-Бабушки, которые идут еще дальше. Суммирование Кэхена-Бабушки отслеживает небольшие ошибки, которые обычно теряются при стандартном сложении. Благодаря поддержанию компенсационной переменной для потерянных младших битов, оно гарантирует, что даже незначительные вклады от крошечных чисел не будут «забыты». Это делает алгоритм Кэхена-Бабушки особенно ценным в ситуациях, требующих высокой точности при длительных последовательностях сложений.

Кроме того, попарное суммирование обычно требует заранее знать все значения, что делает его идеальным для офлайн суммирования, когда доступен полный список. Однако во многих практических сценариях нам нужно суммировать значения по мере их поступления, одно за другим, это называется онлайн суммированием. В таких случаях суммирование Кэхена-Бабушки особенно эффективно.

Вот его реализация:

def kahan_sum(arr):
    s, c = Float8.from_float(0), Float8.from_float(0)
    for v in arr:
        t = s + (v - c)
        c = (t - s) - (v - c)
        s = t
    return s

Переменная s является аккумулятором, а c — переменной компенсации. Полезно заметить, что мы можем приблизительно оценить погрешность s, поскольку она меняется от итерации к итерации. Для этого рассмотрим выражение

Алгебраически это выражение равно нулю. Однако численно это может быть не так. В частности, сумма может быть округлена с точностью до плавающей запятой. Вычитая поочерёдно и , мы получаем приблизительную оценку погрешности приближения . Удаление и из интуитивно переходит от больших порядков величины к меньшим, а не наоборот, и, следовательно, с меньшей вероятностью вызывает ошибку округления, чем вычисление суммы ; это наблюдение объясняет, почему оценка ошибки сама по себе не так подвержена проблемам округления, как исходная операция.

Давайте сравним все реализации на данных, более близких к реальным.
Здесь я хочу сложить 128 случайных чисел из диапазона :

import random
random.seed(1)
arr = [ Float8.from_float(random.uniform(-0.25, .25)) for _ in range(128) ]

print('naive sum:\t',    float(   naive_sum(arr)))
print('pairwise sum:\t', float(pairwise_sum(arr)))
print('Kahan sum:\t',    float(   kahan_sum(arr)))
print('True sum:\t',     float(Float8.from_float(sum([float(v) for v in arr]))))

Вот результаты исполнения:

naive sum:       0.1875
pairwise sum:    -0.03125
Kahan sum:       0.015625
True sum:        0.015625

Как и ожидалось, наивная сумма сильно отличается от истинного результата, попарное суммирование дает разумное приближение, а суммирование Кэхена-Бабушки является лучшим.

Подводим итог

Сложение чисел с плавающей запятой обманчиво кажется простым. Несмотря на то, что математическая операция проста, для ее правильной реализации требуется тщательное внимание к деталям.

Ограниченная точность не является изъяном проектирования, это фундаментальное ограничение представления бесконечного множества вещественных чисел при помощи конечного числа битов. Важно понимать, где возникают ошибки и как их уменьшить:

переписывайте выражения для улучшения численной стабильности,
избегайте вычитания почти равных чисел,
используйте попарное или суммирование Кэхена при суммировании длинных рядов чисел.

С помощью этих инструментов даже крошечный 8-битный формат с плавающей запятой может вести себя предсказуемо и прозрачно, помогая раскрыть структуру реального аппаратного обеспечения с плавающей запятой и численные алгоритмы, которые оно позволяет использовать.

Ликбез о плавающей точке: сложение, катастрофическое сокращение и бабушка Кэхена +13

Float8

Округление

Наивное сложение

Равные показатели степени:

Различные показатели:

GRS биты: правильное округление

Домашнее задание № 1

Домашнее задание № 2

Численные ошибки

Катастрофическое сокращение

Суммирование

Подводим итог

Комментарии (0)