Table-Maker's Dilemma, или почему почти все трансцендентные элементарные функции округляются неправильно

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



Друзья, для вашего удобства статья также доступна в формате презентации на видео (почти 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 ("%.8f\n", z);  // Вывод результата с округлением до 8 знаков после точки.
  z = powf(2.0f, x);  // z = 2**x другим способом
  printf ("%.8f\n", z);  // Такой же вывод.
  return 0;
}

Число x намеренно записано с таким количеством значащих цифр, чтобы оно было точнопредставимым в типе float, то есть чтобы компилятор преобразовал его в бинарный код без округления. Ведь вы прекрасно знаете, что некоторые компиляторы не умеют округлять без ошибок (если не знаете, укажите в комментариях, я напишу отдельную статью с примерами). Далее в программе нам нужно вычислить 2x, но давайте сделаем это двумя способами: функция 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) с некоторыми другими элементарными функциями. Привожу результаты в таблице. У правильных ответов (когда они есть) стоят символы «*» на конце.

Функция Аргумент MS C++ Intel C++ GCC Clang
log10(x) 2.60575359533670695e129 0x40602D4F53729E44 0x40602D4F53729E45* 0x40602D4F53729E44 0x40602D4F53729E44
expm1(x) -1.31267823646623444e-7 0xBE819E53E96DFFA9* 0xBE819E53E96DFFA8 0xBE819E53E96DFFA8 0xBE819E53E96DFFA8
pow(10.0, x) 3.326929759608827789e-15 0x3FF0000000000022 0x3FF0000000000022 0x3FF0000000000022 0x3FF0000000000022
logp1(x) -1.3969831951387235e-9 0xBE17FFFF4017FCFF* 0xBE17FFFF4017FCFE 0xBE17FFFF4017FCFE 0xBE17FFFF4017FCFE


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

$2^x = 2^{[x]}\cdot2^{\{x\}}.$



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

MS C++ Intel C++ GCC Clang
1910726 (0,97%) 90231 (0,05%) 0 0


Из программы ниже понятно, что число проверяемых аргументов x составило 197612997 штук. Получается, что, например, Microsoft С++ неверно вычисляет функцию 2x для почти одного процента из них. Не радуйтесь, уважаемые любители 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%08X\n", 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 можно проверить полным перебором! Так ли сложно это было сделать? Совсем не сложно, и я уже показал пример.

В коде нашего перебора есть «самописная» функция правильного вычисления 2x с помощью аппроксимирующего полинома 10-й степени, и написана она была за несколько минут, потому что такие полиномы выводятся автоматически, например, в системе компьютерной алгебры Maple. Достаточно задать условие, чтобы полином обеспечивал 54 бита точности (именно для этой функции, 2x). Почему 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,0010000000000000000000000000000000000001

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

1,001000000000000000000000000000

Не подозревая, что на 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) то 2x будет на интервале (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, чтобы получить правильными столько битов после округления, сколько нам требуется. А что если их будет сто или тысяча? Мы не знаем заранее!

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

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

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

1,0020604729652405753669743044108123031635398201893943954577320057…

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

1,000000001000011100001001000000000000000000000000000001101111101

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

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

1,00000000100001110000100011111111111111111111111111111?

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

Мы этого не знаем, пока наша точность не достигнет 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)):

p m
24 (binary32) 1865828
53 (binary64) 6017142
113 (binary128) 17570144


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


На практике величина 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...» («мы надеемся...»).

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

image

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

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

Функция x f(x) (обрезанное) 53-й бит и последующие
log2(x) 1.B4EBE40C95A01P0 1.8ADEAC981E00DP-1 10531011...
cosh(x) 1.7FFFFFFFFFFF7P-23 1.0000000000047P0 11890010...
ln(1+x) 1.8000000000003P-50 1.7FFFFFFFFFFFEP-50 10991000...


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

Далее представьте, что вычисляется функция f(x) с бесконечной точностью и у неё отрезают (без округления!) первые 53 бита. Именно эти 53 бита (один из них — до запятой) указываются в столбце f(x). Последующие биты указаны в последнем столбце. Знак «степени» у битовой последовательности в последнем столбце означает число повторений бита, то есть, например, 10531011 означает, что сначала идёт бит, равный 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 Z2”. 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.

Средняя зарплата в IT

111 000 ₽/мес.
Средняя зарплата по всем IT-специализациям на основании 7 247 анкет, за 2-ое пол. 2020 года Узнать свою зарплату
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее
Реклама

