Числа с плавающей запятой

Кувшинов Д.Р.

2015


Общее оглавление


IEEE-754

Введение

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

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

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

Наиболее распространенным стандартом представления и правил работы с числами с плавающей запятой является стандарт IEEE-754 (в последней редакции известный также под обозначением ISO/IEC/IEEE 60559:2011). Его поддерживают большинство современных процессоров, умеющих работать с числами с плавающей запятой “непосредственно”.

Формат

IEEE-754 задаёт следующее двоичное представление числа.

Двоичное представление числа с плавающей запятой
Знак s Экспонента E Мантисса M
1 разряд e разрядов m разрядов
старший бит биты (m + e – 1) … m биты (m – 1) … 0


Экспонента exponent — обычно то же, что и порядок, однако в конкретных форматах порядок (как целое число) закодирован специальным образом (не равен численно беззнаковому представлению E), см. ниже.

Мантисса mantissa — специальным образом закодированное значение множителя степени основания, называемое также показателем significand или коэффициентом coefficient.

В таблице ниже приведены сведения о четырёх основных двоичных (т.е. 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·1038
binary64 двойная double 64 11 52 17 2.22·10–16 4.94·10–324 2.2·10–308 1.8·10308
binary128 четверная quadruple 128 15 112 36 1.93·10–34 6.48·10–4966 3.4·10–4932 1.2·104932

Порядок представлен в сдвинутом коде со сдвигом b = 2e–1 – 1. Максимальное и минимальное значение E, взятого как беззнаковое число, интерпретируются особым образом.

Интерпретация представления числа с плавающей точкой, М — биты мантиссы (двоичная дробь)
Тип значения E M Значение числа с плав. точкой
нуль 0 0 (–1)s·0
денормализованное 0 ≠ 0 (–1)s·0.(M)·21–b
нормализованное 1 … 2 b любое (–1)s·1.(M)·2Eb
бесконечность 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++

Литералы констант с плавающей точкой обязаны содержать точку (целая или дробная часть числа при этом могут быть опущены, если они нулевые) или символ, отделяющий порядок от коэффициента (e, E или p, P), в противном случае они трактуются как целые числа. Например, 1010 как число в плавающей точке будет 1e10, 10’000’000’000. или 10.0e+9 — все эти представления равнозначны.

Для использования точного двоичного представления константы следует воспользоваться префиксом 0x, 0X (C99, C++11). В этом случае коэффициент записывается в виде шестнадцатеричной дроби и позволяет задать мантиссу итогого представления с точностью до бита. Показатель записывается после буквы p или P в десятичной форме, но задаёт степень r (как правило, двойки), например 0xA.BCDp+4 равно (10 + 11·16–1 + 12·16–2 + 13·16–3)·24 = 171.80078125.

Одна из типичных ошибок новичков — использование целочисленных констант там, где необходимы константы в плавающей точке.

double bad_cubic_root(double x)
{
  return pow(x, 1/3); // возводит в степень 1/3?
}

В примере выше 1 и 3 — целые числа, поэтому к ним применяется целочисленное деление (с отбрасыванием остатка), т.е. 1/3 = 0. Поэтому функция bad_cubic_root не вычисляет кубический корень. Правильным вариантом будет следующий код (также с учётом знака, так как стандартная функция возведения в степень не принимает отрицательное основание и нецелый показатель)

double cubic_root(double x)
{
  return x < 0.0?
      -pow(-x, 1./3)
    :  pow( x, 1./3);
}

В стандарте 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, 2x
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(x2 + y2), без промежуточного переполнения
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)


Утилиты
Функция Область определения Область значений Замечание
abs(x), fabs(x) [0, +∞) модуль x
fmod(x, y) ℝ × (ℝ \ {0}) (–|y|, |y|) обобщение остатка от деления: xy·trunc(x/y), имеет знак x
fmax(x, y) 2 2 C++11, максимум, если один аргумент — нечисло, возвращает другой аргумент
fmin(x, y) 2 2 C++11, минимум, если один аргумент — нечисло, возвращает другой аргумент
fdim(x, y) 2 [0, +∞) C++11, максимум из нуля и (xy)
fma(a, b, c) 3 C++11, вычисляет a·b + c без промежуточных округлений
ldexp(x, n) ℝ × ℤ вычисляет x·2n, n — целое число (т.е. добавляет n к показателю)
frexp(x, p) ±[0.5, 1) × ℤ сохраняет целочисленный показатель по указателю p, возвращает норм. мантиссу
modf(x, p) (-1, 1) × ℤ “разбирает” число на дробную часть (возвращает) и целую часть (через указатель p)
ilogb(x) int C++11, извлекает значение показателя
nextafter(x, y) 2 C++11, возвращает следующее после x в направлении y представимое число
copysign(x, y) 2 C++11, возвращает |x|·sgn y (подменяет знак x на знак y)
signbit(x) bool C++11, извлекает бит знака
nanf(s), nan(s) строки нечисла C++11, создаёт нечисло, отвечающее ошибке s (определяется реализацией библиотеки)


