Обратите внимание, что в примерах статьи используется тип данных float32. Простой перевод примеров на float64 может создать иллюзию, что использование float64 решает все проблемы с точностью, чего на самом деле не происходит.
Мы выбрали float32, чтобы примеры и числа оставались компактными, а основные идеи алгоритма были проще для восприятия.
Преимущества:
- Существенно снижает ошибки округления при суммировании
- Простая реализация
Недостатки:
- Добавляет накладные расходы, так как требует дополнительных операций
- Не устраняет ошибки полностью, особенно когда речь идет об очень больших количествах операций
- Код
- Вывод
package main
import (
"fmt"
)
func main() {
var sum float32
var c float32 // компенсация ошибки
for i := 0; i < 100_000; i++ {
x := float32(0.1)
y := x - c
t := sum + y
c = (t - sum) - y
sum = t
if c != 0 && (i+1%10_000) == 0 {
fmt.Printf("Шаг %d: sum: %.23f; компенсация c = %.23f\n", i+1, sum, c)
}
}
fmt.Printf("Итоговая сумма: %.23f, итоговая компенсация c = %.23f\n", sum, c)
}Шаг 10000: sum: 1000.00000000000000000000000; компенсация c = -0.00001490116119384765625
Шаг 20000: sum: 2000.00000000000000000000000; компенсация c = -0.00002980232238769531250
Шаг 30000: sum: 3000.00000000000000000000000; компенсация c = -0.00004470348358154296875
Шаг 40000: sum: 4000.00000000000000000000000; компенсация c = -0.00005960464477539062500
Шаг 50000: sum: 5000.00000000000000000000000; компенсация c = -0.00007450580596923828125
Шаг 60000: sum: 6000.00000000000000000000000; компенсация c = -0.00008940696716308593750
Шаг 70000: sum: 7000.00000000000000000000000; компенсация c = -0.00010430812835693359375
Шаг 80000: sum: 8000.00000000000000000000000; компенсация c = -0.00011920928955078125000
Шаг 90000: sum: 9000.00000000000000000000000; компенсация c = -0.00013411045074462890625
Шаг 100000: sum: 10000.00000000000000000000000; компенсация c = -0.00014901161193847656250
Итоговая сумма: 10000.00000000000000000000000, итоговая компенсация c = -0.00014901161193847656250
Идея алгоритма заключается в том, что на каждом шаге мы корректируем ошибку суммирования. Для этого мы постоянно вычисляем погрешность, которая происходит на каждом шаге сложения и учитываем ее на следующих шагах.
В нашем примере значения выглядят «идеальными», но на самом деле это не так.
Почему значения на шагах, которые кратны 10 выглядят идеально?
Так как каждую итерацию в коде мы прибавляем 0.1, а float32 хранит это значение как степень двойки, то фактически мы имеем значение ~0.10000000149.
Если мы сложим несколько таких числе, то ошибки будут накоплены и компенсированны, в итоге каждые 10 шагов сумма будет округлять к ближайшему числу кратному 1.0.
0.1×5≈0.5
в float32 с компенсацией Кэхана
Важно отметить, что математически «обнуление» будет происходить не каждые 5 шагов, а через N шагов, при котором сумма становится близкой к числу которое кратно единице.
Когда Кэхан может быть недостаточным?
Если разброс наших чисел значительный, то компенсация Кэхана будет давать заметную ошибку, например:
package main
import (
"fmt"
"math"
)
// KahanSum выполняет суммирование Кэхана
func KahanSum(numbers []float64) float64 {
var sum, c float64
for _, x := range numbers {
y := x - c
t := sum + y
c = (t - sum) - y
sum = t
}
return sum
}
func main() {
// Создаём массив с большим разбросом: очень большие и очень маленькие числа
numbers := []float64{
1e16, 1.0, -1e16, 2.0, 3.0,
}
// Наивное суммирование
var naive float64
for _, x := range numbers {
naive += x
}
// Суммирование Кэхана
kahan := KahanSum(numbers)
fmt.Printf("Наивное суммирование: %.16f\n", naive) // 5.00000
fmt.Printf("Суммирование Кэхана: %.16f\n", kahan) // 5.00000
fmt.Printf("Ожидаемый результат: 6.0\n", ) // 6.0
}
Что происходит в примере?
По законам математики мы должны получить 6.0, так как 1.0 + 2.0 + 3.0 = 6.0, а 1e16 и -1e16 должны дать 0.
Рассмотрим подробнее, что происходит, на примере суммирования слайса:
[]float64{1e16, 1.0, -1e16, 2.0, 3.0,}Шаг 1.
- Слагаемое:
1e16 - Сумма:
0 + 1e16 -> 1e16 - Погрешность:
0.0
Так как число помещается в мантису, то сложение на первом шаге происходит без ошибок и погрешностей.
Шаг 2.
- Слагаемое:
1.0(52 бита дробной части) - Сумма:
1e16 + 1 -> 1e16 - Погрешность:
0
Здесь начинается потеря точности. 1e16 уже занимает все значимые биты мантиссы, а мы пытаемся прибавить маленькое число (1.0). В результате единица теряется — её невозможно точно отразить в сумме, и Kahan не может компенсировать эту потерю полностью на этом шаге.
Шаг 3.
- Слагаемое:
-1e16 - Текущая сумма:
1e16 - Погрешность:
0
Так как мы прибавляем то же самое большое число, но со знаком минус, то мы получаем 0. А наша потерянное число 1.0 слишком мало для компенсации, таким образом мы получаем нулевую погрешность.
Потеря точности произошла из-за того, что мантисса хранит только 52 бита, и когда мы пытаемся сложить 1e16 и 1, то получается, что число выходит за границы диапазона.
Большое число: 1e16
Мантисса (52 бита):
[1 0 1 0 1 1 ... 0 0 0] <- все значимые биты заняты
Малое число: 1.0
Мантисса (52 бита):
[0 0 0 0 ... 0 0 1] <- единица в младших битах
Сумма: 1e16 + 1.0
Смещаем маленькое число вправо, чтобы совпало с порядком большого числа:
[1 0 1 0 1 1 ... 0 0 0] <- младшая единица сдвинута за предел мантиссы
Результат: 1e16 <- маленькое число потеряноШаг 4.
- Слагаемое:
2.0 - Текущая сумма:
0 - Погрешность:
0
Теперь у нас нет огромных чисел, поэтому мы без проблем прибавляем 2.
Шаг 5.
- Слагаемое:
3.0 - Текущая сумма:
2.0 - Погрешность:
0
Как и в шаге 4 число увеличивается без погрешности.
На выходе мы получили 5
Для решения проблем больших чисел используются следующие алгоритмы и подходы:
Попарное суммирование (Pairwise summation) позволяет уменьшать накопление ошибок за счет попраного сложения одинаковых чисел с разными знаками.
Различные улучшения алгорима Кэхэна: суммирование Ноймайера (Neumaier summation), компенсирующее суммирование (Compensated summation).
Использованияе библиотек для вычислений с произвольной точностью.
Сортировка перед суммированием – сначала складываем маленькие числа, а потом большие за счет чего уменьшаем потерю точности.
Эти и другие подобные подходы выходят за рамки нашего цикла статей, поэтому мы не будем рассматривать их подробно.
Leave a Reply