Комментарии 83

    0

    Да, с FMA тоже бывают приколы. Суть в том, что это (float)(a * b + c), а не (float)((float)(a * b) + c). То есть не происходит промежуточного округления после умножения. Поэтому результат может отличаться.

      +1
      Совершенно верно, арифметика с плавающей запятой обладает целым рядом особенностей, поэтому такие вот вспомогательные функции типа FMA могут сделать жизнь лучше. Однако есть обратная сторона. Иногда (не всегда) FMA может работать заметно медленнее, чем отдельное умножение и сложение. Примером тому выступает код полного перебора (скрытый текст в статье), где как раз вычисляется полином в точке. Там если написать код через FMA, на моём компьютере становится в разы медленнее.
        0

        Это, скорее, связано с зависимостью операций по данным. Throughput у FMA, может, и хороший, но latency всё перечёркивает. Решение проблемы: одновременное вычисление значений для x и для y (а лучше сразу для векторов — AVX для кого придумали?) с чередованием операций.

      0
      Интересно, будет ли приводить к накоплению ошибок округление не к последнему четному биту (речь идет только о вычислении округления иррациональных чисел, если что), а строго вверх…

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

      Ровно так же последовательность 0.1101111111 может быть просто приближением снизу к точному значению 0.111(0), но тут мы всегда округляем вниз, считая (вполне обоснованно) что минимум один из неизвестных битов равен нулю. Почему бы не делать то же — округлять вверх хвост из нулей, считая, что где-то дальше за ними появится неизвестная сейчас единичка?

      Или я где-то ошибаюсь?
        +2

        Благодарю за вопрос! Вы ошибаетесь в том, что для округления вверх может быть другой пример, когда у нас, например, число 1.11000000001, и вы округлите вверх, получая 1.111, но вы заранее не знаете, этот паровоз из нулей — это уже правильный паровоз, или при последующих итерациях алгоритма на самом деле будет 1.1011111111? И тогда округление вверх должно дать 1.110. Опять ошибка в последнем бите. Короче говоря, для любого варианта округления есть своя ситуация, когда мы находимся очень близко к границе между двух направлений округления, и чтобы точно понять, с какой стороны от этой границы мы находимся, нужно чтобы паровоз из нулей или единиц точно закончился.

          0

          Ещё добавлю к цитате:


          то дальше в ней либо все (до бесконечности) нули, либо где-то встретится единичка и результат будет ближе к верхнему значению округления

          Нет, есть ещё третий вариант: нули превратятся в единички по ходу дальнейших расчётов. Потому что эти нули за пределами 0.5ulp. И вот если это произойдёт, тогда округление вверх даст другой результат.


          UPD: Даже могу ещё точнее пример привести. Допустим, дробная часть мантиссы занимает два бита. Берём два числа, которые отличаются друг от друга меньше чем на 0.5ulp:
          0.11000001
          0.10111111
          Наши 0.5ulp равны 0.125 (ценность последнего — второго — бита равна 0.25). Разница между числами значительно меньше. Но если округлять вверх оба, то будут разные результаты. Поэтому чтобы понять, куда же округлять на самом деле, нам нужно добраться в точности до 8-го бита. Аналогичные примеры можно придумать для любых вариантов округления.

            –1
            Если я правильно понимаю алгебру, то 1.9(9) = 2.0(0), строго равно, математически строго. То есть, 1.99999- это абсолютно такие же _точные_ пять цифр после запятой, что и 2.00000, хотя в этих двух числах и нету ни одной одинаковой цифры вообще. С учетом этого не совсем понятно, о каком «точном» последнем бите может идти речь, если у одного и того же числа в двоичной записи тоже может быть два разных представления: b1.(0) = b0.(1), и как же тогда округлять значение b0.11111111111111111?
              +3

              Вы правильно помните алгебру, но неправильно поняли задачу. У нас не бесконечно много девяток после запятой, а всего шесть. И это именно два разных числа: 1,999999 и 2,0.


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

          +3
          What every computer scientist should know about floating point.jpghtml

          Спасибо за напоминание, полезно.

          ps:
          > #define FAU(x) (*(unsigned int*)(&x))
          > #define DAU(x) (*(unsigned long long*)(&x))

          Не думаю, что это повлияло на результат, но вообще-то в С++ так делать нельзя, это UB (только в С можно). Доступ к памяти некастуемого типа разрешен только через char* (::std::memcpy() в частности или std::bit_cast (since C++20)).
            +1

            Я с вами полностью согласен, это UB, но поясню вам свою позицию. Дело в том, что в моей работе код пишется только под одну архитектуру и если нужно написать под другую, пишется другой код, такое правило. Причём пишется заново. Поэтому те вещи, которые в обычной жизни считаются UB, у меня считаются нормой и всегда работают в точности так, как я себе это представляю. То же касается битовых трюков, которые также нельзя было бы использовать. Ещё я пишу на этакой смеси C и C++, это тоже обусловлено постоянной работой на низком уровне. В общем, не судите строго, благодарю за замечание!


            Здесь мне важно, чтобы читатели поняли суть, а как они это реализуют для себя — это уже пусть смотрят сами.

              +4
              Тут главное помнить, что это всё работает лишь благодаря человеколюбивой милости компиляторописателей)) Чтобы вот кто-то коммитился «мы гарантируем, что это всегда во всех версиях наших компиляторов будет работать именно так» — я такого не видел…
            0
            Эта болезнь не лечится.
            Есть МК51, имеет 100% совпадение с проверкой на бумаге, но иногда не совпадает с GCC, и для этого есть простое объяснение. Этот древний мамонт считает в десятичной системе!!!
            В двоичной системе для округления достаточно отбросить лишнее — без анализа содержимого. Вы уже выбрали необходимую вам точность, остальное туда просто не влезет.
            В десятичной системе так не получится — один младший бит может влиять на старший десятичный разряд, даже без округления (разговор про одинарную и двойную точность).

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

              Да ладно? Такое возможно только в "пограничных" случаях, и такие пограничные случаи в десятичной системе встречаются в 5 раз реже чем в двоичной!

                –1
                Для классической десятичной системы, где каждая десятичная цифра занимает 4 фиксированных бита — «пограничных случаев» не существует в принципе.
                А вот для одинарной и двойной точности в двоичном варианте — огромное количество. Тут я лично споткнулся, когда писал свой вариант printf. Либо получается быстро с уникальными ошибками перевода, либо точно- но ещё медленные чем стандартная функция.
                Кстати, printf для устранения ошибки перевода — почти двух кратно повышает точность. И его перевод считается образцовым, но не идеально точным. Потому как обратное преобразование достаточно часто даёт ошибку.
                  +1

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


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

                    0
                    При чём тут вообще битовое представление? Преобразование между системами счисления — это далеко не единственное место, где может образовываться погрешность.

                    При том что это очень эффектный способ демонстрации погрешности. Именно демонстрации!!! Когда младший бит влияет на старшую десятичную цифру.
                      0

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

                        –1
                        Одинарная двоичная точность 0,1 + 0,2 != 0,3. Число 0,2 представлено как 0x3e4ccccd, последний младший бит влияет на старший знак. Это число всегда округляется, потому как не может быть представлено в десятичном виде без ошибки. Это и есть эффект демонстрации.
                        Десятичная арифметика 0,1 + 0,2 = 0,3. Число 0,2 представлено как 2e-1 в BCD формате, после двойки все остальные числа равны нулю. Как учили в школе.

                        Нужно просто перестать натягивать сову на глобус. Методы округления чисел в двоичном и десятичном варианте — должны выполняться в пределах формата числа, без принудительного перевода формата.
                          0

                          Простите, но кто тут предлагает округлять с принудительным переводом формата?

                            –1
                            Вы похоже не поняли суть статьи, мои высказывания, и общие правила округления. Нельзя подгонять результат вычислений в одинарной или двойной точности под красивый вывод данных через printf.
                              0

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

              0
              Интересно, что если округлить число с КДПВД до 9ти знаков в Wolfram Mathematica, то получается что-то среднее = 1.00206047.
              Код
              image
                +1

                У меня округление выполняется в типе данных float, то есть имеется только 23 бита после запятой. В Wolfram используется иная внутренняя структура хранения чисел, там округление будет более точным, потому что двоичных битов больше.

                  0
                  > двоичных битов

                  Вода мокрая, масло масляное, биты двоичные. Хотя они двоичные по определению.
                    0

                    Благодарю за поправку :)

                +3
                Если я правильно понял автора, то речь идёт о том, что решение должно быть индивидуальным для каждого из алгоритмов вычисления трансценентных чисел, но проблема так же заключается в сложности верификации результата. То есть, по сути, для каких-то алгоритмов могут сработать очень простые эврестики (например, мы знаем на каком шаге с какой стороны мы прближаемся к числу с большей или с меньшей и в противоположную сторону и можем делать округление), но проверить для 128 бит не представляется возможным. Тогда ждём статью именно о способах проверки :-)
                  +2

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


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

                  –2

                  Еще полезно знать, что для чисел с плавающей точкой данные утверждения бывают истины:


                  1. a + b — b != a
                  2. a * b / b != a
                  3. (a + b) + c != a + (b + c)
                    И т.д.
                    +1

                    Благодарю, но это далеко не все особенности, о которых полезно знать. У меня на канале есть по этому поводу информация, если интересно. Да и вообще в интернете много про это пишут.


                    Есть ещё полезная статья David Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic, 1991

                    0
                    Насколько принципиальна разница между 51 и 52 битами точности? Лично мне кажется более важной проблема накопления погрешности при конкретной последовательности вычислений, и она вовсе не обязательно будет коррелировать с точностью вычисления последнего бита отдельно взятой функции.
                      0
                      Ну как-то раз нас спросили почему в 12-ом знаке после запятой в балансе абонента не 0. Пришлось срочно переделывать арифметику.
                        +3
                        Так это скорее вопрос правильного выбора типа данных для денежных расчётов — вряд ли у вас баланс считался с использованием трансцендентных функций (наверное).
                          0

                          Денежные расчёты принято делать в десятичной арифметике с плавающей запятой. Если будет интересно, объясню почему.

                            +1
                            Был уверен, что с фиксированной. Но интересно, расскажите.
                              0

                              Ну можно и с фиксированной, однако есть же Стандарт IEEE-754, который описывает типы данных Decimal64, например как раз для этого дела. Причина проста. В финансовой математике действуют законы государства, прописанные в десятичной системе счисления. Например, если написано по закону, что, скажем, налог составляет 12.4 рубля, значит должно быть ровно столько, а когда мы воспользуемся типом данных binary64, будет (ближайшее число):


                              12.4000000000000003552713678800500929355621337890625


                              А это уже прямое нарушение закона. Когда выполняется множество действий с такими вот числами (видели, наверное, эти ипотечные формулы с возведением в 365-ю степень и прочее), тогда погрешность будет нарастать, и нарушение закона будет нарастать ещё больше. Как-то так, если кратко.

                              0

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

                                0

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

                                  0

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


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

                                    0

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

                          +1

                          Благодарю за комментарий! Вы говорите о совершенно другой проблеме. Вы говорите о том, насколько точным будет итоговый результат после множества операций с числами с плавающей запятой. Эти вещи просчитываются аналитически и описаны в книге [3] (там есть оценки на погрешность при сложениях, умножениях и т. д.). Я же говорю о том, что вычислять значения элементарных функций нужно точно до последнего бита и причина необходимости делать это описана в последних абзацах.


                          А вот для обычных прикладных программистов эта проблема (описанная в статье) не важна. Им без разницы, что там с битами, им важно, чтобы конечный ответ был бы с определённой точностью. Между 51 и 52 битами точности для их целей может не быть разницы вообще.

                            0
                            им важно, чтобы конечный ответ был бы с определённой точностью
                            А разве бывает по-другому? У нас так и так точное значение недоступно, так и так значения округляются — причём в логарифмическом, а не линейном масштабе.
                              0

                              Опять мы говорим о разном. Я же приводил пример с числами 1,999999 и 2,0. Прикладному программисту не важно, что у него получилось "почти два, но не совсем". Он округлит — и дело с концом. Разработчику библиотеки ВАЖНО, чтобы результат был корректно округлён до последнего бита, ему ВАЖНО получить именно 2,0 (в данном примере).

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

                                  Тогда я не понимаю, чего вы хотите от меня, если всё что я хотел сказать, написано в статье. Я описал открытую научную проблему над которой учёные работают активно минимум 15 последних лет. Если вам показалось, что они зря тратят время, то лучше сказать об этом им, но не мне :)

                                    +1
                                    От вас — ничего, мне статья понравилась. Но комментарии с дополнениями статье не повредят — как минимум, добавляют просмотров и дают повод для дискуссий.
                                      0

                                      В таком случае благодарю за полезные комментарии!

                            +2
                            Насколько принципиальна разница между 51 и 52 битами точности?

                            Дело касается переносимости результатов вычислений.

                              0
                              Да, с этим не поспоришь. Тогда очевидно должен быть стандарт.
                                0

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

                                  +1
                                  Но можно застандартизировать алгоритм для вычисления конкретной функции.
                                    –1

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

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

                            Сверился с википедией на которую вы даёте ссылку, что бы проверить ваше «неарест тиес то евен», так вот, ссылка по википедии под таким названием ведёт на «Round half to even», что является банковским округлением, в котором чётная цифра смотрится не после, а до округляемого разряда и никаких сложностей не представляет. И мне понятно о чём вы говорите: если для математического округления ">=" это единый случай и полностью покрыается одним разрядом, то для банковского, где разделяют ">" и "=", и в случае равенства делаются лишние телодвижения, то для трансцендентных чисел само понятия равно в дробной части не имеет смысла, потому что числа не кончаются и после любого количестава нулей может пойти «нормальная» серия. И всё что делается, это долбёжка чисел до упора в поисках любой еденички, что бы мифическое равно отбросить и перейти к чёткому ">". Соотвественно, банковское округление для транцендентых чисел бессмысленно, потому что равенства не бывает. Но использование для банкинга трансцендентных чисел, это вообще говоря, сомнительное занятие, тем более, что банковского округление это частный случай и не решает оно проблему полностью, а является лишь костылями подходящими лишь для определённых операций… По сути, проблема решается тем, что если и появляются где-то транцендентные числа, то для них надо использовать математическое округление и уже только дальше банковское/случайное, если есть какие-то статистически значимые операции.
                              0

                              А разве я говорил где-то, что смотрится чётная цифра ПОСЛЕ округляемой? Вы что-то не то говорите. Цифра, которая в моих примерах подчёркивается — это так называемый round-bit (назовём его b), по его значению проверяем направление округления, есть ещё sticky-bit (s), который является логическим или над всеми остальными битами, идущими после round-bit. Далее смотрится та самая цифра, которая остаётся последней после округления это least-bit (l), именно она должна стать чётной в случае, если b=1 и s=0. Таким образом, имеем три бита lrs информации по которым однозначно определяем округление. Но я намеренно не хотел как-то путать читателей этими тремя битами, а постарался объяснить на словах.


                              Для трансцендентных функций никогда не бывает s=0, вы правы, там всегда будет какая-то единичка дальше. Проблема лишь в том, что мы можем не знать b=0 или b=1. Алгоритм может плясать от одного к другому и обратно, пока не стабилизируется. Думаю, вы правильно меня поняли, но мне просто слово банковское непонятно, не знаком с такой терминологией.

                                +1
                                Округление до ближайшего чётного (в английском языке известно под названием англ. banker's rounding — «округление банкира») — округление для этого случая происходит к ближайшему чётному числу, то есть 2,5 → 2; 3,5 → 4.
                                  0

                                  Спасибо, буду знать, это это банковское округление :) А бывает ещё nearest-ties-to-away, это когда величина посередине округляется дальше от нуля. То есть это то, что у нас в математике принято: 12.5→13; 13.5→14. Именно этот режим прописан в Стандарте для десятичных чисел с плавающей запятой, и я думал, что банковская система должна работать по этому варианту округления, раз она работает с десятичными числами.

                                    0
                                    Округление работает безотносительно системы исчисления, просто равность расстояний проявляется только в «чётных» и только в них имеет смысл играться с тем, куда отнести равноудалённую середину. Если середина не равноудалена, то и вопросов никаких нет.
                                      0

                                      Как это нет вопросов? :) Я же столько написал о том, что мы НЕ ЗНАЕМ, с какой стороны мы от середины :)

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

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


                                          1,01011111111111111111111111111111111111111 → 1,01
                                          1,01100000000000000000000000000000000000001 → 1,10


                                          Между ними разница 2 в степени -40. Оба эти числа были получены в ходе вычислений пошаговым алгоритмом (одно на 100-м шаге, другое — на 101-м), которое из них по-вашему нужно взять в качестве достаточно-точного для того, чтобы сформировать правильное округление?


                                          Как видим, 40 битов точности недостаточно, чтобы понять, с какой стороны мы от серединки. А сколько достаточно? Не знаем! Вот в этом и состоит суть проблемы TMD.

                                            –1
                                            100 ша ни о чём не говорит, для любых распространённых алгоритмов есть гарантии точного знака после запятой. в данном случае достаточно быть уверенным в 3 знаке, он у у вас разный, соотвественно точность не достаточна. Вы сами приводили как пример, вычисление квадратного корня, где количество ТОЧНЫХ знаков удваивается после каждой итерации. Точные — означает, что они не меняются при дальнейших вычислениях.
                                              0

                                              Что-то я не понимаю. Какие такие гарантии точного знака? Не слышал. Метод Ньютона даёт, грубо говоря, удвоение знаков, но это жаргон такой, так просто принято говорить. В действительности он даёт возведение точности в квадрат по сравнению с предыдущим шагом. То есть, скажем, была точность 2 в минус 2, стала 2 в минус 4, потом 2 в минус 16.


                                              Мне кажется, вы не понимаете суть. В моём примере чтобы получить точными ВСЕГО ДВА БИТА, неизвестно как долго ещё считать! Уже получили 41 бит после запятой, а всё ещё не знаем, какими точно будут ВСЕГО ДВА бита. А сколько ещё считать? Никто не знает! Это это незнание характерно для ВСЕХ алгоритмов всех трансцендентных функций, никаких гарантий никто никогда не давал.


                                              Вот вы сами пишите: "соотвественно точность не достаточна".
                                              Вопрос в том КОГДА она будет достаточной? Этого никто не знает.

                                                –1
                                                Дело в том, что то, о чём вы говорите совсем не очевидно. Примеры округлений показываются на точно известных разрядах. Например:
                                                0,12345X, округляя до 0,1235 мы внесём погрешность в 0,00005 максимум и разряд X отличный от нуля, может только уменьшить нашу погрешность. Для итогового же числа 0,1235, было бы уместно записать погрешность в +-0,00005, хотя в реальности она была меньше, просто это покажет, что число округлённое, и 4 знак после запятой не является точным и могло быть как 0,1234x так и 0,1235x.
                                                Собственно говоря, с какой бы стороны мы не приближались к числу, математическое округление внесёт наименьшую погрешность величиной в половину округляемого разряда. В вашем примере, округлив до 2-39 мы полчим одинаковые числа. Понимаю, о чём вы говорите, о том, что сколь угодно малое число может превращать 1цу в набор девяток. Но если мы берём еденицу и +- 0,005 погрешность, то это будут числа от 0,995 до 1,004 и эти числа всегда будут округляться до 1цы, фаткически, пока погрешность не вырастет до 0,5 (а она будет только падать 0,9995 1,0004 и т.п) мы можем быть уверены в этой еденице, никакие младшие разряды 0.5 не переплюнут по величине. Поэтому при повышении точности старшие разряды, после округления в младшем разряде, не будут никак изменяться. Но при этом, если вы не уверены в фактических значениях разрядов вы не можете колдовать с 0.5 значением для десятичной, 0.4 для восмеричной, 0.3 для шестиричной, 0.2 для четырёхричной и 0.1 для двоичной, как в банковском округлении, потому что фаткически этот разряд у вас не настолько точный в любой позиции.

                                                  +1
                                                  0,12345X, округляя до 0,1235 мы внесём погрешность в 0,00005 максимум и разряд X отличный от нуля, может только уменьшить нашу погрешность.

                                                  Ну да, если у нас есть число 0,12345... с гарантированно верными 5 цифрами — оно может быть округлено только до 0,1235. Вопрос в том, как получить эти самые 5 верных цифр; всегда может оказаться что на самом деле у нас 0,123449876... и округлять надо до 0,1234.

                                                0

                                                Да нету такой гарантии! Вот я взял число 0.0199999994**2, и попробовал найти квадратный корень. Первые 10 итераций получились такими:


                                                0: 0.500199999988
                                                1: 0.2504998400339936
                                                2: 0.12604832367255547
                                                3: 0.06461085479793548
                                                4: 0.03540088231324618
                                                5: 0.023350017521141124
                                                6: 0.020240312312009873
                                                7: 0.020001426015725795
                                                8: 0.01999999945087718
                                                9: 0.0199999994

                                                Третий разряд после запятой "стабилизировался" только на 9й итерации! Ну и где тут удвоение числа верных разрядов-то?


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

                                                  +1

                                                  Откройте, пожалуйста, хотя бы википедию и посмотрите там раздел Proof of quadratic convergence for Newton's iterative method.


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


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

                                                    0

                                                    Там уже на 4-5й итерациях число довольно близкое к правильному, но число верных разрядов, в отличии от точности, удваиваться всё равно не хочет.


                                                    Вот в таком виде хорошо видно как удвоение точности, так и то что происходит с "верными" разрядами:


                                                    4: 0.035400882313246180 ± 0.01540088291324
                                                    5: 0.023350017521141124 ± 0.00335001812114
                                                    6: 0.020240312312009873 ± 0.00024031291200
                                                    7: 0.020001426015725795 ± 0.00000142661572
                                                    8: 0.019999999450877180 ± 5.08771809404e-11
                                                      –1

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


                                                      Благодарю за дискуссию!

                                                        0

                                                        Один вопрос — зачем вы дали мне источник? Я и так знаю, что мои формулы в него укладываются, спорю-то я не с вами...

                                                          0

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

                                                            0

                                                            Я спорю не с жарнонизмом, а с вот с этой цитатой:


                                                            Вы сами приводили как пример, вычисление квадратного корня, где количество ТОЧНЫХ знаков удваивается после каждой итерации. Точные — означает, что они не меняются при дальнейших вычислениях.

                                                            Вот и я привёл пример., где вроде как "точный" разряд ещё как меняется при дальнейших вычислениях.

                                                              –1

                                                              Ох, так и быть, отвечу, нарушу своё обещание :)


                                                              Вы цитируете человека, который приписал мне то, чего я не говорил. Вот с ним, пожалуйста, и дискутируйте.


                                                              Хотя может это я не прав, не вижу, кому отвечаете. Тогда прошу прощения, не освоился, видимо :)

                                                                +1

                                                                Да, всё, я увидел, что вы спорите не со мной, ещё раз прошу прощения :) У меня возникли основания полагать, что вы поверили, что я сказал глупость про точные числа в расчёте квадратного корня, хотя я просто не мог такое сказать.

                                                      0
                                                      Спасибо за информацию, не углублялся в этот вопрос и верил, что «удваивается количество точных знаков», как сказано во многих туториалах.
                                                        0

                                                        Да, так пишут потому что это просто удобно так думать, это условно, при этом обычно подразумевается, что есть два обстоятельства:


                                                        • Бывают пограничные случаи, где ни одной правильной цифры не будет, под удвоением подразумевают возведение точности в квадрат.
                                                        • Бывают условия, при которых квадратичная сходимость скатывается до линейной. И тогда одна итерация, грубо говоря (условно), добавляет только одну цифру.
                                      +1

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


                                      Для примера рассмотрю простейшее округление — усечение. Вот мы считали-считали какую-то формулу, и получили в итоге число 1,00 с абсолютной погрешностью не более 0,005, сокращенно 1,00 ± 0,005. Как его усекать — до 0,99 или до 1,00? Хорошо если ответ "1,00 ± 0,005" нас устраивает, но что делать если у нас задача — обязательно получить два корректных разряда?

                                      +3
                                      В компиляторе Intel есть ключи, которые позволяют управлять точностью трансцендентных функций (типа -fimf-accuracy-bits), и по умолчанию стоит -fimf-precision=medium, т.е. не самая точная из имеющихся.

                                      Как правило, разница в скорости исполнения между 1ulp функцией и корректно округлённой функцией — 2-4 раза. А ведь тут мы говорим, про [почти] один распоследний бит ответа, который никто не заметит.
                                        0
                                        Совсем прикладной программист может воспользоваться библиотекой интервальной арифметики и не гадать с ответом.
                                        Для реализации надо иметь только округление вверх и округление вниз.
                                        Будет медленно, но математически точно (диапазон в котором лежит ответ).
                                        Есть принятый стандарт: IEEE 1788-2015.
                                        Отмечу, что для интервальной арифметики и анализа Уильямом Кэхэн так же сделал очень много.
                                          +1
                                          Корректная округлённость на практике нужна, в основном, для воспроизводимости результата на другой платформе. Библиотека с интервальной арифметикой не гарантирует, что интервал на двух разных платформах будет одинаков.
                                            0
                                            Спасибо за пояснение практической цели.
                                            Для интервалов гарантии одинаковости нет, но есть гарантия того, что на любой платформе точный результат f(x)=0 будет внутри вычисленного интервала. Кроме того, пересечение интервалов, полученных от всех платформ, так же будет содержать точный результат. Но как обычно у всего есть цена.
                                          0

                                          Ну, для начала, интересно, на какой платформе вы делали описанные тесты. Версия компилятора может ничего не значить, особенно в Unix мире.


                                          Вот на Ubuntu 18.04/amd64 мне одинаково GCC и Clang выдали два раза 1.00206053 (хотя я мог не проверять на разных — очевидно, что эти функции сидят в glibc — хотя и выставлены отдельно в libm). А вот FreeBSD выдала такие же значения, как у вас в тесте — с разницей в 1 ULP, причём одинаково на 32 и 64 битах. Но что FreeBSDʼшная libm достаточно часто даёт ошибку в 1 ULP, а glibc — нет, известно, наверно, каждому, кто этим занимается серьёзно.


                                          Теперь к этому:


                                          Причина ошибок не всегда кроется в лени или низкой квалификации разработчиков, а в одном фундаментальном обстоятельстве, преодолеть которое современная наука пока не смогла.

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

                                          Вообще-то известно, но решается чуть не так: для всех IEEEʼшных методов округления, и даже для расширений (типа 8 методов в SystemZ) это максимум 3 дополнительных бита (у них даже свои названия есть: guard, round и sticky — вы их знаете, но оставим для остальных) плюс факт ненулевого значения дальше этих битов (точное значение уже не нужно, нужно просто знать, получилось там твёрдое 0 или нет). Всё. Причём сам факт ненулевого хвоста превращается в установку sticky bit в 1, после чего эти три бита вместе с последним значащим подаются в логику округления.
                                          В отдельных реализациях могут делать больше 3 бит, но умеренно — посмотрите например в Berkeley softfloat library — для float32 это 6 бит, но потому, что вообще у них работа идёт с 30-битовыми значениями в 32-битном поле, и удобно делать именно так.


                                          Да, тут решение достигается за счёт того факта "хвост не ноль", который при формальном точном следовании подходу "сколько бит?" мог бы требовать неограниченной точности. Но в том и соль, что это ограничение такими формальными рамками совершенно искусственное. И для большинства этих функций для большинства аргументов, более того, не надо сложных вычислений для "хвост не ноль": есть очень мало аргументов, когда он может стать нулём (типа случай sin 30° == 0.5, тут можно эти 30° в радианах вынести в отдельную проверку), в остальных случаях можно "не ноль" сразу постулировать.


                                          А разрешение ошибки в 1 ULP для трансцедентных функций вызвано, насколько мне известно, не проблемой округления, а тем, что решили, что делать тут полную точность при том, что промежуточные значения могут иметь ту же точность, оказывается слишком дорого. Напомню, что IEEE754-1985 писался по аппаратным реализациям Intel и Motorola, где вся логика вшивалась в сам процессор. Это уже потом поняли, что лучше в процессоре оставить только базовую арифметику и служебные операции, а трансцедентные вынести в софт (как оно сейчас и работает с SSE, Neon и всеми ровесниками и более поздними реализациями). Рядом говорят, что этот один ULP может требовать в 2-4 раза больше времени. Да, возможно. Но сама реализация в софте возможна за подъёмные деньги, в отличие от железа.


                                          UPD: прочёл ваши замечания в другом комментарии про "пограничные" случаи, где микроошибка может перевернуть длинную цепочку бит. Но таких случаев от всех, где реализации дают ошибку в 1 ULP, дай бог чтобы 1 из 1000. Остальные — это просто экономия на последних шагах, которые не просто дороги, а крайне сложны в реализации. Если посмотреть на софтовые реализации трансцедентных функций — видно, какие тонкие игры там применяются, чтобы не потерять этот самый один бит, оставаясь в поле той же точности. Вот это и есть, мне кажется, основной источник проблемы. Если бы можно было легко нарастить точность промежуточных значений при расчёте какого-нибудь log10 рядами — мы бы об этой проблеме ничего не слышали.
                                          (Собственно, 80-битный формат с 64-битной мантиссой FPU Intel и Motorola — как раз попытка решения этой проблемы, насколько известно из исторических документов.)


                                          причём разные компиляторы ошибаются в разных местах, но ни один не будет работать правильно.

                                          Проверьте на линуксе со свежим glibc. Уже интересно, что там найдёте. Они достаточно точно реализовали операции, насколько мне известно, причём достаточно давно. Ну а проверить можно через mpfr, у неё гарантированно своя реализация с управляемо неограниченной точностью.


                                          Во-первых, задача носит фундаментальный характер.

                                          … и давно таки решена? — посмотрите на тот же MPFR. Ну или покажите, что она тоже ошибается в 1 ULP… что-то я такого про неё не слышал.


                                          Замечание по форме: есть такой формат у printf — %a. Лучше пользоваться им, чем переводить состояние в памяти в hex-форму и печатать из неё; и лучше задавать константы в таком контексте тоже лучше в виде 0x1.921fb54442d18p+1, чем в длинной десятичной константе.

                                            0

                                            Благодарю за совет по поводу %a, но в остальном я не могу прокомментировать ваши замечания, потому что из них очевидно, что вы не поняли содержание статьи. Возможно, вы захотите мне возразить, но не стоит, а лучше напишите сразу тем авторам из списка источников, которые более 15 лет уже пытаются решить эту задачу и всё ещё далеки от решения. Возможно, они удивятся, что оказывается всё уже давно решено.


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

                                              –1
                                              но в остальном я не могу прокомментировать ваши замечания, потому что из них очевидно, что вы не поняли содержание статьи.

                                              Как раз из вашего комментария очевидно, что я его таки понял, но не считаю настолько важной акцентуацию именно на тех пограничных случаях с "нечёткими битами" (буду пользоваться таким описательным выражением), о которых вы говорили (причём только в комментариях подчеркнули это явно — до того было непонятно некоторым другим комментаторам).
                                              Вы почему-то сочли, что основная или даже единственная причина (тут уточните сами) того самого разрешения погрешности в 1 ULP — это нечёткие биты на пограничных значениях. Так вот — это не просто сомнительно, это скорее всего неправда, и я встречал объяснение (в документах по истории IEEE754), почему это неправда.


                                              Я внимательно писал статью, а ваша задача — внимательно её прочитать, если намереваетесь дискутировать.

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


                                              Если хотите тут действительно осмысленного рассмотрения и статьи, которая имела бы реально значение повыше простого копирования известных истин — я считаю, что вам следует:


                                              1. Опереться на opensource — его тут достаточно (в смысле количества вариантов), а не на неизвестные "компиляторы" мира Windows, когда вы даже не определили, от какого компонента зависит реализация, и как связаны версионности компилятора и этого компонента (hint: связь неоднозначна, а версия компилятора недостаточна).


                                              2. Провести сравнение типовых реализаций: MSun исходная; MSun осовремененная (в рамках FreeBSD, например); IBM Accurate Mathematical Library исходная; она же доработанная; MPFR. Насколько есть данные, и про x87 FPU (и про его точность — там есть несколько парадоксальных моментов).


                                              3. Взять хотя бы пару функций для примера и проверить утверждение про бо́льшую (или полную) точность glibc. Опровержение меня тут (на какой-то функции, которую я не проверял) было бы на порядки полезнее всего это "вы меня не читаете" ;)


                                              4. Объяснить, с помощью каких же приёмов код той же самой IBM AML ухитряется обеспечить точность до последней цифры (всегда или почти всегда) без расширения точности используемых чисел.



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

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

                                          Самое читаемое