Несколько интересных примеров чисел с плавающей запятой

задняя часть PHP

Эта статья#130Ответы на первые 6 вопросов:

Вопрос 1

int c = 0.58 * 100;
printf("%d\n",c);           //57

int x = 0.58f * 100.0f;
printf("%d\n",x);          //58

Когда оба числа имеют тип float, результатом умножения будет 58. Процесс разгадки тайны включает в себя то, как хранятся числа с плавающей запятой, правила умножения чисел с плавающей запятой и округления.

О методе хранения чисел с плавающей запятой см.Двоичное представление чисел с плавающей запятой и несколько примеров, которые здесь повторяться не будут.

0.58 * 100

Двойное представление числа 0,58:

Двойное представление числа 100 (с плавающей точкой): 0 10000000101 10010000000000000000000000000000000000000000000000000

Давайте посмотрим на процесс расчета 0,58 * 100:

Переход к стандартному закону о науке и технологиях:

0,58 равно 1,0010100011110101110000101000111101011100001010001111 * 2^-1 100 равно 1,100100000000000000000000000000000000000000000000 20^0000

добавить два показателя степени

-1 + 6 = 5

Умножить две мантиссы

1.0010100011110101110000101000111101011100001010001111 * 1.1001000000000000000000000000000000000000000000000000 = 1.11001111111111111111111111111111111111111111111111110000000000000000000000000000000000000000000000000000

нормализовать

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

добавление знакового бита mod2

0 + 0 = 0

включены

Дело здесь в том, что поскольку в двойной точности всего 52 цифры, мы должны отбросить часть мантиссы. IEEE по умолчанию использует округление до четного. Таким образом, результат округления выше: 1.1100111111111111111111111111111111111111111111111

Переполнение суждения

Проверьте, нет ли переполнения в сумме показателей. никто

Итак, конечный результат равен 0 10000000100 110011111111111111111111111111111111111111111111111111111111111111111111111111111111111

Результат: 57.999999999999992894572642399

Поскольку float/double в языке C используется для движения к нулю, окончательный результат равен 0,58 * 100 = 57.

Давайте посмотрим на процесс расчета 0,58f * 100,0f:

0.58f * 100.0f

Представление 0,58 с одинарной точностью: 0 0 01111110 00101000111101011100001

Представление 100 с одинарной точностью (с плавающей точкой): 0 10000101 10010000000000000000000

Конкретный процесс:

Переход к стандартному закону о науке и технологиях

0,58 равно 1,00101000111101011100001 * 2^-1 100 равно 1,1001000000000000000000000 * 2^6

добавить два показателя степени

-1 + 6 = 5

Умножить две мантиссы

1.00101000111101011100001 * 1.10010000000000000000000 = 1.1100111111111111111111110010000000000000000000

нормализовать

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

добавление знакового бита mod2

0 + 0 = 0

включены

Так как в одинарной точности всего 23 цифры, мы должны отбросить часть мантиссы. Приведенное выше округление с двойной точностью округляет все последние нули. Однако при округлении до одинарной точности здесь вам нужно округлить одно место в большую сторону.Это также является причиной того, что результаты двух умножений не равны.Таким образом, результат округления выше:
1.11010000000000000000000

Переполнение суждения

Проверьте, нет ли переполнения в сумме показателей. никто

Таким образом, окончательный результат равен 0 10000100 110100000000000000000000.

Результат ровно 58,0. Итак, 0,58f * 100,0f = 58.

И последнее замечание: округление до четного числа — это не совсем то же самое, что и округление, как мы узнали. могу взглянутьДвоичное представление чисел с плавающей запятой и несколько примеровЧасть округления с плавающей запятой .

вопрос 2

printf("%d",0.1 + 0.2 == 0.3);        //0,not 1

Эта проблема связана с добавлением чисел с плавающей запятой.

Давайте посмотрим на процесс добавления чисел с плавающей запятой:

Двойное представление числа 0,1: 0 0 111 111 10 11 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 0 0 1 1 0 0 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 1 0

Двойное представление 0,2: 0 0 01111111100 1001100110011001100110011001100110011001100110011010

Переход к стандартному закону о науке и технологиях

Во-первых, как и в случае с умножением, мы сначала преобразуем 0,1 и 0,2 в стандартную научную запись:
0,1 равно 1,1001100110011001100110011001100110011001100110011010 * 2^-4
0,2 равно 1,1001100110011001100110011001100110011001100110011010 * 2^-3

выравнивание экспоненты

Представьте себе добавление десятичной научной записи: 1,234567 × 10 ^ 5 + 1,017654 × 10 ^ 2, мы должны изменить последнее на 0,001017654 × 10 ^ 5, чтобы показатели степени были одинаковыми для добавления мантиссы. То же самое верно и для двоичных чисел с плавающей запятой.

