Запятая, Карл!
Как устроены типы float и double

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

Экскурс

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

$$N = [s] \times m \times 2^{e},$$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Более подробно о типах с плавающей запятой можно прочитать в обучающих статьях. Описанный подход к представлению чисел является компромиссом между точностью, диапазоном значений и скоростью выполнения операций с числами с плавающей запятой. Во-первых, число хранится в двоичной форме, но не всякое десятичное число можно представить в виде непериодической дроби в двоичной системе; это одна из причин потери точности (\(0{,}1_{10} = 0{,}000(1100)_{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.11327e{-}005\)), а во втором случае — \(0{,}0000212193\) (\(2.12193e{-}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 \times 10^{38}_{10}\), \(N_{4} = 1{,}7014122 \times 10^{38}_{10}\). Разница между ними равна \(0{,}0000002 \times 10^{38}_{10}\). В рамках этого типа данных нет никакого способа задать некоторое число \(N\), находящееся в интервале между \(N_{3}\) и \(N_{4}\). На секунду, \(0{,}0000002 \times 10^{38}\) — это 20 нониллионов, иначе говоря, двойка и 31 ноль! Маленький шаг для мантиссы и огромный скачок для всего числа.

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

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

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

Заключение

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

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

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

16 июля 2015
© MMXIMMXX. RSS
Светлая тема / тёмная тема