Русский
Русский
English
Статистика
Реклама

Ieee-754

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

17.10.2020 10:23:13 | Автор: admin
Здравствуйте, друзья, как вы думаете, если мы напишем такой код:

s = a+b;z = s-a;t = b-z;

то не кажется ли вам, что в результате его выполнения получится, что t=0? С точки зрения привычной математики действительных чисел это и правда так, а вот с точки зрения арифметики с плавающей запятой в переменной t будет кое-что другое. Там будет то, что спасает нас от потери точности при сложении чисел $a$ и $b$. Кого интересует данная тема, прошу под кат.


По моей традиции текстовая статья дублируется на видео. По содержанию текст и видео идентичны.



Читателю статьи для её восприятия нужно понимать способ представления чисел с плавающей запятой в формате IEEE-754 и понимать, почему, например, 0,1+0,20,3, но если такого навыка у вас по каким-то причинам нет, но вы хотели бы его приобрести, то прошу обратить внимание на список источников конце статьи, на пункты [1] и [3-5].

Итак, у нас следующая проблема. Сумма двух чисел с плавающей запятой c=a+b может вычисляться с некоторой ошибкой. Знаменитый на весь интернет пример 0,3 0,2+0,1 хорошо это показывает. Наша задача в том, чтобы всё-таки складывать числа без этой ошибки. То есть сделать так, чтобы мы могли как-то сложить 0,2+0,1 и хоть в каком-то виде знать точный результат. В каком смысле точный, если даже исходные числа 0,1 и 0,2 не имеют точного представления в формате IEEE-754? Вот сейчас и поясню.

Начнём с того, что чисел 0,1 и 0,2 в двоичной арифметике с плавающей запятой быть не может, а наиболее близкие к ним значения для типа данных double (число удвоенной точности binary64, так его называют в Стандарте IEEE-754) следующие:

$0{,}1 0{,}1000000000000000055511151231257827021181583404541015625 ={}\\ {}=\frac{3602879701896397}{36028797018963968}$


$0{,}2 0{,}200000000000000011102230246251565404236316680908203125 ={}\\ {}=\frac{3602879701896397}{18014398509481984}$


Десятичный результат получен с помощью правильного онлайн-конвертера, а дробь посчитана с помощью библиотеки MPIR и функции mpq_set_d, которая умеет переводить double к рациональной дроби.

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

0,1000000000000000055511151231257827021181583404541015625 +0,2000000000000000111022302462515654042363166809082031250 =0,3000000000000000166533453693773481063544750213623046875.

К сожалению, если просто написать код s=0.1+0.2, то мы получим кое-что другое, а именно

s == 0.3000000000000000444089209850062616169452667236328125

что отличается от правильного ответа ровно на

$t = 2{,}77555756156289135105907917022705078125 \cdot 10^{-17}$


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

Зачем может понадобиться точное сложение двух чисел? Приложений много, но вот одно из них, которое будет понятно всем. Если вы хотите сложить все числа из массива X[N] и сделаете это вот так:

S = 0.0;for (i=0; i<N; ++i)  S = S + X[i];

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

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

У нас есть два числа $a$ и $b$ типа double. Нам не важно, что эти числа $a$ и $b$ уже содержат в себе какую-то погрешность, полученную, например, в результате конвертирования десятичной строки в формат IEEE-754. Важно, что они сейчас перед нами, а их предыстория нас не интересует. Нужно сложить эти два числа c=a+b так, чтобы в результате сложения не возникло погрешности.

Но ведь это невозможно!

Да, это невозможно, спасибо что дочитали, до новых встреч :)

Шучу, конечно. Это невозможно только если мы храним результат сложения в одной переменной того же типа данных. Но вспомните пример выше. У нас была переменная s, и переменная t. Причём s была округлённым результатом сложения a+b, а t была равна разности этого округлённого результата и точного значения суммы. При этом выполняется равенство s+t=a+b.

И вот хорошая новость состоит в том, друзья, что такое представление суммы a+b в виде суммы s+t можно выполнить всегда (если a+b)! Если s=a+b оказывается точно-представимым значением в типе double, то очевидно, что t=0. Если это не так, то t будет равно некоторому очень маленькому (по сравнению с s) числу, которое является точно-представимым. Итак, вот оно, фундаментальное свойство суммы чисел с плавающей запятой: ошибка округления в результате суммирования чисел типа double всегда будет точно-представимым числом типа double! А это означает, что пара чисел (s, t) всегда даёт нам возможность сохранить сумму чисел a+b как бы без погрешности. Да, мы вынуждены хранить две переменные вместо одной, но прелесть в том, что их всего две, а не больше.

Теперь опишем это свойство на математическом языке. Введём обозначение RN(x) это результат приведения произвольного вещественного числа x к типу данных double по правилу округления round-nearest-ties-to-even, то есть это округление к ближайшему, но в случае равного удаления от двух ближайших к тому, у которого последний бит мантиссы равен нулю (чётный). Именно этот режим работает почти везде по умолчанию, то есть если вы не знаете, о чём я сейчас говорю, то в вашем процессоре на 100% работает именно этот режим, так что не беспокойтесь.

Пусть a и b числа типа double. Пусть |a||b| и RN(a+b) . Тогда следующий алгоритм

s = RN (a+b);z = RN (s-a);t = RN (b-z);return (s, t);

вычисляет пару (s, t) таких чисел, для которых:
  1. s+t = a+b (в точности)
  2. s число типа double, наиболее близкое к a+b.


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


Прямоугольник олицетворяет переменную с плавающей запятой фиксированного размера, поэтому прямоугольники, помеченные символами а и b имеют одинаковый размер, но $b$ смещён относительно $a$ вправо, потому что у нас есть условие: |a||b|. Третий прямоугольник Точное а+b это некоторое промежуточное значение, которое может иметь больший размер, чем размер переменных $a$ и $b$, поэтому оно будет округлено, и хвостик, вылезающий за границы нашего типа данных, будет отброшен. Таким образом, мы переходим к четвёртому прямоугольнику Округлённое a+b, это и есть наше значение s, полученное в первой строке алгоритма. Если теперь из s снова вычесть $a$, то останется только b без хвостика. Это делается во второй строке алгоритма. Теперь мы хотим получить хвостик от $b$, и это делается в третьей строке алгоритма, когда из исходной переменной $b$ мы вычитаем b без хвостика. Остаётся хвостик, это и есть переменная t, которая показывает ту самую погрешность, которая возникла в ходе округления. При этом из данных рассуждений видно, что такой хвостик всегда будет умещаться в одну переменную, потому что он не может превышать размера переменной $b$. В крайнем случае, когда смещение $b$ относительно $a$ будет настолько большим, что s=RN(a+b)=a, то в этом случае, очевидно, z=0 и t=b. Скучный рисунок на эту тему я вам предлагаю изобразить самостоятельно.

