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

IEEE-754

Введение

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

• дроби вида m/n: пара целых чисел (m, n);
• числа с фиксированной запятой fixed point: целое число m представляет дробь вида m·r^–p, где r и p — заранее заданные константы;
• числа с плавающей запятой floating point: пара целых чисел (m, p) задает дробь вида m·r^p, где только основание r задано заранее. Таким образом, позиция запятой (порядок) p (в r-ичной системе счисления) указана в самом представлении числа.

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

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

Формат

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

В таблице ниже приведены сведения о четырёх основных двоичных (т.е. r = 2) форматах стандарта IEEE-754. Число ε — наименьшее положительное представимое число такое, что 1 + ε ≠ 1 в заданном формате. Смысл понятий “денорм.” (денормализованные числа) и “норм.” (нормализованные числа) дан ниже. В качестве числа d в таблице указано минимальное число десятичных знаков после запятой, которого достаточно для сохранения двоичного представления числа при преобразовании в десятичную форму и обратно. Необходимо помнить, что многие десятичные дроби являются бесконечными в двоичной системе: например, 0.2 не представимо точно в указанных форматах. Двоичные же дроби хотя и представимы точно конечным числом знаков, могут потребовать сотни десятичных цифр для точной записи.

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

Запись 0.(M) или 1.(M) означает, что m бит мантиссы выписываются после точки в двоичной записи числа. Под нормализацией числа понимается его домножение на такую степень двойки, что слева от точки получается единица. Денормализованные числа требуют недоступных значений экспоненты для нормализации.

Замечания

• Денормализованные числа позволяют постепенно приблизиться к нулю (с постоянным шагом).
• Наличие двух нулей (отрицательного и положительного) иногда полезно: так, при получении очень малых величин, округлившихся к нулю, мы по крайней мере сохраним знак (с какой стороны от нуля).
• Тем не менее оба нуля считаются равными друг другу. Кроме того, полагается sqrt(–0) = –0, log(–0) = –∞, ±0^±0 = 1.
• Бесконечность можно получить, например, поделив ненулевое число на нуль, или любым другим способом, при котором результат вычислений пусть даже и конечен, но слишком велик по абсолютной величине, чтобы быть представленным в используемом формате.
• Так как нулей два, то, например, (–1)/(–0) = +∞.
• Нечисло можно получить, если, например, вычесть бесконечность из бесконечности или попробовать вычислить квадратный корень отрицательного числа (т.е. когда результат операции не определен).
• Арифметические операции и большинство математических функций возвращают нечисла, если хотя бы один аргумент или операнд является нечислом. Поэтому появившееся в одном месте нечисло может породить лавину нечисел в последующих вычислениях. Эта ситуация называется “заражение нечислами” NaN plagueing.
• Для нечисел операторы сравнения всегда возвращают false (нечисло не равно даже само себе).
• Как правило, процессоры гораздо медленнее выполняют арифметические операции с денормализованными числами, бесконечностями и нечислами, чем с нормализованными числами, что может вызывать резкое падение производительности в некоторых случаях.
• Сумма и произведение чисел с плавающей запятой, строго говоря, не являются ассоциативными операциями.

Режимы округления

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

При округлении к ближайшему представимому выбор в пользу округления половинок именно к чётному (как правило, используется именно этот вариант, например, на x86) или нечётному вместо “школьного” округления “от нуля” был сделан с целью балансировки накапливающихся округлений и уменьшению случайной погрешности.

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

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

Литералы констант с плавающей точкой обязаны содержать точку (целая или дробная часть числа при этом могут быть опущены, если они нулевые) или символ, отделяющий порядок от коэффициента (e, E или p, P), в противном случае они трактуются как целые числа. Например, 10¹⁰ как число в плавающей точке будет 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)·2⁴ = 171.80078125.

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

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

Языки C и C++ предоставляют три встроенных типа для работы с числами с плавающей точкой:

• float обычно соответствует IEEE-754 binary32. Литералы этого типа задаются с помощью суффиксов f, F, например 3.14f;
• double обычно соответствует IEEE-754 binary64. Данный тип используется “по умолчанию” для литералов без суффиксов;
• long double на x86 нередко отождествляется с 80-битным форматом сопроцессора Intel 8087 (использование которого на платформе x86-64 объявлено устаревшим), в то время как Visual C++ использует представление, идентичное типу double, а некоторые другие компиляторы могут использовать IEEE-754 binary128. Литералы данного типа задаются с помощью суффиксов l, L, например 1e-1000L.

