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

Экскурс

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

$$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♯‎.