Показатель степени 0,1 равен -4, а показатель степени 0,2 равен -3. Здесь одно правило:Маленькие шаги совпадают с большими.

Таким образом, показатель степени 0,1 становится равным -3, а соответствующая мантисса становится равной 0,11001100110011001100110011001100110011001100110011010.

Добавление мантиссы

0.11001100110011001100110011001100110011001100110011010 +
1.1001100110011001100110011001100110011001100110011010 =
10.01100110011001100110011001100110011001100110011001110

нормализовать

1.001100110011001100110011001100110011001100110011001110 * 2^-2

включены

Мантисса двойной точности имеет только бит 52. Округляя до четного числа, мы получаем 1.0011001100110011001100110011001100110011001100110100 * 2^-2

Итак, окончательный результат равен 0 1111111101 0011001100110011001100110011001100110011001100110100

0,300000000000000044408920985006 в десятичной системе

Таким образом, 0,1 + 0,2 = 0,300000000000000044408920985006 != 0,3

Вопрос 3

printf("%f\n",3.14f + 1e10f - 1e10f);            //0.0,not 3.14
printf("%f\n",3.14f + (1e10f - 1e10f));         //3.14

Этот вопрос не совпадает с предыдущим вопросом, и я надеюсь, что читатели смогут взглянуть на него.

Давайте проанализируем это шаг за шагом.
Представление одинарной точности 3.14f:
0 10000000 10010001111010111000011

Представление 1e10f с одинарной точностью: 0 10100000 00101010000001011111001

Преобразовать в экспоненциальное представление

Представление 3.14f: 1.10010001111010111000011 * 2^1 Представление 1e10f: 1.00101010000001011111001 * 2^33

выравнивание экспоненты

Малый порядок выравнивается с большим порядком 3.14f преобразуется в 0,0000000000000000000000000000000110010001111010111000011* 2^33 (31 последовательный 0 с справа от десятичной точки)

Добавление мантиссы

0.0000000000000000000000000000000110010001111010111000011 + 1.00101010000001011111001 = 1.0010101000000101111100100000000110010001111010111000011

нормализовать

уже является нормализованным числом, преобразование не требуется, т.е.

включены

1.00101010000001011111001 * 2^33

Результат 0 10100000 00101010000001011111001

Этот результат по-прежнему 1e10. Мы обнаружили, что 3,14 с ним и без него не имеет значения. Потому что при сложении порядка и мантиссы фактические значащие цифры 3,14 добавляются в самый конец, и все эти числа округляются.

Эта проблема больше не существует, когда число с плавающей запятой имеет двойную точность, то есть тип double, потому что мантисса типа double имеет 52 бита, что достаточно для того, чтобы гарантировать, что при округлении не будет отброшена мантисса числа 3,14.

Вопрос 4

float d1, d2, d3, d4;

d1 = 194268.02f;
d2 = 194268f;
d4 = 0.02f;
    
d3 = d1 - d2;
if (d3 > d4)
   printf(">0.02\n");
else if (d3 < d4)
   printf("<0.02\n");     //true
else
   printf("=0.02\n");     //false

printf("%f - %f = %f \n", d1,d2,d3);   //194268.015625 - 194268.000000 = 0.015625,not 194268.02 - 194268 = 0.02

Шаги аналогичны сложению, за исключением того, что операция мантиссы — это вычитание.

Представление одинарной точности 194268.02f: 0 10010000 01111011011011100000001

Представление 194268f с одинарной точностью: 0 10010000 01111011011011100000000

Преобразовать в экспоненциальное представление

194268.02f = 1.01111011011011100000001 * 2^17 194268f = 1.01111011011011100000000 * 2^17

Противоположный порядок

Код тот же. обоим по 17

Вычитание мантиссы

Легко видеть, что ответ равен 0,000000000000000000000001

нормализовать

Это довольно специфично.После нормализации мантисса равна 0, то есть 1.000000000000000000000000 * 2^-6

включены

никто

Окончательный результат:
0 01111001 00000000000000000000000
то есть 0,015625

Так почему же printf 194268.02f приводит к 194268.015625?
Ответ очень прост, двоичное число 194268 равно 101111011011011100, убрав первую 1, всего 17 бит, а мантисса числа с плавающей запятой одинарной точности всего 23 бита, всего 6 бит для двоичного представления 0,02, который равен 000001, поэтому:
0 10010000 01111011011011100000001 = 194268.015625

Вопрос 5

gcc compare.c -o compare.o
float p3x = 80838.0f;
float p2y = -2499.0f;
double v321 = p3x * p2y;
printf("%f",v321);          //-202014160,not -202014162