При использовании C или для явного указания того, что аргумент и результат функции “абсолютное значение” имеют тип double, следует использовать функцию fabs вместо abs. См. также о целочисленном варианте abs.

Классификация (все из C++11)
Функция Замечание
fpclassify(x) возвращает описание x: одну из предопределённых констант FP_INFINITE, FP_NAN, FP_NORMAL, FP_SUBNORMAL, FP_ZERO
isfinite(x) возвращает fpclassify(x) ∈ { FP_NORMAL, FP_SUBNORMAL, FP_ZERO }
isinf(x) возвращает fpclassify(x) = FP_INFINITE
isnan(x) возвращает fpclassify(x) = FP_NAN
isnormal(x) возвращает fpclassify(x) = FP_NORMAL
isunordered(x, y) возвращает истину, если x или y — нечисло, ложь в противном случае
islessgreater(x, y) возвращает истину, если x и y не являются нечислами и не равны
isgreater(x, y) возвращает истину, если x > y, иначе ложь, нечисла не вызывают исключительную ситуацию
isgreaterequal(x, y) аналогично isgreater для x >= y
isless(x, y) аналогично isgreater для x < y
islessequal(x, y) аналогично isgreater для x <= y


Некоторые константы, определённые в cmath.
Константа Смысл
INFINITY бесконечность, тип float
NAN “тихое” нечисло, тип float
HUGE_VALF число типа float, вызывающее переполнение (в случае IEEE-754 — бесконечность)
HUGE_VAL то же, что HUGE_VALF, но типа double


Управление режимом округления cfenv

Подключив заголовочный файл cfenv можно управлять текущим режимом округления. Возможен выбор из четырёх стандартных (по стандарту C) режимов:

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


Точность операций

Пусть x — приближённое значение некоторой величины x*.

Абсолютная погрешность x равна |xx*|.

Относительная погрешность x при условии x* ≠ 0 равна |xx*|/|x*|.

Точность выполнения операций часто измеряется в единицах, равных весу самого младшего разряда мантиссы результата, называемых ULP units in the last place. Например, расстояние между 1.0 и (1.0 + ε) составляет как раз 1 ULP. При удвоении числа ULP также удваивается (это относительная, а не абсолютная единица).

При выполнении арифметических операций, вычислении функции fma (слитое умножение-сложение fused multiply–add), а также вычислении квадратного корня результат округляется к ближайшему представимому в текущем формате числу, т.е. имеет точность 0.5 ULP (при использовании соответствующего режима округления). Одним из следствий такого округления является точный результат sqrt от квадрата целого числа (если он представим). Вычисление тригонометрических функций, степеней и логарифмов обычно выполняется не столь точно, давая погрешность в 1–5 ULP и более (зависит от реализации). В результате вычисление формул, содержащих такие функции, с одними и теми же значениями переменных может давать разные результаты на разных платформах.

Благодаря двоичному формату чисел IEEE-754, при интерпретации их как целых (на некоторых платформах может потребоваться обращение порядка байт), модуль их разности в случае совпадения знака и порядка равен расстоянию между ними в ULP.

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

#include <cmath> // далее подразумевается
inline bool abs_equal(double a, double b, double eps)
{
  return abs(a - b) <= eps;
}

Однако для сравнения произвольных чисел подобрать фиксированное значение константы eps часто не представляется возможных. В этом случае следует сравнивать числа по относительной величине погрешности. Здесь переменные a и b в некотором смысле равноценны, поэтому берём не модуль a или b, а максимум из их модулей.

inline bool rel_equal(double a, double b, double eps)
{
  return abs(a - b) <= eps * fmax(abs(a), abs(b)); // или std::max
}

