Представление чисел с плавающей точкой

Эта статья находится в разработке!

Плавающая точка

Определение:

Плавающая точка (floating point) - метод представления действительных чисел, при котором число хранится в виде мантиссы и показателя степени, а значение числа вычисляется по формуле: [math] x = (-1)^{sign} \times mant \times base^{exp} [/math], где [math] x [/math] - число, [math] sign [/math] - бит, отвечающий за знак числа, [math] mant [/math] - мантисса, [math] base [/math] - основание степени, [math] exp [/math] - показатель степени.

Такой метод является компромиссом между точностью и диапазоном представляемых значений. Представление чисел с плавающей точкой рассмотрим на примере чисел двойной точности (double precision). Такие числа занимают в памяти два машинных слова (8 байт на 32-битных системах). Наиболее распространенное представление описано в стандарте IEEE 754.

Кроме чисел двойной точности также используются следующие форматы чисел:

• половинной точности (half precision) (16 бит),
• одинарной точности (single precision) (32 бита),
• четверной точности (quadruple precision) (128 бит),
• расширенной точности (extended precision) (80 бит).

При выборе формата программисты идут на разумный компромисс между точностью вычислений и размером числа.

Нормальная и нормализованная формы

Определение:

Нормальной называется форма представления числа, при которой абсолютное значение мантиссы десятичного числа находится на полуинтервале [math] [0,1) [/math].

Недостатком такой записи является тот факт, что числа нельзя записать однозначно: [math] 0.01 = 0.001 \times 10^1 [/math].

Определение:

Нормализованной называется форма представления числа, при которой абсолютное значение мантиссы десятичного числа лежит на полуинтервале [math] [1, 10) [/math], а двоичного на полуинтервале [math] [1, 2) [/math].

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

Число с плавающей точкой хранится в нормализованной форме и состоит из трех частей (в скобках указано количество бит, отводимых на каждую секцию в формате double):

1. знак
2. экспонента (показатель степени) (в виде целого числа в коде со сдвигом)
3. мантисса (в нормализованной форме)

В качестве базы (основания степени) используется число [math] 2 [/math]. Экспонента хранится со сдвигом [math] -1023 [/math].

Знак
	Экспонента (11 бит)											Мантисса (52+1 бит)
0	0	0	0	0	0	0	0	0	0	0	0	1,	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
	62				52								51																																																0

Утверждение:

Итоговое значение числа вычисляется по формуле: [math] x = (-1)^{sign} \times (1.mant) \times 2^{exp} [/math].

Свойства чисел с плавающей точкой

1. В нормализованном виде любое отличное от нуля число представимо в единственном виде. Недостатком такой записи является тот факт, что невозможно представить число 0.
2. Так как старший бит двоичного числа, записанного в нормализованной форме, всегда равен 1, его можно опустить. Это используется в стандарте IEEE 754.
3. В отличие от целочисленных стандартов (например, integer), имеющих равномерное распределение на всем множестве значений, числа с плавающей точкой (double, например) имеют квазиравномерное распределение.
4. Вследствие свойства 3, числа с плавающей точкой имеют постоянную относительную погрешность (в отличие от целочисленных, которые имеют постоянную абсолютную погрешность).
5. Очевидно, не все действительные числа возможно представить в виде числа с плавающей точкой.
6. Точно в таком формате представимы только числа, являющиеся суммой некоторых обратных степеней двойки (не ниже -53). Остальные числа попадают в некоторый диапазон и округляются до ближайшей его границы. Таким образом, абсолютная погрешность составляет половину величины младшего бита.
7. В формате double представимы числа в диапазоне [math] [2.3 \times 10^{-308}, 1.7 \times 10^{308}] [/math].

Особые значения чисел с плавающей точкой

Ноль (со знаком)

В нормализованной форме невозможно представить ноль. Для его представления в стандарте зарезервированы специальные значения мантиссы и экспоненты.

Знак
	Экспонента					Мантисса
0/1	0	0	0	0	0	1,	0	0	0	0	0	0	0	0	0	0	= [math]\pm0[/math]

Согласно стандарту выполняются следующие свойства:

