LINUX.ORG.RU

Машинное эпсилон, sse и fpu87 и спецификатор для printf для __float128

 , , ,


0

3

Не могу подобрать. %lle и %Le не подходят. Другие спецификации перепробовал тоже. Кое-где %Q упоминается, но у меня не работает

Вычисляю машинное эпсилон и не получается для 128-битной плавучки его правильно отобразить.

Вот пример кода, вычисления машинного эпсилон для разных типов (вообще там занятные особенности всплыли работы gcc, особенно если в ассемблерный листинг глянуть)

#include <stdio.h>

int main()
{

float        e_float,eps_float=1.0;
double       e_double,eps_double=1.0;
long double  e_longdouble,eps_longdouble=1.0;
__float80    e__float80,eps__float80=1.0;
__float128 e__float128,eps__float128=1.0;

 e_float=1+eps_float;
 e_double=1+eps_double;
 e_longdouble=1+eps_longdouble;
 e__float80=1+eps__float80;
 e__float128=1+eps__float128;

 while (e_float!=1) 
 {
   eps_float=(float)(eps_float/2);
   e_float=(float)(eps_float+1);
 }
 
 while (e_double!=1) 
 {
   eps_double=(double)(eps_double/2);
   e_double=(double)(eps_double+1);
 }

 while (e_longdouble!=1) 
 {
   eps_longdouble=(long double)(eps_longdouble/2);
   e_longdouble=(long double)(eps_longdouble+1);
 }

 while (e__float80!=1) 
 {
   eps__float80=(__float80)(eps__float80/2);
   e__float80=(__float80)(eps__float80+1);
 }

 while (e__float128!=1) 
 {
   eps__float128=(__float128)(eps__float128/2);
   e__float128=(__float128)(eps__float128+1);
 }

 printf("size: float = %d  double = %d long double = %d  __float80= %d __float128= %d\n",sizeof(float),sizeof(double),sizeof(long double),sizeof(__float80),sizeof(__float128));
 printf("float eps = %e \ndouble eps = %e \nlong double eps = %lle\n__float80 eps = %lle\n__float128 eps = %Le \n",eps_float,eps_double,eps_longdouble,eps__float80,eps__float128);

return 0;
}

Компилирую

$ gcc -mfpmath=387 gccfloats.c -o gccfloats-387

$ gcc -mfpmath=sse  gccfloats.c -o gccfloats-sse

Запускаю. Оба варианта с 387 и sse дают один и тот же результат

size: float = 4  double = 8 long double = 16  __float80= 16 __float128= 16
float eps = 5.960464e-08 
double eps = 1.110223e-16 
long double eps = 5.421011e-20
__float80 eps = 5.421011e-20
__float128 eps = 0.000000e+00

Понятно, ноль там для __float128 быть не может, это результат какого-то округления при выводе.

Если с sse компилировать, то тоже самое. Пробовал вместо %Le написать %.числоf, так хоть %.50000f то будет 50 тысяч нулей. А это уже какая-то фигня.

Кто виноват и что делать? Как узнать машинное эпсилон для 128 битного типа с плавающей запятой.

Попутно любопытные открытия для себя сделал при просмотре ассемблерного листинга.

1. long double и __float80 - одно и тоже. По крайней мере для gcc 10.2.1 Впрочем, вроде где-то об этом было. Что интереснее, для расчета всегда используется x87-й блок, даже если -mfpmath=sse -msse2

2. Хотя long double - 80-ти битный тип, длинна 16 байт. Видимо из-за выравнивания

3. Для __float128 всегда используется sse2, даже если указано -mfpmath=387

Не понял, как можно заставить использовать avx и avx2

★★★★★

Последнее исправление: praseodim (всего исправлений: 1)

Вычисляю машинное эпсилон и не получается для 128-битной плавучки его правильно отобразить.

glibc printf может и не поддерживать. В мане нету например. А даже если будет поддерживать - то несовместимо с другими libc. Напиши свой (печатаешь цифру слева от точки, вычитаешь её, умножаешь на 10 и так по циклу).

Понятно, ноль там для __float128 быть не может, это результат какого-то округления при выводе.

Скорее всего не округления даже, а не того типа. Он там ждёт double и интерпретирует полученный байты как double, а ему дают __float128 и хорошо если в том же регистре/месте памяти а не в другом.

Что интереснее, для расчета всегда используется x87-й блок, даже если -mfpmath=sse -msse2

Потому что sse не поддерживает этот тип (у х87 есть три типа: одиночной, двойной и расширенной точности).

