Компенсация погрешностей при операциях с числами с плавающей запятой

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

В данной работе примеры приведены на языке програмиирования C.

Представление вещественных чисел

Рассмотрим представление конечного действительного числа в стандарте IEEE 754-2008 в виде выражения, которое характеризуется тремя элеменами: S (0 или 1), мантисса M и порядок E:

v = -1^S * b^{(E — BIAS)} * M

• основанием b (стандартом определено три двоичных формата и два десятичных);
• длиной мантиссы p, определяющей точность представления числа, измеряется в цифрах (двоичных или десятичных);
• максимальным значением порядка emax;
• минимальным значением порядка emin; стандарт требует, чтобы выполнялось условие:
  

   emin = 1-emax 

• BIAS (смещение); порядок записывается в т. н. смещенной форме, т. е.
  реальное значение показателя степени равно E-BIAS.

Ниже в таблице даны параметры стандартных форматов чисел с плавающей запятой. Здесь: w — ширина битового поля для представления порядка, t — ширина битового поля для представления мантиссы, k — полная ширина битовой строки.



В следующией таблице представлены диапазоны изменения и точность стандартных 32 и 64-х разрядных форматов вещественных чисел с плавающей запятой.



Здесь «Точноть, Эпсилон» — наименьшее число, для которого истинно выражение:

1 + EPSILON != 1

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

Рассмотрим представление вещественного числа для типа данных float(Binary32).



В данном примере:

• Число V = 0,15625₁₀
• Знак S = 0, т.е. +
• Порядок (E — BIAS) = 01111100₂ — 01111111₂ = 124₁₀ — 127₁₀ = -3₁₀
• Мантисса M = 1,01000000000000000000000₂

Таким образом, число V = 1,01₂*2₁₀^-3₁₀ = 101₂*2₁₀^-5₁₀ = 5₁₀*2₁₀^-5₁₀ = 0,15625₁₀

Заметим, что при выполнении операций с вещественными числами часто результат не будет помещаться в N-разрядные матиссы, т.е. будет происходить округление. Округление в одной вычислительной операции не превывает порядок EPSILON/2, но когда нам требуется делать много операций, то для повышения точности вычисления результата нам надо научиться находить, насколько именно результат каждой конкретной операции округлён.

Более подробно про округления и нарушение аксиоматики рассказано в файле makarov_float.pdf (ссылка на материал внизу).

Способы нахождения погрешностей округления у чисел с плавающей запятой

Данную проблему исследовали многие специалисты, наиболее известными из них являются: Дэвид Голдберг, Уильям Кэхен, Джонатан Ричард Шевчук.

Ниже мы рассмотрим алгоримы нахождения погрешностей округления, приведенные в работе Шевчука, на примере двух функций:

• TwoSum — функция нахождения погрешности округления при сложении.
• TwoProduct — функция нахождения погрешности округления при умножении.

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

Говоря, что число a является p-битовым числом, имеется ввиду, что длина матиссы числа a представляется p битами.

TwoSum

Теорема: пусть числа a и b являются p-битовыми числами с плавающей запятой, где p >= 3, тогда следуя данному алгоритму мы получим 2 числа: x и y, для которых выполняется условие: a + b = x + y. Причем x является апроксиммацией суммы a и b, а y является погрешностью округления вычисления числа x.

TwoProduct

В действительности, алогритм находжения погрешности огругления при умножении двух вещественных чисел состоит из 2-х функций: Split — вспомогательной функции и функции TwoProduct, где мы находим погрешность.

Рассмотрим алгоритм функции Split.

Split

Теорема: число a — p-битовое число с плавающей запятой, где p >= 3. Выберем точку разрыва s, где p/2 <= s <= p-1. Следуя алгоритму, получим (p — s)-битное число — число a_hi и (s-1)-битное число — a_lo, где |a_hi| >= |a_low| и a = a_hi + a_low.



Теперь перейдем к анализу функции TwoProduct.

TwoProduct

Теорема: пусть числа a и b являются p-битовыми числами с плавающей запятой, где p >= 6. Тогда, выполняя данный алгоритм, мы получим 2 числа: x и y, для которых выполняется условие: a*b = x + y. Причем x является апроксиммацией произведения чисел a и b, а число y является погрешностью округления вычисления числа x.



Видно, что конечный результат представляется парой N-битовых вещественных чисел: результат + погрешность. Причем второе должно иметь порядок не более EPSILON по отношению к первому.

Пример компесации погрешностей округления у чисел с плавающей запятой

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

Погрешность округления суммы

Приведем пример простейшей программы суммирования:

#include <stdio.h>
#include <math.h>

