Double ошибка округления

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

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

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

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

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

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

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

#include <stdio.h>

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.0lfn", a, b);
  return 0;
}

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

1. Вы используете функции, которые созданы специально для double, но подаёте им на вход числа float. Скажем, есть функция sqrt для double. Если передать ей на вход число float, то произойдёт расширение его до double, затем функция вычислит результат с округлением, а затем выполнится второе округление, если вы присваиваете результат обратно во float. К сожалению, каких-то реальных примеров подобной ошибки мне отыскать не удалось, но чисто теоретически любая подмена float на double и обратно может дать ошибку двойного округления. Это будет хорошо видно на последнем, самом любопытном примере.

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

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

1,000000000000000000000010111111111111111111111111111111₍₂₎

Полужирным шрифтом я отметил биты 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 <stdio.h>

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 = %016llxn", *(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

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

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

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

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

When I run the following code:

double x = 4.35;
double y = x * 100;
System.out.println(y);

I get 434.99999999999994 which is expected due to the nature of floating point precision. However when I run:

double x = 4.35;
double y = x * 1000;
System.out.println(y);

I get 4350.0 in the console. Similarly when I replace the 1000 with 10, I get 43.5, which is peculiar, since there is no rounding error for these numbers. It seems that there is only rounding error when I mulptiply the double by 10 to the power of the number of decimal places. Why is it like that?

Also, when I tried running the three blocks of code using 4.25 as my value for x, I got a nice 425.0 in the console. Does rounding error only apply to some decimal values?

EDIT: Dividing 4.35 by 10 also gave me a wacky decimal (0.43499999999999994 to be exact). Also, if I set x to 43.5, there were no rounding errors at all. Multiplying by 10 gave me 435.0.

EDIT 2: I think some people are misinterpreting my question. I’m asking why this happens and not how to avoid running into these errors.

30.12.2019Java

x = 0.7999999999999999
y = 0.8
false

SIGN FRACTION * 2 ^ EXP 

x = 0.7999999999999999
y = 0.8
false

a = 0.09999999999999998

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

Друзья, для вашего удобства статья также доступна в формате презентации на видео (почти 34 минуты), такой формат больше подойдёт тем читателям, которым трудно выстраивать нужные математические образы в голове, так как в презентации много иллюстративного материала. Информация на видео полностью идентична содержанию статьи. Действуйте, пожалуйста, как вам удобно.

Повторюсь, что это не научная, а научно-популярная статья, прочитав которую, вы кратко познакомитесь вот с чем.

• Трансцендентные элементарные функции (exp, sin, log, cosh и другие), работающие с арифметикой с плавающей запятой, округляются некорректно, иногда они допускают ошибку в последнем бите.
• Причина ошибок не всегда кроется в лени или низкой квалификации разработчиков, а в одном фундаментальном обстоятельстве, преодолеть которое современная наука пока не смогла.
• Существуют «костыли», которые позволяют худо-бедно справляться с обсуждаемой проблемой в некоторых случаях.
• Некоторые функции, которые должны делать вроде бы одно и то же, могут выдавать различный результат в одной и той же программе, например, exp2(x) и pow(2.0, x).

Чтобы понять содержание статьи, вам нужно быть знакомым с форматом чисел с плавающей запятой IEEE-754. Будет достаточно, если вы хотя бы просто понимаете что, например, вот это: 0x400921FB54442D18 – число пи в формате с удвоенной точностью (binary64, или double), то есть просто понимаете, что я имею в виду этой записью; я не требую уметь «на лету» делать подобные преобразования. А про режимы округления я вам напомню в этой статье, это важная часть повествования. Ещё желательно знать «программистский» английский, потому что будут термины и цитаты из западной литературы, но можете обойтись и онлайн-переводчиком.

Сначала примеры, чтобы вы сразу поняли, в чём состоит предмет разговора. Сейчас я дам код на C++, но если это не ваш язык, то уверен, вы всё равно без труда поймёте что написано. Посмотрите, пожалуйста, на этот код:

#include <stdio.h>
#include <cmath>

int main() {
  float x = 0.00296957581304013729095458984375f;  // Аргумент, записанный точно.
  float z;
  z = exp2f(x);  // z = 2**x одним способом.
  printf ("%.8fn", z);  // Вывод результата с округлением до 8 знаков после точки.
  z = powf(2.0f, x);  // z = 2**x другим способом
  printf ("%.8fn", z);  // Такой же вывод.
  return 0;
}

