(ст)IEEE-754 fun

[grozny]

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

бумажка с описанием IEEE-754 стоит $42. Жлобы позорные.
http://grouper.ieee.org/groups/754/faq.html

http://www.scri.fsu.edu/~jac/MAD3401/Backgrnd/ieee-expl.html

http://www.cs.nyu.edu/cs/faculty/overton/book/

полезные константы:
http://www.gnu.org/manual/glibc-2.2.3/html_node/libc_674.html


http://research.microsoft.com/~hollasch/cgindex/coding/ieeefloat.html

0. float - число в степенной записи по основанию 2.

1. 1 бит знак, 8 бит поле порядка, 23 бита - поле мантиссы

2. В поле порядка пишется целая часть логарифма числа по основанию 2 плюс 127 (называют bias, power offset etc.)

3. Единица в мантиссе не пишется. Т.е. если у нас есть плавающее IEEE-754 {sign;power;mantis} то само число в десятичной нотации есть: x = (знак)*(1.(мантисса/( 1<<23 ) )*(2^порядок-127))

4. Есть некоторое количество "странных чисел" - ненормальные, не-числа (NaNQ, NaNS), а также два нуля: +0 и -0.

Следствие 1: (ноль он и в Африке ноль...)
+0.0f есть 0х00000000
-0.0f есть 0х80000000
1.0f в записи IEEE есть 3F800000

Следствие 2: простой способ подсчёта логарифма: берём порядок, вычитаем биас 127. Затем берём поле мантиссы, преобразуем в плавающее, из результата читаем порядок, поправляем на биас и добавляем в пред. результат. Повторить операцию с новой мантиссой пока не будет достигнут желаемый уровень точности.

Следствие 3. А если хватит кусочно-линейной аппроксимации, дак просто складываем порядок с мантиссой как целые. Так что грубый логарифм можно вычислить 1 сдвигом (порядок задвинуть на место), и 2 сложениями (убрать биас и сложить порядок с мантиссой).

Следствие 4. Т.к. знак у плавающего в том же разряде, что и у целого, можно очень быстро сказать отличаются или нет знаки произвольного плавающего и целого. Ну и соответственно, положительное или отрицательное данное плавающее (однако бойтесь +0/-0!)

и т.д.

[RXL]

Прям статья получилась Интересно

К пункту 4: ты забыл про +-бесконечность.

Я бы еще добавил это:
https://forum.shelek.ru/index.php/topic,894.0.html
Если один из операндов вычисляемый, то при прямом сравнении возможны ошибки.

[grozny]

Мож и впрямь доработать, да повесить в статьи, чтоб под рукой было...

+- бесконечность тоже добавляет радостей, ага. Но обычно главный гемор вокруг +-0 - всё-таки в обычной арифметике ноль он один.

Про сравнение: для максимально точного сравнения лучше брать FLT_EPSILON = 1.19209290E-07F или для двойной точночти DBL_EPSILON = 2.2204460492503131E-016
 (<float.h>).

Соответственно, точная проверка на равенство есть

Код:

#include <float.h>
float f1,f2;

f1=calc_f1();
f2=calc_f2();

if( fabs(f1-f2)<=FLT_EPSILON) //then f1 is mathematically equal to f2

Надо посмотреть, для f1={binary_rep(-0.0f)}; f2={binary_rep(+0.0f)}; сработает эта проверка или нет. Вроде б должна, ага?

[Гром]

grozny, давай - оформляй - ты давно собирался в авторы прописку оформить....

[RXL]

Надо посмотреть, для f1={binary_rep(-0.0f)}; f2={binary_rep(+0.0f)}; сработает эта проверка или нет. Вроде б должна, ага?

gcc делает сравнени через FPU и 0==-0, а при распечатке printf("%f",a) дает "-0.000000".

Код:

#include <stdio.h>
int main() {
    float a,b;
    unsigned long *p;

    p=(unsigned long*)(void*)&a; *p=0x00000000;
    p=(unsigned long*)(void*)&b; *p=0x80000000;
    if(a==b) printf("a=%f b=%f\n",a,b);
    return 0;
    }

