Что такое double rounding error (ошибка двойного округления)?

#ieee_754

Слышал, что при работе с числами с плавающей запятой возможна так называемая ошибка
двойного (повторного) округления (double rounding error), особенно когда имеют место
переходы от одинарной точности (float) к двойной (double) и наоборот. Как это понимать?
    

Ответы

Ответ 1

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

В чём суть ошибки?

Ошибка двойного округления иногда возникает в тех случаях, когда вместо одного округления
выполняется два последующих друг за другом округления. Рассмотрим пример в десятичной
системе счисления. У нас есть число 1,34999. Если мы сразу округлим его до одного знака
после десятичной запятой, то получим 1,3. Но если сначала округлить до двух знаков,
а потом до одного, то получим 1,35, а затем 1,4. Как мы видим, результаты получились
разными, причём существенно!

Как это связано с моей программой?

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


Вы используете тип данных double размером 64 бита, но на вашем компьютере нет команд,
которые работают с double (или компилятор их не поддерживает), а есть только сопроцессор
x87. Таким образом, любая операция с double выглядит так: происходит расширение типа
с 64-х до 80-ти битов, затем, собственно, необходимые вычисления, в которых результат
округляется до 63-х битов дробной части мантиссы (такова структура 80-битового числа
с плавающей запятой), а затем происходит второе округление до 52-х битов дробной части
мантиссы в double.


Пример программы. Числа для кода я позаимствовал у @Mark Dickinson из англоязычного
SO. Вы знаете, что long double размером 80 битов в настоящий момент не поддерживается
(например, в VC++ точно), а заменяется на обычный double. Поэтому пример работать не
будет, если не указать специальную опцию компилятора, которая заставит его работать
через x87. У каждого компилятора своя такая опция (если вообще есть). 

#include 

int main (void) {
  double a = 1e16 + 2.9999, b;
  long double x = 1e16L, y = 2.9999L;
  b = (double)(x+y);
  printf ("a = %16.0lf, b = %16.0lf\n", a, b);
  return 0;
}


Программа выдаст два разных целых числа: 10000000000000002 и 10000000000000004, показывая,
тем самым, разницу, возникшую в результате двойного округления.


Вы используете функции, которые созданы специально для double, но подаёте им на вход
числа float. Скажем, есть функция sqrt для double. Если передать ей на вход число float,
то произойдёт расширение его до double, затем функция вычислит результат с округлением,
а затем выполнится второе округление, если вы присваиваете результат обратно во float.
К сожалению, каких-то реальных примеров подобной ошибки мне отыскать не удалось, но
чисто теоретически любая подмена float на double и обратно может дать ошибку двойного
округления. Это будет хорошо видно на последнем, самом любопытном примере.
Данное наблюдение я подсмотрел в статье “Double Rounding Errors in Floating-Point
Conversions”, но сделал свой собственный пример. Общий смысл таков, что константа может
быть неверно преобразована в бинарный формат на этапе компиляции, если произойдёт двойное
округление. Это может произойти когда когда вы забыли суффикс f для чисел типа float
или если пользуетесь VC++ любой версии.


Для демонстрации за основу будет взято число 1+2-23+2-24-2-54. Оно равно 1.000000178813934270660723768742172978818416595458984375.
Это число интересно тем, что в двоичном формате оно равно 

1,000000000000000000000010111111111111111111111111111111(2)

Полужирным шрифтом я отметил биты 24 и 53-54. Если читатель не знаком с правилами
округления чисел в формате IEEE-754, то ему придётся сначала их изучить. Например,
здесь. Если мы хотим перевести данное число в формат float, то никаких проблем нет:
24-й бит равен нулю, а значит округление до 23-х бит, необходимых float, произойдёт
вниз без вопросов. Получим число 1,00000000000000000000001, что в формате IEEE-754
будет иметь вид 0x3f800001. Если же мы сначала приводим число к типу double, то округление
до 52-х бит сначала произойдёт вверх, и мы получим 1,0000000000000000000000110000000000000000000000000000,
а затем, из-за того, что 23-й и 24-й биты теперь равны 1, в результате второго округления
до float получим 1,00000000000000000000010 (снова округление произошло вверх), что
в шестнадцатеричной записи будет иметь вид 0x3f800002. Как видно, мы получили два разных
числа.

Рассмотрим демонстрационный код (не плюйтесь, что в коде есть грязные трюки, здесь
я хочу просто и без хитростей показать биты числа):

#include 

typedef unsigned int u32;
typedef unsigned long long u64;

int main (void) {
  float a, b;
  double c;
  a = 1.000000178813934270660723768742172978818416595458984375f;
  b = 1.000000178813934270660723768742172978818416595458984375;
  c = 1.000000178813934270660723768742172978818416595458984375;
  printf ("a = %08x, b = %08x, c = %016llx\n", *(u32*)&a, *(u32*)&b, *(u64*)&c);
  return 0;
}


Здесь переменная a инициализируется правильно, мы используем суффикс f, инициализация
переменной b происходит так, будто программист забыл про суффикс f (к счастью, многие
компиляторы ему об этом напомнят), поэтому здесь произойдёт двойное округление (сначала
константа превратиться в double, а затем будет присвоена к float), а переменная c здесь
просто так, для демонстрации промежуточного округления в double. Функция printf выводит
шестнадцатеричный код всех трёх переменных. Компилятор Tiny C выдаёт верный ответ:

a = 3f800001, b = 3f800002, c = 3ff0000030000000


Такой же точно ответ выдали Intel C 15 и Clang. 

Компилятор VC 2015 выдал

a = 3f800002, b = 3f800002, c = 3ff0000030000000


Это значит, что он ИГНОРИРУЕТ суффикс f, если он установлен, и в любом случае приводит
результат сначала к double, а потом к float, если нужно. По-хорошему, разработчики
компилятора заработали этим сильный удар линейкой по пальцам.

Копилятор GCC из MinGW выдал также неверный ответ

a = 3f800001, b = 3f800001, c = 3ff0000030000000


То есть с точки зрения исходного желания программиста ответ верный (и здесь он молодец,
что угадал это желание), но с точки зрения Стандарта - нет.

Что же теперь делать?

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

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