Число x намеренно записано с таким количеством значащих цифр, чтобы оно было точнопредставимым в типе float, то есть чтобы компилятор преобразовал его в бинарный код без округления. Ведь вы прекрасно знаете, что некоторые компиляторы не умеют округлять без ошибок (если не знаете, укажите в комментариях, я напишу отдельную статью с примерами). Далее в программе нам нужно вычислить 2^x, но давайте сделаем это двумя способами: функция exp2f(x), и явное возведение двойки в степень powf(2.0f, x). Результат, естественно, будет различным, потому что я же сказал выше, что не могут элементарные функции работать правильно во всех случаях, а я специально подобрал пример, чтобы это показать. Вот что будет на выходе:

1.00206053
1.00206041

Эти значения мне выдали четыре компилятора: Microsoft C++ (19.00.23026), Intel C++ 15.0, GCC (6.3.0) и Clang (3.7.0). Они отличаются одним младшим битом. Вот шестнадцатеричный код этих чисел:

0x3F804385  // Правильно
0x3F804384  // Неправильно

Запомните, пожалуйста, этот пример, именно на нём мы рассмотрим суть проблемы чуть ниже, а пока, чтобы у вас сложилось более ясное впечатление, посмотрите, пожалуйста, примеры для типа данных с удвоенной точностью (double, binary64) с некоторыми другими элементарными функциями. Привожу результаты в таблице. У правильных ответов (когда они есть) стоят символы «*» на конце.

Надеюсь, у вас не сложилось впечатления, что я специально взял какие-то прямо совсем уникальные тесты, которые с трудом можно отыскать? Если сложилось, то давайте состряпаем на коленках полный перебор всех возможных дробных аргументов для функции 2^x для типа данных float. Понятно, что нас интересуют только значения x между 0 и 1, потому что другие аргументы будут давать результат, отличающийся только значением в поле экспоненты и интереса не представляют. Сами понимаете:

Написав такую программу (скрытый текст будет ниже), я проверил функцию exp2f и то, сколько ошибочных значений она выдаёт на интервале x от 0 до 1.

Из программы ниже понятно, что число проверяемых аргументов x составило 197612997 штук. Получается, что, например, Microsoft С++ неверно вычисляет функцию 2^x для почти одного процента из них. Не радуйтесь, уважаемые любители GCC и Clang, просто именно эта функция в данных компиляторах реализована правильно, но полно ошибок в других.

#include <stdio.h>
#include <cmath>

    // Этими макросами мы обращаемся к битовому представлению чисел float и double
#define FAU(x) (*(unsigned int*)(&x))
#define DAU(x) (*(unsigned long long*)(&x))

    // Эта функция вычисляет 2**x точно до последнего бита для 0<=x<=1.
    // Страшный вид, возможно, не даёт сразу увидеть, что 
    // здесь вычисляется аппроксимирующий полином 10-й степени.
    // Промежуточные расчёты делаются в double (ошибки двойного округления тут не мешают).
    // Не нужно пытаться оптимизировать этот код через FMA-инструкции, 
    // практика показывает, что будет хуже, но... процессоры бывают разными.
float __fastcall pow2_minimax_poly_double (float x) {
  double a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10;
  DAU(a0) = 0x3ff0000000000001;
  DAU(a1) = 0x3fe62e42fefa3763;
  DAU(a2) = 0x3fcebfbdff845acb;
  DAU(a3) = 0x3fac6b08d6a26a5b;
  DAU(a4) = 0x3f83b2ab7bece641;
  DAU(a5) = 0x3f55d87e23a1a122;
  DAU(a6) = 0x3f2430b9e07cb06c;
  DAU(a7) = 0x3eeff80ef154bd8b;
  DAU(a8) = 0x3eb65836e5af42ac;
  DAU(a9) = 0x3e7952f0d1e6fd6b;
  DAU(a10)= 0x3e457d3d6f4e540e;
  return (float)(a0+(a1+(a2+(a3+(a4+(a5+(a6+(a7+(a8+(a9+a10*x)*x)*x)*x)*x)*x)*x)*x)*x)*x);
} 