Если вы не находите картинки убедительными для себя, то в книге [1, раздел 4.3.1, теорема 4.3] есть строгое математическое доказательство, но оно не умещается в формат научно-популярной статьи. Его краткая суть в том, что в нём показано почему в строках 2 и 3 алгоритма не будет выполняться округления, то есть эти выражения вычисляются точно, а если так, значит t=b-z=b-(s-a) = (a+b)-s в точности, что нам как раз и нужно: t является разностью между реальной суммой a+b и её округлённым значением s.

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

Пример 1.

$\begin{array}{rcl} a &{}={}& 9007199254740991 = 2^{53}-1, \\ b &{}={}& 2, \\ s &{}={}& 9007199254740992, \\ t &{}={}& 1. \end{array}$


Пояснение. Реальное значение a+b=253+1, однако в типе данных double младший бит, равный единице, не влезет в 52 бита дробной части мантиссы, по какой причине будет выброшен при округлении, но переменная t подхватит его и сохранит в качестве погрешности, а переменная s сумеет сохранить только 253 ровно. Легко видеть, что s+t будет равно 253+1.

Пример 2.

$\begin{array}{rcl} a &{}={}& 1152921504606846976 = 2^{60}, \\ b &{}={}& 1073741823 = 2^{30}-1, \\ s &{}={}& 1152921505680588800, \\ t &{}={}& -1. \end{array}$


Пояснение выполнено в двоичном коде.

  a = 1000000000000000000000000000000000000000000000000000000000000  b =                                111111111111111111111111111111a+b = 1000000000000000000000000000000111111111111111111111111111111      -------------------------------------бит округления--^   s = 1000000000000000000000000000001000000000000000000000000000000

Если выравнивание кода не поехало, то вы видите, где у суммы a+b будет бит округления, и что само округление будет выполнено вверх. Чтобы снова получить паровоз из 30 единиц, нужно вычесть единицу из s=260+230. Поэтому t=-1. Как видим, t может быть отрицательным, когда округление выполняется вверх.

Я не говорил этого явно, но алгоритм работает даже когда b<0, то есть фактически происходит вычитание.

Пример 3.

$\begin{array}{rcl} a &{}={}& 9007199254740991 = 2^{53}-1, \\ b &{}={}& -2251799813685247{,}75 = -\frac{2^{53}-1}{4}, \\ s &{}={}& 6755399441055743, \\ t &{}={}& 0{,}25. \end{array}$


В этом примере число $b$ является сдвигом вправо на 2 бита числа $a$, поэтому если $a$ влезает в double впритык, ясно, что вычитание здесь точным получиться не может. Останется небольшой хвостик размером в два бита после запятой.

Из доказательства, приведённого в [1], следует один положительный момент данного алгоритма: в ходе его выполнения не может возникнуть переполнения во второй и третьей строках, если оно не возникло при вычислении s=RN(a+b).

Недостатком является то обстоятельство, что алгоритм работает только для |a||b| (однако если вы посмотрели доказательство, то знаете, что в действительности достаточно более слабого условия, чтобы экспонента $a$ была не меньше экспоненты $b$, но проверить это условие намного сложнее). То есть нам нужно либо делать лишнюю проверку, либо как-то иначе гарантировать правильный порядок чисел на входе алгоритма. Ветвление не всегда обрабатывается достаточно быстро, это зависит от процессора и от того, как на это ветвление посмотрел компилятор (иногда он может его убрать, иногда нет, а иногда он его убрал, но оказалось так плохо, что лучше бы не трогал). Возникает вопрос: а можно сложить числа без ошибок округления, если их порядок друг относительно друга заранее неизвестен?