[grozny]

Да, добавим к списку фактов:
+0.0f есть 0х00000000
-0.0f есть 0х80000000

[Vorlon]

Что-то я не совсем понял:

1.0f есть 3F800000 или 0'011 1111 1'000 0000 0000 0000 0000 0000
------------------------- знак |. порядок .| мантисса
так, что ли? :?  Надо бы детальней разобраться.

А вообще, я хотел спросить о том, какую точность в десятичном представлении обеспечивают float и double, т. е. какое N мне указать в printf("%.NG",FloatOrDoubleVar), что бы ни чего не потерять и не отхватить лишнего.

[grozny]

Нуууу - дык про биас в порядке надо помнить. 127 вычитаем из порядка, получаем 0, 0 в 0-й есть 1, ага?  

ско-ко знаков после точки брать? 7 у плавающего, 16 у двойного (см. FLT_EPS, DBL_EPS). По умолчанию столько и печатают, так что если %G напишешь, то возьмутся сколько есть и с максимальной умеющей смысл точностью.

FLT_EPSILON = 1.19209290E-07F
DBL_EPSILON = 2.2204460492503131E-016
(<float.h>).

[Vorlon]

Не, без указаний выводит и float, и double одинаково - 7 цифр.
У меня результат расчета с переменными типа float выводится таким образом:
1.64473 - по умолчанию
1.644725 - "%.8G"
т. е. по умолчанию происходит округление, а мне его не надо.

Про биас сразу не заметил , но все равно что-то у меня не то:
25.5 = 0x41CC000
знак - 0
порядок - 4
мантисса - 4980736 (здесь я абсолютно не уверен)
и по формуле: (+1)*(1.4980736)*(2^4)=23.9691776 ни как не 25.5
И где я торможу?

[RXL]

Про биас сразу не заметил , но все равно что-то у меня не то:
25.5 = 0x41CC000
знак - 0
порядок - 4
мантисса - 4980736 (здесь я абсолютно не уверен)
и по формуле: (+1)*(1.4980736)*(2^4)=23.9691776 ни как не 25.5
И где я торможу?

Тормозишь   - двоичная арифметика
4980736/0x800000  (23 бита) = 0.59375

[Vorlon]

Ааааа... понял :!:
Благодарю за разъяснения.

