2.1. Числа с плавающей точкой

В математическом анализе мы имеем дело с действительным числами \(\real\). Однако, на практике мы используем числа с плавающей точкой \(\floatset \subset \real\), которые являются подмножеством рациональных чисел.

Эти два множества отличаются.

• Множество действительных чисел неограничено и непрерывно.

• Множество чисел с плавающей точкой ограничено и дискретно.

Существуют разные типы (множества) чисел с плавающей точкой. Например,

• binary32 – числа одинарной точности

• binary64 – числа двойной точности

• extended precision – 80-битные числа

• bfloat16 – используются в машинном обучении

Они задокументированы в виде стандартов, а реализованы в floating-point unit (FPU) процессоров или программно. Реализация может отличаться от документации, но важно одно.

Примечание

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

Цель этого раздела: разобраться с тем, что такое float-числа и как с ними работать.

2.1.1. Множество чисел с плавающей точкой

Все типы чисел с плавающей точкой определяются одинаково.

Определение 2.1 (Число с плавающей точкой)

Множество чисел с плавающей точкой \(\floatset\) состоит из нуля и чисел вида

(2.1)\[\pm (1 + f) \times 2^n,\]

где целое число \(n\in\integer\) называют экспонентой, а \(1+f\) мантиссой, при этом \(f\) имеет вид

(2.2)\[f = \sum_{i=1}^d b_i \, 2^{-i}, \quad b_i\in\{0,1\},\]

где натуральное \(d\) называют (в данном случае) двоичной точностью.

Можно сказать, что \(f\) (2.2) является дробной частью мантиссы.

Используя нотацию из позиционных систем счисления, можно записать в двоичной системе

\[1 + f = 1 + 0.\overline{b_1 b_2 ... b_d} = 1.\overline{b_1 b_2 ... b_d} \in [1,2).\]

Отсюда видно, что ближайшие мантисcы отстоят друг от друга на \(2^{-d}\) (последний бит \(b_d\)). Поэтому, между числами \(2^n\) и \(2^{n+1}\) находится \(2^d\) равномерно распределённых чисел. В свою очередь, отсюда следует, что числа с плавающей точной распределены неравномерно на числовой оси.

Почему мантисса в виде \(1 + f\)?

Попробуем задать число \(1/2\) в формате с плавающей точкой. Принципиально есть два варианта: \(\overline{1.0}_{2} \times 2^{-1}\) и \(\overline{0.1}_{2} \times 2^{0}\). Первый вариант более разумен – более привычней считать, что первый бит мантиссы по умолчанию равняется 1, а не второй бит. В реализации формата этот бит мантиссы вовсе не хранится и диапазон представимых форматом чисел становится шире.

Float-числа с неявной единицей в мантиссе называют нормализованными. Однако вблизи нуля для сохранения «плотности» чисел от правила неявной единицы отказываются, и такие числа называются денормализованными. В случае формата binary64 денормализованные числа заполняют интервал \((0, 2^{-1022})\). Но это уже технические детали.

Пример 2.2

Рассмотрим числа с точностью \(d = 3\). Для \(n=0\) получим \(8\) чисел

\[ 1, \: 1 + \frac{1}{8}, \: 1 + \frac{1}{4}, \: 1 + \frac{1}{4} + \frac{1}{8}, \: 1 + \frac{1}{2}, \: 1 + \frac{1}{2} + \frac{1}{8}, \: 1 + \frac{1}{2} + \frac{1}{4}, \: 1 + \frac{1}{2} + \frac{1}{4} + \frac{1}{8} = 1 + \frac{7}{8}. \]

которые заключены в полуинтервале \([1,2)\) с промежутком \(2^0 2^{-3} = 1/8\).

В свою очередь, для \(n=1\) получим \(8\) чисел в полуинтервале \([2, 4)\) с промежутком \(2^1 2^{-3} = 1/4\)

\[ (1+0)\times 2^1 = 2, \quad \bigg(1 + \frac{1}{8}\bigg) \times 2^1 = 2 + \frac{1}{4} , \quad ..., \quad \bigg(1 + \frac{1}{2} + \frac{1}{4} + \frac{1}{8}\bigg) \times 2^1 = 3 + \frac{3}{4}. \]

Аналогично, для \(n=-1\) получим \(8\) чисел в полуинтервале \([1/2,1)\) с промежуком \(2^{-1}2^{-3}=1/16\).