int main() {
  unsigned int n = 0;  // Счётчик ошибок.
  // Цикл по всем возможным значениям x в интервале (0,1)
  // Старт цикла: 0x33B8AA3B = 0.00000008599132428344091749750077724456787109375
  // Это минимальное значение, для которого 2**x > 1.0f
  // Конец цикла: 0x3F800000 = 1.0 ровно.
  for (unsigned int a=0x33B8AA3B; a<0x3F800000; ++a) {  
   float x;
    FAU(x) = a;
    float z1 = exp2f (x);	// Подопытная функция.
    float z2 = pow2_minimax_poly_double (x);	// Точный ответ.
    if (FAU(z1) != FAU(z2)) {	// Побитовое сравнение.
      // Закомментируйте это, чтобы не выводить все ошибки на экран (их может быть много).
      //fprintf (stderr, "2**(0x%08X) = 0x%08X, but correct is 0x%08Xn", a, FAU(z1), FAU(z2));
      ++n;
    }		
  }
  const unsigned int N = 0x3F800000-0x33B8AA3B;  // Сколько всего аргументов было проверено.
  printf ("%u wrong results of %u arguments (%.2lf%%)n", n, N, (float)n/N*100.0f);
  return 0;
} 

Не буду утомлять читателя этими примерами, здесь главное было показать, что современные реализации трансцендентных функций могут неправильно округлять последний бит, причём разные компиляторы ошибаются в разных местах, но ни один не будет работать правильно. Кстати, Стандарт IEEE-754 эту ошибку в последнем бите разрешает (о чём скажу ниже), но мне всё равно кажется странным вот что: ладно double, это большой тип данных, но ведь float можно проверить полным перебором! Так ли сложно это было сделать? Совсем не сложно, и я уже показал пример.

В коде нашего перебора есть «самописная» функция правильного вычисления 2^x с помощью аппроксимирующего полинома 10-й степени, и написана она была за несколько минут, потому что такие полиномы выводятся автоматически, например, в системе компьютерной алгебры Maple. Достаточно задать условие, чтобы полином обеспечивал 54 бита точности (именно для этой функции, 2^x). Почему 54? А вот скоро узнаете, сразу после того как я расскажу суть проблемы и поведую о том, почему для типа данных учетверённой точности (binary128) в принципе сейчас невозможно создать быстрые и правильные трансцендентные функции, хотя попытки атаковать эту проблему в теории уже есть.

Округление по умолчанию и проблема, с ним связанная

Если вы не погружены в разработку математических библиотек, то нет ничего зазорного, чтобы забыть правило округления чисел с плавающей запятой по Стандарту IEEE-754, используемое по умолчанию. Поэтому я вам его напомню. Если вы всё хорошо помните, всё равно бросьте взгляд хотя бы на конец этого раздела, вас ждёт сюрприз: я покажу ситуацию, когда округлить число может оказаться очень сложным делом.

О том, что такое «округлить вверх» (к плюс бесконечности), «округлить вниз» (к минус бесконечности) или «округлить к нулю» вы легко вспомните по названию (если что, есть википедия). Основные сложности у программистов возникают с округлением «к ближайшему, но в случае равного удаления от ближайших — к тому, у которого последняя цифра чётная». Да, именно так переводится этот режим округления, который западной литературе называют короче: «Round nearest ties to even».

Данный режим округления используется по умолчанию и работает следующим образом. Если в результате расчётов длина мантиссы оказалась больше, чем может вместить в себя результирующий тип данных, округление выполняется к ближайшему из двух возможных значений. Однако может возникнуть ситуация, когда исходное число оказалось ровно посередине между двух ближайших, тогда в качестве результата выбирается то, у которого последний бит (после округления) окажется чётным, то есть равным нулю. Рассмотрим четыре примера, в которых нужно выполнить округление до двух битов после двоичной запятой:

1. Округлить 1,001001. Третий бит после запятой равен 1, но дальше есть ещё 6-й бит, равный 1, значит округление будет вверх, потому что исходное число ближе к 1,01, чем к 1,00.
2. Округлить 1,001000. Теперь округляем вниз, потому что мы находимся ровно посередине между 1,00 и 1,01, но именно у первого варианта последний бит будет чётным.
3. Округлить 1,011000. Мы посередине между 1,01 и 1,10. Придётся округлять вверх, потому что чётный последний бит именно у большего числа.
4. Округлить 1,010111. Округляем вниз, потому что третий бит равен нулю и мы ближе к 1,01, чем к 1,10.

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

