Испытания 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<float>::quiet_NaN(). К слову, в этой библиотеке также не обнаружились реализации для элементарных математических функций (ну, кроме корня) — это ещё один повод заподозрить неладное.

Тест 1. Умножение комплексных чисел

Преобразование Фурье, вычисляемое с использованием комплексной арифметики, используется везде, где только можно. Поначалу я хотел испытать Posit именно на преобразовании Фурье; но поскольку его точность достаточно сильно зависит от реализации, для корректного тестирования придётся рассмотреть все основные алгоритмы, что несколько затратно по времени; поэтому начать можно и с более простой операции — умножения комплексных чисел.

Если взять некоторый вектор и повернуть его 360 раз на 1°, то по итогу мы должны получить тот же самый исходный вектор. По факту результат будет слегка отличаться из-за накопления погрешностей — и чем большее количество поворотов, тем больше будет погрешность. Итак, используя вот такой простой код

complex<T> rot(cos(a), sin(a));
complex<T> 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. Возьмём ещё бо́льшую длину — 10¹⁰, чтобы в полной мере оценить преимущества форматов с плавающей точкой:

Здесь самое интересное в начале, на 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. исходный код для проверки и критики.

P.P.S. Продолжение.