int main() {
    float val = 2.7892;
    printf("%0.7g \n", val);
    val = val/10000000000.0;
    float result = 0.0;
    for (long long i = 0; i < 10000000000; i++) {
        result += val;
    }
    printf("%0.7g \n", result);
    return 0;
}

В результате работы программы, мы должны получить два одинаковых числа: 2,7892 и 2,7892. Но в консоль было выведено: 2,7892 и 0,0078125. Отсюда видно, что погрешность накопилась очень большая.

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

#include <stdio.h>
#include <math.h>

float TwoSum(float a, float b, float& error) {
    float x = a + b;
    float b_virt = x - a;
    float a_virt = x - b_virt;
    float b_roundoff = b - b_virt;
    float a_roudnoff = a - a_virt;
    float y = a_roudnoff + b_roundoff;
    error += y;
    return x;
}

int main() {
    float val = 2.7892;
    printf("%0.7g \n", val);
    val = val/10000000000.0;
    float result = 0.0;
    float error = 0.0;
    for (long long i = 0; i < 10000000000; i++) {
        result = TwoSum(result, val, error);
    }
    result += error;
    printf("%0.7g \n", result);
    return 0;
}

В итоге мы получаем 2 числа: 2,7892 и 0,015625. Результат улучшился, но погрешность все равно дает о себе знать. В данном примере не была учтена погрешость, возникающая в операции сложения, в функции TwoSum():

error += y;

Будем компенсировать результат на каждой итерации цикла и перезаписывать ошибку в переменную, которая накапливает погрешность округления в операции суммирования. Для этого модифицируем функцию TwoSum(): добавим переменную isNull типа bool которая указывает на то, стоит ли нам накапливать погрешность или же стоит ее перезаписать.

В итоге, result будет представляться 2 переменными: result — главная переменная, error1 — погрешность операции result += val.

Код будет выглядеть следующим образом:

#include <stdio.h>
#include <math.h>

float TwoSum(float a, float b, float& error, bool isNull) {
    float x = a + b;
    float b_virt = x - a;
    float a_virt = x - b_virt;
    float b_roundoff = b - b_virt;
    float a_roudnoff = a - a_virt;
    float y = a_roudnoff + b_roundoff;
    if (isNull) {
        error = y;
    } else {
        error += y;
    }
    return x;
}

int main() {
    float val = 2.7892;
    printf("%0.7g \n", val);
    val = val/10000000000.0;
    float result = 0.0;
    float error1 = 0.0;
    for (long long i = 0; i < 10000000000; i++) {
        result = TwoSum(result, val, error1, false);
        result = TwoSum(error1, result, error1, true);
    }
    printf("%0.7g \n", result);
    return 0;
}

Программа выведет числа: 2,7892 и 2,789195.

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

val = val*(1/10000000000.0);

Учтем эту погрешность путем добавления функций учета погрешностей умножения, которые разработал Д.Р. Шевчук. При этом переменная val будет представлена двумя переменными:

val_real = val + errorMult

Таким образом, result будет представляться 3 переменными: result — главная переменная, error1 — погрешность операции result += val, error2 — погрешность операции result += errorMult.

Так же мы будем складывать переменные error1 и error2, а погрешность от этой операции записывать в error2.

В итоге, код:

#include <stdio.h>
#include <math.h>

float TwoSum(float a, float b, float& error, bool isNull) {
//isNull отвечает за то, стоит ли нам накопить возникающую погрешность
// или же стоит перепизаписать ее
    float x = a + b;
    float b_virt = x - a;
    float a_virt = x - b_virt;
    float b_roundoff = b - b_virt;
    float a_roudnoff = a - a_virt;
    float y = a_roudnoff + b_roundoff;
    if (isNull) {
        error = y;
    } else {
        error += y;
    }
    return x;
}

void Split(float a, int s, float& a_hi, float& a_lo) {
    float c = (pow(2, s) + 1)*a;
    float a_big = c - a;
    a_hi = c - a_big;
    a_lo = a - a_hi;
}

float TwoProduct(float a, float b, float& err) {
    float x = a*b;
    float a_hi, a_low, b_hi, b_low;
    Split(a, 12, a_hi, a_low);
    Split(b, 12, b_hi, b_low);
    float err1, err2, err3;
    err1 = x - (a_hi*b_hi);
    err2 = err1 - (a_low*b_hi);
    err3 = err2 - (a_hi*b_low);
    err += ((a_low * b_low) - err3);
    return x;
}

