Нормализованная запись числа

Любое действительное число можно записать в экспоненциальной форме

Экспоненциальная форма

Представление числа в виде \(a=\pm m \times P^q\), где \(P\) – основание системы счисления, \(m\) называется мантиссой числа, а \(q\) – порядком или экспонентой.

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

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

Представление числа в экспоненциальной форме, в которой мантисса числа имеет одну цифру перед запятой, т.е. \(1 \le m < P\).

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

Числа с плавающей запятой хранятся в памяти компьютера в нормализованном виде.

Погрешности

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

Существует два типа погрешностей, которые нас интересуют

Абсолютная погрешность

это модуль разности \(|a - a'|\) между числом \(a\) и его представлением \(a'\)

Относительная погрешность

это отношение абсолютной погрешности к величине \(\left| \frac{a-a'}{a} \right|\)

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

Формат с фиксированной запятой

Почти ничем не отличается от целых цифр. Гарантирует конечную абсолютную погрешность.

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

Формат с плавающей запятой

Гарантирует конечную относительную погрешность.

Существует международный стандарт для представления действительных чисел в двоичном коде, IEEE 754. 754 определяет следующие правила округления:

1. Если следующий за последним (младшим) значащим битом мантиссы разряд = 0, то ничего не делать.
2. Если следующий за последним (младшим) значащим битом мантиссы разряд = 1, но все последующие равны нулю, то:
   1. Если последний значащий бит мантиссы = 1, прибавить к нему 1.
   2. Если последний значащий бит мантиссы = 1, ничего не делать
3. Если следующий за последним (младшим) значащим битом мантиссы разряд = 1, и есть еще хотя бы один последующий разряд = 1, то прибавить к последнему значащему биту мантиссы 1.

В памяти хранится в виде знака, мантиссы и порядка.

Поскольку для хранения используется двоичный код, в нормализованной форме мантисса \(1\le m < 2\), следовательно, цифра перед запятой всегда должна быть \(1\). Поскольку она фиксирована, нет смысла ее хранить, и по соглашению она подразумевается. В памяти хранится, таким образом, только дробная часть мантиссы.

Порядок хранится как целое число со сдвигом \(b = 2^{n-1}-1\), где \(n\) – число бит порядка. То есть, порядок вычисляется как \(q=q' - b\), где \(q'\) – представление числа в памяти (в виде беззнакового целого).

Максимально и минимально представимые значения q получаются естественным образом:

\[q'_ \mathrm{min} = 0\] \[q'_ \mathrm{max} = 2^{n}-1\] \[q^\mathrm{theor}_ \mathrm{min} = 0 - b = 1-2^{n-1}\] \[q^\mathrm{theor}_ \mathrm{max} = 2^{n}-1 -b = 2^n - 2^{n-1} = 2^{n-1} (2-1) = 2^{n-1}\]

На практике, значения \(q=q^\mathrm{theor}_ \mathrm{min}\) и \(q=q^\mathrm{theor}_ \mathrm{max}\) отводятся под особые нужды, о которых ниже. Поэтому фактические пределы на 1 меньше с обеих сторон:

\[q_ \mathrm{min} = 2-2^{n-1}\] \[q_ \mathrm{max} = 2^{n-1}-1\]

Таким образом, максимально представимое оказывается \(2^{n-1}-1\), а минимально представимое – \(2-2^{n-1}\), где \(n\) – число разрядов, отводимое под порядок.

Что касается мантиссы, в нормализованном представлении

\[m_ \mathrm{min} = 1,\] \[m_ \mathrm{max} = 1+\frac{2^r-1}{2^r} = 2 - 2^{-r},\]

где \(r\) – количество бит, отводимых под дробную часть мантиссы.

С учетом денормализованного представления (см. ниже),

\[m_ \mathrm{min}^\mathrm{denorm} = 2^{-r}\]

Минимально и максимально представимые числа можно получить по формулам

\[ a_ \mathrm{max/min} = \pm m_ \mathrm{max} \cdot 2^{q_ \mathrm{max}} = \pm 2^{2^{n-1}} (1-\frac{1}{2^{r+1}}) \]

Минимально представимое положительное число

\[ a_ \mathrm{min}^+ = m_ \mathrm{min}^\mathrm{denorm} \cdot 2^{q_ \mathrm{min}} = 2^{2-2^{n-1}-r}\]

Типичные форматы:

Хранение в памяти

В общем случае, формат с плавающей запятой хранится следующим образом:

Где \(s\) – знак числа (\(0\) для \(+\) и \(1\) для \(-\)), \(q_{n-1} \ldots q_0\) – порядок, \(m_{r-1} \ldots m_0\) – дробная часть мантиссы.

Число \(a\), хранимое таким образом, можно вычислить как \[a = (s) \left(1+\sum_{i=0}^{r-1} m_i 2^{i-r} \right) \cdot 2^{\left(\sum_{j=0}^{n-1} q_i 2^i\right) - 2^{n-1} + 1}\]

Денормализованное представление

Для представления очень малых величин, используется так называемое денормализованное или субнормальное представление.

Если все биты порядка установлены в ноль, то это денормализованное число. Для этих чисел, целая часть мантиссы (цифра перед запятой) считается равной нулю, а порядок считается минимальным (не нулевым!). Это позволяет представлять величины меньшие минимального порядка, включая ноль (при нулевой мантиссе).

Где \(s\) – знак числа (\(0\) для \(+\) и \(1\) для \(-\)), \(m_{r-1} \ldots m_0\) – дробная часть мантиссы.

Число \(a\), хранимое таким образом, можно вычислить как \[a = (s) \left(\sum_{i=0}^{r-1} m_i 2^{i-r} \right) \cdot 2^{2-2^{n-1}}\]

Следует заметить, что “ноль” может иметь знак: “+0” и “-0”. Связано это с тем, что очень малые числа ввиду погрешности округления имеют нулевую мантиссу – в таких случаях полезно по крайней мере сохранять информацию о знаке.

Другие специальные значения включают “бесконечность” и “не число”.

Бесконечность получается, если все биты порядка равны 1, а биты мантиссы – 0.

“Не число” (NaN) получается, если все биты порядка равны 1, а биты мантиссы не все равны 0.

Арифметические операции

Умножение/деление

1. Знаковые биты складываются
2. Порядки складываются/вычитаются
3. Мантиссы умножаются/делятся
4. Результат нормализуется

При сложении знаковых битов мы получаем интуитивно ожидаемый знак результата. Вспомнив, что 0 эквивалентен плюсу, а 1 – минусу,

• При умножении положительных чисел, 0+0 = 0
• При умножении положительного и отрицательного числа, 0+1 = 1
• При умножении отрицательных чисел 1+1 = 10, по правилам целочисленной арифметики в конечном числе разрядов, оставляем 0.

Пример:

Допустим, мы умножаем \(10_2\) и \(1010_2\). В нормализованном виде

\[1\cdot2^1 \times 1,010\cdot 2^3\]

Порядки складываются, а мантиссы умножаются:

\[(1 \times 1,010)\cdot 2^{3+1}\]

\[1,010\cdot 2^4\]

В двоичном коде:

Напрямую складывать сдвинутые представления порядков нельзя, поскольку тогда сдвиг войдет в сумму дважды. Можно, однако, преобразовать представление со сдвигом к дополнительным кодам. Для этого старший разряд инвертируется, и прибавляется 1:

Суммируя и преобразуя обратно (вычитая 1 и инвертируя старший бит):

Теперь перемножаем мантиссы, не забывая про неявный бит:

Результат:

Сложение/вычитание

1. Порядки приводятся к наибольшему
2. Мантиссы складываются/вычитаются
3. Результат нормализуется

Приведение порядков работает следующим образом:

Допустим, мы складываем \(10_2\) и \(1010_2\). В нормализованном виде

\[1\cdot2^1 + 1,010\cdot 2^3\]

Мы не можем сложить мантиссы напрямую, поскольку разряды имеют разный “вес”. Поскольку мы храним только дробную часть, приведем число с меньшим порядком \(1\cdot2^1\) к большему порядку. Для этого, перенесем запятую влево:

\[0,01\cdot2^3 + 1,010\cdot 2^3\] \[(0,01 + 1,010)\cdot 2^3\] \[(1,100)\cdot 2^3\]

Рассмотрим происходящее в двоичном коде. Воспользуемся 16-битным форматом:

При приведении порядков, мантисса числа с меньшим порядком смещается вправо на разность порядков:

Затем мантиссы складываются:

Особые операции

Количество и распределение чисел с плавающей запятой

Понятно, что существует конечное количество представимых в формате с плавающей запятой чисел. В отличие от целых или чисел с фиксированной запятой, распределение представимых чисел оказывается неравномерным.

Общее количество представимых чисел посчитать сравнительно несложно. Это количество возможных значений мантиссы \(2^r\), где \(r\) – битность мантиссы, умноженное на количество возможных значений порядка, за исключением “не числа” и “бесконечностей” \(2^n-1\), умноженное на количество возможных значений знака \(2\) и минус один, поскольку ноль имеет два в общем-то эквивалентных представления.

Итого,

\[ N = 2\cdot 2^r \cdot (2^n - 1) -1 \]

Например, для формата одинарной точности количество представимых чисел равно

\[ N_ \mathrm{float} = 2\cdot2^{23} \cdot (2^8-1) - 1 = 4 278 190 079 \]

Их распределение по числовой оси, однако, неравномерно:

\(2^r\) представимых чисел находятся в диапазоне \([0; 2^{q_ \mathrm{min}})\), \(2^r\) – в диапазоне \([2^{q_ \mathrm{min}}; 2^{q_ \mathrm{min}+1})\), и так далее. При увеличении порядка на 1, диапазон увеличивается вдвое: \(2^{q}-2^{q-1} = 0.5*2^q = 2^{q-1}\), однако количество представимых чисел в диапазоне не меняется. Следовательно, “шаг” между представимыми числами так же увеличивается вдвое, а “плотность”, т.е. количество чисел на диапазон, соответственно, уменьшается вдвое.

Оценка для плотности в области числа \(x\) \[\rho(x) = 2^{r - \log_2 x}\]

Зависимость плотности от величины числа для формата одинарной точности (линейный масштаб)

Зависимость плотности от величины числа для формата одинарной точности (логарифмический масштаб)

Машинный эпсилон

Существует два определения машинного эпсилон.

Машинный эпсилон

это такое минимальное число \(\varepsilon > 0\), что, при вычислениях с плавающей запятой, \(1+\varepsilon \neq 1\)

это такое максимальное число \(\varepsilon>0\), что, при вычислениях с плавающей запятой, \(1+\varepsilon = 1\)

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

Увы, у нас нет неограниченной точности, поэтому на практике \(\varepsilon\), соответствующие этим определениям, несколько отличаются. Мы будем полагаться на первое определение, поскольку оно закреплено в стандарте языка ISO C.

Теоретические значения машинного эпсилон достаточно легко получить – это такое число, при приведении которого к нулевому порядку получается минимально представимая мантисса.

Рассмотрим число с мантиссой \(m=1\) и порядком \(q = -r\). При приведении этого числа к нулевому порядку, необходимо сдвинуть биты мантиссы на \(r\) разрядов вправо. При этом получаем \(m=2^{-r}\). Это минимальная денормализованная мантисса. Рассмотрим теперь число, меньшее этого \(m = m_ \mathrm{max}\), \(q= -r-1\). Сдвигая на \(r+1\) разрядов вправо, получаем \(0\).

Таким образом, получаем, что \(\varepsilon = 2^{-r}\).

Приведем типичные значения \(\varepsilon\):

Получить значение \(\varepsilon\) для типа в C++ можно, используя уже знакомую структуру std::numeric_limits, определенную в заголовке <limits>, используя функцию epsilon():

#include <iostream>
#include <limits>

using namespace std;

int main() {
  cout << numeric_limits<float>::epsilon() << endl;
  cout << numeric_limits<double>::epsilon() << endl;
  cout << numeric_limits<long double>::epsilon() << endl;
}

Так же можно написать несложную функцию для вычисления машинного эпсилон:

template<typename T>
T epsilon() {
    T f = 1;
    while(1+f>1)
        f/=2;
    return f*2;
}

Машинный эпсилон интересен тем, что относительная погрешность округления ограничена им сверху:

\[ \left| \frac{a-a'}{a} \right| < \varepsilon \]

Или, если используется корректное округление

\[ \left| \frac{a-a'}{a} \right| < \frac{1}{2}\varepsilon \]

Правила округления чисел с плавающей запятой

IEEE 754 определяет следующие правила округления:

1. Если следующий за последним (младшим) значащим битом мантиссы разряд = 0, то ничего не делать.
2. Если следующий за последним (младшим) значащим битом мантиссы разряд = 1, но все последующие равны нулю, то:
   1. Если последний значащий бит мантиссы = 1, прибавить к нему 1.
   2. Иначе, ничего не делать
3. Если следующий за последним (младшим) значащим битом мантиссы разряд = 1, и есть еще хотя бы один последующий разряд = 1, то прибавить к последнему значащему биту мантиссы 1.

По сути, правило (2) означает, что “половина” округляется к ближайшему чётному.

Например, для 10-битной мантиссы:

\(1,0000000001(00000\ldots)\) округляется до \(1,0000000001\)

\(1,0000000001(10000\ldots)\) округляется до \(1,0000000010\)

\(1,0000000000(10000\ldots)\) округляется до \(1,0000000000\)

\(1,0000000010(10000\ldots)\) округляется до \(1,0000000010\)

\(1,0000000010(11000\ldots)\) округляется до \(1,0000000011\)

\(1,0000000010(10000\ldots 1\ldots)\) округляется до \(1,0000000011\)

\(1,0000000001(10000\ldots 1\ldots)\) округляется до \(1,0000000010\)

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

Погрешность округления

Речь идет о погрешностях, связанных с конечностью представимых чисел. Так, например, число \(0.0001\) непредставимо в виде конечной двоичной дроби, и суммирование 10000 раз \(0.0001\) в итоге не даст ровно 1.

По расширению, так же связано с машинным эпсилон: если сложить очень большое и очень маленькое числа, последнее при приведении порядков обратится в ноль, и в результат не изменится.

Вычитание близких значений

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

Переполнение

Происходит при попытке представить слишком большое или слишком малое число. В первом случае получается ∞, а во втором – 0.

По сути, нужно всегда помнить, что любое представление \(a'\) имеет бесконечно много соответствующих ему действительных чисел, и по сути представляет некий интервал \([a'(1-\varepsilon/2); a'(1+\varepsilon/2)]\). В результате вычислений, этот интервал имеет тенденцию расширятся.

Элементы теории погрешностей

Пусть \(y(x)\) – функция одного аргумента, и аргумент \(x\) задан приближенным значением \(x'\) с абсолютной погрешностью \(\Delta x\). Рассматривая интервал значений \([y-\Delta y/2;y+\Delta y/2]\), получаемый в результате вычисления с приближенным значением \(x'\), и заменяя приращение функции \(\frac{\Delta y}{\Delta x}\) на ее производную в точке (теорема Лагранжа), получаем \[\Delta y = \Delta x \left|\left.\frac{\mathrm{d}f}{\mathrm{d}x}\right|_ {x=x'}\right|\]

Это легко обобщается на случай функции многих переменных:

\[ \Delta y = \sum_i \Delta x_i \left|\left.\frac{\partial f}{\partial x_i}\right|_ {x_1=x_1',\,x_2=x_2',\,\ldots}\right|\]

Так, например, для суммы или разности \(f(x_1,x_2) = x_1+x_2\), абсолютная погрешность будет равна \[\Delta (x_1+x_2) = \Delta x_1 + \Delta x_2 \]

Относительная погрешность \[\delta (x_1+x_2) = \left|\frac{\Delta x_1 + \Delta x_2}{x_1+x_2}\right| = \left|\frac{|x_1| \delta x_1 + |x_2| \delta x_2}{x_1+x_2}\right|\]

Для суммы, если относительные погрешности слагаемых равны \(\delta x_1 = \delta x_2 = \varepsilon\), то

\[\delta (x_1+x_2) = \varepsilon\]

Вывод 1

при суммировании, относительная погрешность сохраняется

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

Погрешность умножения \(f(x_1,x_2) = x_1 x_2\), \[\Delta (x_1 x_2) = |x_2| \Delta x_1 + |x_1| \Delta x_2 \]

\[\delta (x_1 x_2) = \left|\frac{|x_1| |x_2| \delta x_1 + |x_2| |x_1| \delta x_2}{x_1x_2}\right| = \delta x_1 + \delta x_2\]

Вывод 2

При умножении, относительные погрешности складываются. При одинаковых относительных погрешностях \(\forall i :\:\delta x_i = \varepsilon\), каждая операция умножения увеличивает погрешность на \(\varepsilon\).

Погрешность деления \(f(x_1,x_2) = x_1/ x_2\), \[\Delta (x_1 / x_2) = \left|\frac{1}{x_2}\right| \Delta x_1 + \left|\frac{x_1}{x_2^2}\right| \Delta x_2 \]

\[\delta (x_1 / x_2) = \left|\frac{\left|\frac{1}{x_2}\right| |x_1| \delta x_1 + \left|\frac{x_1}{x_2^2}\right| |x_2| \delta x_2}{x_1/x_2}\right| = \delta x_1 + \delta x_2\]

Вывод 3

При делении, относительные погрешности складываются. При одинаковых относительных погрешностях \(\forall i :\:\delta x_i = \varepsilon\), каждая операция деления увеличивает погрешность на \(\varepsilon\).