gcc компилируется с параметром -mfpmath:

gcc -mfpmath=387 compare.c -o compare.o 
float p3x = 80838.0f;
float p2y = -2499.0f;
double v321 = p3x * p2y;
printf("%f",v321);          //-202014162,not -202014160

Почему результаты двух прогонов различаются? И 80838,0 * -2499 без учета точности ответ равен -202014162, а не -202014160.

Давайте сначала посмотрим на значение параметра gcc -mfpmath:

Generate floating point arithmetics for selected unit unit.  The choices for unit are:

387 Use the standard 387 floating point coprocessor present majority of chips and emulated otherwise.  Code
    compiled with this option will run almost everywhere.  The temporary results are computed in 80bit
    precision instead of precision specified by the type resulting in slightly different results compared to
    most of other chips.  See -ffloat-store for more detailed description.

    This is the default choice for i386 compiler.

sse Use scalar floating point instructions present in the SSE instruction set.  This instruction set is
    supported by Pentium3 and newer chips, in the AMD line by Athlon-4, Athlon-xp and Athlon-mp chips.  The
    earlier version of SSE instruction set supports only single precision arithmetics, thus the double and 
    extended precision arithmetics is still done using 387.  Later version, present only in Pentium4 and the
    future AMD x86-64 chips supports double precision arithmetics too.

    For the i386 compiler, you need to use -march=cpu-type, -msse or -msse2 switches to enable SSE 
    extensions and make this option effective.  For the x86-64 compiler, these extensions are enabled by 
    default.   

    The resulting code should be considerably faster in the majority of cases and avoid the numerical
    instability problems of 387 code, but may break some existing code that expects temporaries to be 80bit.

    This is the default choice for the x86-64 compiler.
скопировать код

Опция 387 сохранит временный результат операции как 80 бит, а затем усекет его до 32 бит или 64 бита в зависимости от целевого типа (с плавающей запятой/двойной). Выполнение этого с помощью FPU уменьшит проблемы, связанные с округлением. (В этом примере тип FPU вычисляет «правильный» результат).

Но опция sse будет напрямую усекать результат до 32-битной одинарной точности, а затем назначать его целевому типу 64-бит.

О процессе работы с числами с плавающей запятой уже говорилось ранее, и здесь мы не будем вычислять пошагово, только нормализованная десятичная дробь для умножения мантиссы равна 1.10000001010011111011101001. Однако, поскольку целевой тип — double, FPU не усекает, а напрямую присваивает ему значение double, и в результате получается -202014162. Однако в режиме SSE будет выполнено усечение, и 001 в конце будет усечено, чтобы гарантировать, что мантисса имеет только 23 бита, поэтому окончательный результат будет -202014160.

Так что меняем тип цели на float, и других результатов не будет. Или просто измените оба операнда на тип double.

В конце разбора посмотрим на ассемблерный код, сгенерированный режимами FPU и SSE:

  • FPU
.cfi_startproc
pushq   %rbp
.cfi_def_cfa_offset 16
.cfi_offset 6, -16
movq    %rsp, %rbp
.cfi_def_cfa_register 6
subq    $16, %rsp
movl    $0x479de300, %eax
movl    %eax, -16(%rbp)
movl    $0xc51c3000, %eax
movl    %eax, -12(%rbp)
flds    -16(%rbp)
fmuls   -12(%rbp)
fstpl   -8(%rbp)
movl    $.LC2, %eax
movsd   -8(%rbp), %xmm0
movq    %rax, %rdi
movl    $1, %eax
call    printf  
leave
.cfi_def_cfa 7, 8
ret
.cfi_endproc
скопировать код
  • SSE
.cfi_startproc
pushq   %rbp
.cfi_def_cfa_offset 16
.cfi_offset 6, -16
movq    %rsp, %rbp
.cfi_def_cfa_register 6
subq    $16, %rsp
movl    $0x479de300, %eax
movl    %eax, -16(%rbp)
movl    $0xc51c3000, %eax
movl    %eax, -12(%rbp)
movss   -16(%rbp), %xmm0
mulss   -12(%rbp), %xmm0
unpcklps        %xmm0, %xmm0
cvtps2pd        %xmm0, %xmm0
movsd   %xmm0, -8(%rbp)
movl    $.LC2, %eax
movsd   -8(%rbp), %xmm0
movq    %rax, %rdi
movl    $1, %eax
call    printf
leave
.cfi_def_cfa 7, 8
ret
.cfi_endproc
скопировать код

Очевидно, что одна из них — серия flds/fmuls/fstpl, а другая — серия movss/mulss/movsd.

