Русская Википедия:Алгоритм Кэхэна
В вычислительной математике алгоритм Кэхэна (также известный как компенсационное суммирование) — это алгоритм вычисления суммы последовательности чисел c плавающей запятой, который значительно уменьшает Шаблон:Нп3 по сравнению с наивным подходом. Уменьшение погрешности достигается введением дополнительной переменной для хранения нарастающей суммы погрешностей.
В частности, простое суммирование <math>n</math> чисел в худшем случае имеет погрешность, которая растёт пропорционально <math>n</math> и при суммировании случайных чисел имеет среднее квадратичное отклонение, пропорциональное <math>\sqrt{n}</math> (ошибки округления вызывают случайное блуждание)[1]. При компенсационном суммировании погрешность даже в худшем случае не зависит от <math>n</math>, так что большое число слагаемых могут быть просуммированы с погрешностью, зависящей только от точности числа с плавающей запятой[1].
Авторство алгоритма приписывают Уильяму Кэхэну[2].
Алгоритм
В псевдокоде алгоритм можно записать так:
function kahan_sum(input)
var sum = 0.0
var c = 0.0 // Сумма погрешностей.
for i = 1 to len(input) do
y = input[i] - c // Пока всё хорошо: c - ноль.
t = sum + y // Увы, sum велика, y мало, так что младшие разряды y потеряны.
c = (t - sum) - y // (t - sum) восстанавливает старшие разряды y; вычитание y восстанавливает -(младшие разряды y)
sum = t // Алгебраически, c всегда должна равняться нулю. Берегитесь слишком оптимизирующих компиляторов!
// В следующий раз потерянные младшие разряды будут заново прибавлены к y.
return sum
Пример исполнения
В данном примере будем использовать десятичные дроби. Компьютеры обычно используют двоичную арифметику, но иллюстрируемый алгоритм от этого не меняется. Представим что используется шестиразрядная арифметика с плавающей точкой, sum было присвоено значение 10000,0, и следующие два значения input(i) равны 3,14159 и 2,71828. Точный результат равен 10005,85987, что округляется до 10005,9. При простом суммировании порядок каждого входящего значения был бы выровнен с sum, и много младших разрядов было бы потеряно в результате округления. Первый результат после округления был бы 10003,1. Второй результат был бы 10005,81828 до округления, и 10005,8 — после, то есть в последний знак результата была бы внесена ошибка.
При компенсационном суммировании мы бы получили правильный округлённый результат 10005,9.
Предположим, что начальное значение c — 0.
y = 3.14159 - 0 y = input[i] - c t = 10000.0 + 3.14159 = 10003.1 Много разрядов потеряно! c = (10003.1 - 10000.0) - 3.14159 Это нужно вычислять как записано! = 3.10000 - 3.14159 Восстановлена не вошедшая в t часть y , а не всё исходное y. = -0.0415900 Завершающие нули показаны, потому что это шестиразрядная арифметика. sum = 10003.1 Таким образом не все разряды из input(i) включены в sum.
Сумма настолько велика, что сохраняются только старшие разряды слагаемого. К счастью, на следующем шаге c хранит погрешность.
y = 2.71828 - -0.0415900 Учитывается погрешность с предыдущего шага. = 2.75987 Её порядок не слишком отличается от y. Большинство разрядов учтены. t = 10003.1 + 2.75987 Но только немногие разряды попадают в sum. = 10005.85987, округляется до 10005.9 c = (10005.9 - 10003.1) - 2.75987 Здесь извлекается то, что пришло = 2.80000 - 2.75987 В данном случае слишком много. = 0.040130 Так или иначе излишек будет вычтен в следующий раз. sum = 10005.9 Точный результат: 10005.85987, что корректно округлено до 6 знаков.
Таким образом сложение происходит в двух переменных: sum хранит сумму, и c хранит часть суммы, которая не попала в sum, чтобы быть учтённой в sum на следующей итерации. Хотя суммировать, храня младшие разряды в c лучше, чем не храня их нигде, это всё же не так точно, как производить вычисление, используя ввод двойной точности. Тем не менее, просто увеличивать точность вычислений в целом не практично; если input уже имеет двойную точность, немногие системы могут предоставить учетверённую точность вычислений, и, если бы могли, ввод тоже мог бы иметь учетверённую точность.
Недостатки
К сожалению, алгоритм Кэхэна не гарантирует защиту от потери точности, связанной с наличием в сумме слагаемых с существенно различными порядками. Это может происходить, если последовательность суммирования выбрана неудачно. Чтобы убедиться в этом, достаточно в рассмотренном выше примере вместо числа 2,71828 взять −10000. На последнем шаге алгоритма имеем следующее:
y = -10000.0 - -0.0415900 = -10000.0 Результат округляется до 6 значащих цифр t = 10003.1 + -10000.0 = 3.10000 c = (3.10000 - 10003.1) - -10000.0 = -10000.0 + 10000.0 = 0 sum = 3.10000
Точный результат: 3,14159, то есть произошла потеря точности.
Надо заметить, что если предварительно упорядочить слагаемые по убыванию их абсолютной величины, то в результате применения алгоритма получим результат 3,14159, где все знаки верны.
Примечания
Ссылки
- Evan Manning, Floating-point Summation Dr. Dobb’s Journal September, September 01, 1996