3. Для __float128 всегда используется sse2, даже если указано -mfpmath=387

Аналогично, х87 не поддерживает его. Насчёт того как он вообще работает не знаю, мне казалось что предел точности sse это double.

Не понял, как можно заставить использовать avx и avx2

-mfpmath=avx какое-нить?

firkax ★★★★★
()
Последнее исправление: firkax (всего исправлений: 1)
Ответ на: комментарий от MKuznetsov

а покажи-ка свой uname

Да под Linux я 5.16 ядро, Debian 11, впрочем не думаю, чтобы ядро имело значение

long double вполне печатается.

praseodim ★★★★★
() автор топика
Последнее исправление: praseodim (всего исправлений: 1)
Ответ на: комментарий от firkax

glibc printf может и не поддерживать. В мане нету например. А даже если будет поддерживать - то несовместимо с другими libc. Напиши свой (печатаешь цифру слева от точки, вычитаешь её, умножаешь на 10 и так по циклу).

Неужели нет какой-то относительно стандартной библиотеки? Не в glibc, так просто что-то из поставки gcc?

Потому что sse не поддерживает этот тип (у х87 есть три типа: одиночной, двойной и расширенной точности).

Если у меня не ложные воспоминания, то есть еще 64-х битный целый. Был смысл на 16-ти битных x86 и 32-битных i386. Но думал, что если задано sse / sse2, то полагалось бы хоть через отдельную библиотеку и очень медленно, но считать таки 80-ти битный тип без использования x87

Аналогично, х87 не поддерживает его. Насчёт того как он вообще работает не знаю, мне казалось что предел точности sse это double.

sse2 тоже не поддерживает. Вызываются отдельные библиотечные функции для работы с 128-битным типом. И тут sse2 и 387 в равных условиях - для обоих это не аппаратный тип и для обоих это удвоенный double. Я еще понимаю 80-битный тип сильно не стандартный.

-mfpmath=avx какое-нить?

Допускает только 387 или sse или комбинацию 387 и sse

praseodim ★★★★★
() автор топика
Последнее исправление: praseodim (всего исправлений: 3)
Ответ на: комментарий от ox55ff

Бесполезная константа

?????

Используется при математических расчетах. Если разница двух чисел по модулю меньше, их следует считать равными в расчете. Также и при задании относительной погрешности.

Выведи значение регистра в hex и на бумажке по стандарту вычисли значение.

Все же хочется как-то узнать есть ли что-то более-менее стандартное. А то как-то вообще странно. Числовой тип данных есть, но вывод его не предусмотрен.

praseodim ★★★★★
() автор топика
Последнее исправление: praseodim (всего исправлений: 1)
Ответ на: комментарий от firkax

glibc printf может и не поддерживать. В мане нету например. А даже если будет поддерживать - то несовместимо с другими libc. Напиши свой (печатаешь цифру слева от точки, вычитаешь её, умножаешь на 10 и так по циклу).

Порылся, погуглил. Есть вокруг gcc библиотека quadmath, там есть функция quadmath_snprintf https://gcc.gnu.org/onlinedocs/libquadmath/quadmath_005fsnprintf.html

В общем, добавляется что-то вроде

#include <quadmath.h>

...

char buf[128]; 
int n = quadmath_snprintf (buf, sizeof buf, "%+-#*.20Qe", 46, eps__float128);
if ((size_t) n < sizeof buf)  
    printf ("%s\n", buf);

Компилировать с - lquadmath - напечатало +9.62964972193617926528e-35

Несколько напрягает, что вроде просто для вычислений она не обязательна, почти лишняя сущность получается.

Может у кого лучшее решение есть? Кроме колхозинга собственной функции.

praseodim ★★★★★
() автор топика
Последнее исправление: praseodim (всего исправлений: 2)
Ответ на: комментарий от praseodim

Неужели нет какой-то относительно стандартной библиотеки? Не в glibc, так просто что-то из поставки gcc?

__float128 сам по себе не стандартный, а так не знаю, может и есть, но зачем возиться с поисками когда тривиальная печаталка дробей которые тебе нужны пишется за 5 мин?

Если у меня не ложные воспоминания, то есть еще 64-х битный целый.

А ещё есть BCD-формат. Но всё это нативно мейнстримными Си не поддерживается, я писал про дробные форматы.

Но думал, что если задано sse / sse2, то полагалось бы хоть через отдельную библиотеку и очень медленно, но считать таки 80-ти битный тип без использования x87

Зачем использовать медленный эмулятор когда можно быстро без эмулятора? Процов с sse но без х87 не бывает вроде как.