1,00

1

0000000000000000000000000000000000001

Вам сейчас очевидно, что округление должно быть вверх, то есть к числу 1,01. Однако вы смотрите на число, в котором после запятой 40 битов. А что если ваш алгоритм не смог обеспечить 40 битов точности и достигает, скажем, только 30 битов? Тогда он выдаст другое число:

1,00

1

000000000000000000000000000

Не подозревая, что на 40-й позиции (которую алгоритм рассчитать не в состоянии) будет заветная единичка, вы округлите это число книзу и получите 1,00, что неправильно. Вы неправильно округлили последний бит — в этом и состоит предмет нашего обсуждения. Из сказанного выходит, что для того чтобы получить всего лишь 2-й бит правильным, вам придётся вычислять функцию до 40 битов! Вот это да! А если «паровоз» из нулей окажется ещё длиннее? Вот об этом и поговорим в следующем разделе.

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

Суть проблемы округления последнего значащего бита

Проблема проявляется по двум причинам. Первая — намеренный отказ от трудоёмкого вычисления в пользу скорости. В этом случае лишь бы заданная точностью соблюдалась, а какие там будут биты в ответе — дело второстепенное. Вторая причина — Table Maker’s Dilemma, которая является основным предметом нашего разговора. Рассмотрим обе причины подробнее.

Первая причина

Вы, конечно, понимаете, что вычисление трансцендентных функций реализовано какими-то приближёнными методами, например, методом аппроксимирующих полиномов или даже (редко) разложением в ряд. Чтобы вычисления происходили как можно быстрее, разработчики соглашаются выполнить как можно меньшее количество итераций численного метода (или взять полином как можно меньшей степени), лишь бы алгоритм допускал погрешность не превосходящую половину ценности последнего бита мантиссы. В литературе это пишется как 0.5ulp (ulp = unit in the last place).

Например, если речь идёт о числе x типа float на интервале (0,5; 1) величина ulp = 2^-23. На интервале (1;2) ulp = 2^-22. Иными словами, если x находится на интервале (0;1) то 2^x будет на интервале (1,2), и чтобы обеспечить точность 0.5ulp, нужно, грубо говоря, выбрать EPS = 2^-23 (так мы будем обозначать константу «эпсилон», в простонародье именуемую «погрешность», или «точность», кому как нравится, не придирайтесь, пожалуйста).

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

Для тех кто не понимает, приведу пример в десятичной системе счисления. Перед вами два числа: 1,999999 и 2,0. Допустим, что первое — это то, что получил программист, а второе — это эталон, что должно было бы получиться, будь у нас безграничные возможности. Разница между ними всего лишь одна миллионная, то есть ответ рассчитан с погрешностью EPS=10^-6. Однако ни одной правильной цифры в этом ответе нет. Плохо ли это? Нет, с точки зрения прикладной программы это фиолетово, программист округлит ответ, скажем, до двух знаков после запятой и получит 2,00 (например, речь шла о валюте, $2,00), ему больше и не надо, а то, что он заложил в свою программу EPS=10^-6, то молодец, взял запас на погрешность промежуточных вычислений и решил задачу правильно.

Иными словами, не путайте: точность и число правильных битов (или цифр) — это разные вещи. Кому нужна точность (это почти 100% программистов), тех обсуждаемая проблема совершенно не касается. Кому нужно чтобы битовая последовательность совпадала с корректно округлённым эталоном, того эта проблема волнует очень сильно, например, разработчиков библиотек элементарных функций. Тем не менее, для общего развития об этом полезно знать всем.

Напомню, это было первое направление проблемы: последние биты ответа могут быть неправильными потому, что это намеренное решение. Главное — оставить точность 0.5ulp (или выше). Поэтому численный алгоритм подбирается только из этого условия, лишь бы он работал предельно быстро. При этом Стандарт разрешает реализацию элементарных функций без корректного округления последнего бита. Цитирую [1, раздел 12.1] (англ.):