В итоге, формула для перевода типа float:
(знак)*(1. (мантисса/( 1<<23 ) )*(2^(порядок-127))


Но может кто-нить подскажет, как вывести абсолютно все значащие цифры, не вдаваясь в самостоятельный перевод.

[RXL]

Сколько знаков выводить? - grozny уже писал.
Приблизительный подсчет: размер мантисы делить на 10, умножить на 3 и округлить в большую сторону. Т.е. на 10 двоичных порядков приходится примерно 3 десятичных. Чем больше разрядов, тем менее точная эта формула.

[grozny]

В итоге, формула для перевода типа float:
(знак)*(1. (мантисса/( 1<<23 ) )*(2^(порядок-127))

Да, вроде правильная формула и более понятная. Поправил в оригинал.

[Sighted]

Ребята, вы тут все умные, помогите пожалуйста. А то я уже туплю.
Задача такая. У меня есть текстовый массив, представляющий число дабл во внутреннем формате (к примеру "0x3ff0000000000000" или "0x4000000000000000"<-- это тескст). Моя задача получить из него число типа дабл.
Подскажите алгоритм, пожалуйста. Заранее благодарен.

[RXL]

Sighted, ни какого чпециального алгоритма не нужно - просто переведи текст в 64-х битное целое и используй его как double.

Пример:

Код:

union aaa {
    double dbl;
    unsigned long long ull;
    };

union aaa a;

scanf("%ll",&a.ull);
printf("Double: %f\n",a.dbl);

Только тип 'long long' - это GNU расширение. Поищи в коде к своему компилятору подходящий тип. Пример условный и как шаблон использовать не рекомендую - это только идея.

[Sighted]

Спасибо огромное за подсказку, но я, если честно говорить, не очень хорошо понял, как это осуществить. У меня, впрочем, своя специфика. Я пишу на С под... пусть будет линух. Это не совсем так, потому как библиотеки под С не стандартные и большинства функций нету, приходится писать самому. Использую кросс-компилятор. long long он не поддерживает. Причем, похоже, он в принципе не поддерживает 64-х разрядное целое. Впрочем, это в данном случае не важно.
Меня именно сильно интересует как осуществить это: "просто переведи текст в 64-х битное целое".
В общем-то я догадываюсь, что туплю по страшному, но всё же задам тупой вопрос:
"0x3ff0"=(?)0011 1111 1111 0000
Это верно?
Заранее благодарен.

[RXL]

Sighted, вроде верно. Не забудь про порядок байт.

[Sighted]

Я извиняюсь, что наглею, но... Я всё так же туплю...
1 - целое - это 0...0 0001
Тогда как здесь  0011 1111 1111 0000.
2 - целое - это 0...0 0010
Тогда как здесь (0х4000) - 0100 0000 0000 0000.
Но вот 3 - (0...0 0011) здесь уже (0х4008) - 0100 0000 0000 1000....
4 - (0..0 0100) - (0х4010) - 0100 0000 0001 0000
5 - (0..0 0101) - (0х4014) - 0100 0000 0001 0100
и т.д
Собственно... Как мне переводить такую последовательность битов (0100 0000 0001 0100) в 0101?
"Не забудь про порядок байт." <-- Не мог ты напомнить, что это за порядок... А то я новичок в "дабле". IEEE 754 читал, но, видимо, не очень понял...

[RXL]

Sighted, ты сравниваешь целое и плавающее значение в их двоичных представлениях - это некоректно. В начале этой темы есть описание плавающего формата.
Тебя мой пример смутил? union, если не знаешь, заставляет объявленные внутри нее переменные использовать общкю памят, но более ни каких отношений между ними нет. Тем более преобразований.

Код:

union aaa {
    double dbl;
    unsigned long long ull;
    };

union aaa a;

a.ull=0x3ff0000000000000;
printf("Double: %f\n",a.dbl);

Double: 1.000000

[Sighted]

Да нет, что такое union я вроде знаю...
Дело вот в чём...
В этой строке:
a.ull=0x3ff0000000000000;
ты целочисленной переменной присваиваешь число же (правильно?). Может я конечно не догоняю чего-то конкретно, но:
Конкретно проблема такова:
У меня есть массив, скажем (я работаю именно с текстовым массивом поступающих чисел в виде текста, соответственно):
char txt[100]="0x3ff0000000000000";
То есть наборчик аскии-чисел.
Как мне из этого наборчика сделать число типа, ну пусть unsigned long long?

[RXL]

В линуксе - sscanf(txt,"%ll",&a.ull);

[Sighted]

RXL, спасибо большое за подсказки. Это и в самом деле довольно очевидный вариант.

[Sighted]

x=(знак)*(1. (мантисса/( 1<<23 ) )*(2^(порядок-127))

Хотел спросить, а для двойной точности подходит та же формула только с поправками на порядок и мантиссу, то есть:
х=(знак)*(1. (мантисса/( 1<<52 ) )*(2^(порядок-1023))
?
Или здесь требуются какие-либо ещё видоизменения?

[RXL]

Sighted, логика та же. Загляни в доки intel-а - там описны форматы длиной 32, 64 и 80 бит.

[Sighted]

Ок, спасибо большое. Вроде всё получилось.

[Uncle_Mo]

Привет всем! Класную тему вы ребята подняли! Все это здорово и все понятно!
Я вот вопросик подкину! А как обратно то сделать?
Вот например:
a = 1;
Как превратить в мантису и порядок?

[RXL]

Uncle_Mo, в теории: сдвигом преобразовать число в вид 1.ххх. Из числа потребовавшихся сдвигов вычислить порядок, а у числа 1.ххх убрать 1 - вот и мантиса. И про знак не забудь.