Плавающая запятая даёт много преимуществ, но «заплыв» этот полон моментов, требующих понимания и тщательного подхода. Числа с плавающей запятой используются слишком часто, чтобы пренебрегать их глубоким пониманием. Поэтому нужно хорошо разбираться, как представлены числа с плавающей запятой и как работают вычисления с ними.

Экскурс

Число с плавающей запятой представлено в следующем виде:

$$N=(-1)^{s}\cdot m \cdot 2^{e},$$

где $m$ — мантисса (23 бита — дробная часть мантиссы и неявно заданный ведущий бит, всегда равный единице, поскольку мантисса хранится в нормализованном виде), $e$ — смещённая экспонента/порядок (8 бит). Один бит отводится под знак ($s$, бит равен нулю — число положительное, единице — отрицательное).

Тип double полностью аналогичен float, он лишь содержит больше битов мантиссы и порядка.

Рассмотрим какой-нибудь пример. Число $3{,}5$ будет представлено в следующем виде:

Бит знака равен нулю, это говорит о том, что число положительное.

Смещённая экспонента равна $128_{10}$ ($10000000_{2}$), смещена она на $127$, следовательно, в действительности получаем:

$$e = e_{\mathrm{shifted}} - 127 = 128 - 127 = 1.$$

Нормализованная мантисса равна $[1{,}]11000000000000000000000_{2}$, то есть $1{,}75_{10}$:

$$N = m \cdot 2^{e} = 1{,}75 \cdot 2^{1} = 3{,}5.$$

Инвертируем бит знака:

Теперь число стало отрицательным: $-3{,}5$.

Рассмотрим несколько исключений. Как представить число $0$, если мантисса всегда хранится в нормализованном виде, то есть больше либо равна $1$? Решение такое — условимся обозначать $0$, приравнивая все биты смещённого порядка и мантиссы нулю:

Ноль бывает и отрицательным:

Плюс и минус бесконечность обозначаются приравниванием всех битов смещённого порядка единице, а битов мантиссы — нулю:

В некоторых случаях результат не определён, для этого есть метка $\mathrm{NaN}$, что значит not a number (не число). Приравняем все биты смещённого порядка и мантиссы единице (по правде говоря, достаточно того, чтобы хотя бы один бит мантиссы равнялся единице):

Подводим небольшой итог. Возьмём все возможные сочетания значений всех 32 битов, уберём те варианты, которые зарезервированы для $\pm\infty$ и метки $\mathrm{NaN}$, и получим все числа, которые можно хранить в переменной типа float. Их количество очень велико, но пусть это никого не обольщает; из-за особенностей представления, описанных выше, не все числа точно представимы и не все операции выполняются так, как этого можно было бы ожидать.

Более подробно о типах с плавающей запятой можно прочитать в обучающих статьях. Описанный подход к представлению чисел является компромиссом между точностью, диапазоном значений и скоростью выполнения операций с числами с плавающей запятой. Во-первых, число хранится в двоичной форме, но не всякое десятичное число можно представить в виде непериодической дроби в двоичной системе; это одна из причин потери точности ($0{,}1_{10} = 0{,}0(0011)_{2}$). Во-вторых, значащих цифр (цифр мантиссы) не так уж и много: 24 двоичных цифры, или 7,2 десятичных, поэтому придётся много округлять. Но есть и менее очевидные моменты.

Нарушение свойства ассоциативности

Вычислим $\sin(x)$, разложив эту функцию в ряд Тейлора:

$$\sin(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!} - ...$$

Напишем небольшую функцию, которая будет реализовывать эти вычисления. Тип float вместо double выбран для того, чтобы показать, как быстро накапливается погрешность; никаких дополнительных действий не производится, код исключительно демонстрационный. Целью не является показать, что семи членов ряда недостаточно, цель — показать нарушение свойства ассоциативности операций:

$$a+(b+c)\neq(a+b)+c,$$

хотя свойство коммутативности сохраняется:

$$a + b = b + a.$$

float sine(float x)

{

    const int n = 13;

    float sine = x;

    bool isNegative = true;

    for (int i = 3; i <= n; i += 2)

    {

        sine += (isNegative ?

            - (pow(x, i) / Factorial(i))

            : (pow(x, i) / Factorial(i)));

        isNegative = !isNegative;

    }

    return sine;

}

Теперь напишем аналогичную функцию, но сделаем так, чтобы те же элементы ряда суммировались в обратном порядке.

float sineReversed(float x)

{

    const int n = 13;

    float sine = 0;

    bool isNegative = false;

    for (int i = n; i >= 3; i -= 2)

    {

        sine += (isNegative ?

            - (pow(x, i) / Factorial(i))

            : (pow(x, i) / Factorial(i)));

        isNegative = !isNegative;

    }

    sine += x;

    return sine;

}

Запустим.

#include <iostream>

#include <math.h>

 

int main()

{

    float sine1 = sine(3.1415926f);

    float sine2 = sineReversed(3.1415926f);

 

    std::cout << sine1 << std::endl << sine2 << std::endl;

}

Посмотрим на это с точки зрения машины:

$\sin(\pi) \approx 3{,}141592 - 5{,}16771$ $+\ 2{,}55016$ $-\ 0{,}599265$ $+\ 0{,}0821459$ $-\ 0{,}00737043$ $+\ 0{,}000466303,$

$\sin(\pi) \approx 0{,}000466303 - 0{,}00737043$ $+\ 0{,}0821459$ $-\ 0{,}599265$ $+\ 2{,}55016$ $-\ 5{,}16771$ $+\ 3{,}141592.$

