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.

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

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

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).