Суммирование Кэхэна

Миллион покупок по 149,90 тенге — 149 900 000 тенге ровно. Сложим их обычным циклом с числами float: выйдет 149 900 000,0028. Почти ничего, но откуда-то взялись лишние тиыны, и с каждым миллионом их становится больше. В float32, которым считают видеокарты и многие модели, та же сумма недосчитывается больше ста двадцати тысяч тенге. Ниже — сумма растёт прямо в вашем браузере тремя способами, и график показывает, насколько каждый отстал от правды.

Синяя линия, float32, уходит вверх быстрее всех: у него всего около семи значащих цифр. Оранжевая, float64, растёт медленнее, но тоже растёт. Зелёная — та же float64, но с поправкой Кэхэна: её ошибка стоит на месте, сколько бы слагаемых ни было. Переключите на гармонический ряд: в float32 он перестаёт расти на 2 097 152-м слагаемом и навсегда замирает на 15,403683, хотя в математике растёт до бесконечности. Почему 0,1 + 0,2 не равно 0,3 и как устроены числа с плавающей запятой, рассказано в главах курса о переменных и о битах; здесь — что бывает, когда таких сложений миллионы.

Идея

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

Вот весь алгоритм. Переменная c — поправка, сначала ноль.

y = x - c          # слагаемое вместе с тем, что потеряли в прошлый раз
t = s + y          # новая сумма: здесь младшие биты y теряются
c = (t - s) - y    # (t - s) — сколько на самом деле прибавилось; минус y — сколько потерялось
s = t

Третья строка выглядит как ноль: (s + y − s) − y. В точной арифметике так и есть, а в машинной t − s — это то, что реально прибавилось к сумме, и разница с y — ровно потерянные биты, причём сама эта разница вычисляется без ошибки. Уильям Кэхэн опубликовал этот приём в 1965 году. Для простого цикла ошибка суммы n слагаемых в худшем случае пропорциональна n; с поправкой Кэхэна — ограничена двумя единицами младшего разряда, умноженными на сумму модулей слагаемых, и от n почти не зависит.

Где Кэхэн бессилен

Поправка — одно число, и она хороша, когда к большой сумме много раз прибавляют маленькие числа. Но если маленькое число зажато между огромными разного знака — 10¹⁰⁰, потом 1, потом −10¹⁰⁰, — поправка сама теряется на фоне огромной суммы, и результат — ноль вместо единицы. В 1974 году Арнольд Ноймайер предложил вариант, который сравнивает модули суммы и слагаемого и запоминает младшие биты того, кто меньше; такой вариант с этим примером справляется, и именно он с версии 3.12 работает внутри встроенного sum() в Python.

Есть ещё два способа. Попарное суммирование складывает числа как дерево: половину с половиной, четверть с четвертью. Каждое число проходит не через n сложений, а через log₂ n, и ошибка растёт как логарифм; так складывает numpy.sum. Точное суммирование не округляет ни разу до конца: math.fsum хранит сумму как несколько чисел float, которые вместе представляют её без потерь (алгоритм Джонатана Шевчука), а в Go то же можно сделать длинным big.Float.

А деньги дробными числами не считают вовсе. 149,90 тенге — это 14 990 тиынов, целое число, и сумма целых точна, пока помещается в 64 бита, то есть до 92 квадриллионов тенге. Базы данных для денег держат тип DECIMAL, в Python есть модуль decimal.

Код

Пять способов на двух опытах и гармонический ряд в float32:

Поправки Кэхэна и Ноймайера на этих опытах дают точный ответ. Попарное суммирование 149,90 ошиблось на одну единицу младшего разряда — гораздо лучше простого цикла, но не точно. Точная сумма — та, что напечатана в десятичной строке, с одной оговоркой: в первом опыте слагаемое — не одна десятая, а ближайшее к ней двоичное число, чуть большее, и точная сумма миллиона таких чисел округляется к 100 000 только потому, что разница меньше половины шага.

Сколько это стоит

Кэхэн и Ноймайер — O(n) и O(1) памяти, но четыре арифметические операции на слагаемое вместо одной. На практике разница во времени часто меньше: длинный цикл упирается в чтение памяти, а не в сложения. Попарное — те же O(n) сложений плюс рекурсия, которую реализации разворачивают в блоки. Точное — O(n) в среднем, но с длинными числами в Go заметно медленнее.

Когда что брать. Сумма сотни чисел — простой цикл. Миллионы слагаемых в модели, накопление времени или координат шаг за шагом — поправка Кэхэна или Ноймайера. Итог, который должен совпасть до последнего бита, — точная сумма. Деньги — целые числа в минимальных единицах.

Задачи

Задача — посчитать и сумму, которую даёт обычный цикл, и точную. Среди тестов — суточный дождь в Алматы за 85 лет и числа, которые гасят друг друга.