До этого момента мнение компьютера и человека, казалось бы, сходится. Но результатом в первом случае окажется $0{,}0000211327$ ($2.11327\mathrm{E}{-}005$), а во втором случае — $0{,}0000212193$ ($2.12193\mathrm{E}{-}005$), совпадают лишь две значащие цифры!

Разгадка проста: у чисел представленного ряда шесть различных (двоичных) порядков: $1$, $2$, $1$, $-1$, $-4$, $-8$, $-12$. Когда складываются (вычитаются) два числа одного порядка или близких порядков, потери точности, как правило, небольшие. Если бы мы складывали огромное число и много маленьких чисел одного порядка, то мы заметили бы, что лучше в плане точности сперва сложить все маленькие числа, а затем уже прибавить большое. Рассмотрим обратный сценарий: сложим большое и первое маленькое число; поскольку порядки значительно различаются, маленькое число будет (фигурально выражаясь) «раздавлено» большим из-за приведения порядков; получилось новое большое число, не очень точное, но пока ещё достаточно близкое к точному результату; к получившемуся большому числу прибавляем второе маленькое, порядки снова значительно различаются, снова маленькое число оказывается раздавленным, уже две «жертвы». И так далее. Погрешность накопилась достаточно большая.

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

Большие числа

Множество чисел, представимых с помощью float (и double), конечно. Из этого, а также из того, что количество разрядов мантиссы весьма ограничено, сразу же следует, что очень большие числа представить не получится; придётся изобретать пути обхода переполнения. Менее очевидным следствием из способа задания числа в виде, описанном выше, является то, что пространство таких чисел не является равномерным. Что это значит?

Возьмём порядок $e = 0$ и рассмотрим два стоящих подряд числа $N_{1}$ и $N_{2}$: у первого мантисса $m_{1}=[1{,}]00000000000000000000001_{2}$, у второго — $m_{2}=[1{,}]00000000000000000000010_{2}$.

Воспользуемся конвертером: $N_{1}=1{,}0000001_{10}$, $N_{2} = 1{,}0000002_{10}$. Разница между ними равна $0{,}0000001_{10}$.

Теперь возьмём порядок $e = 127$ и снова рассмотрим два стоящих подряд числа: у первого мантисса $m_{3}=[1{,}]00000000000000000000001_{2}$, у второго — $m_{4}=[1{,}]00000000000000000000010_{2}$.

$N_{3}=1{,}701412 \cdot 10^{38}_{10}$, $N_{4} = 1{,}7014122 \cdot 10^{38}_{10}$. Разница между ними равна $0{,}0000002 \cdot 10^{38}_{10}$. В рамках этого типа данных нет никакого способа задать некоторое число $N$, находящееся в интервале между $N_{3}$ и $N_{4}$. На секунду, $0{,}0000002 \cdot 10^{38}$ — это 20 нониллионов, иначе говоря, двойка и 31 ноль! Маленький шаг для мантиссы и огромный скачок для всего числа.

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

• порядок $0$, значения: $1{,}0...1{,}9999999$,
• порядок $1$ значения: $2{,}0...3{,}9999998$,
• порядок $2$, значения: $4{,}0...7{,}9999995$,
• порядок $3$, значения: $8{,}0...15{,}999999$,
• порядок $4$, значения: $16{,}0...31{,}999998$,
• ...
• порядок $126$, значения: $8{,}507059 \cdot 10^{37}...1{,}7014117 \cdot 10^{38}$,
• порядок $127$, значения: $1{,}7014118 \cdot 10^{38}...3{,}4028235 \cdot 10^{38}$.

Каждый из этих диапазонов разбивается на равное количество интервалов. Следовательно, от $1{,}0$ до $1{,}9999999$, от $16{,}0$ до $31{,}999998$, от $1{,}7014118 \cdot 10^{38}$ до $3{,}4028235 \cdot 10^{38}$ находится одинаковое количество доступных значений — $(2^{23} - 2)$ (поскольку мантисса, за исключением скрытого бита, имеет $23$ бита; минус сами границы).

Всё сказанное в равной степени касается отрицательных чисел и чисел с отрицательными порядками.

Заключение

Когда в программе используются числа с плавающей запятой, необходимо обращать внимание на операции сложения и вычитания из-за нарушения свойства ассоциативности вследствие ограниченной точности типа float (и double). Особенно в том случае, когда порядки чисел сильно различаются. Другой проблемой является большой шаг между ближайшими представимыми числами больши́х порядков. Что и было показано.

Здесь были освещены только эти два аспекта. Также следует помнить о том, что нельзя непосредственно сравнивать два числа с плавающей запятой (if (x == y)); что у диапазона есть ограничения; что следует использовать логарифм, когда происходят вычисления с использованием огромных чисел. Ну и совсем не следует забывать о бесконечностях и метке $\mathrm{NaN}$. Из экзотики — флаги компилятора, касающиеся точности промежуточных результатов, из простых вещей (хотя и не менее коварных, чем что-либо упомянутое выше) — особенности приведения типов. В общем, в выражении «вагон и маленькая тележка» эта заметка — даже не тележка, а маленькая горстка. Маленькая, но достаточно важная.

И последнее. Когда требуется точность вычислений (например, при финансовых расчётах), следует использовать собственные или сторонние классы чисел с фиксированной запятой. Если речь идёт о языке, поддерживающем такие типы, следует использовать их. Таковым, например, является тип decimal в языке C♯‎.