The 1985 version of the IEEE 754 Standard for Floating-Point Arithmetic did not specify anything concerning the elementary function. This was because it has been believed for years that correctly rounded functions would be much too slow at least for some input arguments. The situation changed since then and the 2008 version of the standard recommends (yet does not require) that some functions be correctly rounded.

Далее перечисляются эти функции, которые рекомендуется, но не требуется округлять корректно:

Вторая причина

Наконец-то мы перешли к теме разговора: Table Maker’s Dilemma (сокращённо TMD). Её название я не смог адекватно перевести на русский язык, оно было введено Уильямом Кэхэном (отцом-основателем IEEE-754) в статье [2]. Возможно, если вы прочитаете статью, то поймёте, почему название именно такое. Если кратко, то суть дилеммы в том, нам нужно получить абсолютно точное округление функции z=f(x), как если бы у нас в распоряжении была бесконечная битовая запись идеально посчитанного результата z. Но всем ясно, что бесконечную последовательность мы не можем получить. А сколько битов тогда взять? Выше я показал пример, когда нам нужно увидеть 40 битов результата, чтобы получить хотя бы 2 корректных бита после округления. А суть проблемы TMD в том, что мы заранее не знаем, до скольких битов рассчитать величину z, чтобы получить правильными столько битов после округления, сколько нам требуется. А что если их будет сто или тысяча? Мы не знаем заранее!

Например, как я сказал, для функции 2^x, для типа данных float, где дробная часть мантиссы имеет всего 23 бита, нам надо выполнить расчёт с точностью 2^-54, чтобы округление произошло правильно для всех без исключения возможных аргументов x. Эту оценку нетрудно получить полным перебором, но для большинства других функций, особенно для типа double или long double (ставь «класс», если знаешь что это), подобные оценки неизвестны.

Давайте уже разбираться с тем, почему так происходит. Я намеренно привёл самый первый пример в этой статье с типом данных float и просил вас его запомнить, потому что в этом типе всего 32 бита и проще будет на него смотреть, в остальных типах данных ситуация аналогична.

Мы начали с числа x = 0.00296957581304013729095458984375, это точнопредставимое число в типе данных float, то есть оно записано так, чтобы его конвертирование в двоичную систему float выполнялось без округления. Мы вычисляем 2^x, и если бы у нас был калькулятор с бесконечной точностью, то мы должны были бы получить (Чтобы вы могли проверять меня, расчёт выполнен в онлайн-системе WolframAlpha):

1,0020604729652405753669743044108123031635398201893943954577320057…

Переведём это число в двоичную систему, скажем, 64 бита будет достаточно:

1,00000000100001110000100

1

000000000000000000000000000001101111101

Бит округления (24-й бит после запятой) подчёркнут. Вопрос: куда округлять? Вверх или вниз? Ясное дело, вверх, вы знаете это, потому что видите достаточно битов и можете принять решение. Но посмотрите внимательно…

После бита округления у нас 29 нулей. Это означает, что мы очень-очень близко находимся к середине между двумя ближайшими числами и достаточно только чуть-чуть сдвинуться вниз, как направление округления изменится. Но вот вопрос: куда будет этот сдвиг? Численный алгоритм может последовательно шаг за шагом подбираться к точному значению с разных сторон и до тех пор, пока мы не пройдём все эти 29 нулей и не достигнем точности, которая превысит ценность самого последнего нуля в этом «паровозе», мы не будем знать направление округления. А вдруг в действительности правильный ответ должен быть таким:

1,00000000100001110000100

0

11111111111111111111111111111?

Тогда округление будет вниз.

Мы этого не знаем, пока наша точность не достигнет 54-го бита после запятой. Когда 54-й бит будет известен точно, мы будем знать точно, к какому из двух ближайших чисел мы на самом деле ближе. Подобные числа называются hardest-to-round-points [1, раздел 12.3] (критические для округления точки), а число 54 называется hardness-to-round (трудоёмкость округления) и в цитируемой книге обозначается буквой m.

Трудоёмкость округления (m) — это число битов, минимально необходимое для того, чтобы для всех без исключения аргументов некоторой функции f(x) и для заранее выбранного диапазона функция f(x) округлялась корректно до последнего бита (для разных режимов округления может быть разное значение m). Иными словами, для типа данных float и для аргумента x из диапазона (0;1) для режима округления к «ближайшему чётному» трудоёмкость округления m=54. Это значит что абсолютно для всех x из интервала (0;1) мы можем заложить в алгоритм одну и ту же точность ESP=2^-54, и все результаты будут округлены корректно до 23-х битов после двоичной запятой.