Вопрос 6

Сравните время работы следующих двух программ: test.c

float x=1.1;
float z=1.123;
float y=x;
for(int j=0;j<90000000;j++)
{
    y*=x;
    y/=z;
    y+=0.1f;
    y-=0.1f;
}

test1.c

float x=1.1;
float z=1.123;
float y=x;
for(int j=0;j<90000000;j++)
{
    y*=x;
    y/=z;
    y+=0;        //diffenent
    y-=0;        //different
}

Время выполнения второй программы примерно в 10 раз больше, чем первой, почему такой разрыв?

Мы печатаем j при j=50000000 и в конце цикла: Результат первой программы:

0.0000001751516833792265970259904861450195312500000000000000000000
0.0000001788139343261718750000000000000000000000000000000000000000
скопировать код

Результат второй программы:

0.0000000000000000000000000000000000000000000630584308946167681916
0.0000000000000000000000000000000000000000000630584308946167681916
скопировать код

Почему возникает такой результат? (Напоминаем, что двоичное представление 0.1 равно 1.10011001100110011001101 * 2^-4, вы можете оглянуться на ранее проанализированное 3.14f + 1e10f - 1e10f != 0.0f. Оставьте пустое место для размышления).
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
^-^
Поскольку показатель степени 0,1 равен -4, при сложении и вычитании чисел с плавающей запятой показатель степени будет выравниваться от меньшего порядка к большему. Поэтому, когда программа запустится до определенной степени, появится форма, похожая на следующую:

0.00000000...(很多0)..1yyy   * 2^-4   +
1.xxxxxxxxx                  * 2^-4
скопировать код

Из соображений точности и округления окончательный результат может быть 1.xxxx...1yy * 2^-4, слагаемое может быть отброшено с некоторыми значащими цифрами, но поскольку показатель степени по-прежнему равен -4, результат по-прежнему является нормализованным десятичный . -0.1f то же самое, результат по-прежнему нормированный десятичный.

Но вторая программа, поскольку нет нормализованного десятичного показателя степени ограничения 0,1, окончательный результат будет становиться все меньше и меньше и, наконец, выродится в денормализованное десятичное число: 0 00000000 xxxxxxxxxxxx.

Когда компьютер вычисляет денормализованные десятичные дроби, он примерно в 10-100 раз медленнее, чем вычисление нормализованных десятичных дробей.Почему это происходит? У автора мало знаний, а знания в железе и схемах не систематизированы и не всегда понятны. (Знание аппаратного обеспечения и схем всегда было для меня проблемой, часто мешавшей мне глубже вникнуть в проблему.)

Чтобы узнать, как избежать денормализованных десятичных знаков, ознакомьтесь со следующими двумя статьями:

На моем компьютере используйте

#include <xmmintrin.h>
_mm_setcsr( _mm_getcsr() | 0x8040 );
скопировать код

Выполнение второй программы можно сделать таким же, как и первой программы. Но точность намного хуже.После добавления двух приведенных выше предложений текущий результат программы равен 0,0, а результат разрешения денормализованной десятичной операции равен 0,0000000000000000000000000000000000000000000000630584308946167681916.

Кратко опишите, почему две приведенные выше строки избегают денормализованных десятичных вычислений: В наборах инструкций SSE и SSE2 есть два режима FTZ (Flush-to-Zero) и DAZ (Denormals-Are-Zero) Скорость работы с нормализованными десятичными знаками. Где FTZ означает, что когда результат операции дает денормализованное десятичное число, режим FTZ устанавливает результат операции равным 0. И DAZ означает, что когда операнд имеет денормализованное десятичное число, он сначала устанавливается в 0, а затем обрабатывается другими операндами. Можно также сказать, что FTZ влияет на результат вывода, а DAZ влияет на результат ввода.

image

Чтобы режим FTZ вступил в силу, необходимо выполнить два условия:

  • Бит 15 регистра MXCSR (бит FTZ) равен 1.
  • Бит исключения потери значимости, то есть 15-й бит равен 1 (исключение потери значимости).

Чтобы режим DAZ вступил в силу, необходимо выполнить два условия:

  • Бит 6 регистра MXCSR (бит DAZ) равен 1.
  • Процессор должен поддерживать режим DAZ.

Приведенный выше код сначала использует _mm_getcsr() для получения значения регистра MXCSR, и мы печатаем результат, чтобы увидеть:

printf("%d",_mm_getcsr());    //8064  二进制1111110000000

Бит 11 равен 0. Но бит FTZ и бит DAZ равны 0, поэтому мы или 0x8040 (1000000001000000) можем заставить действовать режимы FTZ и DAZ.

(Заканчивать)

Использованная литература:

Категории