Tag Archives: исчезновение

Эксперименты с hypot II

Небольшое дополнение к предыдущему посту. Ввиду низкого быстродействия frexp/ldexp может быть разумно заменить их на выбор множителя (степени двойки) банальными if-else.

inline float hypot_v2b(float x, float y)
{
  x = std::fabs(x);
  y = std::fabs(y);

  float
    min_ = x < y? x: y,
    max_ = x < y? y: x,
    mul, invmul;

  if (max_ <= 0x1p-126f)
  {
    invmul = 0x1p-126f;
    mul = 0x1p+126f;
  }
  else if (max_ <= 0x1p-63f)
  {
    invmul = 0x1p-63f;
    mul = 0x1p+63f;
  }
  else if (max_ <= 1.f)
  {
    invmul = 1.f;
    mul = 1.f;
  }
  else if (max_ <= 0x1p+63f)
  {
    invmul = 0x1p+63f;
    mul = 0x1p-63f;
  }
  else
  {
    invmul = 0x1p+126f;
    mul = 0x1p-126f;
  }

  max_ *= mul;
  min_ *= mul;

  return invmul * std::sqrt(max_* max_ + min_* min_);
}

Данный код предотвращает некорректное переполнение и исчезновение и обеспечивает погрешность 1 ULP. Конечно, гребёнка if-ов сказывается на скорости работы — у меня получилось, что hypot_v2b примерно втрое медленнее hypot_v1 (но всё же в десять раз быстрее hypot_v2a). Эта оценка довольно условна, так как скорость исполнения ветвлений зависит от множества факторов.

Аналогичный способ можно применить для замены frexp/ldexp при вычислении произведений. Например, если произведение получилось меньше по абсолютной величине наименьшего положительного нормального числа, то пересчитать его, предварительно домножив на 0x1p+63f; если произведение стало бесконечным (переполнение) или приблизилось по величине к наибольшему представимому числу, то пересчитать его, домножив на 0x1p-63f.

Эксперименты с hypot

В стандартной библиотеке математических функций C++ (унаследованной от C) имеется функция hypot(x, y), возвращающая евклидову длину вектора (x, y).

Зачем особая функция, если можно просто определить «в лоб»?

inline float hypot_v1(float x, float y)
{
  return std::sqrt(x * x + y * y);
}

(Здесь я использую float и только float по причине, о которой будет сказано ниже.)

Причин может быть несколько:

  1. Гарантировать отсутствие переполнения (из-за возведения в квадрат большого числа или сложения больших чисел) в случае, когда конечный результат (длина) представим.
  2. Гарантировать отсутствие влияния исчезновения (из-за возведения в квадрат слишком маленького числа) на результат. В частности, в случае, когда один из аргументов x, y есть нуль, а другой — число, hypot(x, y) должен возвращать модуль другого аргумента, даже если это субнормальное число.
  3. Обеспечивать погрешность лучше 1 ULP.

Легко показать, что hypot_v1 не удовлетворяет всем этим пунктам. Экспериментально установлено, что погрешность hypot_v1 на основной части диапазона, где не возникает ситуация вредоносного переполнения или исчезновения составляет 1 ULP (проверка для y из набора { 0.f, 1e-40f, 1e-30f, 1e-20f, 1e-15f, 1e-6f, 1e-2f, 1.f, 1e+2f, 1e+6f, 1e+15f, 1e+20f, 1e+30f } и для всех возможных значений x).

В активе hypot_v1 сравнительно большая скорость вычисления (у меня получилось, что она примерно в 10 раз быстрее, чем стандартная hypot на g++ -O3).

Можно ли получить какой-то иной компромисс в плане скорости и качества? В Википедии можно наблюдать следующий вариант (на момент написания). Ниже дан исправленный код:

inline float hypot_v2(float x, float y)
{
  x = std::fabs(x);
  y = std::fabs(y);

  const float
    min_ = x < y? x: y,
    max_ = x < y? y: x;

  if (max_ == 0.f || std::isinf(max_))
    return max_;

  const float q = min_ / max_;
  return max_ * std::sqrt(1.f + q*q);
}

Хорош ли этот код? Он страхует нас от патологического переполнения и исчезновения. Что можно сказать о его скорости? В моём эксперименте получилось, что hypot_v2 примерно вдвое медленнее hypot_v1, что можно считать приемлемым. Что можно сказать о погрешности? Экспериментально получено максимальное значение 2 ULP.

В стремлении убрать операции, вносящие погрешность, я написал иной вариант, использующий вынесение степеней двойки:

inline float hypot_v2a(float x, float y)
{
  x = std::fabs(x);
  y = std::fabs(y);

  if (x == 0.f)
    return y;
  if (y == 0.f)
    return x;

  float
    min_ = x < y? x: y,
    max_ = x < y? y: x;

  int min_e, max_e;
  min_ = std::frexp(min_, &min_e);
  max_ = std::frexp(max_, &max_e);

  min_ *= min_;
  max_ *= max_;

  return std::ldexp(
    std::sqrt(std::ldexp(min_, 2*(min_e - max_e)) + max_),
    max_e);
}

Этот вариант действительно обеспечивает погрешность <= 1 ULP (экспериментально). К сожалению, этот код в моём случае работал втрое медленнее стандартной hypot. (Что может говорить о потенциально низком быстродействии кода предложенного мной для перемножения массивов чисел.)