В действительности, некоторые алгоритмы способны обеспечить точный результат и на основе 53-х и даже 52-х битов, полный перебор это показывает, но теоретически, нужно именно 54. Если бы не было возможности провернуть полный перебор, мы бы не могли «схитрить» и сэкономить пару битов, как это сделал я в программе перебора, что была дана выше. Я взял полином степенью ниже, чем следовало, но он всё равно работает, просто потому что повезло.

Итак, независимо от режима округления, у нас возможны две ситуации: либо в области округления возникает «паровоз» из нулей, либо «паровоз» из единиц. Задача правильного алгоритма расчёта трансцендентной функции f(x) в том, чтобы уточнять значение этой функции до тех пор, пока точность не превысит ценность последнего бита этого «паровоза», и пока не станет точно понятно, что в результате последующих флуктуаций численного алгоритма расчёта f(x) нули не превратятся в единицы, или наоборот. Как только всё стабилизировалось, и алгоритм «вышел» на такую точность, которая находится за пределами «паровоза», тогда мы можем округлять так, как если бы у нас было бесконечное число битов. И это округление будет с корректным последним битом. Но как этого добиться?

«Костыли»

Как было сказано, основная проблема — заставить алгоритм преодолеть паровоз из нулей или единиц, которые идут сразу после бита округления. Когда паровоз преодолён и мы видим его целиком, тогда это эквивалентно тому, что эти нули или единицы уже рассчитаны точно, и мы уже точно знаем, в какую сторону теперь произойдёт округление. Но если мы не знаем длину паровоза, то как же тогда проектировать алгоритм?

Первый «костыль»

Читателю может показаться, что ответ очевиден: взять арифметику с бесконечной точностью и заложить заведомо избыточное число битов, а если будет мало, то заложить ещё и пересчитать. В общем-то правильно. Так и делается, когда скорость и ресурсы компьютера не играют особой роли. У этого подхода есть имя: «многоуровневая стратегия Зива» (Ziv’s multilevel strategy) [1, раздел 12.3]. Суть её предельно проста. Алгоритм должен поддерживать расчёты на нескольких уровнях: быстрый предварительный расчёт (в большинстве случаев он же оказывается финальным), более медленный, но более точный расчёт (спасает в большинстве критических случаев), ещё более медленный, но ещё более точный расчёт (когда совсем «худо» пришлось) и так далее.

В подавляющем большинстве случаев достаточно взять точность чуть выше чем 0.5ulp, но если появился «паровоз», то увеличиваем её. Пока «паровоз» сохраняется, наращиваем точность до тех пор, пока не будет строго понятно, что дальнейшие флуктуации численного метода на этот «паровоз» не повлияют. Так, скажем, в нашем случае если мы достигли ESP=2^-54, то на 54-й позиции появляется единица, которая как бы «охраняет» паровоз из нулей и гарантирует, что уже не произойдёт вычитания величины, больше либо равной 2^-53 и нули не превратятся в единицы, утащив за собой бит округления в ноль.

Это было научно-популярное изложение, всё то же самое с теоремой Зива (Ziv’s rounding test), где показано как быстро, одним действием проверять достигли ли мы желаемой точности, можно прочитать в [1, глава 12], либо в [3, раздел 10.5].

Проблема этого подхода очевидна. Нужно проектировать алгоритм вычисления каждой трансцендентной функции f(x) так, чтобы по ходу пьесы можно было увеличивать точность расчётов. Для программной реализации это ещё не так страшно, например, метод Ньютона позволяет, грубо говоря, удваивать число точных битов после запятой на каждой итерации. Можно удваивать до тех пор, пока не станет «достаточно», хотя это довольно трудоёмкий процесс, надо признаться и не везде метод Ньютона оправдан, потому что требует вычисления обратной функции f^-1(x), что в некоторых случаях может оказаться ничуть не проще вычисления самой f(x). Для аппаратной реализации «стратегия Зива» совершенно не годится. Алгоритм, зашитый в процессоре, должен выполнить ряд действий с уже предустановленным числом битов и это достаточно проблематично реализовать, если мы этого числа заранее не знаем. Взять запас? А сколько?