Для чисел с плавающей точкой верно равенство \(1 + \macheps / 2 = 1\).

Из примера видно, что наименьшее число с плавающей точкой, большее \(1\) это \(1 + 2^{-d}\). Разницу между этими числами \(2^{-d}\) называют машинным нулём (машинным эпсилон, машинной точностью) \(\macheps\).

Поскольку числа с плавающей точкой имеют конечную точность, сделаем оценку лучшего приближения. Пусть \(\float(x)\in\floatset\) ближайшее к действительному \(x\in\real\) число с плавающей точкой. Расстояние для чисел в промежутке \([2^n, 2^{n+1})\) равняется \(2^n \macheps = 2^{n-d}\). Поэтому число \(x\) находится не далее, чем на \(2^{n-d} / 2 = 2^{n-d-1}\) от \(\float(x)\). Относительная ошибка приближения тогда

\[\frac{|\float(x)-x|}{|x|} \le \frac{2^{n-d-1}}{2^n} \le \frac{1}{2}\macheps.\]

В зависимости от того, какая требуется точность арифметики для решения задачи, можно выбирать те или иные мантиссу и экспоненту. Например, у стандартов binary16 и bfloat16 одинаковое общее число бит \(n + d + 1\) (один бит на знак), но распределение бит на мантиссу \(d\) и экспоненту \(n\) разные. В некоторых задачах хватает точности чисел binary32 и их можно использовать вместо binary64, а это даёт прирост быстродействия примерно в два раза (см. SIMD).

2.1.2. Числа двойной точности

Последняя версия стандарта – IEEE 754-2019.

Поскольку в двоичном представлении первая цифра мантиссы нормализованного числа всегда 1, то её смысла хранить нет.

Число с плавающей точкой двойной точности Float64 определил стандарт IEEE 754-1985. Официальное название – IEEE 754 binary64.

Для Float64 отводится 64 бита: один под знак числа, 11 под экспоненту и 52 под дробную часть мантиссы \(f\).

  s eeeeeeeeeee mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm
  ↑  экспонента                        мантисса
знак
1 бит  11 бит                           52 бита

Таким образом для Float64 \(d=52\) (2.2), а машинная точность \(\macheps\) составляет

(2.3)\[\macheps = 2^{-52} \approx 2.2 \times 10^{-16}.\]

Экспонента \(n\) (2.1) принимает значения \(-1022 \le n \le 1023\). Максимальное значение в Float64 составляет почти \(2^{1024} \approx 2 \times 10^{308}\). Числа крупнее принимают специальное значение-бесконечность Inf. Наименьшим положительным числом является \(2^{-1022} \approx 2 \times 10^{-308}\), меньшие его становятся нулём. Отрицательные аналоги работают схожим образом.

Стоит отметить, что нулей в Float64, как и бесконечностей, два: +0 и -0.

Да, на ноль делить можно.

Кроме того, существует ещё одно специальное значение Not a Number, NaN. Оно появляется в результате неопределённых операций: ±0/±0, ±Inf/±Inf и операций, в которых участвует NaN, что приводит к распространению NaN в вычислениях.

Демонстрация 2.3

Float64 в Julia

Битовое представление Float64

@show typeof(2.784);

typeof(2.784) = Float64

bitstring(2.784)

"0100000000000110010001011010000111001010110000001000001100010010"

Функции nextfloat и prevfloat возвращают ближайшие следующее и предыдущее числа

[bitstring(prevfloat(1.0)), bitstring(1.0), bitstring(nextfloat(1.0))]

3-element Vector{String}:
 "0011111111101111111111111111111111111111111111111111111111111111"
 "0011111111110000000000000000000000000000000000000000000000000000"
 "0011111111110000000000000000000000000000000000000000000000000001"

[prevfloat(1.0), 1.0, nextfloat(1.0)]

3-element Vector{Float64}:
 0.9999999999999999
 1.0
 1.0000000000000002

Определённые в (2.1) экспоненту и мантиссу также можно получить

x = 2.784
@show sign(x), exponent(x), significand(x);

(sign(x), exponent(x), significand(x)) = (1.0, 1, 1.392)

Расстояние между числами в диапазоне \([2^n, 2^{n+1})\) можно получить с помощью функции eps(x). По умолчанию, eps() возвращает значение для единицы 1.0, что является машинным нулём \(\macheps\)