sse2 тоже не поддерживает. Вызываются отдельные библиотечные функции для работы с 128-битным типом. И тут sse2 и 387 в равных условиях - для обоих это не аппаратный тип и для обоих это удвоенный double.

Да нет думаю там всё хуже - вручную парсится формат и всё делается эмулятором через целочисленные операции. Не представляю, как можно сэмулировать float128 через какие бы то ни было операции с double. А целочисленные быстрее делать sse чем обычной целой арифметикой, видимо.

Я еще понимаю 80-битный тип сильно не стандартный.

Он как раз максимально стандартный, уже около 45 лет как, был с первых fpu.

firkax ★★★★★
()
Ответ на: комментарий от praseodim

Используется при математических расчетах. Если разница двух чисел по модулю меньше, их следует считать равными в расчете. Также и при задании относительной погрешности.

Эм, нет. Я думал ты для развлечения или тренировки эту прогу проверял и не стал писать о практической стороне вопроса.

Во-первых, это «эпсилон» само по себе относительное, то есть это не константа которую прибавляют а коэфициент. Ты его с единицей складывал, поэтому умножать ничего не надо было. А вот если б ты прибавлял его к 16, то обнаружил бы что и «эпсилон» стало в 16 раз больше.

Во-вторых, не «если разница меньше то следует считать равными», а разницу меньше тебе проц и не покажет (ты ж сам сравниваешь строгое равенство 1+eps==1 которое оказывается истиной). То есть это погрешность, которая уже случилась и ты можешь только оценить её величину (но надо учитывать, что она применяется к каждой математической операции, а не только к итоговому результату формулы, так что накопиться может и больше).

В-третьих, эту величину можно безо всяких опытов легко найти их формата float/double итд. Например у float 24 значащих двоичных разряда, у double - 53, у формата расширенной точности (80бит) - 64. За исключением очень маленьких чисел на границе допустимых порядков - там точность уменьшается на 1 бит с каждым порядком, пока не дойдёт до нуля где и закончится предел возможностей типа.

firkax ★★★★★
()
Ответ на: комментарий от firkax

Во-первых, это «эпсилон» само по себе относительное, то есть это не константа которую прибавляют а коэфициент. Ты его с единицей складывал, поэтому умножать ничего не надо было. А вот если б ты прибавлял его к 16, то обнаружил бы что и «эпсилон» стало в 16 раз больше.

Ну да, коэффициент должен быть. Но я не согласился с тем, что бесполезная константа, в математических методах она вполне учитывается. Если я не сильно отстал от жизни конечно.

В-третьих, эту величину можно безо всяких опытов легко найти их формата float/double итд.

Можно конечно, но если пишется программа лучше все-таки вычислить прямо в ней. Чтобы недоразумений не было, а то где-нибудь под что-то нестандартное ее скомпилируют.

Он как раз максимально стандартный, уже около 45 лет как, был с первых fpu.

Я имел ввиду стандарт IEE754 для вычислений с плавающей точкой. Насколько я помню 80-битный у x87 формально ему не соответствует.

praseodim ★★★★★
() автор топика
Последнее исправление: praseodim (всего исправлений: 1)
Ответ на: комментарий от praseodim

просто для вычислений она не обязательна, почти лишняя

в ней реализованы как-раз таки вычислительные функции. И описан вожделенный вами epsilon __FLT128_EPSILON. Посмотрите quadmath.h

выводить вручную epsilon лучше всё-таки на этапе configure/cmake (с целью сравнить с указанным) а отнюдь не в рантайм при старте программы. Вот какими макросами, тайными словами и шаманскими жестами - не знаю

Может у кого лучшее решение есть? Кроме колхозинга собственной функции.

вместо колхозинга - внедрятинг quadmath_snprintf() в обычный printf. Можно ведь добавлять свои спецификаторы или переопределять/дополнять существующие

MKuznetsov ★★★★★
()
Ответ на: комментарий от MKuznetsov

в ней реализованы как-раз таки вычислительные функции. И описан вожделенный вами epsilon __FLT128_EPSILON. Посмотрите quadmath.h

Хорошо =)

вместо колхозинга - внедрятинг quadmath_snprintf() в обычный printf. Можно ведь добавлять свои спецификаторы или переопределять/дополнять существующие

Тоже немного колхозинг, но я понял, стандартного варианта нет. Хотя все же странно. Реализация какая-то неконсистентная получается. Могли бы все же, если сделали расширение, добавив тип, добавить и работу с ним всяких printf. Между прочим, надо будет еще проверить работу вычислительных функций, а то подозреваю, что для них тоже libquadmath нужен.