Вероятностный подход к решению проблемы [1, раздел 12.6] позволяет оценить величину m (напомню, это число битов, которого достаточно для корректного округления). Оказывается, что длина «паровоза» в вероятностном смысле чуть больше длины мантиссы числа. Таким образом, в большинстве случаев будет достаточно брать m чуть более чем вдвое больше величины мантиссы и только в очень редких случаях придётся брать ещё больше. Цитирую авторов указанной работы: «we deduce that in practice, m must be slightly greater than 2p» (у них p — длина мантиссы вместе с целой частью, то есть p=24 для float). Далее в тексте они показывают, что вероятность ошибки при такой стратегии близка к нулю, но всё-таки положительна, и подтверждают это экспериментами.

Тем не менее, всё равно остаются случаи, когда величину m приходится брать ещё больше, и худший случай заранее неизвестен. Теоретические оценки худшего случая существуют [1, раздел 12.7.2], но они дают немыслимые миллионы битов, это никуда не годится. Вот таблица из цитируемой работы (это для функции exp(x) на интервале от -ln(2) до ln(2)):

Второй «костыль»

На практике величина m не будет такой ужасно-большой. И чтобы определить худший случай применяется второй «костыль», который называется «исчерпывающий предподсчёт». Для типа данных float (32 бита), если функция f имеет один аргумент (x), то мы можем запросто «прогнать» все возможные значения x. Проблема возникнет только с функциями, у которых больше одного аргумента (среди них pow(x, y)), для них ничего такого придумать не удалось. Проверив все возможные значения x, мы вычислим свою константу m для каждой функции f(x) и для каждого режима округления. Затем алгоритмы расчёта, которые нужно реализовать аппаратно, проектируются так, чтобы обеспечить точность 2^-m. Тогда округление f(x) будет гарантированно-корректным во всех случаях.

Для типа double (64 бита) простой перебор практически невозможен. Однако ведь перебирают! Но как? Ответ даётся в [4]. Расскажу о нём очень кратко.

Область определения функции f(x) разбивается на очень маленькие сегменты так, чтобы внутри каждого сегмента можно было заменить f(x) на линейную функцию вида b-ax (коэффициенты a и b, конечно же, разные для разных сегментов). Размер этих сегментов вычисляется аналитически, чтобы такая линейная функция действительно была бы почти неотличима от исходной в каждом сегменте.

Далее после некоторых операций масштабирования и сдвига мы приходим к следующей задаче: может ли прямая линия b-ax пройти «достаточно близко» от целочисленной точки?

Оказывается, что можно относительно просто дать ответ в форме «да» или «нет». То есть «да» — если потенциально опасные точки будут близки к прямой линии, и «нет» — если ни одна такая точка в принципе не может подойти к линии близко. Прелесть этого метода в том, что ответ «нет» на практике получается в подавляющем большинстве случаев, а ответ «да», который получается редко, вынуждает пройтись по сегменту полным перебором, чтобы определить какие конкретно точки оказались критическими.

Тем не менее, перебор аргументов функции f(x) сокращается во много-много раз и делает возможным обнаруживать критические точки для чисел типа double (binary64) и long double (80 битов!). Делается это на суперкомпьютерах и, конечно же, видеокартах… в свободное от майнинга время. Тем не менее, что делать с типом данных binary128 пока никто не знает. Напомню, что дробная часть мантиссы таких чисел занимает 112 битов. Поэтому в иностранной литературе по данному поводу пока можно отыскать только полуфилософские рассуждения, начинающиеся с «we hope…» («мы надеемся…»).

Подробности метода, который позволяет быстро определить прохождение линии вблизи от целочисленных точек, здесь излагать неуместно. Желающим познать процесс более тщательно рекомендую смотреть в сторону проблемы поиска расстояния между прямой и Z², например, в статье [5]. В ней описан усовершенствованный алгоритм, который по ходу построения напоминает знаменитый алгоритм Евклида для нахождения наибольшего общего делителя. Приведу одну и ту же картинку из [4] и [5], где изображена дальнейшая трансформация задачи:

Существуют обширные таблицы, содержащие худшие случаи округления на разных интервалах для каждой трансцендентной функции. Они есть в [1 раздел 12.8.4] и в [3, раздел 10.5.3.2], а также в отдельных статьях, например, в [6].