int main() {
    float val = 2.7892;
    printf("%0.7g \n", val);
    float errorMult = 0;//погрешность умножения
    val = TwoProduct(val, 1.0/10000000000.0, errorMult);
    float result = 0.0;
    float error1 = 0.0;
    float error2 = 0.0;
    for (long long i = 0; i < 10000000000; i++) {
        result = TwoSum(result, val, error1, false);
        result = TwoSum(result, errorMult, error2, false);
        error1 = TwoSum(error2, error1, error2, true);
        result = TwoSum(error1, result, error1, true);
    }
    printf("%0.7g \n", result);
    return 0;
}

В консоль были выведены следующие числа: 2,7892 и 2,789195.

Это говорит о том, что погрешность округления умножения достаточно мала, чтоб проявиться даже на 10 миллиардах итераций. Данный результат является максимально приближенным к исходному числу, если учитывать погрешности в операциях сложения и умножения. Для получение более точного результата можно ввести дополнительные переменные, учитывающие погрешности. Скажем, добавить переменную, учитывающую погрешность в операции накопления основной погрешности в функции TwoSum(). Тогда эта погрешность будет иметь порядок EPSILON², по отношению к главному результату (первая погрешность будет иметь порядок EPSILON).

Погрешность округления умножения

Посчитаем число 1.0012¹⁰¹ в цикле, т.е. сделаем следующие:

#include <stdio.h>
#include <math.h>

int main() {
    float val = 1.0012;
    float result = 1.0012;
    for (long long i = 0; i < 100; i++) {
        result *= val;
    }
    printf("%0.15g \n", result);
    return 0;
}

Заметим, что точный результат, с точностью до пятнадцатого знака после запятой, равен 1.128768638496750. Мы получим: 1.12876391410828. Как видно, погрешность оказалось достаточно большой.

Выведем переменную val, приведя ее к типу данных double, и посмотрим, что в нее на самом деле записалось:

printf("%0.15g \n", (double)val);

Мы получим число 1.00119996070862. Это говорит о том, что в программировании даже самая точная константа не является ни надежной, ни константой. Поэтому, наш реальный точный результат будет равен 1.128764164435784, с точностью до пятнадцатого знака после запятой.

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

Код:

#include <stdio.h>
#include <math.h>

float TwoSum(float a, float b, float& error, bool isNull) {
//isNull отвечает за то, стоит ли нам накопить возникающую погрешность
// или же стоит перепизаписать ее
    float x = a + b;
    float b_virt = x - a;
    float a_virt = x - b_virt;
    float b_roundoff = b - b_virt;
    float a_roudnoff = a - a_virt;
    float y = a_roudnoff + b_roundoff;
    if (isNull) {
        error = y;
    } else {
        error += y;
    }
    return x;
}

void Split(float a, int s, float& a_hi, float& a_lo) {
    float c = (pow(2, s) + 1)*a;
    float a_big = c - a;
    a_hi = c - a_big;
    a_lo = a - a_hi;
}

float TwoProduct(float a, float b, float& err) {
    float x = a*b;
    float a_hi, a_low, b_hi, b_low;
    Split(a, 12, a_hi, a_low);
    Split(b, 12, b_hi, b_low);
    float err1, err2, err3;
    err1 = x - (a_hi*b_hi);
    err2 = err1 - (a_low*b_hi);
    err3 = err2 - (a_hi*b_low);
    err += ((a_low * b_low) - err3);
    return x;
}

int main() {
    float val = 1.0012;
    float result = 1.0012;
    float errorMain = 0.0;
    for (long long i = 0; i < 100; i++) {
        result = TwoProduct(result, val, errorMain);
        result = TwoSum(errorMain, result, errorMain, true);
    }
    printf("%0.15g \n", result);
    return 0;
}

Программа выводит следующее число: 1.12876415252686. Мы получили погрешность 1.0e-008, что меньше чем EPSILON/2 для типа данных float. Таким образом, данный результат можно считать достаточно хорошим.

Итоги

1) В данной работе было рассмотрено представление чисел с плавающей запятой в формате стандарта IEEE 754-2008.
2) Был показан способ нахождения погрешностей округления в операциях сложения и умножения у чисел с плавающей запятой.
3) Были рассмотрены простые примеры компенсации погрешностей округления у чисел с плавающей запятой.

Работу выполнил Виктор Фадеев.
Консультировал Макаров А.В.

P.S. Спасибо за найденные ошибки пользователям:
xeioex, Albom.

Использованная литература

• Работы Макарова Андрея Владимировича на тему «Представление вещественных чисел». В файле makarov_float .pdf показано, как возникают погрешности округления у чисел с плавающей запятой.
• Википедия. Число одинарной точности.
• Adaptive Precision Floating-Point Arithmetic
  and Fast Robust Geometric Predicates
  Jonathan Richard Shewchuk