В качестве eps можно брать кратное DBL_EPSILON — оценивать разность в “эпсилонах”.

inline bool rel_epsilon_equal(double a, double b, double epsilons)
{
  return abs(a - b) <= (DBL_EPSILON * epsilons) * fmax(abs(a), abs(b)); // или std::max
}

Наконец, можно оценивать погрешность непосредственно в ULP.

В отличие от предыдущих вариантов, ulp_equal вернет true для бесконечностей одного знака. В то же время расстояние в ULP от DBL_MAX до +∞ равно всего лишь единице. Нечисла также могут оказаться “почти равными”. Напротив, расстояние между –0 и +0 получается очень большим (числа разного знака таким способом сравнивать в общем-то бессмысленно, и определённая ниже функция возвращает ложь для неравных чисел разного знака).

#include <cstdint> // int64_t
#include <cstdlib> // abs (в cmath abs для плавающей точки)
#include <cstring> // memcpy

bool ulp_equal(double a, double b, int ulps)
{
  using namespace std;

  // Предполагаем формат binary64 для double, проверим хотя бы размер во время компиляции.
  static_assert(sizeof(int64_t) == sizeof(double), "double is not 64 bits on this platform");

  if (a == b)
    return true;

  int64_t abits, bbits; // двоичное представление в дополнительном коде
  memcpy(&abits, &a, sizeof(a));
  memcpy(&bbits, &b, sizeof(b));

  // разные знаки? -- всегда не равны
  if ((abits ^ bbits) < 0)
    return false;

  // разность в ULP
  return abs(abits - bbits) <= ulps;
}

Функция memcpy копирует кусок памяти.

Продолжая тему оперирования непосредственно двоичными представлениями чисел с плавающей запятой формата IEEE-754, можно обратить внимание на следующее. Предположим, что a ≥ +0 (поведение отрицательных чисел аналогично), тогда при увеличении a, интерпретируемого как целое, на единицу получаем:

Простым циклом можно перечислить все представимые неотрицательные числа с плавающей точкой формата binary32:

// на примере float, сколько их всего? а в случае double?
void print_all_nonnegative_floats()
{
  using namespace std;
  float f;
  for (uint32_t i = 0; i <= 0x7F800000u; i++)
  {
    memcpy(&f, &i, sizeof(float));
    cout << f << '\n';
  }
}

Стандартным средством получения следующего представимого числа является функция nextafter (см. в таблице, приведённой в разделе, описывающем стандартный заголовочный файл cmath):

// на примере float
void print_all_nonnegative_floats()
{
  using namespace std;
  for (float f = 0.f; f <= FLT_MAX; f = nextafter(f, INFINITY))
    cout << f << '\n';
}


Примеры

Параметр линейной интерполяции

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

// t -- время, t1 -- следующий узел, h -- шаг
inline float lerp_param_bad(float t, float t1, float h)
{
  const float t0 = t1 - h; // предыдущий узел
  return (t - t0) / h;
}

Lerp — распространённая аббревиатура английского linear interpolation — “линейная интерполяция”.

Однако округление при вычитании существенно различных по величине t1 и h может приводить к выходу конечного результата за пределы отрезка [0, 1]. Исправленный вариант может выглядеть так:

inline float lerp_param(float t, float t1, float h)
{
  const float dt = t1 - t; // сколько осталось
  return (h - dt) / h;
}  


Вычисление значений многочленов

Метод Горнера

Пусть дан многочлен вида Σai xi, где ai — коэффициенты, а x — свободная переменная. Широко известен метод эффективного вычисления значения многочлена для заданного x, называемый методом Горнера. Метод заключается в циклическом выполнении добавления очередного коэффициента и домножении суммы на x.

// Передаём коэффициенты многочлена указателем на массив.
double poly_horner(double x, double a[], size_t n)
{
  assert(n != 0);
  auto s = a[--n];
  while (n != 0)
    s = s * x + a[--n];
  return s;
}

Воспользовавшись стандартной функцией fma можно попробовать увеличить точность функции poly_horner (сократив половину промежуточных округлений).

// fma(a, b, c) --> (a * b + c) с одним округлением (C99, C++11)
double poly_horner_fma(double x, double a[], size_t n)
{
  assert(n != 0);
  auto s = a[--n];
  while (n != 0)
    s = fma(s, x, a[--n]);
  return s;
}

