Испытания 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:

$\frac{-\frac{481959816488503 x^{11}}{363275871831577908403200}+\frac{23704595077729 x^9}{42339845201815607040}-\frac{2933434889971 x^7}{33603051747472704}+\frac{3605886663403 x^5}{617703157122660}-\frac{109061004303 x^3}{722459832892}+x}{\frac{37291724011 x^{10}}{11008359752472057830400}+\frac{3924840709 x^8}{2016183104848362240}+\frac{101555058991 x^6}{168015258737363520}+\frac{1679739379 x^4}{13726736824948}+\frac{34046903537 x^2}{2167379498676}+1}$


Непосредственно для вычислений обычно используют схему Горнера, чтобы сэкономить на вычислениях. В нашем случае (с помощью 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. исходный код для проверки и критики.

© Habrahabr.ru