praseodim ★★★★★
() автор топика
Ответ на: комментарий от praseodim

но я понял, стандартного варианта нет.

и видимо не будет :-( стандарты это хорошо, но есть реальность. А она например в том что crt windows обошла стороной long double и соотв. vc тоже (там видимо отдельно свои нашлёпки, не скажу,студию снёс года 3-4 назад и не жалею). А MS это порядка 50% вообще всего C/C++. Поэтому все тоже подзабили болт

на прущих в гору ARM неизвестно как там с длинными дабл. Как там clang, что в эмбедщине творится, как в embarcadero реализовано..

пока все в реальности не договорятся, стандарт это бумага хуже туалетной

надо будет еще проверить работу вычислительных функций

успехов.. в итоге прочувствуете отчего продают и покупают мат.библиотеки и почему дедушка фортран это вещь

MKuznetsov ★★★★★
()
Ответ на: комментарий от praseodim

Используется при математических расчетах. Если разница двух чисел по модулю меньше, их следует считать равными в расчете.

Машинное эпсилон это минимальная ненулевая разница двух чисел. По определению. Меньше неё только ноль. Так что если твоя хотелка состоит в том, что искать пары, отличающиеся меньше, чем на эпсилон, можешь просто сравнивать напрямую (==).

i-rinat ★★★★★
()
Ответ на: комментарий от i-rinat

Нет. По определению машинный эпсилон это максимальное такое число, что 1+eps=1.

$ cat 1.c
#include <stdio.h>

int main() {
        double eps = 1./(1L<<53);
        double not_eps = 1./(1L<<52);
        double a = 1e-17;
        double b = 2e-17;
        printf("%e %e %d %d %d %d %e %e\n", eps, not_eps, (1.+eps==1.), (1.+not_eps==1.), (a==b), (a<eps), a, b);
}
$ ./a.out
1.110223e-16 2.220446e-16 1 0 0 1 1.000000e-17 2.000000e-17

Покурите на досуге представление чисел с плавающей точкой в памяти.

Reset ★★★★★
()
Последнее исправление: Reset (всего исправлений: 1)
Ответ на: комментарий от Reset

Если не трудно, добавь на страницу https://ru.wikipedia.org/wiki/Машинный_ноль ссылку на твоё определение. А то для приведённого мной варианта там есть ссылка на литературу, а для твоего любимого варианта нет.

Раз уж ты выучил именно такой вариант, наверняка не составит труда вспомнить учебники, по которым ты учился. Полезное дело сделаешь.

i-rinat ★★★★★
()
Ответ на: комментарий от i-rinat

Английская википедия фактически с этого определения начинается. Меня устроит любое определение, которое написано во втором абзаце. Там есть и со ссылками. Вот это «Машинный ноль (Машинный нуль) — числовое значение с таким отрицательным порядком, которое воспринимается машиной как ноль» просто безграмотно. Так как как ноль может восприниматься и число с положительным порядком (внезапно), а может как ноль не восприниматься число меньшее 1e-16 (сюрприз!). У Вержбицкого, кстати, такого нет, специально скачал проверил. Там все правильно написано.

Reset ★★★★★
()
Ответ на: комментарий от Reset

Английская википедия фактически с этого определения начинается.

И другим продолжается. Но ты его в упор не увидишь, это я уже понял.

Вот это «Машинный ноль (Машинный нуль) — числовое значение с таким отрицательным порядком, которое воспринимается машиной как ноль» просто безграмотно. Так как как ноль может восприниматься и число с положительным порядком (внезапно), а может как ноль не восприниматься число меньшее 1e-16 (сюрприз!).

лицоладонь

Там же есть ссылка на книгу. Причём текст доступен онлайн, просто кликнуть и прочитать. Меньше десяти строк.

Или ты считаешь «машинный ноль» и «машинное эпсилон» тождественными понятиями? Вряд ли, конечно, но всё возможно.

i-rinat ★★★★★
()
Ответ на: комментарий от Reset

Вержбицкий какого года? Машинный ноль — это из книг 70-80х годов прошлого века, когда стандартно было только нормализованное представление и превращение в ноль всего, что туда не влазило. И получалось, что область вблизи нуля значительно больше, чем расстояние между двумя соседними ненулевыми числами (underflow gap).

просто безграмотно

Не надо фокусироваться только на IEEE 754.

mky ★★★★★
()
Последнее исправление: mky (всего исправлений: 1)
Ответ на: комментарий от fsb4000

На С тоже работает: https://gcc.godbolt.org/z/4rW147sGo

size: float = 4  double = 8 long double = 16  __float80= 16 __float128= 16
float eps = 5.960464e-08 
double eps = 1.110223e-16 
long double eps = 5.421011e-20
__float80 eps = 5.421011e-20
__float128 eps = 9.6296497219361792652798897129246365926905082411e-35

quadmath_snprintf

https://gcc.gnu.org/onlinedocs/gcc-13.1.0/libquadmath/quadmath_005fsnprintf.html

fsb4000 ★★★★★
()
Ответ на: комментарий от praseodim

Вот что я нашёл на С без сторонних библиотек:

Принятый пропозал в стандарт: https://open-std.org/JTC1/SC22/WG14/www/docs/n2601.pdf (X.12.1String from floating и X.12.2String to floating)

Видимо они забили на printf/scanf.

strfromf128

и

strtof128

https://gcc.godbolt.org/z/rqj489YMh

fsb4000 ★★★★★
()
Последнее исправление: fsb4000 (всего исправлений: 3)
Ответ на: комментарий от i-rinat

Я не понимаю что ты пытаешься оспорить? Вот с этим «Машинное эпсилон это минимальная ненулевая разница двух чисел.» ты обосрался, думаю корона с тебя не свалится, если ты просто это признаешь.

Reset ★★★★★
()
Ответ на: комментарий от firkax

(печатаешь цифру слева от точки, вычитаешь её, умножаешь на 10 и так по циклу)

Кстати, это 100% неправильный алгоритм.

Я не углублялся как делать правильно. Но умножение/деление на 10 для float типов 100% будет печатать неверно.

Если есть желание то можно поизучать: https://github.com/ulfjack/ryu

fsb4000 ★★★★★
()
Ответ на: комментарий от fsb4000

Не надо обобщений. Деление - да, часто потеряет точность (впрочем, где-то в последних знаках, который всё равно обычно округляют). Умножать можно, это практически точная операция в данном случае (ну, есть чуть исключений, но они настолько незначительны, что для демо-задачи автора на них точно следует забить).

firkax ★★★★★
()
Ответ на: комментарий от Reset

И.Ю. Алибеков. Численные методы, У/П. — МГИУ, 2008-01-01. — 221 с. — ISBN 9785276014623

стр. 23.

Важной характеристикой является число ε, называемое машинный эпсилон. Эта характеристика определяется как расстояние между единицей и ближайшим следующим за ней числом системы машинных чисел с плавающей точкой.

i-rinat ★★★★★
()
Ответ на: комментарий от fsb4000

Но умножение/деление на 10 для float типов 100% будет печатать неверно.

Если есть желание то можно поизучать

там как-раз таки сказано что ОБА АЛГОРИТМА ДАЮТ ВЕРНЫЙ РЕЗУЛЬТАТ.(и не один внутри имплементаций не делит float на 10.0f).. и было-бы странно если алг.с делением на 10 работал неверно - он просто переводит из одной системы счисления в другую. Ryu ищет (вычисляет) ближайшую десятичную дробь. Разница вроде невелика, но она есть.

MKuznetsov ★★★★★
()
Последнее исправление: MKuznetsov (всего исправлений: 1)
Ответ на: комментарий от MKuznetsov

Нет. Там сказано что предыдущий алгоритм неверный.

Вот, например, неправильные результаты которые выдает Double.ToString(). (или выдавал в 2020, может в OpenJDK исправили за 3 года)

  in: -2.147483648E32
out: -2.1474836479999998E32
  C: -2.147483648e+32
Java implementation is incorrect (the original string is shortest)

 in: -3.1802768007018752E17
out: -3.180276800701875E17
  C: -3.180276800701875e+17
Not a bug; Java and C are provably correct

 in: -6.4432225704600013E18
out: -6.443222570460001E18
  C: -6.443222570460001e+18
Not a bug; Java and C are provably correct

 in: 1.3867957848030121E18
out: 1.386795784803012E18
  C: 1.386795784803012e+18
Not a bug; Java and C are provably correct

 in: 1.41027712162200013E18
out: 1.4102771216220001E18
  C: 1.410277121622e+18
Java implementation is incorrect (the original string had unnecessary digits)

 in: 8.8243617212999997E18
out: 8.8243617213E18
  C: 8.8243617213e+18
Not a bug; Java and C are provably correct

Therefore, there does appear to be incorrect behavior here, but only 2 of these cases are affected. Hope this helps! 

https://github.com/ulfjack/ryu/issues/187

fsb4000 ★★★★★
()