@show eps();
@show log2(eps());

eps() = 2.220446049250313e-16
log2(eps()) = -52.0

Пример вызова для для числа, отличного от 1.0

@show eps(3.0);
@show prevfloat(3.0);
@show nextfloat(3.0);

eps(3.0) = 4.440892098500626e-16
prevfloat(3.0) = 2.9999999999999996
nextfloat(3.0) = 3.0000000000000004

Однако, стоит отметить, что eps(x) определяется как

eps(x) == max(x-prevfloat(x), nextfloat(x)-x)

например, для числа 2.0 проявится то, что оно находится на границе диапазонов \([1, 2), [2, 4)\)

@show eps(2.0);
@show prevfloat(2.0);
@show nextfloat(2.0);

eps(2.0) = 4.440892098500626e-16
prevfloat(2.0) = 1.9999999999999998
nextfloat(2.0) = 2.0000000000000004

Наименьшее и наибольшее числа (но не Inf или 0) получить можно так

@show floatmin(), floatmax();

(floatmin(), floatmax()) = (2.2250738585072014e-308, 1.7976931348623157e308)

Специальные значения можно получить с помощью констант Inf и NaN

[bitstring(Inf), bitstring(-Inf), bitstring(NaN)]

3-element Vector{String}:
 "0111111111110000000000000000000000000000000000000000000000000000"
 "1111111111110000000000000000000000000000000000000000000000000000"
 "0111111111111000000000000000000000000000000000000000000000000000"

2.1.3. Особенности арифметики чисел с плавающей точкой

Конечная точность чисел с плавающей точкой приводит к ошибкам округления. Рассмотрим на примерах ниже.

Пример 2.4

Oшибки округления

Пусть требуется сложить \(2.25 + 1.125\) c двоичной точностью \(d=3\) Тогда сложение представится в виде

\[ \float(2.25) + \float(1.125) = \bigg(1+\frac{1}{8}\bigg)\times 2^1 + \bigg(1+\frac{1}{8}\bigg)\times 2^0. \]

Чтобы его выполнить, необходимо сначала выравнять экспоненты слагаемых

\[ \require{cancel} \bigg(1+\frac{1}{8}\bigg)\times 2^1 + \bigg(\frac{1}{2}+\cancel{\frac{1}{16}}\bigg)\times 2^1. \]

На этом этапе происходит обрезание \(1/16\), поскольку \(1/16\) меньше точности \(1/8\), с которой можно задать мантиссу.

Когда экспонентны выравнены, можно сложить мантиссы \((1 + 1/8) + (1/2)\). Получаем результат

\[ \bigg(1+ \frac{1}{2}+\frac{1}{8}\bigg)\times 2^1 = 3.25, \]

который не равен сложению в действительныйх числах \(2.25 + 1.125 = 3.375\).

В случае Float64 демонстрация ошибок округления может быть такой

@show 2.0 + nextfloat(1.0);
@show nextfloat(1.0);

2.0 + nextfloat(1.0) = 3.0
nextfloat(1.0) = 1.0000000000000002

Ошибки округления копятся, например, можно вычислить разницу

\[ \frac{k}{10} - \sum_{i=1}^{k} 0.1. \]

При работе с действительными числами ожидался бы ноль, но для float-чисел это не так.

ks = [10.0, 100.0, 1000.0, 10000.0]
@show Σ = [sum(0.1 for _ in 1:k) for k in ks]
@show exact = 0.1 .* ks;

Σ = [sum((0.1 for _ = 1:k)) for k = ks] = [0.9999999999999999, 9.99999999999998, 99.9999999999986, 1000.0000000001588]
exact = 0.1 .* ks = [1.0, 10.0, 100.0, 1000.0]

При этом относительная ошибка копится

(exact .- Σ) ./ exact

4-element Vector{Float64}:
  1.1102230246251565e-16
  1.9539925233402757e-15
  1.4068746168049984e-14
 -1.588205122970976e-13

Ошибка возникает из-за того, что при сложении переменная-аккумулятор acc растёт и меняет «диапазон точности» \([2^{n-1}, 2^{n})\), а инкремент 0.1 не меняется. В какой-то момент можно вовсе прийти к ситуации acc + 0.1 == acc, это наступит, когда \(0.1/\text{acc} < \macheps/2\).