В случае, если целевой процессор имеет встроенные команды для выполнения совмещённого умножения-сложения, и компилятор поддерживает соответствующую оптимизацию, то данный код может оказаться не только точнее, но и быстрее первого варианта (скорость вычисления FMA может совпадать со скоростью простого умножения).

Метод Эстрина

Впрочем, можно попробовать ещё ускорить данную операцию, если у нас в распоряжении имеется процессор, способный исполнять одновременно более одной команды (в данном случае имеется в виду не “многоядерность”, а наличие таких аппаратных особенностей как глубокий конвейер и/или суперскалярность и/или SIMD-команды). Для этого воспользуемся ассоциативностью суммы (впрочем, в плавающей точке ассоциативность теряется, поэтому результат может не совпадать с результатом функции poly_horner_fma). Кроме того, чтобы игра стоила свеч, нужен многочлен высокой степени (с большим количеством членов).

Один шаг: разбить по квадратам x.

// a3*x^3 + a2*x^2 + a1*x + a0 = (a3*x + a2)*x^2 + (a1*x + a0)
// sum(ai*x^i, i=0:(n-1)) = 
//  (...((a(n-1)*x + a(n-2))*x^2 + (a(n-3)*x + a(n-4)))*x^2 + ...)*x^2 + (a1*x + a0)
// -- метод Горнера по x^2
double poly_estrin2_fma(double x, double a[], size_t n)
{
  assert(n != 0);
  if (n == 1)
    return a[0];
 
  const auto x2 = x * x;
  auto s = fma(a[n-1], x, a[n-2]);
  for (--n; n > 2; n -= 2)
  {
    const auto p = fma(a[n-2], x, a[n-3]);
    s = fma(s, x2, p);
  }
 
  // если был нечётный размер массива, то n == 2, иначе 1
  return n != 2? s: fma(s, x, a[0]);
}

C++, HTML

Второй шаг: разбить по квадратам квадратов x.

Обобщение метода Горнера — метод Эстрина: дерево суммирования по квадратам
Обобщение метода Горнера — метод Эстрина: дерево суммирования по квадратам
// Вспомогательная функция для линейного двучлена.
inline double poly_1(double x, double a[])
{
  return fma(a[1], x, a[0]);
}

// Вспомогательная функция для квадратного трёхчлена.
inline double poly_2(double x, double a[])
{
  return fma(fma(a[2], x, a[1]), x, a[0]);
}

// Вычисление произвольного многочлена.
double poly_estrin4_fma(double x, double a[], size_t n)
{
  assert(n != 0);
  double s, x2 = x * x;
  
  // начало считаем по разным формулам в зависимости от остатка
  const size_t remainder = n % 4;
  switch (remainder)
  {
  case 1: s = a[n-1]; break;
  case 2: s = poly_1(x, a + (n-2)); break;
  case 3: s = poly_2(x, a + (n-3)); break;
  default:
    s = fma(fma(a[n-1], x, a[n-2]), x2, fma(a[n-3], x, a[n-4]));
    n -= 4;
  }

  if (n -= remainder)
  {
    auto x4 = x2 * x2;
    do
      s = fma(s, x4,
        fma(fma(a[n-1], x, a[n-2]), x2, fma(a[n-3], x, a[n-4])));
    while (n -= 4);
  }

  return s;
}


Вычисление произведения последовательности без промежуточного переполнения или исчезновения

Иногда случается ситуация, когда требуется вычислить произведение элементов некоторой конечной последовательности чисел в плавающей точке, хранящихся в массиве, читаемых с потока или вычисляемых по формуле. Так или иначе, в процессе такого вычисления может возникнуть неприятная ситуация: нехватка диапазона экспоненты для представления величины промежуточного результата, в то время как правильно вычисленных конечный результат представим. Например, для простоты представим себе последовательность, состоящую из двух тысяч двоек и двух тысяч “половинок” (0.5). Произведение равно единице. Однако, даже при использовании формата binary64 ещё на первой половине произойдёт переполнение, так как 22000 слишком велико, и общий результат наивного последовательного перемножения получится равным бесконечности. Если перемножать с конца к началу, то произойдёт исчезновение, так как 2–2000 слишком мало, и общий результат получится равным нулю.

