На Хабре уже было несколько статей (раз, два, два с половиной), посвящённых новому формату чисел с плавающей запятой 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(). К слову, в этой библиотеке также не обнаружились реализации для элементарных математических функций (ну, кроме корня) — это ещё один повод заподозрить неладное.
Преобразование Фурье, вычисляемое с использованием комплексной арифметики, используется везде, где только можно. Поначалу я хотел испытать Posit именно на преобразовании Фурье; но поскольку его точность достаточно сильно зависит от реализации, для корректного тестирования придётся рассмотреть все основные алгоритмы, что несколько затратно по времени; поэтому начать можно и с более простой операции — умножения комплексных чисел.
Если взять некоторый вектор и повернуть его 360 раз на 1°, то по итогу мы должны получить тот же самый исходный вектор. По факту результат будет слегка отличаться из-за накопления погрешностей — и чем большее количество поворотов, тем больше будет погрешность. Итак, используя вот такой простой код
повращаем вектор с разными типами данных, а ошибку будем считать как среднее квадратическое отклонение результирующего вектора от исходного (что также можно интерпретировать как длину разностного вектора).
Для начала возьмём единичный вектор, как наиболее благосклонный к Posit:
Здесь пока не видно явного лидера — преимущество в два раз то у одного, то у другого. Увеличим длину вращаемого вектора до 1000:
Значения погрешностей практически сравнялись. Идём дальше — 1000000:
Вот тут Posit уже уверенно отстаёт, а погрешность double начинает заползать во float. Возьмём ещё бо?льшую длину — 1010, чтобы в полной мере оценить преимущества форматов с плавающей точкой:
Здесь самое интересное в начале, на 4-х итерациях — когда float даёт погрешность соизмеримую с double, а у Posit-а уже полностью некорректный результат.
Так как в оригинальной библиотеке реализации математических функций не оказалось, попробуем сделать что-нибудь самостоятельно. Многие функции плохо приближаются разложением в ряд Тейлора, и их удобнее вычислять через приближение рациональным полиномом. Приближение это можно получить разными способами, в том числе и чисто аналитически — через аппроксимацию Паде. Его мы и возьмем для испытаний, причём с достаточно большими коэффициентами, чтобы они также подверглись округлению перед вычислением.
Используя Wolfram Mathematica и команду PadeApproximant[Sin[x], {x, 0, {11, 11}}]
получим вот такой вот рациональный полином для аппроксимации синуса, обеспечивающий точность double в диапазоне примерно от -2 до 2:
Непосредственно для вычислений обычно используют схему Горнера, чтобы сэкономить на вычислениях. В нашем случае (с помощью HornerForm) это будет выглядеть как
Посмотрим:
Как видно, дела у Posit здесь выглядят плачевно — еле-еле две значащих цифры набирается.
К сожалению, чуда не произошло и революция отменяется. Преимущество Posit, демонстрируемое на единичных вычислениях — не более чем фокус, расплатой за которое является катастрофическое снижение точности в «тяжёлых» реальных вычислениях. Единственная причина, по которой имеет смысл использовать Posit вместо IEEE 754 float или fixed point — религиозная. Использование волшебного формата, точность которого обеспечивает святая вера его создателей — может привнести немало чудес в ваши программы!
P.S. исходный код для проверки и критики.
P.P.S. Продолжение.
Введение
Достоинства нового формата демонстрируются на простых примерах со сложением/умножением чисел схожего порядка, в результате которых получается точность на один-два знака выше. Однако в реальных вычислениях плюс/минус один разряд погрешности единичной операции не имеет значения, поскольку мы так и так вычисляем с ограниченной точностью. Имеет значение накопление погрешности в результате последовательности операций, когда спустя некоторое время погрешности в младших разрядах начинают приводить к погрешности в старших. Это мы и попробуем испытать.
Подготовка
Для испытаний я взял реализацию 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:
Непосредственно для вычислений обычно используют схему Горнера, чтобы сэкономить на вычислениях. В нашем случае (с помощью 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. Продолжение.
0serg
Оба примера — хорошая демонстрация недостатков posit, но оба они искусственно натянуты. Что, к примеру мешает в апроксимации Паде посчитать posit(a/b) вместо того чтобы использовать posit(a)/posit(b) где a и b — числа порядка 10^20?
skrimafonolog
«Что мешает использовать вместо варианта 1 вариант 2» — это тоже натягивание.
0serg
Автор берет вариант где фигурируют 3 операции округления вместо 1, две из которых искусственно подобраны так чтобы быть максимально проблемными для posit. А ведь в исходном варианте посчитанном матлабом подобных проблем не было — простая подстановка коэффициентов прямо по формуле которую выдал матлаб решила бы проблему.
А так да, то что надо остерегаться некоторых граничных случаев — это проблема posit
Pand5461
Если вендоры вдруг поймут всю красоту позитов и уберут аппаратную поддержку IEEE 754, то серьёзный вопрос — а
a / bкаким образом считать?skrimafonolog
Во первых:
«Вдруг» поддержку IEEE 754 никто не уберет.
Вон еще даже COBOL жив (и вычислительные его возможности этому поспособствовали).
Во вторых:
В отсутствии аппаратных возможностей — все прекрасно реализуется через программную эмуляцию. Это широко практикуется еще со времен, когда математический сопроцессор покупался за отдельные деньги (причем включение эмулятора было реализовано незаметным для программы образом, программа думала, что на математическом сопроцессоре работала). И реализовано даже для таких объемных по коду вещей как 3D графика (Mesa), что уж говорить про несколько базовых математических операций.
Собственно автор статьи эту эмуляцию и использовал, только он эмулировал наоборот — новый стандарт.
В общем, описанной вами гипотетической проблемы нет. И не будет.
Pand5461
Ясно, что проблемы не будет, потому что не пойдёт posit в широкие массы (хотя, если присыпать его несколькими триллионами $$...), слишком большой аккуратности требует, хоть авторы формата и утверждают обратное.
ads83
А в чем тут неправильность? Все вычисления проводятся с числами в Posit, разве это плохо?
Я не знаю С++, но кажется все значения приводятся к типу Т. Также смог найти что оператор деления переопределен для Posit.
0serg
posit(a)/posit(b) — это три операции округления (a, b и результата posit(a)/posit(b)).
posit(a/b) — одна.
В данном же случае a и b специально подобраны так чтобы округление для posit(a) было очень неточным.
Refridgerator Автор
Я ничего специально не подбирал, просто взял побольше членов в полиноме.
0serg
Простите, но прямая подстановка коэффициентов из матлабовской формулы дала бы другой результат.
Refridgerator Автор
И вы хотите сказать — просто не используйте схему Горнера при вычислениях с posit? Ну хорошо — раз нельзя использовать схему Горнера, значит возможно, что-то ещё нельзя. Где список всего того, что нельзя — с учётом того, что авторы недостатки своего формата отрицают?
0serg
Проблем со схемой Горнера никаких, там коэффициенты ровно те же самые.
Проблемы у posit никто не скрывает — они начинаются при использовании коэффициентов больше 10^6. Что да, тоже бывает, но которых в данном случае в исходной формуле не было. Домножением числителя и знаменателя в формуле на число с величиной более 10^20 Вы элегантно решили эту проблему. Однако то что в результате вышла формула в которой будут явные проблемы с posit видно невооруженным взглядом :).
Такие полиномы осмыслены только для x << 1. Вводим y = x/10^5 и считаем формулу для y. Результат к слову будет для float/double тоже точнее.
Refridgerator Автор
«Числа posit имеют большую точность, больший динамический диапазон и большее покрытие. Они могут использоваться для получения лучших результатов, чем float той же разрядности, или (что может быть ещё большим конкурентным преимуществом), тех же результатов при меньшей разрядности… Так как они работают как float, а не как интервальная система, они могут рассматриваться как прямая замена float»
Вот и я прямо заменил float на posit.
0serg
Posit МОГУТ давать большую точность? Могут. В ваших же экспериментах это видно. При прямой замене float на posit. ВСЕГДА ЛИ они будут давать большую точность? Очевидно что нет, не всегда и это вроде никто и нигде не отрицает. Какой сценарий чаще встречается для реальных задач? Авторы считают что их.
Berkof
Ой, только не надо этой демагогии тут. Будете вы менять свой бензиновый автомобиль на водородный, если водородный может иногда экономить деньги (если вы ездите в районе, где есть водородные заправки и очень большие налоги на неэкологичные ТС, ездите много и часто въезжаете/выезжаете из города)? Да, может быть кто-то в единичных случаях и будет. А для 99.99% людей хватит стандартного бензинового авто, ну ок, они оплатят катализатор на выхлопе (double вместо float).
Чтобы что-то внедрить быстро — нужно чтобы у него были принципиальные и важные преимущества (как у сотового телефона), чтобы внедрить медленно — нужны заметные, но не критичные преимущества, желательно без недостатков (или с мизерными недостатками) (пример VVTi и прочие фазовращалки в ДВС). А если что-то не лучше существующего и уже распространённого варианта — его внедрят единицы, которым вот так уж получилось именно это сочетание плюсов/минусов поможет… Готов поверить, что этот posit внедрят в какие-то математические библиотеки для отдельных очень тяжёлых алгоритмов рассчёта.
0serg
Имхо posit хорошо пойдет для всяких акселлераторов включая gpu. Для настольных cpu он возможно пойдет как дополнительный набор инструкций. Прямо вот так совсем и везде заменять 754 он конечно не будет. Но широкое применение волна может найти.
Refridgerator Автор
Я попробовал, результаты даже ещё интереснее получились. Пока значения аргумента возле единицы, точность на 1-2 разряда выше float. Как только значение аргумента меняется на порядок — точность резко снижается. При x=10 только 3 правильных цифры, при x=100 только две.
x=1
-------------------------
0.841470984807896 sin(x)
0.841470984807896 double
0.841470956802368 float
0.841470986604690 posit
x=10
-------------------------
-0.642458796040545 double
-0.642459392547607 float
-0.642386730760335 posit
x=100
-------------------------
-35435.826218351600000 double
-35435.824218750000000 float
-35367.474609375000000 posit
x=1000
-------------------------
-391247.500187484000000 double
-391247.500000000000000 float
-265990.968750000000000 posit
0serg
Вот это уже ближе к жизненным проблемам (там x^11 используется что для x>10 дает числа более 10^11)
ads83
Простите, но ведь у вас разложение в полином сделано для интервала [-2..2]? В таком случае нельзя использовать sin(10) для оценки погрешности вычислений: мы не знаем, насколько велика аналитическая погрешность(отклонение sin(10) от pader(10) ).
Вы используете double как точное значение, но ведь у него тоже есть погрешность. Для чистоты эксперимента было бы круто вычислить точное значение pader(10). Правда, не представляю как это сделать :-) Интуиция подсказывает, что скорей всего это не сильно изменит дело, но она часто подводит...
Refridgerator Автор
В данном случае я рассматриваю полином сам по себе, без привязки к синусу. У double безусловно тоже есть погрешность, поэтому предварительно я сверился с Mathematica убедиться, что он выдаёт корректный результат.
ads83
Понятно, спасибо что прояснили.
ads83
На мой дилетантский взгляд, вы не совсем правы.
Во-первых, мы считаем не posit(a)/posit(b), а posit(a)/b_as_posit. У нас одна операция преобразования типа и одно деление. Отдельный вопрос, теряется ли точность при приведении конкретно этих чисел в posit. Возможно, Refridgerator сможет ответить на этот вопрос (равны ли 363275871831577908403200.0 и posit(363275871831577908403200.0), ну и все остальные коэффициенты). То есть не три операции с потерей точности, а одна или две.
Во-вторых, операция деления переопределена для posit-ов. Мы проверяем потерю точности при операциях с такими числами. Если же мы считаем posit(a/b), то узнаем потерю точности перевода в другой формат, не более.
В-третьих, я склонен согласиться с вами: деление двух многочленов высокой степени приведет к гораздо большей потере точности, чем расчет по формуле аппроксимации Паде. Независимо от того, какой формат чисел используется.
0serg
b_as_posit подразумевает округление b к ближайшему posit. Почему мы не должны учитывать эту операцию округления? Она даёт весьма заметный вклад в результат.
То что posit не всегда лучше довольно очевидно. Вопрос в том лучше ли он для практических задач и несколько сложно обойти его проблемы там где он работает хуже.
ads83
Да, вы правы, b_as_posit != b_equal. Но я считаю, что лучше когда мы получили b_as_posit не как округление b_as_double к ближайшему posit, а как выполнение операций над posit-ами. Так как заявленная цель — посчитать накопление ошибки в операциях над posit-ами, то на каждой итерации "округление" posit(b_as_double) нам будет мешать.
0serg
Если мы можем тривиально посчитать posit(a/b), причём ни a, ни b нас не интересуют, то расчёт через posit(a) /posit(b) требующий три операции округления на мой взгляд просто искусственная пессимизация
ads83
На мой взгляд, расчет posit(a/b) привносит ошибки другого формата чисел: на самом деле считается posit(aAsDoubleWithErr/bAsDoubleWithErr) == posit(resultAsDoubleWithMOREErr). И трудно сказать, как связаны ошибки деления double-ов, ошибки конвертации и ошибки posit-ов.
Цель же сравнить одни и те же операции над разными форматами. Поэтому пусть эти вычисления не оптимальны, но они отличаются только типом, мы оставили только один тип ошибок и изучаем только его.
kasyachitche
И какой вывод?
Нельзя делить posit на posit?
0serg
Нельзя использовать posit32 с величиной более ~10^10 для вычислений в которых финальный результат ожидается в районе 1. Желательно при использовании posit масштабировать вычисления к числам в диапазоне 10^-6...10^6.
Refridgerator Автор
Ну почему нельзя-то? Это 257 нельзя в байт запихать, а posit и бо?льшие порядки держит — это же всё-таки формат с плавающей запятой.
0serg
Потому что смысла нет. Гвозди микроскопом тоже можно забивать, но предпочительно использовать инструмент сообразный задаче
Pand5461
На практике это означает, что, если сложная физическая задача разбивается на несколько подзадач, для каждой из них нужно подгонять свою систему измерения, а потом их сопрягать. Казалось бы, при чём тут Mars Climate Orbiter?
picul
Кстати, я когда читал недавнию статью про Posit, тоже был в небольшом недоумении от следующего предложения:
В том смысле, что ничуть не реже такие числа встречаются, а может даже и чаще.Вообще тут можно было бы сказать, что на каждую задачу свой инструмент, только вот IEE 754 реализован в железе практически везде, так что сложно представить задачи, в которых использовать Posit будет выгоднее. Не говоря уже о том, что сам новый формат сильно сложнее, чем классический (ну или мне так показалось), что скажется на скорости аппаратной реализации (если она будет) и, возможно, на простоте использования и отладки.
Regis
picul
А зачем десятичная арифметика в двоичном процессоре? Возможно, оно решит некоторые проблемы вроде 0.300000001, но ведь и диапазон значений очень сильно урезается (насколько я понимаю, предлагается хранить десятичный разряд в 4 битах). Разве что в банковской сфере поможет, или как?
Regis
Все финансовые расчеты, весь бухгалтерский учет по-хорошему должен делаться в десятичной арифметике с использованием правил десятичного округления, а не двоичного. Сейчас ради этого при выполнении финансовых расчетов идут на всяческие ухищрения (используют высокоуровневые эмуляции того или иного рода).
Вообще говоря, двоичные вычисления как раз почти никому не нужны. Просто в подавляющем большинстве задач детали округления и количество значащих цифр в результате не очень существенны, поэтому особой разницы в том, была ли использована двоичная или десятичная арифметика нет и поэтому все привыкли относиться к двоичным вычислениям, которые нам обычно предлагают по-умолчанию, как к чему-то самом собой разумеющемуся или даже "единственно правильному".
Нет, не сильно. IEE 754 decimal64 хранит каждые 3 десятичных разряда в блоках размером 10 бит. 10^3=1000, 2^10=1024, т. е. потеря тут небольшая. См. вики.
agmt
По-хорошему — это в целочисленной арифметике и всячески избегая округлений.
Regis
Как только оказывается, что fixed point вам не подходит (так как в одном типе могут быть и очень большие и очень маленькие значения), то сразу появляется необходимость постоянно контролировать положение точки, следить за переполнением и т.д. В конечном счете получаете в том или ином виде эмуляцию десятичной арифметики с плавающей запятой.
Пример проблемной ситуации: вам хочется в одном типе данных хранить суммы в Биткоинах (в том числе очень маленькие значения, 0.00000001 и меньше) и в Зимбабвийских долларах (очень большие значения, триллиарды и больше). Fixed point с целым числом вам не подойдет — ему не хватит значащих цифр (19 десятичных значащих цифр в эквиваленте для int64).
netch80
> Пример проблемной ситуации: вам хочется в одном типе данных хранить суммы в Биткоинах (в том числе очень маленькие значения, 0.00000001 и меньше) и в Зимбабвийских долларах (очень большие значения, триллиарды и больше).
Тут уже идея хранить в одном типе данных такие разные значения является той исходной диверсией, которую надо избегать всеми средствами. Если нормальный язык, то в нём не должно быть проблемы сделать что-то в духе
template <class CurrencyName, int DigitsAfterDot> class Currency<CurrencyName, DigitsAfterDot> { private: intmax_t mValue; ... }; struct ZimbabweDollar : Currency { }; using ZimbabweCurrency = Currency<ZimbabweDollar, -6>;и не допускать конверсию даже между разными представлениями одной валюты без явного оператора конверсии с заданным округлением и с контролем переполнения.
> следить за переполнением
Лучше и проще следить за явным переполнением, чем за молчаливой потерей точности (если мы всё ещё о финансах, где в бухгалтерии каждая потерянная копейка это основание выть на несхождение баланса).
> и в Зимбабвийских долларах (очень большие значения, триллиарды и больше)
Что, там до 10^21 дошли? :)
netch80
> Вообще говоря, двоичные вычисления как раз почти никому не нужны.
Они нужны ровно там, откуда возникли — из математических расчётов физического (или близкого к нему) происхождения.
Двоичный порядок на длинных расчётах значительно выгоднее десятичного из-за меньших потерь при округлении.
На коротких расчётах при десятичном вводе/выводе это, конечно, незаметно, за счёт проблем переквантования. Но так как, например, АЦП сейчас чаще генерируют двоичный выход, чем какой-то другой, десятичный транспорт всё равно лишний.
Для бухгалтерского учёта плавающая точка тоже плоха — если вы, например, в decimal32 подсчитаете 9999999+4 и молча (ну, inexact flag подымет, если умеет) получите 10000000 вместо 10000003 (потому что в нём ровно 7 десятичных цифр)… лучше бы вылетело грубое переполнение, на нём было бы лучше заметно, что что-то не так. Поэтому там наиболее желательный вариант это фиксированная точка (копейки, сотые доли копеек, где как лучше).
В итоге, десятичная плавучка остаётся только для задачи «нам нужно точно поддерживать decimal(n,m) из SQL, PL/1, etc., но нам лень самим это делать».
> IEEE 754 decimal64 хранит каждые 3 десятичных разряда в блоках размером 10 бит. 10^3=1000, 2^10=1024, т. е. потеря тут небольшая.
Нет, IEEE754 даёт два равноправных варианта рядом — с двоичной и десятичной мантиссой (последняя — через Chen-Ho BCD). Это и по вашей ссылке сказано.
Реализация с Chen-Ho BCD это у IBM? Я их не щупал. Но она не единственная, и для железа IMHO практически всегда удобнее использовать не BCD, а с двоичной мантиссой. Есть какие-то примеры, когда BCD действительно удобнее?
VolCh
Спросите у математиков и физиков, нужны ли им именно двоичные расчёты при прочих равных? Если десятичные или, например, троичные будут обеспечивать ту же точность, скорость, стоимость не будет ли им всё равно что "под капотом"?
netch80
«Если» были они такое обеспечили — да, было бы всё равно.
Проблема в том, что этого «если» не будет. То, что чем меньше основание для порядка, тем меньше ошибки округления, это банальный доказанный факт. Если бы было основание меньше 2, пытались бы его применить, но меньше не бывает (дробные не считаем, реализация усложняется просто чудовищно).
BD9
Двоичные расчёты эффективнее десятичных.
VolCh
На двоичной аппаратной базе.
red_andr
Начнём с того, что decimal32/64/128 новые форматы и появились только в стандарте 2008 года. И эти экзотические процессоры, включая IBM System z9, на самом деле очень распространённы там, где нужны эти типы, а именно в банках и страховых конторах. В остальном мире достаточно программной эммуляции.
Chaos_Optima
picul
Как еще один формат хранения данных в памяти — да, можно, хотя вряд ли оно того стоит.
0serg
float, double и т.п. это тоже форматы хранения данных в памяти
picul
float, double — это еще и формат хранения данных в регистрах. Я имел в виду, что posit подойдет именно для записи в память после серии обычных floating-вычислений. Но делать регистры в posit формате, добавляя к лейтенси практически любой операции сложную распаковку-упаковку — ИМХО так себе затея.
0serg
Распаковка там не особенно сложная а в регистрах, вообще говоря, вообще не обязательно float / double хранить: в x87 например регистр всегда 80-битный
VolCh
Вроде как предполагается, что это именно хранение в памяти, регистры в другом формате, и они 512 бит. Заточен на цепочки/циклы регистровых вычислений.
0serg
Для posit32 «близкие к единице» — это от 1e-6 до 1e+6. (2^-20...2^+20).
Реализация собственно FP-блока даже проще, там просто добавляется, грубо говоря «распаковка» из posit в float, затем собственно вычисления и «упаковка» результата обратно. На фоне сложности самого FP-блока эта распаковка-запаковка минимальна.
Вот отличная статья на тему: hal.inria.fr/hal-01959581v3/document
DrSmile
Собственно, первое, что стоит сделать для приближения формата к железу — это выкинуть дополнение до двух для отрицательных чисел. Сразу серьезно уменьшит latency.
0serg
У FP-блока много транзисторов потому что много стадий и latency тоже большая потому что много стадий. Pack и unpack добавят не так много на этом фоне, да и конкретно для FP величина latency вообще в большинстве случаев не критична. Ширина — да, будет больше, но это компенсируется за счет отказа работы с NaN / denormals и возможностью использовать более эффективное представление чем 754 для внутренней реализации логики.
picul
Эмм, так а куда тогда пристроить повышеную точность Posit'а, если операции будут выполняться с обычными float'ми? Размерность операций повышать что ли?
0serg
Грубо говоря любой posit32 может быть точно представлен в виде double. Т.е. внутри блока вычисления производятся с точностью double, но в памяти они занимают только 32 бита.
old_bear
Это будет верно только если любой double может быть точно представлен в виде posit32, а не наоборот.
0serg
Большинство вычислений с double производят результат который не имеет точного представления в double. Не вижу никакого смысла в свете этого рассуждать о «точном представлении». Так или иначе результат вычислений округляется.
ads83
В соседней статье приводится рекуррентное соотношение Мюллера
В этой же статье скрипт на Питоне для float дает ошибку на 13 итерации.
С одной стороны, эта формула искусственная в том смысле, что создана для демонстрации проблем с округлением. Возможно, она нацелена именно на текущий формат хранения чисел.
С другой, ряд должен сходиться к числу 5 (небольшое число), используются только базовые арифметические операции, и ошибка проявляет себя довольно быстро.
Было бы интересно увидеть, как ведет себя Posit на этом примере. Уважаемый автор, вас не затруднит провести тесты и рассказать о результатах?
P.S. Согласно статье, числа с фиксированной точкой на 23 шаге выдают неверный результат. Не уверен, имеет ли в данном случае значение язык (не знаю "внутренности" ни Python, ни С), но примерное представление получить будет возможно.
Refridgerator Автор
У posit результаты как у float без какого-либо намёка на преимущество.
4.64475
4.77073
4.85972
4.99343
6.59878
30.0151
88.4203
99.3479
99.9673
99.9984
99.9999
100
100
100
ads83
Спасибо!
BD9
Нет, надо правильно формализовать задачу. Решение может быть численным.
ads83
Я не знаю С/С++/С#, проверяльщик из меня так себе. Но для интереса глянул код.
1.1. Кажется, корень не используется в ваших проверках, но стоит ли доверять библиотеке? Особенно в свете того, что для запуска вам пришло что-то там исправлять. Корректно ли там реализованы другие операции?
2.1. В качестве примера возьмем результат для N=4, lenght=10^10. Если я все верно понял, при N=4 имеем угол 90°, то есть sin=1, cos=0. Умножение-сложение чисел 0, 1 и 10^10 друг на друга 4 раза у вас выдает 9216 — что очевидно неправильно.
2.2. Я осознаю, что в библиотеке нет тригонометрических функций, но можно выбрать такие углы поворота, что sin и cos будут рациональными числами. Тогда их можно будет задать явно. Думаю, sin=0, cos=1 не покажет ничего, а вот sin=1/2, cos=1/2 (угол 45°) было бы интересно посмотреть.
2.3. Цель — не облегчить испытания для Posit, а избежать ситуации Posit(doubleWithError). Поэтому задаем sin=cos=0.5 для всех трех "конкурсантов"
2.4. Как по мне, удлинение вектора переносит ошибку с младших разрядов в старшие (err=1^-38*lenght=1^10), но мало говорит о том, как ошибка накапливается. Было бы интересно посмотреть на сотню-тысячу вращений на 45°.
BD9
Что?!
Теорема Пифагора не освоена?
ads83
Точно, писал поздним вечером, ошибся.
ok8
В 2.2, наверное, не хватает квадратного корня в знаменателе у двойки.
Refridgerator Автор
Это тоже вопрос к авторам. Что-то пришлось исправлять только потому, что компилятор другой.
Да, но у меня не было другого варианта — там нет ни своих синусов, ни конструктора из символьного представления числа.
0serg
Любой posit32 имеет точное предоставление в double, так что здесь это не существенно
ads83
На мой взгляд, проблема в том, что этот double получен в результате вычислений и поэтому содержит некоторую ошибку. Пример с поворотом на 90 градусов в этом смысле очень показателен. То есть мы сначала задали некоторое произвольное отклонение, а потом смотрим на ее поведение. Для меня это отличается от задачи как быстро ошибки вычислений/округлений в pure posit приведут к погрешности, большей чем ошибки вычислений в pure float.
Возможно, вычисление сходящегося бесконечного ряда будет более чистым экспериментом. Во-первых, процесс итеративен. Во-вторых, мы можем отличить погрешность расчета от погрешности округлений. Например, 1/2+1/4+..1/n^2 сходится к единице. Ошибка расчета не должна превышать последний член. Зададим n=100,1000 и сравним отклонения. Имеет ли этот тест недостатки?
Другой отличный тест — специальная рекуррентная формула, в которой posit отказывает так же быстро, как и float.
Refridgerator Автор
ads83
Сегодня утром осознал, что кроме неинтересных случаев с (0;1) ничего не получится. Мотивация использовать рациональные числа в том, чтобы вычисления начинались с posit, а не posit(double_with_error). Мне не интересно вытащить posit в лидеры, мне интересно продумать такие тесты, которые будут сравнивать два формата
в идеальных условияхв чистом виде и на более-менее реальной задаче.Более реальная задача, на мой взгляд — это "крутить" не один круг, а десятки-сотни раз, но с небольшой длиной вектора. Это лучше покажет скорость нарастания ошибки со временем. Хотя не избавляет от проблемы, что вращать мы будем по неидельному кругу, т.к. sin даже в double посчитаем приблизительно.
Более идеальные условия — вращать вектор не по кругу. А умножать на десяток (заранее заданных) комплексных чисел с рациональными коэффициентами. Правда, тут вопросы к реальности задачи.
Вычисление сходящихся рядов, мне кажется, довольно неплохой сценарий. "Хитрая" формула Мюллера (еще раз спасибо! за тест), ряды типа знакопеременного гармонического или ряда обратных квадратов позволяют сравнить отклонение от теоретического значения, но используют числа, меньшие 1.
Еще была идея с итеративным извлечением корня. Этот процесс сходится к единице, можно сравнивать число итераций когда next_sqrt < epsilon. Но в используемой вами библиотеке функция sqrt какая-то сомнительная, и может оказаться что мы измерим расхождения в алгоритмах, а не форматах числа.
DrSmile
Refridgerator Автор
Да, действительно, был неправ. Но вектор вернуться в единицу уже не сможет.
ads83
А если подобрать количество "оборотов" так, чтобы вернуться? Пусть даже не в единицу, а другое, вычислимое аналитически, состояние вроде (1;1) или (0.5;2) ?
Refridgerator Автор
Ну так это просто, достаточно просто к комплексному виду перейти. После n-вращений получим вектор с рациональными коэффициентами, и тенденция к упрощению у них не наблюдается.
lain8dono
Теоретически алгоритмы нужно оптимизировать под семантику этого posit. Грубо — вы взяли алгоритмы, которые уже пофиксили для IEEE 754, а для posit нужны другие хаки (которые ещё и найти нужно). И при сравнениях надо лезть в битовое представление этого всего.
Даже если результат известен, то всё равно было бы интересно посмотреть на это всё.