Испытания Posit по-взрослому
На Хабре уже было несколько статей (раз, два, два с половиной), посвящённых новому формату чисел с плавающей запятой Posit, авторы которого преподносят его его как превосходящий стандартный IEEE 754 float по всем параметрам. У нового формата нашлись и критики (раз, два) утверждающих, что недостатки Posit перевешивают его достоинства. Но что, если у нас действительно появился новый революционный формат, а критика просто вызвана завистью и некомпетентностью критикующих? Что же, лучший способ выяснить это — взять и повычислять самостоятельно.
Введение
Достоинства нового формата демонстрируются на простых примерах со сложением/умножением чисел схожего порядка, в результате которых получается выше точность на один-два знака выше. Однако в реальных вычислениях плюс/минус один разряд погрешности единичной операции не имеет значения, поскольку мы так и так вычисляем с ограниченной точностью. Имеет значение накопление погрешности в результате последовательности операций, когда спустя некоторое время погрешности в младших разрядах начинают приводить к погрешности в старших. Это мы и попробуем испытать.
Подготовка
Для испытаний я взял реализацию Posit с пафосным названием отсюда. Чтобы она откомпилировалась в Visual Studio, пришлось в файле util.h добавить строчку #define CLZ (n) __lzcnt (n) и в файле posit.cpp заменить 0.f / 0.f на std: numeric_limits: quiet_NaN (). К слову, в этой библиотеке также не обнаружились реализации для элементарных математических функций (ну, кроме корня) — это ещё один повод заподозрить неладное.
Тест 1. Умножение комплексных чисел
Преобразование Фурье, вычисляемое с использованием комплексной арифметики, используется везде, где только можно. Поначалу я хотел испытать Posit именно на преобразовании Фурье;, но поскольку его точность достаточно сильно зависит от реализации, для корректного тестирования придётся рассмотреть все основные алгоритмы, что несколько затратно по времени; поэтому начать можно и с более простой операции — умножения комплексных чисел.
Если взять некоторый вектор и повернуть его 360 раз на 1°, то по итогу мы должны получить тот же самый исходный вектор. По факту результат будет слегка отличаться из-за накопления погрешностей — и чем большее количество поворотов, тем больше будет погрешность. Итак, используя вот такой простой код
complex rot(cos(a), sin(a));
complex vec(length, 0);
for (long i = 0; i < count; i++)
{
vec *= rot;
}
cout << "error: " << stdev(vec.real() - length, vec.imag()) << endl;
повращаем вектор с разными типами данных, а ошибку будем считать как среднее квадратическое отклонение результирующего вектора от исходного (что также можно интерпретировать как длину разностного вектора).
Для начала возьмём единичный вектор, как наиболее благосклонный к Posit:
Здесь пока не видно явного лидера — преимущество в два раз то у одного, то у другого. Увеличим длину вращаемого вектора до 1000:
Значения погрешностей практически сравнялись. Идём дальше — 1000000:
Вот тут Posit уже уверенно отстаёт, а погрешность double начинает заползать во float. Возьмём ещё бо́льшую длину — 1010, чтобы в полной мере оценить преимущества форматов с плавающей точкой:
Здесь самое интересное в начале, на 4-х итерациях — когда float даёт погрешность соизмеримую с double, а у Posit-а уже полностью некорректный результат.
Тест 2. Вычисление рационального полинома
Так как в оригинальной библиотеке реализации математических функций не оказалось, попробуем сделать что-нибудь самостоятельно. Многие функции плохо приближаются разложением в ряд Тейлора, и их удобнее вычислять через приближение рациональным полиномом. Приближение это можно получить разными способами, в том числе и чисто аналитически — через аппроксимацию Паде. Его мы и возьмем для испытаний, причём с достаточно большими коэффициентами, чтобы они также подверглись округлению перед вычислением.
Используя Wolfram Mathematica и команду PadeApproximant[Sin[x], {x, 0, {11, 11}}]
получим вот такой вот рациональный полином для аппроксимации синуса, обеспечивающий точность double в диапазоне примерно от -2 до 2:
Непосредственно для вычислений обычно используют схему Горнера, чтобы сэкономить на вычислениях. В нашем случае (с помощью HornerForm) это будет выглядеть как
template< typename T >
T padesin(T x) {
T xx = x*x;
return
(x*(T(363275871831577908403200.) +
xx*(-T(54839355237791393068800.) +
xx*(T(2120649063015013090560.) +
xx*(-T(31712777908498486800.) +
xx*(T(203385425766914820.) -
T(481959816488503.) * xx)))))) /
(T(363275871831577908403200.) +
xx*(T(5706623400804924998400.) +
xx*(T(44454031219351353600.) + xx*
(T(219578286347980560.) +
xx*(T(707177798947620.) +
T(1230626892363.) * xx)))));
}
Посмотрим:
Как видно, дела у Posit здесь выглядят плачевно — еле-еле две значащих цифры набирается.
Заключение
К сожалению, чуда не произошло и революция отменяется. Преимущество Posit, демонстрируемое на единичных вычислениях — не более чем фокус, расплатой за которое является катастрофическое снижение точности в «тяжёлых» реальных вычислениях. Единственная причина, по которой имеет смысл использовать Posit вместо IEEE 754 float или fixed point — религиозная. Использование волшебного формата, точность которого обеспечивает святая вера его создателей — может привнести немало чудес в ваши программы!
P.S. исходный код для проверки и критики.