Запятая, Карл!

Как устроены типы float и double

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

Экскурс

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

N = [s] × m × 2e,

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

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

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

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

Смещённая экспонента равна 12810 (100000002), смещена она на 127, следовательно, в действительности получаем:

e = eshifted − 127 = 128 − 127 = 1.

Нормализованная мантисса равна [1,]110000000000000000000002, то есть 1,7510:

N = m × 2e = 1,75 × 21 = 3,5.

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

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

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

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

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

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

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

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

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

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

sin x = x
x3
3!
+
x5
5!
x7
7!
+
x9
9!
x11
11!
+
x13
13!
− ...

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

a + (b + c) ≠ (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 π ≈ 3,141592 − 5,16771 + 2,55016 − 0,599265 + 0,0821459 − 0,00737043 + 0,000466303,

sin π ≈ 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 и рассмотрим два стоящих подряд числа N1 и N2: у первого мантисса m1 = [1,]000000000000000000000012, у второго — m2 = [1,]000000000000000000000102.

Воспользуемся конвертером: N1 = 1,000000110, N2 = 1,000000210. Разница между ними равна 0,000000110.

Теперь возьмём порядок e = 127 и снова рассмотрим два стоящих подряд числа: у первого мантисса m3 = [1,]000000000000000000000012, у второго — m4 = [1,]000000000000000000000102.

N3 = 1,701412×103810, N4 = 1,7014122×103810. Разница между ними равна 0,0000002×103810. В рамках этого типа данных нет никакого способа задать некоторое число N, находящееся в интервале между N3 и N4. На секунду, 0,0000002×1038 — это 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×1037...1,7014117×1038,
порядок 127, значения: 1,7014118×1038...3,4028235×1038.

Каждый из этих диапазонов разбивается на равное количество интервалов. Следовательно, от 1,0 до 1,9999999, от 16,0 до 31,999998, от 1,7014118×1038 до 3,4028235×1038 находится одинаковое количество доступных значений — 223 − 2 (поскольку мантисса, за исключением скрытого бита, имеет 23 бита; минус сами границы).

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

Заключение

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

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

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

15 июля 2015 · float