Я приведу несколько примеров, взяв случайные строки из таких таблиц. Подчёркиваю, это не худшие случаи для всех x, а только для каких-то небольших интервалов, смотрите первоисточник, если заинтересовались.

Как читать таблицу? Величина x указана в шестнадцатеричной нотации числа с плавающей запятой double. Сначала, как положено, идёт лидирующая единица, затем 52 бита дробной части мантиссы и буква P. Эта буква означает «умножить на два в степени» далее идёт степень. Например, P-23 означает, что указанную мантиссу нужно умножить на 2^-23.

Далее представьте, что вычисляется функция f(x) с бесконечной точностью и у неё отрезают (без округления!) первые 53 бита. Именно эти 53 бита (один из них — до запятой) указываются в столбце f(x). Последующие биты указаны в последнем столбце. Знак «степени» у битовой последовательности в последнем столбце означает число повторений бита, то есть, например, 10⁵³1011 означает, что сначала идёт бит, равный 1, затем 53 нуля и далее 1011. Потом троеточие, которое означает, что остальные биты нам, в общем-то, и не нужны вовсе.

Далее дело техники — мы знаем худшие случаи для каждого интервала отдельно взятой функции и можем для этого интервала подобрать такую аппроксимацию, чтобы она своей точностью накрывала худший случай. Всего лишь годы расчётов на суперкомпьютерах дают возможность создавать быстрые и точные аппаратные реализации элементарных функций. Дело за малым: осталось научить хотя бы разработчиков компиляторов пользоваться этими таблицами.

Зачем это нужно?

Прекрасный вопрос! Ведь я выше неоднократно высказался о том, что почти 100% программистам не нужно знать элементарную функцию с точностью до корректно-округлённого последнего бита (часто им и половина битов не нужна), зачем учёные гоняют суперкомпьютеры и составляют таблицы ради решения «бесполезной» задачи?

Во-первых, задача носит фундаментальный характер. Интерес представляет скорее не то, чтобы получить точное округление ради точного округления, а чтобы в принципе понять, как можно было бы решить эту интересную задачу, какие тайны вычислительной математики откроет нам её решение? Как можно было бы эти тайны использовать в других задачах? Фундаментальные науки — они такие, можно десятки лет заниматься какой-то «ерундой», а потом через сто лет благодаря этой «ерунде» происходит научный провыв в какой-то другой области.

Во-вторых, вопрос переносимости кода. Если функция может позволить себе как угодно обращаться с последними битами результата, то значит на различных платформах и на различных компиляторах могут получаться чуть-чуть различные результаты (пусть даже в пределах заданной погрешности). В каких-то случаях это не важно, а в каких-то может быть существенным, особенно когда в программе есть ошибка, которая проявляется на одной, но не проявляется на другой платформе как раз именно из-за различных битов результата. Но зачем я расписываю вам прекрасно известную головную боль, связанную с различным поведением программы? Вы и без меня всё это знаете. Было бы здорово иметь математическую систему, которая работает абсолютно одинаково на всех платформах, чем бы её не компилировали. Вот для этого и нужно корректно округлять последний бит.

Список источников

[1] Jean-Michel Muller, “Elementary Functions: Algorithms and Implementation”, 2016

[2] William Kahan, “A Logarithm Too Clever by Half”, 2004

[3] Jean-Michel Muller, “Handbook of floating-point arithmetic”, 2018

[4] Vincent Lefèvre, Jean-Michel Muller, “Toward Correctly Rounded Transcendentals”, IEEE TRANSACTIONS ON COMPUTERS, VOL. 47, NO. 11, NOVEMBER 1998. pp. 1235-1243

[5] Vincent Lefèvre. “New Results on the Distance Between a Segment and Z²”. Application to the Exact Rounding. 17th IEEE Symposium on Computer Arithmetic — Arith’17, Jun 2005, Cape Cod, MA,
United States. pp.68-75

[6] Vincent Lefèvre, Jean-Michel Muller, “Worst Cases for Correct Rounding of the Elementary Functions in Double Precision”, Rapport de recherche (INSTITUT NA TIONAL DE RECHERCHE EN INFORMA TIQUE ET EN AUTOMA TIQUE) n˚4044 — Novembre 2000 — 19 pages.