В этом разделе рассматриваются несколько известных и/или интересных примеров алгоритмов вычислений с плавающей точкой для демонстрации того, каким образом может быть применено понимание чисел с плавающей точкой.
Не все из этих алгоритмов на данный момент практичны для своей изначальной задачи, но они имеют высокую иллюстративную ценность, и, кроме того, рассматриваемые подходы могут быть полезны для решения других задач.
Алгоритм суммирования Кэхена
Достаточно известный в вычислительной математике алгоритм, авторство которого обычно приписывают Уильяму Кэхену, ссылаясь на статью Pracniques: further remarks on reducing truncation errors.
Цель алгоритма - вычисление суммы последовательности с более высокой точностью, чем в результате прямолинейного алгоритма.
На C++ для float алгоритм может выглядеть примерно так:
float KahanSum(constfloat *x,unsignedint n)
{
float sum=0.0f;
float y,t; // Temporary values.float c=0.0f; // A running compensation for lost low-order bits.for(unsignedint i=0;i<n;++i)
{
y=x[i]-c; // So far, so good: c is zero.
t=sum+y; // Alas, sum is big, y small, so low-order digits of y are lost.
c=(t-sum)-y; // (t - sum) recovers the high-order part of y; subtracting y recovers -(low part of y)
sum=t; // Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
}
return sum;
}
По сути, сумма разбита на 2 части - старшую (sum) и младшую (c). Если sum существенно больше y по абсолютной величине, то t-sum - точная операция (вычитаются близкие числа), как и (t-sum)-y, что и позволяет восстановить младшие биты y.
Если относительная погрешность округления не превышает \(\varepsilon\), то максимальная ошибка для прямолинейного суммирования имеет вид
Абсолютная: \(| E_n | \leq O (n \varepsilon) \sum_i | x_i |\)
Относительная: \(\frac{| E_n |}{| S_n |} \leq O (n \varepsilon) \cdot \frac{\sum_i | x_i |}{| \sum_i x_i |}\)
В среднем их значения:
Абсолютная: \(| E_n | \sim O (\sqrt{n} \varepsilon) \sum_i | x_i |\)
Относительная: \(\frac{| E_n |}{| S_n |} \sim O (\sqrt{n} \varepsilon) \cdot \frac{\sum_i | x_i |}{| \sum_i x_i |}\)
Для достаточно малых n (меньше \(\frac{1}{\varepsilon}\), что для float (т. е. binary32) составляет ~107) слагаемым \(O (n \varepsilon^2)\) обычно можно пренебречь, и таким образом, относительная погрешность метода константа порядка \(\varepsilon\).
Та же идея находит более широкое применение.
Double-double
Родственным приёмом является использование 2-х (или более) чисел базового типа для хранения значения с повышенной точностью. См. напр. https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_… le_arithmetic . Это существенно повышает точность (~106 бит для double-double), но диапазон экспонент остаётся тем же, что и у базового типа.
Double-double иногда используется вместо чисел четверной точности (quadruple-precision, binary128), в случае отсутствия их аппаратной поддержки (и наличия поддержки double).
Также это довольно популярный способ хранения констант повышенной точности в математических библиотеках.
Range reduction
Приведение аргумента функции к диапазону (range reduction) - частая операция при вычислении математических функций, как периодических (напр. sin), так и других (напр. exp).
Вообще говоря, аккуратная range reduction - обширная тема (см. напр. http://www.imada.sdu.dk/~kornerup/papers/RR2.pdf ). Примеры здесь подобраны смыслом напоминающими суммирование Кэхэна.
Данный метод был достаточно популярен (по крайней мере среди разработчиков игр) некоторое время назад, но сейчас менее известен.
Постановка задачи: дано число с плавающей точкой в относительно небольшом диапазоне. Требуется быстро (быстрее, чем встроенное преобразование) привести его к целому (не слишком заботясь о направлении округления).
Хорошее изложение данного метода можно найти в http://kunaifusu.livejournal.com/221098.html .
Соответствующий стандарту код имеет вид:
Этот код - для little-endian архитектуры, его применимость для big-endian - отдельный вопрос.
Исходная идея метода несложная - прибавить к числу такую константу, чтобы 1 ULP суммы стал равен единице, тогда младшие биты окажутся в нужных (тех же, что и в целом числе) позициях, а старшие можно занулить.
Для положительных чисел естественным выбором константы является 223. При этом для всех целых исходных в диапазоне [0;223) хвостовые биты мантиссы суммы будут в точности битами искомого целого. Старшие биты (хранящие экспоненту) можно занулить, например, битовой маской.
Однако, для отрицательных чисел эта константа оказывается неподходящей - сумма будет меньше 223, экспонента - меньше 23 и тогда 1 ULP окажется меньше 1. Из этого есть чрезвычайно элегантный выход - в качестве константы выбирается 1.5·223, а вместо зануления маской используется вычитание целого числа с тем же битовым представлением, что и используемая константа. Легко видеть, что если оба числа (константа и сумма) имеют ULP равный единице, то разность их битовых представлений (интерпретированных как целые; разность понимается в смысле 2-complement) будет в точности равна их разности (как чисел с плавающей точкой), т. е. округленному значению исходного числа. Диапазон значений, на которых эта функция выдаёт правильный ответ: [−222;+222] (другими словами [−4194304;+4194304]).
Естественно, метод с вычитанием можно использовать и для положительной (беззнаковой) версии функции. Благодаря счастливому стечению обстоятельств, его использование позволяет расширить диапазон с [0;223) до [0;223] (т. е. добавляет 223).
Стоит отметить, что (благодаря вычитанию) константа 1.5·223 не является, в каком-то смысле особенной. Любое число с плавающей точкой в диапазоне [223; 224) также подходит, для достаточно малых по абсолютной величине входных значений. Данная конкретная константа обеспечивает максимальный симметричный диапазон.
В каком же направлении округляет предложенный метод? Сумма округляется в выбранном направлении округления, но сумма (для правильного входного диапазона) - всегда положительна. Поэтому, фактически имеют место 3 возможных режима округления - вниз (floor), вверх (ceil) и к ближайшему в смысле roundTiesToEven (режим по умолчанию).
Стоит отметить, что для изначальной задачи метод не особенно актуален, в силу наличия достаточно быстрой инструкции cvttss2si, которую компилятор радостно использует (в т. ч. и в режиме компиляции x87, т. к инструкция (в соответствующей версии) не требует XMM-регистра).
Однако сам подход всё-равно может представлять ценность; в частности - он может быть использован для решения близких по духу задач (например приведение к fixed-point, или приведение вещественного числа диапазона [0;1) (или [0;1] ) к целому диапазона [0;255]).
Как видно, код не содержит ветвлений, и прекрасно поддаётся векторизации.
Краткая версия: оригинальным автором кода является, вероятно, Грег Уолш (Greg Walsh), а популярность этот метод получил благодаря использованию его в исходном коде Quake3.
Благодаря Джону Кармаку, код Quake3 находится в открытом доступе, и читатель может самостоятельно прочитать соответствующий отрывок.
Интересующая функция находится в файле quake3-1.32b\code\game\q_math.c . Её код, в слегка изменённом виде, приводится ниже:
float Q_rsqrt(float number )
{
long i;
float x2, y;
constfloat threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = * (long * ) &y;
i = 0x5f3759df - ( i >>1);
y = * (float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ));
// y = y * ( threehalfs - ( x2 * y * y ) );return y;
}
Код, как можно видеть, приводит к undefined behavior (нарушает strict-aliasing), но достаточно прямолинейный компилятор производит ожидаемый ассемблерный вывод:
+ Показать
− Скрыть
Примечание: изначально Quake3 компилировался в Visual Studio.
Ниже приводится результат для gcc ( x86 gcc 5.2.0 (g++ (GCC-explorer-build) 5.2.0), -std=c++11 -O2 -m32 ):
У кода есть и другие проблемы, кроме strict-aliasing. Например на системах, где sizeof(long)==8 (напр. 64-битные версии Linux) он может выдавать некорректные результаты:
http://rextester.com/SBTJM71863 (x87)
http://rextester.com/LKA47580 (SSE)
Метод состоит в получении начального приближения с последующим его улучшением методом Ньютона. Метод Ньютона (также называемый методом Ньютона — Рафсона) хорошо известен в вычислительной математике и не является предметом данной статьи. Желающие могут прочитать статью на Википедии по приведённой выше ссылке.
Таким образом, интерес представляет выбор начального приближения.
Здесь нам на руку играют несколько счастливых стечений обстоятельств.
Рассмотрим, что происходит с положительным нормализованным числом с плавающей точкой:
\(x=2^{E-bias} \cdot 1.M=2^{a} \cdot (1+b),\; 0 \leq b < 1\)
если его биты сдвинуть вправо на 1.
При этом знак остаётся положительным и ведущий (неявный) бит мантиссы остаётся равным 1 (если только результат не оказался денормализованным, каковой случай мы пока что не рассматриваем). Результатом будет число
\(x'=2^{E'-bias} \cdot 1.M'=2^{a'} \cdot (1+b'),\; 0 \leq b' < 1\)
где
\(E' = \left \lfloor \frac{E}{2} \right \rfloor\) \(M' = 2^{p-1} \cdot (E \:\mathrm{mod}\: 2) + \left \lfloor \frac{M}{2} \right \rfloor\)
таким образом:
\(a' = \left \lfloor \frac{a+bias}{2} \right \rfloor - bias\) \(b' = 0.5 \cdot (E \:\mathrm{mod}\: 2) + \frac{b}{2} + \epsilon, | \epsilon | < 2^{-p}\)
Здесь \(\epsilon\) соответствует младшему биту \(M\).
Явно рассматривая чётные и нечётные \(a\) (нечётные и чётные \(E\) соответственно, т. к. \(bias\) - нечётный):
\(x'=
\left\{\begin{matrix}
2^{a/2+bias/2-1/2} \cdot (1 + 0.5 + b/2 + \epsilon), \, a = 2 n,
\\
2^{a/2+bias/2} \cdot (1 + b/2 + \epsilon), \, a = 2 n + 1.
\end{matrix}\right.\)
(\(bias/2\) - не является целым числом).
Отсюда легко видеть, что \(2^{-bias/2} \cdot x'\) - является довольно неплохим грубым приближением \(\sqrt{x}\). Действительно, выражение \((2^{-bias/2} \cdot x')^2\) равно:
для чётных \(a\):
\(2^{a} \cdot \frac{1}{2} (1.5 + b/2 + \epsilon)^2=2^{a} \cdot (1.125 + 0.75 b + 0.5 \frac{b^2}{4} + O(\epsilon))=x+2^{a} \cdot (0.125 - 0.25 b + 0.5 \frac{b^2}{4} + O(\epsilon))\),
для нечётных \(a\):
\(2^{a} \cdot (1 + b/2 + \epsilon)^2=2^{a} \cdot (1 + b + \frac{b^2}{4} + O(\epsilon))=x+2^{a} \cdot (\frac{b^2}{4} + O(\epsilon))\).
Здесь стоит отметить несколько удачных стечения обстоятельств:
1. Ведущая единица в мантиссе не "портится", т. к. явно не хранится.
2. Сдвиг вправо на 1 бит соответствует делению на 2, и таким образом, для показателя степени - грубо соответствует извлечению квадратного корня.
3. Младший бит экспоненты "вдвигается" в хвостовые биты мантиссы, что (так уж совпало) идёт точности на пользу.
Нас интересует, однако \(1 / \sqrt{x}\).
Т. к. вычитание экспонент соответствует делению, то возникает идея вычесть битовое представление (сдвинутое вправо на 1 бит) из некоторой константы.
Очевидным выбором является константа с экспонентой, соответствующей \(2^{\lfloor bias/2 \rfloor}=2^{bias/2-1/2}\) (численно: E=190) и хвостом мантиссы из всех 1. Для неё при вычитании в мантиссе не происходит переноса, а вычитание экспонент обеспечивает деление с довольно близким к требуемому результатом.
Возникает вопрос: а нельзя ли подобрать константу таким образом, чтобы получить более высокую точность? Как видно из указанного выше - перенос между мантиссой и экспонентой может идти точности не во вред, а на пользу.
Это оказывается возможным. График максимальной (на всех допустимых аргументах) относительной погрешности такого приближения, как функция от константы получается довольно плавным, с единственным минимумом.
Примечание: здесь предполагается использование максимального относительного отклонения в качестве метрики. При желании, можно выбрать и другую: например, среднеквадратичное отклонение (математически, обе эти метрики можно считать частным случаем L-нормы: L∞ и L2 соответственно). Однако обычно интересующей метрикой является именно L∞.
Можно, однако, оптимизировать не начальное приближение, а конечный результат функции. Оптимальная константа в этом случае получится несколько другая.
Графики относительной ошибки функции и её (десятичного) логарифма, в зависимости от константы, приводятся ниже:
Стоит отметить, что константа 0x5f3759df, используемая в коде функции не является оптимальной (ни в смысле начального приближения, ни в смысле результата функции), хотя и близка к оптимальной для результата функции (0x5f375a86).
Данный метод легко модифицируется для других типов(напр. double).
Рассмотрев механизм работы функции, перейдём к её практическим свойствам.
Максимальная относительная погрешность функции на всех допустимых (подробнее - см. ниже) входных значениях не превышает 0.00176, что соответствует ~11 битам точности.
Функция корректно возвращает NaN для NaN аргумента.
Для отрицательных чисел функция выдаёт либо +∞ либо малые конечные отрицательные числа (математически \(1 / \sqrt{x}\) для этих аргументов не определён (в действительных числах), поэтому разумным ответом был бы NaN).
Q_rsqrt(1.0f)=0.998307168f , при истинном ответе 1.0, т. е. единица не сохраняется (тем же недостатком обладает и инструкция rsqrtss).
Q_rsqrt(+0.0f)=1.98177537e+19 ("штатное" значение +∞). Это свойство, впрочем может быть удобным, т. к. позволяет вычислять \(\sqrt{x}\) как \(x \cdot \frac{1}{\sqrt{x}}\) без проблем (в случае \(\frac{1}{\sqrt{x}}=+\infty\) второе выражение дало бы NaN). Это имеет смысл, т. к. вычисление функции \(1 / \sqrt{x}\) часто быстрее (для схожей точности), чем \(\sqrt{x}\). Это связано со свойствами функций \(\sqrt{x}\) и \(1 / \sqrt{x}\), важными для алгоритма Ньютона.
Q_rsqrt(-0.0f)=-5.82391438e-20. Штатный способ вычисления (1.0f/sqrtf(-0.0f)) даст −∞ (т. к. по стандарту \(\sqrt{-0}=-0\)), что может в некоторых ситуациях приводить к неожиданным неприятностям.
Q_rsqrt(+INFINITY)=-INFINITY (штатный ответ: +0.0) что сложно назвать удачным.
Q_rsqrt(-INFINITY)=+INFINITY (штатный ответ: NaN) что тоже имеет сомнительную ценность.
Для денормализованных чисел результат может отличаться от правильного во много раз (код всё-равно предполагает ведущий бит мантиссы единичным), но хотя бы примерно лежит в разумной области.
Для всех же конечных нормализованных чисел с плавающей точкой функция обеспечивает указанную выше точность.
При наличии (и возможности использовать) SSE ценность данного метода для изначальной задачи сомнительна - в SSE существует встроенная инструкция для быстрого приблизительного вычисления обратного квадратного корня: rsqrtss (и её векторная форма rsqrtps). Её заявленная точность составляет 12 бит, и некоторые процессоры могут вычислять её и с более высокой точностью. По скорости она в несколько раз превышает Q_rsqrt. В качестве примера ниже приводится без особых комментариев таблица из http://assemblyrequired.crashworks.org/timing-square-root/ :
Method
Total time
Time per float
Avg Error
Carmack’s Magic Number rsqrt
59.378ms
3.54ns
0.0990%
SSE rsqrtss
14.202ms
0.85ns
0.0094%
SSE rsqrtss with one NR step
45.952ms
2.74ns
0.0000%
SSE позволяет также вычислять несколько обратных корней одновременно (rsqrtps). Впрочем, Q_rsqrt тоже можно преобразовать к векторной форме.
Стоит отметить, что т. к. rsqrtss является приближённой, то по умолчанию компилятор не использует её для обычного кода (1.0f/sqrtf(x)). При этом, например, x*Q_rsqrt(x) всё ещё заметно быстрее, чем ssqrts (штатная (правильно округленная) инструкция вычисления квадратного корня). Для доступа к rsqrtss/rsqrtps можно использовать, например, ассемблерные вставки или интринсики.
Читатель, желающий проверить своё умение обращаться с числами с плавающей точкой, но не собирающийся придумывать задачи на свою голову, может воспользоваться задачами, предложенными ниже.
Задача 1. Реализовать функцию, проверяющую, что для всех числел типа binary32 две заданне функции дают одинаковые результаты, со следующим заголовком:
bool equal(float(*f0)(float),float(*f1)(float));
Как зависит ответ от определения понятия "одинаковые"?
Задача 2. (автор: Уильям Кэхэн). Дано число x базового типа (для определённости - binary64) в "разумном" диапазоне (для определённости - \([1;\,1000]\)). Написать функцию, возвращающую 2 числа, a и b, того же типа, такие, что:
1. a - близко к x (в пределах нескольких ulp, в идеале - ближайшее возможное).
2. b в точности (математически) равно 10*a.
Задача 3. Реализовать функцию ulp с заголовком:
float ulp(float);
и поведением, идентичным функции Math.ulp(float) математической библиотеки языка Java.