Так и поступают для вычисления рядов, например, \(e^x = \sum_{n=0}^\infty x^n /n!\).

Также, см. алгоритм Кэхэна.

Ошибку округления можно уменьшить, если изменить порядок сложения

@show 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1;
@show (0.1 + 0.1) + (0.1 + 0.1) + (0.1 + 0.1) + (0.1 + 0.1) + (0.1 + 0.1);

0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 = 0.9999999999999999
(0.1 + 0.1) + (0.1 + 0.1) + (0.1 + 0.1) + (0.1 + 0.1) + (0.1 + 0.1) = 1.0

Другими словами, сложение для чисел с плавающей точкой неассоциативно.

Более прикладным примером является работа с числами разных порядков

ϵ = 1e-6
@show (1.0 + ϵ)^2;
@show 1.0 + 2.0*ϵ + ϵ^2;
@show 1.0 + (2.0*ϵ + ϵ^2);

(1.0 + ϵ) ^ 2 = 1.000002000001
1.0 + 2.0ϵ + ϵ ^ 2 = 1.0000020000010001
1.0 + (2.0ϵ + ϵ ^ 2) = 1.000002000001

Машинный ноль (2.3) на практике означает, что при сложении чисел двойной точности, отличающихся примерно в \(10^{16}\) раз или более, одно из слагаемых потеряется

@show 1.0 + 1e-15;
@show 1.0 + 1e-16;
@show 1.0 + 1e15;
@show 1.0 + 1e16;

1.0 + 1.0e-15 = 1.000000000000001
1.0 + 1.0e-16 = 1.0
1.0 + 1.0e15 = 1.000000000000001e15
1.0 + 1.0e16 = 1.0e16

Ещё один пример неассоциативности. Он возникает на границе \(2^n\) диапазонов \([2^{n-1}, 2^n)\) и \([2^n,2^{n+1})\)

В англоязычной литературе см. subtractive cancellation.

@show (1.0 + eps()/2) - 1.0;
@show 1.0 + (eps()/2 - 1.0);

(1.0 + eps() / 2) - 1.0 = 0.0
1.0 + (eps() / 2 - 1.0) = 1.1102230246251565e-16

2.1.4. Меры предосторожности и практика работы с float-числами

Точность представления чисел

julia> 0.1 + 0.2  # Float64
0.30000000000000004

julia> 0.1f0 + 0.2f0  # Float32
0.3f0

Не все числа представимы в float точно, в реальности программист имеет дело с ближайшим к действительному float-числу. Сложение float-чисел может быть выполнено точно, например, результат 0.30000000000000004 верный, просто сложение происходило для приближений, а не исходных чисел. В разных представлениях могут быть разные результаты, так для Float32 те же числа складываются иначе.

Сравнение чисел

Когда необходимо сравнивать числа, используйте неточное сравнение isapprox, при этом лучше использовать относительную погрешность, а не абсолютную.

Будьте аккуратнее в точных сравнениях. Например, вам необходимо найти отрицательную энергию, тогда сравнение вида E < 0 иногда не подходит, поскольку E может быть порядка \(-\macheps\): считать ли такое значение значение физически отрицательным – решать вам.

Числа сильно разных точностей

Работая с Float64, старайтесь не допускать сложений чисел, отличающихся примерно в \(10^{16}\) раз или более.

Сложение неассоциативно

Малое складывать с малым, а крупное с крупным.

Осторожнее с

• С суммами большого числа слагаемых (ряды, циклы);

• Возведением в степень.

Subtractive cancellation

Осторожнее с вычитанием чисел близких по значению. Примеры из Википедии. Операции происходят с конечной точностью, поэтому иногда ошибка вычитания приближений может сильно превышать вычитание точных чисел. Алгебраические преобразования с вычисляемым выражением могут помочь снизить ошибку.

2.1.5. Какие ещё бывают представления дробных чисел?

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

В случае Julia можете взглянуть на следующие библиотеки чисел.

• BigFloat (стандартная библиотека) – числа с плавающей точкой произвольной точности, обёртка библиотеки GNU MPFR

• Quadmath – реализация стандарта IEEE binary12, обёртка библиотеки GCC libquadmath

• Extended precision числа

  • DoubleFloats

  • MultiFloats

2.1.6. Дополнительные материалы

Детальную справку о числах с плавающей точкой можно получить в работе Давида Голдберга (David Goldberg) [Gol91].