Стандарт гарантирует лишь отношение float ⊆ double ⊆ long double.

Часть характеристик форматов с плавающей запятой, используемых конкретным компилятором, можно получить, подключив заголовочный файл cfloat, где определен ряд макросов. Ниже * соответствует FLT для типа float, DBL для double и LDBL для long double.

• FLT_RADIX база системы счисления r (обычно 2);
• FLT_ROUNDS режим округления по умолчанию;
• *_MANT_DIG разряды полной мантиссы (m + 1);
• *_DIG достаточное число десятичных знаков после запятой d;
• *_MIN_EXP минимальная смещенная экспонента для нормализованных;
• *_MIN_10_EXP минимальная десятичная экспонента для нормализованных;
• *_MAX_EXP максимальная смещенная экспонента для нормализованных;
• *_MAX_10_EXP максимальная десятичная экспонента;
• *_MIN наименьшее положительное нормализованное число;
• *_TRUE_MIN наименьшее положительное (денормализованное) число;
• *_MAX наибольшее конечное представимое число;
• *_EPSILON наименьшее положительное представимое число ε такое, что 1 + ε ≠ 1 в заданном формате.

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

Заголовочный файл cmath

Заголовочный файл cmath предоставляет ряд определений математических функций (вычисляемых приближённо в плавающей точке, тем не менее, для простоты в обозначениях используются действительные числа вплоть до всей числовой оси). Здесь приводятся сведения по определениям для C++, принимающим разнотипные значения (float, double и long double). Если предел той или иной функции при стремлении к нулю (с той или иной стороны) или бесконечностям определён (включая значения ±∞), то он берётся в качестве значения этой функции в соответствующих предельных точках (поэтому многие интервалы и полуинтервалы в колонке “область значений” могут считаться отрезками).

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

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

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

• FE_DOWNWARD — округление к минус-бесконечности;
• FE_TONEAREST — округление к ближайшему (стандарт не указывает, какой вариант округления “половинок” будет использоваться при использовании данного режима);
• FE_TOWARDZERO — округление к нулю;
• FE_UPWARD — округление к плюс-бесконечности.

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

• int fegetround() — функция возвращает текущий режим округления.
• int fesetround(int mode) — функция устанавливает новый режим округления. Возвращает 0 в случае успеха и не 0 в случае ошибки.

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

Точность выполнения операций часто измеряется в единицах, равных весу самого младшего разряда мантиссы результата, называемых 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, интерпретируемого как целое, на единицу получаем:

• следующее представимое число b = (1 + ε) a;
• либо +∞, если a — уже наибольшее нормализованное число;
• либо нечисло, если 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;
}

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

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

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

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

Пусть дан многочлен вида Σa_i xⁱ, где a_i — коэффициенты, а 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 ещё на первой половине произойдёт переполнение, так как 2²⁰⁰⁰ слишком велико, и общий результат наивного последовательного перемножения получится равным бесконечности. Если перемножать с конца к началу, то произойдёт исчезновение, так как 2^–2000 слишком мало, и общий результат получится равным нулю.

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

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 даст 2²⁴, что примерно 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) = x² – 2, a = 0, b = 2.

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

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

Однако, при использовании физических компьютеров, абсолютно точная арифметика недоступна. При использовании чисел с плавающей запятой через конечное число шагов алгоритма станет истинно (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) = x³/18 – 1, a = 1, b = 3.

Легко заметить, что, с одной стороны, метод секущих уже не гарантирует нам сколь угодно малый интервал — в нашем примере бесконечное число шагов даст 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] функции f : ℝ → ℝ.

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

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

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

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

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

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

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

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

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

c = b – (b – a)/φ,

d = a + (b – a)/φ.

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

c’ = d – (d – a)/φ,

d’ = a + (d – a)/φ = a + (a + (b – a)/φ – a)/φ = a + (b – a)/φ² = a + (b – a)/(φ + 1) = (a(φ + 1) + b – a)/(φ + 1) = b – (b – a)φ/(φ + 1) = b – (b – a)/φ = c.

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

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

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

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

C++, HTML

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

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