Светлая тема /
тёмная тема
Запятая, Карл!
Как устроены типы float и double
Плавающая запятая даёт много преимуществ, но «заплыв» этот полон моментов, требующих понимания и тщательного подхода. Числа с плавающей запятой используются слишком часто, чтобы пренебрегать их глубоким пониманием. Поэтому нужно хорошо разбираться, как представлены числа с плавающей запятой и как работают вычисления с ними.
Экскурс
Число с плавающей запятой представлено в следующем виде:
$$N = [s] \cdot m \cdot 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. Каждому порядку соответствует одинаковое количество значений мантиссы, следовательно — одинаковое количество представимых чисел:
- порядок \(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 \times 10^{37}...1{,}7014117 \times 10^{38}\),
- порядок \(127\), значения: \(1{,}7014118 \times 10^{38}...3{,}4028235 \times 10^{38}\).
Каждый из этих диапазонов разбивается на равное количество интервалов. Следовательно, от \(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♯.