Можно. Для этого есть второй алгоритм, но в нём вместо 3-х операций их уже 6.

  s = RN (a+b);  b' = RN (s-a);  a' = RN (s-b');  d_b = RN (b-b');  d_a = RN (a-a');  t = RN (d_a+d_b);  return (s, t)

Как это работает? В случае когда |a||b|, алгоритм работает по сути так же как предыдущий, и строгое доказательство для этого случая, которое вы можете отыскать в [2, раздел 2.3, теорема 7], сводится к тому, чтобы показать, что лишние три операции не портят этого результата. А вот основная сложность доказательства в том, чтобы рассмотреть случай |a|<|b|, он значительно труднее и тоже выходит за рамки лёгкого научно-популярного изложения. Схематичная картинка общей ситуации при |a|<|b| даётся ниже с пояснением.


Первые два прямоугольника показывают соотношение между $a$ и $b$, на этот раз $a$ смещено относительно $b$ вправо. Следующие три прямоугольника показывают процедуру сложения и отбрасывания хвостика. Дальше начинается основная сложность. Вычисляем b=s-a, то есть мы вычитаем из округлённого значения s число $a$ и выброшенный xвостик, а это значит, что b будет чуть меньше, чем $b$, и этот недостаток отмечен красным прямоугольничком. Другое дело, когда мы вычисляем a=s-b, то есть вычитаем из суммы величину b с недостатком, а значит получим а без хвостика с избытком, и этот избыток обозначается зелёным прямоугольничком. Если мы теперь вычтем из $b$ величину b (которая с недостатком), то получим избыток db. Далее вычитаем из $a$ величину a (которая с избытком), и получаем хвостик с недостатком da. Если теперь сложить хвостик с недостатком и избыток, мы получим чистый хвостик.

У читателя может возникнуть вопрос: а в каких случаях можно с уверенностью сказать, что t=0? Иными словами, про какие пары чисел (a, b) точно известно, что их сумма будет точно-представимым (без округления) числом в том же типе данных?

Мне известны только три теоремы на этот счёт, которые также приводятся в [1] и [2].

Теорема 1 [книга 1, раздел 4.2.1, теорема 4.2]. Если сумма RN(a+b) оказалась денормализованным числом, то сумма a+b является точно-представимой.

Теорема 2 [статья 2, раздел 2.2, лемма 4] Если |a+b||a| и |a+b||b|, то число a+b является точно-представимым (аналогичное верно и для разности).

Теорема 3 [книга 1, раздел 4.2.1, лемма 4.1] Пусть a0 и a/2 b 2a, тогда разность a-b будет точно представимой.

Бдительный читатель сразу вспомнит про так называемую ошибку критического сокращения старших значащих битов, называемую в литературе Catastrophic Cancellation и спросит: Как же так, ведь давно известно, что когда мы вычитаем очень близкие друг к другу числа, там возникает такая погрешность, что ой-ой-ой! Как ты говоришь, что в этом случае разность будет точно-представимой.

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

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

Благодарю за внимание!

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


[1] Jean-Michel Muller, Handbook of floating-point arithmetic, 2018.
[2] Jonathan Richard Shewchuk, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates, Discrete & Computational Geometry 18(3), 1997, pp. 305363.
[3] На Хабре: Что нужно знать про арифметику с плавающей запятой.
[4] На Хабре: Наглядное объяснение чисел с плавающей запятой.
[5] Учебный видео-курс для самых маленьких, предельно наглядное разъяснение чисел с плавающей запятой в 8-ми уроках. Первый урок на на YouTube.
Подробнее..

Наиболее точное скалярное произведение векторов типа double. Вычисление значения полинома

28.11.2020 08:12:33 | Автор: admin
В прошлой статье речь шла о том, как можно сложить массив из чисел типа double наиболее точно, то есть получить такую сумму, как если бы мы выполняли вычисления в рамках арифметики с бесконечной точностью, а затем один раз округлили бы результат. Был показан алгоритм, который эквивалентен применению типа данных double-double, в котором сложение происходит сразу в двух переменных: основная сумма и хвостик-погрешность. Опытные читатели сразу догадались, что сложение хвостиков-погрешностей также допускает по отношению к себе рекурсивное применение того же алгоритма, что приводит не к удвоенной, а к утроенной точности, и вообще, можно организовать каскад сложений произвольного размера, получая любую наперёд заданную точность расчётов, поэтому фактически в прошлой статье была показана предпосылка к так называемой дробной длинной арифметике. Опытный программист без труда разберётся как её реализовать, ну а я обещал дать аналогичные фундаментальные основы для скалярного произведения и вычисления полинома в точке. Поскольку все базовые вводные слова уже были сказаны в двух предшествующих статьях, в этой будет меньше воды и лишних, по мнению опытных математиков, сведений. Прошу под кат.





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



Умножение двух чисел без потери точности (ну почти)


О сложении мы знаем уже достаточно: с помощью функции two_sum можно сложить N чисел наиболее точно, дело осталось за малым: научиться перемножать два числа типа double так, чтобы не терять отбрасываемый при округлении хвостик, сохраняя его в отдельной переменной. Если в случае сложения сделать это можно очень просто при любом раскладе, лишь бы промежуточная сумма $a+b$ не переполнялась (кто не читал: Сложение двух чисел с плавающей запятой без потери точности), то в случае умножения есть дополнительное ограничение.

Обозначим через $e_{min}$ минимальную экспоненту вашей системы чисел с плавающей запятой, в рамках которой выполняется умножение (для float: 126, для double: 1022), буквой $k$ обозначим количество битов мантиссы с учётом скрытой единицы (24 для float и 53 для double). Пусть $e_a$ и $e_b$ экспоненты наших исходных чисел $a$ и $b$, тогда требуется, чтобы было выполнено условие:

$e_a+e_b \geq e_{min} +k-1.$


Если данное условие не выполняется, то могут возникнуть случаи, когда $ab$ нельзя представить в виде суммы $p+t$, где $p = RN(ab)$ и $t=ab-p$. Указанное неравенство является сложным способом сказать простую фразу: хвостик $t$ не должен быть денормализованным. И это вполне понятно: когда мы перемножаем два числа, получаем удвоенное количество битов, которые гарантированно влезают в две переменные того же типа данных. Если же одна из этих переменных попадает в денормализованный диапазон, происходит битовый сдвиг мантиссы с потерей младших битов. И кроме переполнения это единственная ситуация, при которой произведение двух чисел нельзя точно представить в форме суммы двух чисел того же типа.

Так ли существенно это ограничение? Для float нужно чтобы сумма экспонент оказалась не меньше 126+241=103, то есть, грубо говоря, если каждый из множителей больше 251, то всё будет хорошо. Спросите себя: часто ли вы применяете числа такого порядка? Если да, то можете выполнить масштабирование и перейти к числам побольше. Можете применить double, и работать с числами больше 2485. Но даже если плюнуть на всё и применить алгоритм в неподходящих условиях, в простых прикладных задачах он всё равно отработает правильно чуть чаще чем всегда.

Если указанное условие выполнено, то с помощью инструкции FMA (Fused Multiply-Add) точное умножение двух чисел выполняется следующим алгоритмом:

double __fastcall two_prod (double &t, double a, double b) {  double p = a*b;  t = fma (a, b, -p);  // t = a*b-p, не забудьте подключить <cmath>  return p;}

Результатом работы функции будет пара чисел $(p, t)$, такая, что $p$ является ближайшим к $ab$ числом, а $t$ хвостиком-погрешностью, причём $p+t=ab$в точности. Разумеется, по ходу расчётов не должно произойти переполнения, за этим придётся следить отдельно: нужно оставаться в рамках условия $RN(ab)\neq\infty$.

Если операция FMA для вас недоступна, то ситуация сильно усложняется. Желаемое произведение можно получить путём разбиения исходных переменных $a$ и $b$ на две части примерно посередине так, чтобы старшая часть операндов (ah и bh) содержала половину старших битов, а младшая (al и bl) половину младших. Далее нужные операции выполняются именно с этими половинками, в общей сложности алгоритм потребует 17 операций сложения и умножения. Данный трюк в чём-то аналогичен известному алгоритмическому трюку для вычисления произведения двух целых чисел с сохранением результата в удвоенном типе данных (кто не знал, читайте Г. Уоррен Алгоритмические трюки для программистов, в издании 2014 года это раздел 8.2). Нижеследующий алгоритм написан мною по мотивам книги Jean-Michel Muller, Handbook of floating-point arithmetic, также как и остальная часть этой статьи.

    // Разбиение x на старшую и младшую частьvoid __fastcall split (double &hi, double &lo, double x) {  double gamma = 0x8000001*x; // gamma=(2**27+1)*x; 27 = ceil(53/2), где 53 - размер мантиссы double (со скрытым битом).  double delta = x - gamma;  hi = gamma + delta;  lo = x - hi;}    // Произведение с разбиением, результат идентичен функции two_proddouble __fastcall two_prod_s (double &t, double a, double b) {  double ah, al, bh, bl, p;  split (ah, al, a);  split (bh, bl, b);  p = a*b;  t = (ah*bh-p) + ah*bl + al*bh + al*bl;  return p;}

Данный алгоритм для двоичных типов данных из Стандарта IEEE-754 работает всегда, для десятичных типов он будет работать только если количество цифр мантиссы чётно, то для есть для decimal32 он может давать ошибочный результат. Если же вы изобретаете свой тип данных с плавающей запятой, читайте полное условие теоремы 4.8 в упомянутой выше книге в разделе 4.4.2.2 Dekkers multiplication algorithm. Там же, но разделом ранее разъясняется как подбирать магическую константу для правильного разбиения в функции split. В двоичной системе она подбирается равной 2s+1, при этом младшая часть lo будет содержать s битов, а старшая все остальные.

Вычисление скалярного произведения


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

    // Функция точного сложения двух чисел из предыдущей статьиdouble __fastcall two_sum (double &t, double a, double b) {  double s = a+b;  double bs = s-a;  double as = s-bs;  t = (b-bs) + (a-as);  return s;}using dvect_t = vector<double>;double __fastcall ogita (const dvect_t &X, const dvect_t &Y) {  double s=0.0, c=0.0, p, pi, t;  for (size_t i=0; i<X.size(); ++i) {    p = two_prod (pi, X[i], Y[i]);// p = X[i]*Y[i], но в pi остался хвостик.    s = two_sum (t, s, p);// s += p, но в t остался хвостик.    c += pi+t;// Собираем хвостики  }  return s+c;// Общая сумма и хвостик}

Понятно, что two_prod нужно заменить на two_prod_s, если вам недоступна FMA.

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

Помимо этой версии скалярного произведения есть ещё две наивные функции, разъяснение которых в тексте статьи будет моветоном.

double __fastcall naive (const dvect_t &X, const dvect_t &Y) {  double s=0.0;  for (size_t i=0; i<X.size(); ++i)    s += X[i]*Y[i];  return s;}double __fastcall naive_fma (const dvect_t &X, const dvect_t &Y) {  double s=0.0;  for (size_t i=0; i<X.size(); ++i)    s = fma (X[i], Y[i], s);  return s;}

Тестирование выполнялось точно также как в статье про сумму N чисел. Напомню только, что верхнее число в ячейке средняя погрешность в ULP на 100 тестах, а нижнее максимальная погрешность на них же. Размер каждого теста N=106. В первом столбце указано выбранное распределение случайных чисел.

Naive Naive FMA OgitaRumpOishi
U[1, 2) 95.31
275
95.36
276
0.00
0
U[1, 2) 441.61
6878
432.99
7489
0.00
0
U[1e-10, 1e10) 2478.61
2740
2478.55
2739
0.00
0
U[1e-10, 1e10) 392.50
11893
394.86
12126
0.00
0
exp[2] 183.70
513
183.65
512
0.00
0
exp[2] 329.21
5962
347.31
7251
0.00
0
N(0, 1) 1554.07
109700
1408.30
95896
0.00
0

Ожидаемо, что алгоритм с применением double-double сработал безошибочно, и так будет происходить на всех типичных бытовых задачах. Таблица не требует особых комментариев кроме одного: наивный алгоритм с FMA может давать менее точный результат, чем ещё более наивный алгоритм с двумя округлениями на каждом шаге. При этом ДА, я выполнил дизассемблирование и убедился, что в этом наивном методе компилятор не всунул FMA в качестве оптимизации. Вообще говоря, он не имеет на это права без моего разрешения, но всё же осторожность не помешает.

Вычисление значения полинома


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

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

double __fastcall naive (const dvect_t &A, double x) {  double s=0.0;  for (auto a: A) {    s *= x;    s += a;  }  return s;}double __fastcall naive_fma (const dvect_t &A, double x) {  double s=0.0;  for (auto a: A)    s = fma (s, x, a);  return s;}    // Алгоритм по именам трёх авторов: GraillatLangloisLouvetdouble __fastcall graillat (const dvect_t &A, double x) {  double s=0.0, c=0.0, p, pi, t;  for (auto a: A) {    // Полином из основных значений    p = two_prod (pi, s, x);    s = two_sum (t, p, a);    // Полином из хвостиков    c *= x;    c += pi+t;  // Можете применить здесь FMA, но не вижу особого смысла  }  return s+c;// Фактически, здесь сумма значений двух полиномов.}

Поскольку значение полинома на случайных тестах очень быстро возрастает и создаёт трудности для точного расчёта (эталонное значение ответа, как вы знаете, я вычисляю через дроби), для тестов взято ограничение N=100 элементов в массиве. Комментарии снова излишни, кроме одного: здесь применение FMA для наивного расчёта становится оправданным точность заметно повышается.

Naive Naive FMA GraillatLangloisLouvet
U[1, 2) 3.15
11
2.67
9
0.00
0
U[1, 2) 3.96
29
3.03
27
0.00
0
U[1/10, 10) 1.86
11
1.47
7
0.00
0
U[1/10, 10) 2.20
61
1.29
10
0.00
0
exp[2] 0.59
5
0.46
4
0.00
0
exp[2] 0.90
30
0.71
6
0.00
0
N(0, 1) 0.90
30
0.71
6
0.00
0

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

Одна из лучших наград автору успешное применение описанных им знаний на практике для всеобщей пользы :) Благодарю за внимание!

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


[1] Разделы 4.4.2.2 и 5.4-5.5 книгиJean-Michel Muller, Handbook of floating-point arithmetic, 2018.

[2] Раздел 8.2 книги (по поводу умножения целых чисел) Генри С. Уоррен мл. Алгоритмические трюки для программистов, 2014.
Подробнее..

Перевод числа в строку с помощью FPU

17.01.2021 06:18:51 | Автор: admin

Введение

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

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

Постановка задачи

Формализуем задачу. Необходимо перевести восьмибайтовое число из формата IEEE-754 в текстовую строку, где указан знак числа, мантисса заданной длины и десятичный порядок с признаком E и своим знаком.

Так, если задано число 1234.567890, то строка с мантиссой, например, в 16 знаков должна выглядеть как 1.23456788999999E+003. Вместо плюса у мантиссы нужно писать пробел, а сама мантисса должна быть не меньше единицы. Кстати, данный пример иллюстрирует дискретность и приближенность представления чисел в формате IEEE-754: приведенное исходное число не может быть точно представлено как 1.23456789000000E+003, что, может быть, интуитивно ожидалось.

Использование команд FPU для преобразования

На первый взгляд, решение выглядит простым. В устройстве FPU (Floating Point Unit) процессора даже имеются команды, явно предназначенные и для перевода чисел из формата IEEE-754 в текст. Это, во-первых, команда EXTRACT разделения числа на мантиссу и порядок, а во-вторых, команда FBSTP выдающая мантиссу сразу в виде двоично-десятичного значения. Однако при ближайшем рассмотрении эти команды не дают нужного результата.

Например, если применить команду FBSTP к указанному выше числу, то получится 10 байт со значением 35 12 00 00 00 00 00 00 00 00, поскольку нецелочисленные значения сначала округляются согласно полю RC управляющего слова FPU. Это двухбитовое поле RC может принимать 4 значения, но среди них нет режима отключить округление. Поэтому часть мантиссы просто теряется, единственное чего можно добиться режимами округления это еще получить значение 34 12 при округлении к меньшему.

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

Итак, при использовании возможностей FPU придется решать две задачи.

Во-первых, умножением (или делением) исходного числа на десять нужно перегнать всю мантиссу в целочисленную часть числа, тогда команда FBSTP выдаст, наконец, в двоично-десятичном виде все цифры мантиссы.

Во-вторых, нужно определить такой десятичный порядок (опять-таки умножением или делением на десять), при котором результат попадет в диапазон между 1 и 10. Это нужно для представления числа в виде мантиссы с одной цифрой перед точкой и найденным порядком. Увы, совместить эти две задачи в едином цикле умножения/деления невозможно.

Причем есть и подводный камень в виде значения числа, максимально приближенного к единице, но не равного единице. Циклом деления или умножения легко можно ошибиться в показателе степени и вместо требуемого в данном случае 9.999E-001 получить неправильное значение типа 9.999E+000. К сожалению, при всем богатстве команд FPU обойтись без циклов деления и умножения на десять не удается.

Алгоритм преобразования

При формировании мантиссы для числа типа double умножение или деление на 10 продолжаются до тех пор, пока двоичный порядок числа не достигнет (или не станет больше/меньше) 53, поскольку под мантиссу выделено именно 53 двоичных разряда.

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

Пример реализации преобразования

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

Пара комментариев к тексту. Подпрограмма написана на языке ассемблера RASM, который имеет некоторые особенности.

В частности здесь используются временные метки в виде символа @. Алгоритм их реализации следующий: транслятор имеет внутренний счетчик меток. Когда метка @ встречается в команде перехода, к ней автоматически дописывается значение этого счетчика (т.е. реально создаются метки @0000, @0001 и т.д.). Когда встречается символ @ с двоеточием, к нему также автоматически приписывается значение счетчика и после этого счетчик увеличивается на 1.

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

Кроме этого, для команд FPU в RASM не анализируется размер операнда, он имеется прямо в названиях команд в виде значений 16, 32, 64 или 80.

Но вернемся к задаче. Подпрограмме по адресу в EBX передается исходное восьмибайтовое число в формате IEEE-754 и требуемая длина текстовой строки-результата в AL.

В результате работы число записывается в стек в виде текста в научной нотации и в AL сохраняется длина этой строки. Другие регистры не сохраняются.

С целью разгрузить конвейеры процессора некоторые команды рассыпаны по тексту и не образуют непрерывных логических цепочек, что несколько затрудняет чтение текста. Но таких мест немного. Программа реентерабельна и может быть использована в параллельных вычислениях. При начале работы программы управляющее слово FPU CW=0360.

  1. Сначала в стеке выделяется место для ответа и заполняется нулевым шаблоном, т.е. значением 0.000 E+000.

  2. Затем проверяется знак числа и в зависимости от него формируется мантисса умножением или делением числа на десять.

  3. Командой FBSTP мантисса переписывается в память из FPU в двоично-десятичном виде и ее часть (заданной длины) переносится в ответ.

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

      CSEGPUBLIC ?QFCMW:;---- ПОДГОТОВКА МЕСТА В СТЕКЕ ----      MOVZX     ECX,AL          ;ЗАДАННАЯ ДЛИНА СТРОКИ      POP       EAX             ;ДОСТАЛИ АДРЕС ВОЗВРАТА      SUB       ESP,ECX         ;ВДЕЛИЛИ МЕСТО В СТЕКЕ      MOV       EDI,ESP         ;ОНО ЖЕ НАЧАЛО СТРОКИ-РЕЗУЛЬТАТА      MOV       ESI,ESP         ;ЗАПОМНИЛИ НАЧАЛО      PUSH      EAX             ;АДРЕС ВОЗВРАТА ВЕРНУЛИ НА МЕСТО      PUSH      ECX             ;ЗАПОМНИЛИ ЗАДАННУЮ ДЛИНУ ОТВЕТА;---- СНАЧАЛА ЗАПОЛНЕНИЕ СТРОКИ-РЕЗУЛЬТАТА НУЛЕВМ ШАБЛОНОМ ----      MOV       EAX,'0.0 '      ;ЗНАК И ПЕРВОЕ ЧИСЛО С ТОЧКОЙ      STOSD      DEC       EDI             ;ОСТАВЛЯЕМ ПЕРВЕ ТРИ СИМВОЛА ШАБЛОНА      SUB       CL,8            ;ВЧЛИ СЛУЖЕБНЕ ЗНАКИ И ПОРЯДОК      MOV       AL,'0'      JBE       @               ;ОШИБКА, НЕТ МЕСТА ПОД МАНТИССУ      REP       STOSB           ;ЗАПОЛНИЛИ МАНТИССУ НУЛЯМИ@:    MOV       EAX,'00+E'      ;ПОКА ДВА НУЛЯ ПОРЯДКА      STOSD                     ;ЗАПИСАЛИ НУЛЕВОЙ ПОРЯДОК      LEA       EBP,[EDI]-2     ;ЗАПОМНИЛИ АДРЕС ПОРЯДКА      FLD64     [EBX]           ;ЗАГРУЗИЛИ ЗАДАННОЕ ЧИСЛО      MOV       B PTR [EDI],'0' ;ПОКА ТРЕТИЙ НУЛЬ ПОРЯДКА;---- ПРОВЕРКА ЗАДАННОГО ЧИСЛА НА ПЛЮС/МИНУС/НОЛЬ ----      FTST                      ;СРАВНИЛИ С НУЛЕМ      FNSTSW    AX      SAHF      JNZ       @               ;ЗАДАННОЕ ЧИСЛО НЕ НОЛЬ;---- СРАЗУ ВХОД С ЧИСТМ НУЛЕМ ----      POP       EAX             ;ВХОД С ЗАДАННОЙ ДЛИНОЙ И НУЛЕМ      FSTP      ST              ;ОЧИСТИЛИ FPU ОТ ИСХОДНОГО ЧИСЛА      RET;---- ОБРАБОТКА ЗНАКА ЗАДАННОГО ЧИСЛА ----@:    JNB       @               ;ЕСЛИ ЧИСЛО ПОЛОЖИТЕЛЬНО      MOV       B PTR [ESI],'-' ;ОТМЕТИЛИ ЗНАК ОТРИЦАТЕЛЬНОГО      FABS                      ;УБРАЛИ ЗНАК В САМОМ ЧИСЛЕ@:    MOV       EDX,OFFSET ДЕСЯТЬ ;ДЕСЯТИЧНАЯ СИСТЕМА;---- ПРОВЕРКА ВЕЛИЧИН ПОРЯДКА ИСХОДНОГО ЧИСЛА ----      FLD       ST              ;РАЗМНОЖИЛИ ПОЛОЖИТЕЛЬНОЕ ЧИСЛО      POP       ECX             ;ОПЯТЬ ДОСТАЛИ ДЛИНУ СТРОКИ      FXTRACT                   ;РАЗДЕЛИЛИ МАНТИССУ И ПОРЯДОК      PUSH      ECX             ;ОПЯТЬ СОХРАНИЛИ ДЛИНУ СТРОКИ      FSTP      ST              ;ВКИНУЛИ МАНТИССУ      SUB       ESP,11          ;ВДЕЛИЛИ МЕСТО ПОД МАНТИССУ КАК BCD      FIST32P   [ESP]           ;ЗАПИСАЛИ ПОРЯДОК      CMP       D PTR [ESP],53  ;ОН УЖЕ 53 РАЗРЯДА ?      JE        M0003           ;ТОГДА МАНТИССА УЖЕ ГОТОВА      JL        M0002           ;ИНАЧЕ ЧИСЛО ОЧЕНЬ МАЛО;---- УМЕНЬШЕНИЕ ПОРЯДКА ЧИСЛА ОТ ОЧЕНЬ БОЛЬШОГО ----M0001:FLD       ST              ;РАЗМНОЖИЛИ ЧИСЛО      FXTRACT                   ;РАЗДЕЛИЛИ МАНТИССУ И ПОРЯДОК      FSTP      ST              ;ВКИНУЛИ МАНТИССУ      FIST32P   [ESP]           ;ДОСТАЛИ ПОРЯДОК      CMP       D PTR [ESP],53  ;УЖЕ 53 РАЗРЯДА ?      JLE       M0003           ;ТОГДА МАНТИССА УЖЕ ГОТОВА      FIDIV16   [EDX]           ;РАЗДЕЛИЛИ ЧИСЛО НА ДЕСЯТЬ      JMPS      M0001           ;ПРОВЕРЯЕМ НОВЙ ПОРЯДОК;---- УВЕЛИЧЕНИЕ ПОРЯДКА ЧИСЛА ОТ ОЧЕНЬ МАЛОГО ----M0002:FLD       ST              ;РАЗМНОЖИЛИ ЧИСЛО      FXTRACT                   ;РАЗДЕЛИЛИ МАНТИССУ И ПОРЯДОК      FSTP      ST              ;ВКИНУЛИ МАНТИССУ      FIST32P   [ESP]           ;ДОСТАЛИ ПОРЯДОК      CMP       D PTR [ESP],53  ;УЖЕ 53 РАЗРЯДА ?      JGE       M0003           ;ТОГДА МАНТИССА УЖЕ ГОТОВА      FIMUL16   [EDX]           ;УМНОЖИЛИ ЧИСЛО НА 10      JMPS      M0002           ;ПРОВЕРЯЕМ НОВЙ ПОРЯДОК;---- ВДАЧА МАНТИСС В ДВОИЧНО-ДЕСЯТИЧНОМ ВИДЕ ----M0003:FBSTP     [ESP]+1         ;ЗАПОМНИЛИ МАНТИССУ КАК BCD-ФОРМАТ;---- ФОРМИРОВАНИЕ ТЕКСТА ИЗ BCD-МАНТИСС ----      LEA       EDI,[ESI]+1     ;АДРЕС ОТВЕТА ПОСЛЕ ЗНАКА      FLD1                      ;ЗАГРУЗИЛИ КОНСТАНТУ 1E0      SUB       CL,7            ;ЗАДАННАЯ ДЛИНА МАНТИСС В ОТВЕТЕ      FLD64     [EBX]           ;ОПЯТЬ ЗАГРУЗИЛИ ИСХОДНОЕ ЧИСЛО      LEA       ESI,[ESP]+10    ;АДРЕС ПЕРВХ ЦИФР МАНТИСС      STD      MOV       DH,0            ;ПОКА НУЛИ НЕ ПИШЕМ      FABS                      ;ЗНАК ИСХОДНОГО ЧИСЛА БОЛЬШЕ НЕ НУЖЕН      MOV       AH,0            ;НАЧИНАЕМ С ЧЕТНОЙ ТЕТРАД      FCOM                      ;ЗАРАНЕЕ СРАВНИВАЕМ ЧИСЛО С 1E0      MOV       BL,1            ;ПОКА ОНО МОЖЕТ БТЬ ВБЛИЗИ +1/-1;---- ЦИКЛ ПЕРЕПИСИ ПО ВСЕМ ТЕТРАДАМ BCD-МАНТИСС ----M0004:XOR       AH,1            ;ОЧЕРЕДНАЯ ТЕТРАДА      JZ        @               ;ЕСЛИ НЕЧЕТНАЯ ТЕТРАДА      LODSB                     ;ДОСТАЛИ БАЙТ МАНТИСС      CMP       ESI,ESP         ;ЗА ПРЕДЕЛАМИ МАНТИСС ?      MOV       DL,AL           ;ДВЕ ТЕТРАД МАНТИСС      JB        M0007           ;ЗАКОНЧИЛИ ВВОД;---- ЧЕТНАЯ ТЕТРАДА ----      SHR       AL,4            ;ЧЕТНАЯ ТЕТРАДА BCD      JMPS      M0005;---- НЕЧЕТНАЯ ТЕТРАДА ----@:    MOV       AL,DL      AND       AL,0FH          ;НЕЧЕТНАЯ ТЕТРАДА BCD;---- ПРОПУСК ЛИДИРУЮЩИХ НУЛЕЙ МАНТИСС ----M0005:OR        DH,AL           ;ЕЩЕ ИДУТ НУЛИ МАНТИСС ?      JNZ       @               ;УЖЕ ИДУТ ЦИФР МАНТИСС      INC       ECX             ;НЕ УЧИТВАЕМ ЭТОТ НОЛЬ В МАНТИССЕ      JMPS      M0006           ;ПРОПУСКАЕМ НЕЗНАЧАЩИЙ НОЛЬ;---- ПРОВЕРКА НА ВСЕ ДЕВЯТКИ (Т.Е. НА ЧИСЛО ВБЛИЗИ +1/-1) ----@:    CMP       AL,9            ;ИДУТ СПЛОШНЕ ДЕВЯТКИ ?      JZ        @               ;ДА, ПРОДОЛЖАЮТСЯ ДЕВЯТКИ      MOV       BL,0            ;УЖЕ НЕ ВБЛИЗИ +1/-1 (НЕ ДЕВЯТКА);---- ПРОПУСК ТОЧКИ, ЗАРАНЕЕ ЗАПИСАННОЙ В ОТВЕТ ----@:    CMP       B PTR [EDI],'.' ;ЗДЕСЬ В ШАБЛОНЕ ТОЧКА ?      JNZ       @      INC       EDI             ;ПРОПУСКАЕМ ТОЧКУ;---- ЗАПИСЬ ОЧЕРЕДНОЙ ЦИФР МАНТИСС КАК ТЕКСТА ----@:    ADD       [EDI],AL        ;ПИШЕМ ОЧЕРЕДНУЮ ЦИФРУ В ОТВЕТ      INC       EDI             ;СЛЕДУЮЩИЙ АДРЕСM0006:LOOP     M0004            ;ЗА СЛЕДУЮЩЕЙ ТЕТРАДОЙM0007:CLD;---- ФОРМИРОВАНИЕ ВЕЛИЧИН ПОРЯДКА ----      MOV       ESI,OFFSET ДЕСЯТЬ      FNSTSW    AX      XOR       EDX,EDX         ;ПОКА ПОРЯДОК НОЛЬ      SAHF      JZ        M0011           ;ЧИСЛО СТРОГО РАВНО 1 - ПОРЯДОК НОЛЬ      JA        M0009           ;ЧИСЛО БОЛЬШЕ 1 - ПОРЯДОК ПОЛОЖИТЕЛЕН      MOV       B PTR [EBP]-1,'-' ;ОТМЕТИЛИ ОТРИЦАТЕЛЬНЙ ПОРЯДОК;---- УВЕЛИЧЕНИЕ ПОРЯДКА ДО ЧИСЛА БОЛЬШЕ 1 ----M0008:FIMUL16   [ESI]           ;УВЕЛИЧИЛИ ЧИСЛО В 10 РАЗ      INC       EDX             ;УВЕЛИЧИЛИ ПОРЯДОК      FCOM                      ;СРАВНИВАЕМ С 1      FNSTSW    AX              ;ЗАПОМНИЛИ РЕЗУЛЬТАТ СРАВНЕНИЯ      SAHF      JZ        M0011           ;СТРОГО РАВНО 1 - НАШЛИ ПОРЯДОК      FLD       ST              ;РАЗМНОЖИЛИ ЧИСЛО      FSUB      ST,ST2          ;РАЗНИЦА С 1      FXTRACT                   ;ДОСТАЛИ МАНТИССУ И ПОРЯДОК      FSTP      ST              ;ВБРОСИЛИ МАНТИССУ      FIST32P   [ESP]           ;ПОРЯДОК РАЗНИЦ      CMP       D PTR [ESP],-53 ;УЖЕ ВПЛОТНУЮ К 1 ?      JLE       M0010           ;ДА, ВПЛОТНУЮ, БЛИЖЕ НЕ БУДЕТ      SAHF                      ;ОПЯТЬ ДОСТАЛИ ФЛАГИ СРАВНЕНИЯ С 1      JA        M0011           ;УЖЕ БОЛЬШЕ - НАШЛИ ПОРЯДОК      JMPS      M0008           ;ПРОДОЛЖАЕМ УМНОЖАТЬ НА 10;---- УМЕНЬШЕНИЕ ПОРЯДКА ДО ЧИСЛА МЕНЬШЕ 1 ----M0009:FIDIV16   [ESI]           ;УМЕНЬШИЛИ ЧИСЛО В 10 РАЗ      INC       EDX             ;УМЕНЬШИЛИ ПОРЯДОК      FCOM                      ;СРАВНИВАЕМ С 1      FNSTSW    AX              ;ЗАПОМНИЛИ РЕЗУЛЬТАТ СРАВНЕНИЯ      SAHF      JZ        M0011           ;СТРОГО РАВНО 1 - НАШЛИ ПОРЯДОК      FLD       ST              ;РАЗМНОЖИЛИ ЧИСЛО      FSUB      ST,ST2          ;РАЗНИЦА С 1      FXTRACT                   ;ДОСТАЛИ МАНТИССУ И ПОРЯДОК      FSTP      ST              ;ВБРОСИЛИ МАНТИССУ      FIST32P   [ESP]           ;ПОРЯДОК РАЗНИЦ      CMP       D PTR [ESP],-53 ;УЖЕ ВПЛОТНУЮ К 1 ?      JG        @               ;ЕЩЕ НЕ ВПЛОТНУЮM0010:INC       EBX             ;ПРИЗНАК НАХОЖДЕНИЯ ВБЛИЗИ 1      JMPS      M0011           ;ПЕРЕСТАЕМ ИСКАТЬ ПОРЯДОК@:    SAHF                      ;ОПЯТЬ ЗАГРУЗИЛИ ФЛАГИ СРАВНЕНИЯ С 1      JNB       M0009           ;ПРОДОЛЖАЕМ ДЕЛИТЬ НА 10      DEC       EDX             ;ЧИСЛО В ОТВЕТЕ ДОЛЖНО БТЬ БОЛЬШЕ 1M0011:ADD       ESP,11          ;ОСВОБОДИЛИ СТЕК ОТ BCD-МАНТИСС;---- КОРРЕКТИРОВКА ПОРЯДКА ВБЛИЗИ +/-1 ----      CMP       BL,2            ;БЛИ ВСЕ ДЕВЯТКИ И ПОЧТИ 1 ?      XCHG      EAX,EDX         ;ДОСТАЛИ ЗНАЧЕНИЕ ПОРЯДКА      JNZ       @               ;НЕТ, ЧИСЛО НЕ ВБЛИЗИ 1      DEC       EAX             ;ВБЛИЗИ 1 СВЕРХУ, ДЕЛАЕМ 0.999...E+000      CMP       B PTR [EBP]-1,'-' ;ЧИСЛО МЕНЬШЕ 1E0 ?      JNZ       @               ;НЕТ, БОЛЬШЕ      INC       EAX             ;ВЕРНУЛИ ПОРЯДОК      INC       EAX             ;ВБЛИЗИ 1 СНИЗУ,  ДЕЛАЕМ 9.999...E-001;---- ЗАПИСЬ СТАРШЕЙ ЦИФР ПОРЯДКА ----@:    PUSH      100      XOR       EDX,EDX      POP       EBX             ;ДЕЛИМ НА КОНСТАНТУ 100      FSTP      ST              ;ОЧИЩАЕМ FPU ОТ ПОИСКА ПОРЯДКА      DIV       EBX             ;ПОЛУЧАЕМ ЧАСТНОЕ - ПЕРВУЮ ЦИФРУ      ADD       [EBP],AL        ;СТАРШАЯ ЦИФРА ПОРЯДКА;---- ЗАПИСЬ ДВУХ МЛАДШИХ ЦИФР ПОРЯДКА ----      MOV       BL,10           ;ДЕЛИМ НА КОНСТАНТУ 10      XCHG      EAX,EDX         ;ОСТАТОК ОТ ДЕЛЕНИЯ НА 100      DIV       BL              ;ЧАСТНОЕ И ОСТАТОК - ДВЕ ЦИФР ПОРЯДКА      FSTP      ST              ;ВБРОСИЛИ КОНСТАНТУ 1 ИЗ FPU      ADD       [EBP]+1,AX      ;ДВЕ ОСТАЛЬНЕ ЦИФР ПОРЯДКА;---- ВХОД С ОТВЕТОМ В СТЕКЕ И ДЛИНОЙ СТРОКИ В AL ----      POP       EAX             ;ДЛИНА СТРОКИ ОТВЕТА В СТЕКЕ      RET      DSEGДЕСЯТЬ DW 10                    ;БАЗА ДЛЯ ПЕРЕВОДА В ДЕСЯТИЧНУЮ СИСТЕМУ

Заключение

Использование команд FPU сделало данную реализацию довольно компактной (326 байт команд) и удовлетворительной по скорости. Например, на компьютере с процессором Intel Core Solo 1.33 GHz сто миллионов преобразований числа 1234.567890 в текст заняли 89 секунд.

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

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

Подробнее..

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

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



Друзья, для вашего удобства статья также доступна в формате презентации на видео (почти 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.002060531.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, потому что другие аргументы будут давать результат, отличающийся только значением в поле экспоненты и интереса не представляют. Сами понимаете:

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



Написав такую программу (скрытый текст будет ниже), я проверил функцию 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 Makers 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) нули не превратятся в единицы, или наоборот. Как только всё стабилизировалось, и алгоритм вышел на такую точность, которая находится за пределами паровоза, тогда мы можем округлять так, как если бы у нас было бесконечное число битов. И это округление будет с корректным последним битом. Но как этого добиться?

Костыли


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

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


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

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

Это было научно-популярное изложение, всё то же самое с теоремой Зива (Zivs 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 Lefvre, Jean-Michel Muller, Toward Correctly Rounded Transcendentals, IEEE TRANSACTIONS ON COMPUTERS, VOL. 47, NO. 11, NOVEMBER 1998. pp. 1235-1243

[5] Vincent Lefvre. New Results on the Distance Between a Segment and Z2. Application to the Exact Rounding. 17th IEEE Symposium on Computer Arithmetic Arith17, Jun 2005, Cape Cod, MA,
United States. pp.68-75

[6] Vincent Lefvre, 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) n4044 Novembre 2000 19 pages.
Подробнее..

Категории

Последние комментарии

  • Имя: Макс
    24.08.2022 | 11:28
    Я разраб в IT компании, работаю на арбитражную команду. Мы работаем с приламы и сайтами, при работе замечаются постоянные баны и лаги. Пацаны посоветовали сервис по анализу исходного кода,https://app Подробнее..
  • Имя: 9055410337
    20.08.2022 | 17:41
    поможем пишите в телеграм Подробнее..
  • Имя: sabbat
    17.08.2022 | 20:42
    Охренеть.. это просто шикарная статья, феноменально круто. Большое спасибо за разбор! Надеюсь как-нибудь с тобой связаться для обсуждений чего-либо) Подробнее..
  • Имя: Мария
    09.08.2022 | 14:44
    Добрый день. Если обладаете такой информацией, то подскажите, пожалуйста, где можно найти много-много материала по Yggdrasil и его уязвимостях для написания диплома? Благодарю. Подробнее..
© 2006-2024, personeltest.ru