Рассмотрим простой алгоритм, который позволяет избежать этой неприятной ситуации. Заметим, что нормальные числа в формате binary32 (например), находятся в полуинтервале [2–126, 2128), поэтому перемножение любой пары чисел из отрезка [2–63, 263] даст нормальное число. Это означает, что если текущее произведение и сомножители “загонять” в отрезок [2–63, 263] (домножая на 2–63 или на 263 в случае нормальных чисел), то переполнения или исчезновение не случится. При этом необходимо запоминать множитель — соответствующую степень двойки.

float cproduct(float a[], size_t n)
{
  float p = 1.f; // текущее произведение.
  int e = 0;     // сколько раз надо p домножить на 2 в 63.
  int s = 0;     // знак произведения.
  
  // Инвариант цикла:
  // p принадлежит [0x1p-63, 0x1p+63].  
  for (size_t i = 0; i < n; ++i)
  {
    float x = a[i]; // следующий сомножитель.
    if (signbit(x)) // поменять знак?
    {
      x = -x;
      s = 1 - s;
    }
    
    if (x < 0x1p-126f) // субнормальное или нуль?
    {
      if (x == 0) // произведение обнулилось?
        return s? -0.f: +0.f;
      
      e -= 2;
      x *= 0x1p+126f;
    }
    else if (x < 0x1p-63f) // маленькое?
    {
      --e;
      x *= 0x1p+63f;
    }
    else if (0x1p+63f < x) // большое?
    {
      if (isinf(x)) // произведение равно бесконечности.
        return s? -INFINITY: +INFINITY;
      
      ++e;
      x *= 0x1p-63f;
    }
    
    p *= x; // домножить.
    
    // Обеспечить выполнение инварианта цикла.
    if (p < 0x1p-63f) // маленькое?
    {
      --e;
      p *= 0x1p+63f;
    }
    else if (0x1p+63f < p) // большое?
    {
      ++e;
      p *= 0x1p-63f;
    }
  }
  
  // Результат не может быть представим в случае
  // e < -3 или e > 3.
  if (e < -3)
    return s? -0.f: +0.f;
  if (3 < e)
    return s? -INFINITY: +INFINITY;
  
  // Домножим на нужное количество 2 в 63 (или 2 в -63).
  if (e < 0)
  {
    do p *= 0x1p-63f;
    while (++e != 0);
  }
  else if (0 < e)
  {
    do p *= 0x1p+63f;
    while (--e != 0);
  }
  
  // Результат.
  return s? -p: p;      
}


Суммирование с компенсацией погрешности округлений

Предположим, требуется вычислить арифметическое среднее большого набора чисел. Простое последовательное суммирование может приводить к очень большой погрешности. Например, последовательное суммирование миллиарда единиц в формате binary32 даст 224, что примерно 17 миллионов, далеко не миллиард (почему?).