• [math] +0 = -0 [/math]
• [math] \frac{-0}{ \left| x \right| } = -0\,\![/math] (если [math]x\ne0[/math])
• [math] (-0) \cdot (-0) = +0\,\![/math]
• [math] \left| x \right| \cdot (-0) = -0\,\![/math]
• [math] x + (\pm 0) = x\,\![/math]
• [math] (-0) + (-0) = -0\,\![/math]
• [math] (+0) + (+0) = +0\,\![/math]
• [math] \frac{-0}{-\infty} = +0\,\![/math]
• [math] \frac{\left|x\right|}{-0} = -\infty\,\![/math] (если [math]x\ne0[/math])

Бесконечность (со знаком)

Для приближения ответа к правильному при переполнении, в double можно записать бесконечное значение. Так же, как и в случае с нолем, для этого используются специальные значение мантиссы и экспоненты.

Знак
	Экспонента					Мантисса
0/1	1	1	1	1	1	1,	0	0	0	0	0	0	0	0	0	0	= [math]\pm\infty[/math]

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

Неопределенность

В математике встречается понятие неопределенности. В стандарте double предусмотрено псевдочисло, которое арифметическая операция может вернуть даже в случае ошибки.

Знак
	Экспонента					Мантисса
0/1	1	1	1	1	1	1,	0/1	0/1	0/1	0/1	0/1	0/1	0/1	0/1	0/1	0/1	= [math]NaN[/math]

Неопределенность можно получить в нескольких случаях. Приведем некоторые из них:

• [math] f(NaN) = NaN [/math], где [math] f [/math] - любая арифметическая операция
• [math] \infty + (-\infty) = NaN [/math]
• [math] 0 \times \infty = NaN [/math]
• [math] \frac{\pm0}{\pm0} = \frac{\pm \infty}{\pm \infty} = NaN [/math]
• [math] \sqrt{x} = NaN [/math], где [math] x \lt 0 [/math]

Денормализованные числа

Денормализованные (denormalized numbers) - способ увеличить количество представимых числе в окрестности нуля. Каждое такое число по модулю меньше самого маленького нормализованного.< Согласно стандарту числа с плавающей точкой можно представить в следующем виде:

• [math] (-1)^s \times 1.M \times 2^E [/math], в нормализованном виде если [math] E_{min} \leq E \leq E_{max} [/math],
• [math] (-1)^s \times 0.M \times 2^E_{min} [/math], в денормализованном виде если [math] E = E_{min} - 1 [/math],

где [math] E_{min} [/math] - минимальное значение порядка, используемое для записи чисел (единичный сдвиг), [math] E_{min} - 1 [/math] - минимальное значение порядка, которое он может принимать - все биты нули, нулевой сдвиг.

Ввиду сложности, денормализованные числа обычно реализуют на программном уровне, а не на аппаратном. Из-за этого резко возрастает время работы с ними. Это недопустимо в областях, где требуется большая скорость вычислений (например, видеокарты). Так как денормализованные числа представляют числа мало отличные от нуля и мало влияют на результат, зачастую они игнорируются (что резко повышает скорость). При этом используются две концепции:

1. Flush To Zero (FTZ) - в качестве результата возвращается нуль, как только становится понятно, что результат будет представляться в денормализованном виде.
2. Denormals Are Zero (DAZ) - денормализованные числа, поступающие на вход, рассматриваются как нули.

Начиная с версии стандарта IEEE 754 2008 года денормализованные числа называются "субнормальными" (subnormal numbers), то есть числа, меньшие "нормальных".

Машинная эпсилон

Определение:

Машинная эпсилон - наибольшее положительное число [math] \varepsilon_m [/math], такое что, [math] 1 \oplus \varepsilon_m = 1 [/math], где [math] \oplus [/math] - машинное сложение.

Утверждение:

Таким образом, компьютер не различает числа [math] x [/math] и [math] y [/math], если [math] 1 \lt \frac{x}{y} \lt 1 + \varepsilon_m [/math].

Утверждение:

Из свойств чисел двойной точности следует, что для них [math] \varepsilon_m = 2^{-54}[/math].

Unit in the last place (Unit of least precision)

Мера единичной точности используется для оценки точности вычислений.

Определение:

Пусть [math] a [/math] - число с плавающей точкой, мантисса которого имеет длину [math] m [/math] бит, а экспонента - [math] e [/math] бит. Тогда [math] ulp(a) = 2^{e - m} [/math].

Приведем пример кода на Python, который показывает, при каком значении числа [math] x [/math] компьютер не различает числа [math] x [/math] и [math] x + 1 [/math].