Теперь следует прояснить вопрос об экспериментальном вычислении погрешности. Так как используются числа float (IEEE-754 binary32, 23 двоичных разряда после запятой), и доступны числа double (IEEE-754 binary64, 52 двоичных разряда после запятой), то нетрудно организовать вычисление ближайшего к истинной величине значения в формате float, попросту выполнив вычисление с точностью double и округлив затем результат до float:

inline float hypot_v3(float x, float y)
{
  const double x_ = x, y_ = y;
  return (float)std::sqrt(x_ * x_ + y_ * y_);
}

Кстати, этот вариант также не страдает от переполнения и исчезновения и вычисляется лишь вдвое медленнее hypot_v1, обеспечивая погрешность <= 0.5 ULP. Он плох тем, что не масштабируется на double (хотя можно задействовать софтовое вычисление с точностью IEEE-754 binary128, заплатив за это резким падением производительности), зато очень удобен для определения качества других вариантов, включая стандартную hypot.

В случае g++ результат hypot везде совпал с результатом hypot_v3 (т.е. точность <= 0.5 ULP стандартной hypot обеспечена). Интересно, что для её вычисления используется древний (1993г.) код от Sun Microsystems, довольно зубодробительного вида.

В принципе, альтернативой hypot и hypot_v1 может стать их гибрид, использующий hypot_v1 там, где она работает хорошо (с погрешностью 1 ULP), а в остальных случаях (предположительно, редких) вызывающий hypot:

inline float hypot_v4(float x, float y)
{
  x = std::fabs(x);
  y = std::fabs(y);

  if (1e-19f <= x && x <= 1e+19f &&
      1e-19f <= y && y <= 1e+19f)
    return std::sqrt(x*x + y*y);

  return std::hypot(x, y);
}

Такой код способен работать со скоростью, в среднем близкой к скорости hypot_v1.

Борьба с переполнением и исчезновением при вычислении произведений

Пусть дан массив чисел с плавающей запятой, требуется найти произведение.

double product(const double x[], size_t n)
{
  double p = 1.;
  for (size_t i = 0; i < n; ++i)
    p *= x[i];
  return p;
}

Какие недостатки у приведённого выше очевидного решения?

Из трёх компиляторов (g++ 7.1, clang 4.0, icc 17.0) векторизацию данного цикла (параметры -O3 -mavx2) выполнил только icc, что лишь подтверждает известный тезис о том, что icc порой действует чересчур агрессивно: операция умножения в плавающей точке не ассоциативна.

Для любителей обобщённого кода можно записать то же самое, но для «любых чисел»:

#include <iterator>
#include <functional>
#include <numeric>

template <class It>
auto product(It from, It to)
{
  return std::accumulate(from, to,
    typename std::iterator_traits<It>::value_type(1),
    std::multiplies<>{});
}

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

Одно из возможных решений — складывать логарифмы, затем сумму возвести в степень:

\prod\limits_i{x_i} = \exp\sum\limits_i{\log{x_i}}.

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

В случае плавающей запятой на основе средств стандартной библиотеки C доступен способ получше: комбинация frexp и ldexp.

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

#include <cmath>
#include <iterator>

template <class It, class NT>
NT product(It from, It to, NT first)
{
  using std::frexp;
  using std::ldexp;
  int e = 0;
  NT p = frexp(first, &e);
  for (int e2 = 0; from != to; ++from)
  {
    p = frexp(p * *from, &e2);
    e += e2;
  }
  return ldexp(p, e);
}

template <class It> inline
auto product(It from, It to)
{
  return product(from, to,
    typename std::iterator_traits<It>::value_type(1));
}

Или то же самое в виде функтора:

template <class NT, class ET = long>
class Frexp_ldexp_product
{
  NT p;
  ET e;

public:
  using Number = NT;
  using Exponent = ET;

  explicit Frexp_ldexp_product(Number first = Number(1))
  {
    using std::frexp;
    int first_e = 0;
    p = frexp(first, &first_e);
    e = first_e;
  }

  Number operator()() const
  {
    using std::ldexp;
    return ldexp(p, e);
  }

  void operator()(Number next)
  {
    using std::frexp;
    int next_e = 0;
    p = frexp(p * next, &next_e);
    e += next_e;
  }

  Number significand() const
  {
    return p;
  }

  Exponent exponent() const
  {
    return e;
  }
};

// Пример использования: при обработке диапазона.
template <class It, class NT>
NT product(It from, It to, NT first)
{
  return std::for_each(from, to,
    Frexp_ldexp_product<NT>{first})();
}

Функтор удобнее, если значения в произведении генерируются на ходу, а не лежат в памяти. Кроме того, есть возможность получить отдельно мантиссу ∈ (0.5, 1] и порядок, даже если результирующее произведение не помещается в диапазон используемого типа чисел с плавающей точкой.

Напоследок хочется отметить, что если проблема переполнения/исчезновения снята, то рост накопленной погрешности остановить не так просто. Хотя результат произведения зависит от порядка, оценка относительной погрешности от порядка выполнения операций не зависит:

(x(1+n\varepsilon))(y(1+m\varepsilon)) \rightarrow xy(1+n\varepsilon)(1+m\varepsilon)(1+\varepsilon) \rightarrow xy(1+(n+m+1)\varepsilon).

В каком порядке не перемножай, всё равно получится относительная погрешность ne, где n — число сомножителей, e — половина «эпсилона» (ULP).