Алгоритм, называемый суммированием Кэхэна или компенсационным суммированием (У.Кэхэн), позволяет получить гораздо более точный результат, суммируя отброшенные при округлении основной суммы части в отдельной переменной — компенсации (в примере реализации — переменная comp). Относительная погрешность этого алгоритма ограничена сверху значением (цитир. по: Higham N. The accuracy of floating point summation // SIAM J. Sci. Comput. 1993. Vol. 14, No. 4. P.783–799)

Погрешность компенсационного суммирования
Погрешность компенсационного суммирования

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

float sum_kahan(float x[], size_t n)
{
  float sum = 0.f, comp = 0.f;
  for (size_t i = 0; i < n; ++i)
  {
    // В случае включенной агрессивной оптимизации операций с плавающей точкой
    // компилятор может сократить этот код до простого суммирования.
    const float y = x[i] - comp, t = sum + y;
    comp = (t - sum) - y;
    sum = t;
  }
  return sum;
}

Улучшенный вариант компенсационного суммирования Ноймайера имеет следующий вид:

float sum_neumaier(float x[], size_t n)
{
  float sum = 0.f, comp = 0.f;
  for (size_t i = 0; i < n; ++i)
  {
    const float xi = x[i], t = sum + xi;
    if (abs(sum) >= abs(xi))
      comp += (sum - t) + xi;
    else
      comp += (xi - t) + sum;
    sum = t;
  }
  return sum + comp;
}


Численное решение уравнений

Пусть дана непрерывная функция f : ℝ → ℝ и два числа a < b такие, что f(a) f(b) < 0, тогда по теореме о промежуточном значении найдётся c ∈ [a, b] такое, что f(c) = 0. Иными словами, уравнение f(x) = 0 имеет решение на отрезке [a, b]. Естественно поставить вопрос о том, как отыскать такое решение.

Метод деления отрезка пополам

Рассмотрим следующий алгоритм (метод бисекции или деления отрезка пополам):

  1. Положим m = (a + b) / 2.
  2. Если f(m) = 0, то c = m, закончить.
  3. Если f(a) f(m) < 0, то положить b = m.
  4. Иначе положить a = m.
  5. Продолжить с п.1.

На каждом шаге алгоритма отрезок [a, b] “уполовинивается”, но продолжает содержать некоторое решение уравнения. Впрочем, записанный выше алгоритм завершится только при чрезвычайно удачном стечении обстоятельств: если на каком-то шаге m совпадёт с решением. Легко привести пример, когда этого не произойдёт (пока рассуждаем в рамках действительных чисел):

f(x) = x2 – 2, a = 0, b = 2.

В указанном примере на каждом шаге, очевидно, истинно m ∈ ℚ, в то же время, решения уравнения f(x) = 0 не принадлежат ℚ. Таким образом, наш алгоритм за n шагов может гарантировать лишь сужение исходного интервала, содержащего корень, в 2n раз, но не обнаружение точного значения корня уравнения. Если достаточно знать корень с погрешностью ε > 0, то достаточно выполнить n = ⌈log2((ba)/ε)⌉ шагов.

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

Однако, при использовании физических компьютеров, абсолютно точная арифметика недоступна. При использовании чисел с плавающей запятой через конечное число шагов алгоритма станет истинно (a = m)∨(b = m) (почему?). Поэтому бисекцию можно записать следующим образом:

// Некая функция, которая задаёт правую часть уравнения.
float f(float);

// Решаем уравнение f(x) == 0 методом бисекции.
// Точность: 1 ULP.
float bisection(float a, float b)
{ 
  float fa = f(a), fb = f(b);
  assert(signbit(fa) != signbit(fb));
  
  while (true)
  {
    if (fa == 0.f)
      return a;
    if (fb == 0.f)
      return b;
    
    const auto m = 0.5f * (a + b);
    if (m == a || m == b) // <-- дальше "делить" отрезок бессмысленно
      return m;
      
    const auto fm = f(m);
    if (fm == 0.f)
      return m;
    
    if (signbit(fa) != signbit(fm))
    {
      b = m;
      fb = fm;
    }
    else
    {
      assert(signbit(fb) != signbit(fm));
      a = m;
      fa = fm;
    }
  }
}


Вилочный метод секущих

Рассмотрим ещё другой алгоритм — вилочный метод секущих (другое название: метод ложного положения regula falsi method), напоминающий метод бисекции в том плане, что по мере приближения корень всегда остаётся в некотором интервале. Всё отличие от предыдущего алгоритма будет состоять в способе вычисления очередной точки m. Вместо того, чтобы брать середину текущего отрезка, возьмём точку пересечения оси абсцисс секущей — прямой, проведённой через координаты (a, f(a)) и (b, f(b)). Делая так, мы “надеемся”, что на интервале [a, b] функция достаточно похожа на прямую. Это действительно так, если интервал достаточно мал, а функция f дифференцируема на нём.

Например, положим f(x) = x3/18 – 1, a = 1, b = 3.

Шаг 1: m = 30/13 ≈ 2.3077
Шаг 1: m = 30/13 ≈ 2.3077
Шаг 2: m = 1028/399 ≈ 2.57644
Шаг 2: m = 1028/399 ≈ 2.57644
Шаг 3: m ≈ 2.61485
Шаг 3: m ≈ 2.61485
Шаг 4: m ≈ 2.61996
Шаг 4: m ≈ 2.61996
первые шаги метода секущих (значения с индексом b — результат бисекции для сравнения)
Шаг m bm 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 намного быстрее подходит к корню, чем в методе бисекции.

// Решаем уравнение f(x) == 0 методом секущих.
// Точность: 1 ULP.
float bracket_secant(float a, float b)
{ 
  float fa = f(a), fb = f(b);
  assert(signbit(fa) != signbit(fb));
  
  while (true)
  {
    if (fa == 0.f)
      return a;
    if (fb == 0.f)
      return b;
    
    const auto m = (a*fb - b*fa) / (fb - fa);
    if (m == a || m == b) // <-- дальше "делить" отрезок бессмысленно
      return m;
      
    const auto fm = f(m);
    if (fm == 0.f)
      return m;
    
    if (signbit(fa) != signbit(fm))
    {
      b = m;
      fb = fm;
    }
    else
    {
      assert(signbit(fb) != signbit(fm));
      a = m;
      fa = fm;
    }
  }
}

Реализации этих двух методов, а также некоторых других методов, см. в примере.


Поиск минимума функции

Функция унимодальна на отрезке [a, b], если существует m ∈ [a, b] такая, что функция строго убывает на [a, m] и строго возрастает на [m, b]. Таким образом, функция достигает единственного на [a, b] минимума в точке m.

Рассмотрим задачу поиска минимума унимодальной на заданном отрезке [a, b] функции f : ℝ → ℝ.

Троичный поиск

Троичный поиск ternary search — простейший метод, решающий с заданной точностью ε поставленную задачу. Принцип его работы можно описать следующим алгоритмом:

  1. Если (ba) ≤ ε, то закончить.
  2. Положить m1 = a + (ba)/3.
  3. Положить m2 = b – (ba)/3.
  4. Если (a = m1)∨(b = m2)∨(m2m1), то закончить (исчерпана точность представления чисел).
  5. Если f(m1) < f(m2), то положить b = m2.
  6. Иначе если f(m2) < f(m1), то положить a = m1.
  7. Иначе положить a = m1, b = m2.
  8. Повторить с п.1.

Примечание. Если снять условие строгой монотонности функции на участках “до минимума” и/или “после минимума”, то троичный поиск может потерять минимум (случай f(m1) = f(m2)).

Пример троичного поиска

Шаг 1: min = –0.223108
Шаг 1: min = –0.223108
Шаг 2: min = –0.0200079
Шаг 2: min = –0.0200079
Шаг 3: min = –0.222849
Шаг 3: min = –0.222849
Шаг 4: min = –0.210824
Шаг 4: min = –0.210824
Шаг 5: min = –0.230368
Шаг 5: min = –0.230368

Легко заметить, что приближение минимума (для данного примера min = –0.230445) происходит немонотонно.

Метод золотого сечения

Существует метод, основанный на том же принципе, что и троичный поиск, но более экономный — позволяющий вычислять одно новое значение функции на каждом шаге, а не два, как троичный поиск. Это метод носит название метод золотого сечения. Сократить вычисления можно в том случае, если одна из предыдущих позиций, выбранных в [a, b] будет сохраняться на каждом шаге. Этого можно добиться, если делить отрезок не на трети, а в отношении 1:φ, где φ — “золотое сечение”. Для любых x > y > 0 таких, что φ = x/y, выполняется также (x + y)/x = x/y = φ. Отсюда легко получить тождество φ = 1 + 1/φ, позволяющее вычислить φ (положительный корень квадратного уравнения).

Например, возьмём

c = b – (ba)/φ,

d = a + (ba)/φ.

Теперь передвинем точку b в точку d и пересчитаем новые c‘и d’:

c’ = d – (da)/φ,

d’ = a + (da)/φ = a + (a + (ba)/φa)/φ = a + (ba)/φ2 = a + (ba)/(φ + 1) = (a(φ + 1) + ba)/(φ + 1) = b – (ba)φ/(φ + 1) = b – (ba)/φ = c.

Итак, алгоритм можно записать следующим образом:

  1. Положить c = b – (ba)/φ.
  2. Положить d = a + (ba)/φ.
  3. Если ((dc) ≤ ε)∨(ca)∨(bd), то закончить.
  4. Если f(c) < f(d), то положить b = d, d = c, c = b – (ba)/φ.
  5. Иначе если f(d) < f(c), то положить a = c, c = d, d = a + (ba)/φ.
  6. Иначе положить a = c, b = d, c = b – (ba)/φ, d = a + (ba)/φ.
  7. Повторить с п.3.

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

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

Шаг 1: min = –0.147106
Шаг 1: min = –0.147106
Шаг 2: min = –0.147106
Шаг 2: min = –0.147106
Шаг 3: min = –0.147106
Шаг 3: min = –0.147106
Шаг 4: min = –0.227573
Шаг 4: min = –0.227573
Шаг 5: min = –0.227573
Шаг 5: min = –0.227573

C++, HTML


Общее оглавление

Кувшинов Д.Р. © 2015