Такой метод является компромиссом между точностью и диапазоном представляемых значений. Представление чисел с плавающей точкой рассмотрим на примере чисел двойной точности (double precision). Такие числа занимают в памяти два машинных слова (8 байт на 32-битных системах). Наиболее распространенное представление описано в стандарте IEEE 754.
Кроме чисел двойной точности также используются следующие форматы чисел:
При выборе формата программисты идут на разумный компромисс между точностью вычислений и размером числа.
Нормальная и нормализованная формы [ править ]
Числа двойной точности [ править ]
Число с плавающей точкой хранится в нормализованной форме и состоит из трех частей (в скобках указано количество бит, отводимых на каждую секцию в формате double):
Знак
Экспонента (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
Утверждение:
Свойства чисел с плавающей точкой [ править ]
Особые значения чисел с плавающей точкой [ править ]
Ноль (со знаком) [ править ]
В нормализованной форме невозможно представить ноль. Для его представления в стандарте зарезервированы специальные значения мантиссы и экспоненты.
Знак
Экспонента
Мантисса
0 /1
0
0
0
0
0
1,
0
0
0
0
0
0
0
0
0
0
= [math]\pm0[/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]
Неопределенность можно получить в нескольких случаях. Приведем некоторые из них:
Денормализованные числа [ править ]
Ввиду сложности, денормализованные числа обычно реализуют на программном уровне, а не на аппаратном. Из-за этого резко возрастает время работы с ними. Это недопустимо в областях, где требуется большая скорость вычислений (например, видеокарты). Так как денормализованные числа представляют числа мало отличные от нуля и мало влияют на результат, зачастую они игнорируются (что резко повышает скорость). При этом используются две концепции:
Начиная с версии стандарта IEEE 754 2008 года денормализованные числа называются «субнормальными» (subnormal numbers), то есть числа, меньшие «нормальных».
Машинная эпсилон [ править ]
Unit in the last place (Unit of least precision) [ править ]
Мера единичной точности используется для оценки точности вычислений.
Заметим, что это довольно грубая оценка. Вполне можно было бы написать [math] \tilde <\epsilon>\lt 4.25 \varepsilon_m \tilde[/math] или [math] \tilde <\epsilon>\lt 4.5 \varepsilon_m \tilde.[/math]
Для работы с вещественными числами могут использоваться различные способы покрытия вещественной оси (конечным) набором представимых чисел. Как правило, они ограничиваются определенным подмножеством рациональных чисел:
Вместо слова запятая часто используется слово точка, так как в англоязычных странах для разделения целой и дробной части числа традиционно используется точка. Поэтому оба этих варианта можно считать синонимами.
Сравнительно с числами с фиксированной запятой числа с плавающей запятой, имеющие тот же размер представления, позволяют существенно расширить покрываемый диапазон вещественной оси, однако при этом теряется однородность покрытия: чем дальше от нуля, тем больше шаг между соседними представимыми числами.
Наиболее распространенным стандартом представления и правил работы с числами с плавающей запятой является стандарт IEEE-754 (в последней редакции известный также под обозначением ISO/IEC/IEEE 60559:2011). Его поддерживают большинство современных процессоров, умеющих работать с числами с плавающей запятой “непосредственно”.
Формат
IEEE-754 задаёт следующее двоичное представление числа.
Двоичное представление числа с плавающей запятой
Знак s
Экспонента E
Мантисса M
1 разряд
e разрядов
m разрядов
старший бит
биты (m + e – 1) … m
биты (m – 1) … 0
Экспонента exponent — обычно то же, что и порядок, однако в конкретных форматах порядок (как целое число) закодирован специальным образом (не равен численно беззнаковому представлению E), см. ниже.
В таблице ниже приведены сведения о четырёх основных двоичных (т.е. r = 2) форматах стандарта IEEE-754. Число ε — наименьшее положительное представимое число такое, что 1 + ε ≠ 1 в заданном формате. Смысл понятий “денорм.” (денормализованные числа) и “норм.” (нормализованные числа) дан ниже. В качестве числа d в таблице указано минимальное число десятичных знаков после запятой, которого достаточно для сохранения двоичного представления числа при преобразовании в десятичную форму и обратно. Необходимо помнить, что многие десятичные дроби являются бесконечными в двоичной системе: например, 0.2 не представимо точно в указанных форматах. Двоичные же дроби хотя и представимы точно конечным числом знаков, могут потребовать сотни десятичных цифр для точной записи.
Четыре стандартных формата
Формат
“Точность” precision
всего бит
e
m
d
ε
min денорм.
min норм.
max норм.
binary16
половинная half
16
5
10
5
9.77·10 –4
5.96·10 –8
6.1·10 –5
65504
binary32
одинарная single
32
8
23
9
1.19·10 –7
1.4·10 –45
1.18·10 –38
3.4·10 38
binary64
двойная double
64
11
52
17
2.22·10 –16
4.94·10 –324
2.2·10 –308
1.8·10 308
binary128
четверная quadruple
128
15
112
36
1.93·10 –34
6.48·10 –4966
3.4·10 –4932
1.2·10 4932
Порядок представлен в сдвинутом коде со сдвигом b = 2 e–1 – 1. Максимальное и минимальное значение E, взятого как беззнаковое число, интерпретируются особым образом.
Интерпретация представления числа с плавающей точкой, М — биты мантиссы (двоичная дробь)
Тип значения
E
M
Значение числа с плав. точкой
нуль
0
0
(–1) s ·0
денормализованное
0
≠ 0
(–1) s ·0.(M)·2 1–b
нормализованное
1 … 2 b
любое
(–1) s ·1.(M)·2 E–b
бесконечность
2 b + 1
0
(–1) s ·∞
нечисло (“NaN”)
2 b + 1
≠ 0
код ошибки
Запись 0.(M) или 1.(M) означает, что m бит мантиссы выписываются после точки в двоичной записи числа. Под нормализацией числа понимается его домножение на такую степень двойки, что слева от точки получается единица. Денормализованные числа требуют недоступных значений экспоненты для нормализации.
Замечания
Режимы округления
Стандарт IEEE-754 предусматривает наличие пяти режимов округления, которые можно проиллюстрировать с помощью следующей таблицы.
Режимы округления IEEE-754
Режим (округление к целому)
–2.7
–2.5
–2.2
–1.5
1.5
2.2
2.5
2.7
к ближайшему, половинки к чётному to nearest, ties to even
–3.0
–2.0
–2.0
–2.0
2.0
2.0
2.0
3.0
к ближайшему, половинки от нуля to nearest, away from zero
–3.0
–3.0
–2.0
–2.0
2.0
2.0
3.0
3.0
к нулю toward zero, truncation
–2.0
–2.0
–2.0
–1.0
1.0
2.0
2.0
2.0
к плюс-бесконечности to positive infinity, upwards, ceiling
–2.0
–2.0
–2.0
–1.0
2.0
3.0
3.0
3.0
к минус-бесконечности to negative infinity, downwards, floor
–3.0
–3.0
–3.0
–2.0
1.0
2.0
2.0
2.0
При округлении к ближайшему представимому выбор в пользу округления половинок именно к чётному (как правило, используется именно этот вариант, например, на x86) или нечётному вместо “школьного” округления “от нуля” был сделан с целью балансировки накапливающихся округлений и уменьшению случайной погрешности.
Режим округления применяется не только при приведении к целым, но и при округлении результата для “вписывания” в мантиссу, т.е. вычислении результирующего младшего бита мантиссы.
Числа с плавающей запятой в C и C++
Одна из типичных ошибок новичков — использование целочисленных констант там, где необходимы константы в плавающей точке.
В примере выше 1 и 3 — целые числа, поэтому к ним применяется целочисленное деление (с отбрасыванием остатка), т.е. 1/3 = 0. Поэтому функция bad_cubic_root не вычисляет кубический корень. Правильным вариантом будет следующий код (также с учётом знака, так как стандартная функция возведения в степень не принимает отрицательное основание и нецелый показатель)
В стандарте C++11, впрочем, определена функция для вычисления кубического корня (cbrt). При её наличии следует использовать её вместо варианта на основе pow.
Языки C и C++ предоставляют три встроенных типа для работы с числами с плавающей точкой:
Стандарт гарантирует лишь отношение float ⊆ double ⊆ long double.
Часть характеристик форматов с плавающей запятой, используемых конкретным компилятором, можно получить, подключив заголовочный файл cfloat, где определен ряд макросов. Ниже * соответствует FLT для типа float, DBL для double и LDBL для long double.
Стандартные средства управления режимом округления и реакцией на исключительные ситуации (ошибки, приводящие к появлению нечисел, переполнения и округления со значительной потерей точности) предоставлены в заголовочном файле cfenv.
Заголовочный файл cmath
Заголовочный файл cmath предоставляет ряд определений математических функций (вычисляемых приближённо в плавающей точке, тем не менее, для простоты в обозначениях используются действительные числа вплоть до всей числовой оси). Здесь приводятся сведения по определениям для C++, принимающим разнотипные значения (float, double и long double). Если предел той или иной функции при стремлении к нулю (с той или иной стороны) или бесконечностям определён (включая значения ±∞), то он берётся в качестве значения этой функции в соответствующих предельных точках (поэтому многие интервалы и полуинтервалы в колонке “область значений” могут считаться отрезками).
Тригонометрические и гиперболические функции
Функция
Область определения
Область значений
Замечание
sin(x)
ℝ
[–1, 1]
x в радианах
cos(x)
ℝ
[–1, 1]
x в радианах
tan(x)
x ∉ πℤ + π/2
ℝ
неограниченный рост погрешности около точек разрыва
asin(x)
[–1, 1]
[–π/2, π/2]
acos(x)
[–1, 1]
[0, π]
atan(x)
ℝ
[–π/2, π/2]
atan2(y, x)
ℝ 2
[–π, π]
угол против часовой стрелки между векторами (1, 0) и (x, y)
sinh(x)
ℝ
ℝ
cosh(x)
ℝ
[1, +∞)
tanh(x)
ℝ
(–1, 1)
asinh(x)
ℝ
ℝ
C++11
acosh(x)
[1, +∞)
ℝ
C++11
atanh(x)
(–1, 1)
ℝ
C++11
Степени и логарифмы
Функция
Область определения
Область значений
Замечание
exp(x)
ℝ
(0, +∞)
константа Эйлера e (основание натурального логарифма) в степени x
expm1(x)
ℝ
(–1, +∞)
C++11, формально то же, что exp(x) – 1, но сохраняет высокую точность при x ≈ 0
exp2(x)
ℝ
(0, +∞)
C++11, 2 x
log(x)
(0, +∞)
ℝ
натуральный логарифм
log10(x)
(0, +∞)
ℝ
десятичный логарифм
log1p(x)
(–1, +∞)
ℝ
C++11, формально то же, что log(x + 1), но сохраняет высокую точность при x ≈ 0
log2(x)
(0, +∞)
ℝ
C++11, двоичный логарифм
sqrt(x)
[0, +∞)
[0, +∞)
квадратный корень
cbrt(x)
ℝ
ℝ
C++11, кубический корень
hypot(x, y)
ℝ 2
[0, +∞)
C++11, sqrt(x 2 + y 2 ), без промежуточного переполнения
pow(x, i)
ℝ × ℤ
ℝ
возведение x в целочисленную степень i
pow(x, y)
(0, +∞) × ℝ
(0, +∞)
возведение x в степень с показателем y, заданным в плавающей точке
Округление
Функция
Замечание
ceil(x)
округление наверх: ⌈x⌉, результат в формате плавающей точки
floor(x)
округление вниз: ⌊x⌋, результат в формате плавающей точки
trunc(x)
C++11, округление к нулю, результат в формате плавающей точки
round(x)
C++11, округление к ближайшему целому (режим “от нуля” для половинок), результат в формате плавающей точки
lround(x)
C++11, округление к ближайшему целому, результат в формате long int
llround(x)
C++11, округление к ближайшему целому, результат в формате long long int
nearbyint(x)
C++11, округляет до целых, используя текущий режим округления (который можно выставить с помощью fesetround)
При использовании C или для явного указания того, что аргумент и результат функции “абсолютное значение” имеют тип double, следует использовать функцию fabs вместо abs. См. также о целочисленном варианте abs.
Управление режимом округления cfenv
Подключив заголовочный файл cfenv можно управлять текущим режимом округления. Возможен выбор из четырёх стандартных (по стандарту C) режимов:
Данные макросы раскрываются в целочисленные константы и определены, если платформа поддерживает соответствующие режимы. Платформа может также определять дополнительные режимы.
Точность операций
Абсолютная погрешностьx равна |x – x * |.
Относительная погрешностьx при условии x * ≠ 0 равна |x – x * |/|x * |.
При выполнении арифметических операций, вычислении функции fma (слитое умножение-сложение fused multiply–add ), а также вычислении квадратного корня результат округляется к ближайшему представимому в текущем формате числу, т.е. имеет точность 0.5 ULP (при использовании соответствующего режима округления). Одним из следствий такого округления является точный результат sqrt от квадрата целого числа (если он представим). Вычисление тригонометрических функций, степеней и логарифмов обычно выполняется не столь точно, давая погрешность в 1–5 ULP и более (зависит от реализации). В результате вычисление формул, содержащих такие функции, с одними и теми же значениями переменных может давать разные результаты на разных платформах.
Благодаря двоичному формату чисел IEEE-754, при интерпретации их как целых (на некоторых платформах может потребоваться обращение порядка байт), модуль их разности в случае совпадения знака и порядка равен расстоянию между ними в ULP.
Так как из-за округлений при вычислениях накапливается погрешность, вместо прямого сравнения чисел с плавающей точкой на равенство или неравенство, обычно разумнее выполнять оценку расстояния между ними. В простейшем случае может использоваться сравнение по абсолютной величине погрешности. Этот вариант лучше всего подходит, когда надо проверять числа на равенство нулю.
В качестве eps можно брать кратное DBL_EPSILON — оценивать разность в “эпсилонах”.
Наконец, можно оценивать погрешность непосредственно в ULP.
В отличие от предыдущих вариантов, ulp_equal вернет true для бесконечностей одного знака. В то же время расстояние в ULP от DBL_MAX до +∞ равно всего лишь единице. Нечисла также могут оказаться “почти равными”. Напротив, расстояние между –0 и +0 получается очень большим (числа разного знака таким способом сравнивать в общем-то бессмысленно, и определённая ниже функция возвращает ложь для неравных чисел разного знака).
Функция memcpy копирует кусок памяти.
Продолжая тему оперирования непосредственно двоичными представлениями чисел с плавающей запятой формата IEEE-754, можно обратить внимание на следующее. Предположим, что a ≥ +0 (поведение отрицательных чисел аналогично), тогда при увеличении a, интерпретируемого как целое, на единицу получаем:
Простым циклом можно перечислить все представимые неотрицательные числа с плавающей точкой формата binary32:
Стандартным средством получения следующего представимого числа является функция nextafter (см. в таблице, приведённой в разделе, описывающем стандартный заголовочный файл cmath):
Примеры
Параметр линейной интерполяции
Пусть некоторая функция времени задана в конечном числе узлов и требуется ее приблизить ломаной. Чтобы получить параметр линейной интерполяции на отрезке ломаной, требуется преобразовать текущее время. Это можно попробовать сделать следующим образом:
Lerp — распространённая аббревиатура английского linear interpolation — “линейная интерполяция”.
Однако округление при вычитании существенно различных по величине t1 и h может приводить к выходу конечного результата за пределы отрезка [0, 1]. Исправленный вариант может выглядеть так:
Вычисление значений многочленов
Метод Горнера
Воспользовавшись стандартной функцией fma можно попробовать увеличить точность функции poly_horner (сократив половину промежуточных округлений).
В случае, если целевой процессор имеет встроенные команды для выполнения совмещённого умножения-сложения, и компилятор поддерживает соответствующую оптимизацию, то данный код может оказаться не только точнее, но и быстрее первого варианта (скорость вычисления FMA может совпадать со скоростью простого умножения).
Метод Эстрина
Впрочем, можно попробовать ещё ускорить данную операцию, если у нас в распоряжении имеется процессор, способный исполнять одновременно более одной команды (в данном случае имеется в виду не “многоядерность”, а наличие таких аппаратных особенностей как глубокий конвейер и/или суперскалярность и/или SIMD-команды). Для этого воспользуемся ассоциативностью суммы (впрочем, в плавающей точке ассоциативность теряется, поэтому результат может не совпадать с результатом функции poly_horner_fma). Кроме того, чтобы игра стоила свеч, нужен многочлен высокой степени (с большим количеством членов).
Один шаг: разбить по квадратам x.
Второй шаг: разбить по квадратам квадратов x.
Обобщение метода Горнера — метод Эстрина: дерево суммирования по квадратам
Вычисление произведения последовательности без промежуточного переполнения или исчезновения
Иногда случается ситуация, когда требуется вычислить произведение элементов некоторой конечной последовательности чисел в плавающей точке, хранящихся в массиве, читаемых с потока или вычисляемых по формуле. Так или иначе, в процессе такого вычисления может возникнуть неприятная ситуация: нехватка диапазона экспоненты для представления величины промежуточного результата, в то время как правильно вычисленных конечный результат представим. Например, для простоты представим себе последовательность, состоящую из двух тысяч двоек и двух тысяч “половинок” (0.5). Произведение равно единице. Однако, даже при использовании формата binary64 ещё на первой половине произойдёт переполнение, так как 2 2000 слишком велико, и общий результат наивного последовательного перемножения получится равным бесконечности. Если перемножать с конца к началу, то произойдёт исчезновение, так как 2 –2000 слишком мало, и общий результат получится равным нулю.
Суммирование с компенсацией погрешности округлений
Алгоритм, называемый суммированием Кэхэна или компенсационным суммированием (У.Кэхэн), позволяет получить гораздо более точный результат, суммируя отброшенные при округлении основной суммы части в отдельной переменной — компенсации (в примере реализации — переменная comp ). Относительная погрешность этого алгоритма ограничена сверху значением (цитир. по: Higham N. The accuracy of floating point summation // SIAM J. Sci. Comput. 1993. Vol. 14, No. 4. P.783–799)
Погрешность компенсационного суммирования
Следуя этой формуле, можно сделать вывод, что если в суммируемой последовательности есть числа с разным знаком, то их надо суммировать отдельно (положительные и отрицательные), затем сложить полученные результаты.
Улучшенный вариант компенсационного суммирования Ноймайера имеет следующий вид:
Численное решение уравнений
Пусть дана непрерывная функция f : ℝ → ℝ и два числа a 2 – 2, a = 0, b = 2.
В указанном примере на каждом шаге, очевидно, истинно m ∈ ℚ, в то же время, решения уравнения f(x) = 0 не принадлежат ℚ. Таким образом, наш алгоритм за n шагов может гарантировать лишь сужение исходного интервала, содержащего корень, в 2 n раз, но не обнаружение точного значения корня уравнения. Если достаточно знать корень с погрешностью ε > 0, то достаточно выполнить n = ⌈log2((b – a)/ε)⌉ шагов.
Это типичная ситуация в вычислительной математике: даже при условии абсолютно точной арифметики рассматривается нахождение решений с некоторой ненулевой погрешностью, которую, впрочем, можно сделать сколь угодно малой.
Однако, при использовании физических компьютеров, абсолютно точная арифметика недоступна. При использовании чисел с плавающей запятой через конечное число шагов алгоритма станет истинно (a = m)∨(b = m) (почему?). Поэтому бисекцию можно записать следующим образом:
Вилочный метод секущих
Рассмотрим ещё другой алгоритм — вилочный метод секущих (другое название: метод ложного положения regula falsi method ), напоминающий метод бисекции в том плане, что по мере приближения корень всегда остаётся в некотором интервале. Всё отличие от предыдущего алгоритма будет состоять в способе вычисления очередной точки m. Вместо того, чтобы брать середину текущего отрезка, возьмём точку пересечения оси абсцисс секущей — прямой, проведённой через координаты (a, f(a)) и (b, f(b)). Делая так, мы “надеемся”, что на интервале [a, b] функция достаточно похожа на прямую. Это действительно так, если интервал достаточно мал, а функция f дифференцируема на нём.
Шаг 1: m = 30/13 ≈ 2.3077 Шаг 2: m = 1028/399 ≈ 2.57644 Шаг 3: m ≈ 2.61485 Шаг 4: m ≈ 2.61996
первые шаги метода секущих (значения с индексом b — результат бисекции для сравнения)
Шаг
m
b – m
f(m)
mb
f(mb)
1
2.30769
0.692308
–0.317251
2
–0.(5)
2
2.57644
0.423559
–0.0498588
2.5
–0.131944
3
2.61485
0.385152
–0.00673158
2.75
0.155382
4
2.61996
0.380036
–0.000889564
2.625
0.00488281
5
2.620639
0.379361
–0.000117219
2.5625
–0.0651991
6
2.620728
0.379272
–0.0000154404
2.59375
–0.0305803
Легко заметить, что, с одной стороны, метод секущих уже не гарантирует нам сколь угодно малый интервал — в нашем примере бесконечное число шагов даст a = m = корень ≈ 2.62074, в то время как b так и останется в точке 3. С другой стороны, m намного быстрее подходит к корню, чем в методе бисекции.
Реализации этих двух методов, а также некоторых других методов, см. в примере.
Поиск минимума функции
Функция унимодальна на отрезке [a, b], если существует m ∈ [a, b] такая, что функция строго убывает на [a, m] и строго возрастает на [m, b]. Таким образом, функция достигает единственного на [a, b] минимума в точке m.
Рассмотрим задачу поиска минимума унимодальной на заданном отрезке [a, b] функции f : ℝ → ℝ.
Троичный поиск
Троичный поиск ternary search — простейший метод, решающий с заданной точностью ε поставленную задачу. Принцип его работы можно описать следующим алгоритмом:
Легко заметить, что приближение минимума (для данного примера min = –0.230445) происходит немонотонно.
Метод золотого сечения
Существует метод, основанный на том же принципе, что и троичный поиск, но более экономный — позволяющий вычислять одно новое значение функции на каждом шаге, а не два, как троичный поиск. Это метод носит название метод золотого сечения. Сократить вычисления можно в том случае, если одна из предыдущих позиций, выбранных в [a, b] будет сохраняться на каждом шаге. Этого можно добиться, если делить отрезок не на трети, а в отношении 1:φ, где φ — “золотое сечение”. Для любых x > y > 0 таких, что φ = x/y, выполняется также (x + y)/x = x/y = φ. Отсюда легко получить тождество φ = 1 + 1/φ, позволяющее вычислить φ (положительный корень квадратного уравнения).
Теперь передвинем точку b в точку d и пересчитаем новые c‘и d’:
Обобщение метода Горнера — метод Эстрина: дерево суммирования по квадратам
Погрешность компенсационного суммирования
Шаг 1: m = 30/13 ≈ 2.3077 
Шаг 2: m = 1028/399 ≈ 2.57644 
Шаг 3: m ≈ 2.61485 
Шаг 4: m ≈ 2.61996