2.1. Числа с плавающей точкой#
В математическом анализе мы имеем дело с действительным числами \(\real\). Однако, на практике мы используем числа с плавающей точкой \(\floatset \subset \real\), которые являются подмножеством рациональных чисел.
Эти два множества отличаются.
Множество действительных чисел неограничено и непрерывно.
Множество чисел с плавающей точкой ограничено и дискретно.
Существуют разные типы (множества) чисел с плавающей точкой. Например,
binary32 – числа одинарной точности
binary64 – числа двойной точности
extended precision – 80-битные числа
bfloat16 – используются в машинном обучении
Они задокументированы в виде стандартов, а реализованы в floating-point unit (FPU) процессоров или программно. Реализация может отличаться от документации, но важно одно.
Примечание
Числа с плавающей точкой имеют конечную точность, их арифметика отличается от арифметики действительных чисел. Чтобы избежать ошибок при программировании, требуется соблюдать меры предосторожности.
Цель этого раздела: разобраться с тем, что такое float-числа и как с ними работать.
2.1.1. Множество чисел с плавающей точкой#
Все типы чисел с плавающей точкой определяются одинаково.
Множество чисел с плавающей точкой \(\floatset\) состоит из нуля и чисел вида
где целое число \(n\in\integer\) называют экспонентой, а \(1+f\) мантиссой, при этом \(f\) имеет вид
где натуральное \(d\) называют (в данном случае) двоичной точностью.
Можно сказать, что \(f\) (2.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})\). Но это уже технические детали.
Рассмотрим числа с точностью \(d = 3\). Для \(n=0\) получим \(8\) чисел
которые заключены в полуинтервале \([1,2)\) с промежутком \(2^0 2^{-3} = 1/8\).
В свою очередь, для \(n=1\) получим \(8\) чисел в полуинтервале \([2, 4)\) с промежутком \(2^1 2^{-3} = 1/4\)
Аналогично, для \(n=-1\) получим \(8\) чисел в полуинтервале \([1/2,1)\) с промежуком \(2^{-1}2^{-3}=1/16\).
Из примера видно, что наименьшее число с плавающей точкой, большее \(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)\). Относительная ошибка приближения тогда
В зависимости от того, какая требуется точность арифметики для решения задачи, можно выбирать те или иные мантиссу и экспоненту. Например, у стандартов binary16 и bfloat16 одинаковое общее число бит \(n + d + 1\) (один бит на знак), но распределение бит на мантиссу \(d\) и экспоненту \(n\) разные. В некоторых задачах хватает точности чисел binary32 и их можно использовать вместо binary64, а это даёт прирост быстродействия примерно в два раза (см. SIMD).
2.1.2. Числа двойной точности#
Число с плавающей точкой двойной точности Float64
определил стандарт IEEE 754-1985.
Официальное название – IEEE 754 binary64.
Для Float64
отводится 64 бита: один под знак числа, 11 под экспоненту и 52 под дробную часть мантиссы \(f\).
s eeeeeeeeeee mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm
↑ экспонента мантисса
знак
1 бит 11 бит 52 бита
Таким образом для Float64
\(d=52\) (2.2), а машинная точность \(\macheps\) составляет
Экспонента \(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
в вычислениях.
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. Особенности арифметики чисел с плавающей точкой#
Конечная точность чисел с плавающей точкой приводит к ошибкам округления. Рассмотрим на примерах ниже.
Oшибки округления
Пусть требуется сложить \(2.25 + 1.125\) c двоичной точностью \(d=3\) Тогда сложение представится в виде
Чтобы его выполнить, необходимо сначала выравнять экспоненты слагаемых
На этом этапе происходит обрезание \(1/16\), поскольку \(1/16\) меньше точности \(1/8\), с которой можно задать мантиссу.
Когда экспонентны выравнены, можно сложить мантиссы \((1 + 1/8) + (1/2)\). Получаем результат
который не равен сложению в действительныйх числах \(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
Ошибки округления копятся, например, можно вычислить разницу
При работе с действительными числами ожидался бы ноль, но для 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\).
Ошибку округления можно уменьшить, если изменить порядок сложения
@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})\)
@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 числа
2.1.6. Дополнительные материалы#
Детальную справку о числах с плавающей точкой можно получить в работе Давида Голдберга (David Goldberg) [Gol91].