Алгоритм Кэхэна

В вычислительной математике алгоритм Кэхэна (также известный как компенсационное суммирование) — это алгоритм вычисления суммы последовательности чисел c плавающей запятой, который значительно уменьшает вычислительную погрешность по сравнению с наивным подходом. Уменьшение погрешности достигается введением дополнительной переменной для хранения нарастающей суммы погрешностей.

В частности, простое суммирование n чисел в худшем случае имеет погрешность, которая растет пропорционально n и при суммировании случайных чисел имеет среднее квадратичное отклонение равное $\sqrt{n}$ (ошибки округления вызывают случайное блуждание)^[1]. При компенсационном суммировании погрешность даже в худшем случае не зависит от n, так что большое число слагаемых могут быть просуммированы с погрешностью, зависящей только от точности числа с плавающей запятой^[1].

Авторство над алгоритмом приписывают Уильяму Кэхэну^[2].

Алгоритм

В псевдокоде алгоритм можно записать так:

function KahanSum(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.
    = -.0415900                     Завершаюшие нули показаны, потому что это шестиразрядная арифметика.
sum = 10003.1                       Таким образом не все разряды из input(i) включены в sum.

Сумма настолько велика, что сохраняются только старшие разряды слагаемого. К счастью, на следующем шаге c хранит погрешность.

  y = 2.71828 - -.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             В данном случае слишком много.
    = .040130                       Так или иначе излишек будет вычтен в следующий раз.
sum = 10005.9                       Точный результат: 10005.85987, что корректно округлено до 6 знаков.

Таким образом сложение происходит в двух переменных: sum хранит сумму, и c хранит часть суммы, которая не попала в sum, чтобы быть учтенной в sum на следующей итерации. Хотя суммировать, храня младшие разряды в c лучше, чем не храня их нигде, это все же не так точно, как производить вычисление, используя ввод двойной точности. Тем не менее, просто увеличивать точность вычислений в целом не практично; если input уже имеет двойную точность, немногие системы могут предоставить учетверенную точность вычислений, и, если бы могли, ввод тоже мог бы иметь учетверённую точность!

Примечания

1. ↑ ^{1 2} Higham, Nicholas J. (1993), "«The accuracy of floating point summation»", SIAM Journal on Scientific Computing Т. 14 (4): 783–799, DOI 10.1137/0914050 
2. ↑ Kahan, William (January 1965), "«Further remarks on reducing truncation errors»", Communications of the ACM Т. 8 (1): 40, DOI 10.1145/363707.363723 

Ссылки

• Floating-point Summation, Dr. Dobb’s Journal September, 1996

Для улучшения этой статьи желательно?: Переработать оформление в соответствии с правилами написания статей.