Программирование
C++
Математика
17 сентября

Испытания 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:
итераций 4 100 1000 10000 100000
double 0 0 0 0 0
float 0 0,00000036 0,00001038 0,0001858 0,0001961
posit 0 0,00000073 0,00000534 0,0000411 0,0004468

Здесь пока не видно явного лидера — преимущество в два раз то у одного, то у другого. Увеличим длину вращаемого вектора до 1000:
итераций 4 100 1000 10000 100000
double 0 0 0 0 0
float 0 0,00028 0,0103 0,18 0,19
posit 0 0,00213 0,0188 0,16 2,45

Значения погрешностей практически сравнялись. Идём дальше — 1000000:
итераций 4 100 1000 10000 100000
double 0 0 0,00000002 0,00000042 0,0000036
float 0 0,33 12,0 185,8 198,1
posit 0 8,12 71,0 769,2 10706,8

Вот тут Posit уже уверенно отстаёт, а погрешность double начинает заползать во float. Возьмём ещё бо́льшую длину — 1010, чтобы в полной мере оценить преимущества форматов с плавающей точкой:
итераций 4 100 1000 10000 100000
double 0,00000245 0,00001536 0,0002041 0,0040951 0,03621497
float 0,00000245 6003,8 88111,8 1836254,0 1965083,0
posit 9216,0 1287208,7 14443543,7 202630144,4 1784050328,2

Здесь самое интересное в начале, на 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)))));
}


Посмотрим:
x = 0.5 x = 1 x = 2
sin(x) 0,479425538604203 0,8414709848078965 0,9092974268256817
double 0,479425538604203 0,8414709848078965 0,9092974268256816
float 0,4794255495071411 0,8414709568023682 0,9092974066734314
posit 0,4788961037993431 0,8424437269568443 0,9110429435968399

x = 3 x = 4 x = 5
sin(x) 0,1411200080598672 -0,7568024953079282 -0,9589242746631385
double 0,1411200080585958 -0,7568024960833886 -0,9589243758030122
float 0,1411200165748596 -0,7568024396896362 -0,9589243531227112
posit 0,1444759201258421 -0,7614213190972805 -0,9691629931330681

Как видно, дела у Posit здесь выглядят плачевно — еле-еле две значащих цифры набирается.

Заключение


К сожалению, чуда не произошло и революция отменяется. Преимущество Posit, демонстрируемое на единичных вычислениях — не более чем фокус, расплатой за которое является катастрофическое снижение точности в «тяжёлых» реальных вычислениях. Единственная причина, по которой имеет смысл использовать Posit вместо IEEE 754 float или fixed point — религиозная. Использование волшебного формата, точность которого обеспечивает святая вера его создателей — может привнести немало чудес в ваши программы!

P.S. исходный код для проверки и критики.

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

+89
10,5k 53
Комментарии 83