>>> from math import *
>>> x = 1.0
>>> while (x != x + 1):
...   x *= 2
... 
>>> x
9007199254740992.0
>>> log(x) / log(2)
53.0

То есть [math] x = 2^{53} [/math], так как мантисса числа двойной точности содержит 53 бита (в памяти хранятся 52). В C++ для расчета расстояния между двумя числами двойной точности можно воспользоваться функцией [math] \mathrm{boost::math::float\_distance(a, b)} [/math].

Погрешность предиката "левый поворот"

Определения

Утверждение:

Пусть [math] D [/math] - множество всех чисел с плавающей точкой с операциями [math] \oplus, \ominus, \otimes. \forall a, b \in D: [/math] [math] a \oplus b = (a + b) (1 + \delta), |\delta| \leq \varepsilon_m [/math], [math] a \ominus b = (a - b) (1 + \delta), |\delta| \leq \varepsilon_m [/math], [math] a \otimes b = ab (1 + \delta), |\delta| \leq \varepsilon_m [/math].

Утверждение:

[math] \forall a, b, c \in D^2, \tilde{v} = (b_x \ominus a_x) \otimes (c_y \ominus a_y) \ominus (b_y \ominus a_y) \otimes (c_x \ominus a_x) [/math] [math] \exists \tilde{\epsilon} \in D: [/math] [math] \tilde{v} \gt \tilde{\epsilon} \Rightarrow (b - a) \times (c - a) \gt 0 [/math] [math] \tilde{v} \lt -\tilde{\epsilon} \Rightarrow (b - a) \times (c - a) \lt 0 [/math]

Расчет [math] \tilde{\epsilon} [/math]

Обозначим [math] v = (b - a) \times (c - a) = (b_x - a_x) (c_y - a_y) - (b_y - a_y) (c_x - a_x)[/math].

Теперь распишем это выражение в дабловой арифметике.

[math]\tilde{v} = (b_x \ominus a_x) \otimes (c_y \ominus a_y) \ominus (b_y \ominus a_y) \otimes (c_x \ominus a_x) = \\ = [ (b_x - a_x) (c_y - a_y) (1 + \delta_1) (1 + \delta_2) (1 + \delta_3) - \\ - (b_y - a_y) (c_x - a_x) (1 + \delta_4) (1 + \delta_5) (1 + \delta_6) ] (1 + \delta_7),[/math]

[math] |\delta_i| \leq \varepsilon_m [/math]

Заметим, что [math] v \approx \tilde{v} [/math]

Теперь оценим абсолютную погрешность [math] \epsilon = |v - \tilde{v}|. [/math]

[math] |v - \tilde{v}| = |(b_x - a_x) (c_y - a_y) - (b_y - a_y) (c_x - a_x) - \\ - (b_x - a_x) (c_y - a_y) (1 + \delta_1) (1 + \delta_2) (1 + \delta_3) (1 + \delta_7) + \\ + (b_y - a_y) (c_x - a_x) (1 + \delta_4) (1 + \delta_5) (1 + \delta_6) (1 + \delta_7)| = \\ = |(b_x - a_x) (c_y - a_y) (1 - (1 + \delta_1) (1 + \delta_2) (1 + \delta_3) (1 + \delta_7)) - \\ - (b_y - a_y) (c_x - a_x) (1 - (1 + \delta_4) (1 + \delta_5) (1 + \delta_6) (1 + \delta_7))| \leq \\ \leq |(b_x - a_x) (c_y - a_y) (1 - (1 + \delta_1) (1 + \delta_2) (1 + \delta_3) (1 + \delta_7))| + \\ + |(b_y - a_y) (c_x - a_x) (1 - (1 + \delta_4) (1 + \delta_5) (1 + \delta_6) (1 + \delta_7))| = \\ = |(b_x - a_x) (c_y - a_y)| \cdot |((1 + \delta_1) (1 + \delta_2) (1 + \delta_3) (1 + \delta_7) - 1)| + \\ + |(b_y - a_y) (c_x - a_x)| \cdot |((1 + \delta_4) (1 + \delta_5) (1 + \delta_6) (1 + \delta_7) - 1)| = \\ = |(b_x - a_x) (c_y - a_y)| \cdot |\delta_1 + \delta_2 + \delta_3 + \delta_7 + \delta_1 \delta_2 \ldots| + \\ + |(b_y - a_y) (c_x - a_x)| \cdot |\delta_4 + \delta_5 + \delta_6 + \delta_7 + \delta_4 \delta_5 \ldots| \leq \\ \leq |(b_x - a_x) (c_y - a_y)| \cdot (|\delta_1| + |\delta_2| + |\delta_3| + |\delta_7| + |\delta_1 \delta_2| \ldots) + \\ + |(b_y - a_y) (c_x - a_x)| \cdot (|\delta_4| + |\delta_5| + |\delta_6| + |\delta_7| + |\delta_4 \delta_5| \ldots) \leq \\ \leq |(b_x - a_x) (c_y - a_y)| \cdot (4 \varepsilon_m + 6 \varepsilon_m^2 + 4 \varepsilon_m^3 + \varepsilon_m^4) + \\ + |(b_y - a_y) (c_x - a_x)| \cdot (4 \varepsilon_m + 6 \varepsilon_m^2 + 4 \varepsilon_m^3 + \varepsilon_m^4) = \\ = (|(b_x - a_x) (c_y - a_y)| + |(b_y - a_y) (c_x - a_x)|)(4 \varepsilon_m + 6 \varepsilon_m^2 + 4 \varepsilon_m^3 + \varepsilon_m^4)[/math]

Пусть [math] t = (|(b_x - a_x) (c_y - a_y)| + |(b_y - a_y) (c_x - a_x)|).[/math] Получаем, что

[math] \epsilon = |v - \tilde{v}| \leq t \cdot (4 \varepsilon_m + 6 \varepsilon_m^2 + 4 \varepsilon_m^3 + \varepsilon_m^4). [/math]

[math]\tilde {t} = (|(b_x - a_x) (c_y - a_y) (1 + \delta_1) (1 + \delta_2) (1 + \delta_3)| + \\ + |(b_y - a_y) (c_x - a_x) (1 + \delta_4) (1 + \delta_5) (1 + \delta_6)|) (1 + \delta_7) \geq \\ \geq |(b_x - a_x) (c_y - a_y) (1 - \varepsilon_m)^3)|(1 - \varepsilon_m) + \\ + |(b_y - a_y) (c_x - a_x) (1 - \varepsilon_m)^3)|(1 - \varepsilon_m) = \\ = |(b_x - a_x) (c_y - a_y)| (1 - \varepsilon_m)^4 + |(b_y - a_y) (c_x - a_x)| (1 - \varepsilon_m)^4 = \\ = (|(b_x - a_x) (c_y - a_y)| + |(b_y - a_y) (c_x - a_x)|) (1 - \varepsilon_m)^4 = t \cdot (1 - \varepsilon_m)^4[/math]

Итого:

[math] t \leq \tilde{t} \frac{1}{(1 - \varepsilon_m)^4} = \tilde{t} (1 + 4 \varepsilon_m + 10 \varepsilon_m^2 + 20 \varepsilon_m^3 + \cdots) [/math]

[math] \epsilon = |v - \tilde{v}| \leq \tilde{\epsilon} \leq \tilde{t} (1 + 4 \varepsilon_m + 10 \varepsilon_m^2 + 20 \varepsilon_m^3 + \cdots) (4 \varepsilon_m + 6 \varepsilon_m^2 + 4 \varepsilon_m^3 + \varepsilon_m^4) [/math]

Ответ

[math] \tilde{\epsilon} \lt 8 \varepsilon_m \tilde{t} [/math]

Заметим, что это довольно грубая оценка. Вполне можно было бы написать [math] \tilde{\epsilon} \lt 4.25 \varepsilon_m \tilde{t} [/math] или [math] \tilde{\epsilon} \lt 4.5 \varepsilon_m \tilde{t}.[/math]

Зачем Антон сделал так? Да просто чтобы не возиться с суммой ряда.

Ссылки

en.wikipedia.org Floating point
en.wikipedia.org Half-precision floating point format
en.wikipedia.org Single precision floating point format
en.wikipedia.org Double precision floating point format
en.wikipedia.org Extended precision floating point format
en.wikipedia.org Quadruple precision floating point format
Goldberg, D. 1991 What every computer scientist should know about floating-point arithmetic
ieee.org IEEE 754
en.wikipedia.org Unit in the last place neerc.ifmo.ru/mediawiki Предикат "левый поворот"