Tag Archives: FMA

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

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

Итак, пусть дан массив чисел, представляющий набор коэффициентов многочлена по целым степеням x. Я буду обозначать его a, а его размер – n.

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

double poly(double x, double a[], int n)
{
  double s = 0;
  for (int i = 0; i < n; ++i)
    s += a[i] * pow(x, i);
  return s;
}

Выдав такое, например, на собеседовании, кандидат рискует быть отправленным куда подальше.

Запретив вызывать дорогую функцию pow, мы получим, скорее всего, следующий вариант

double poly(double x, double a[], int n)
{
  double s = 0, xi = 1;
  for (int i = 0; i < n; ++i)
  {
    s += a[i] * xi;
    xi *= x;
  }
  return s;
}

Впрочем, может оказаться, что студент учился неплохо и вспомнит тот стандартный метод вычисления многочленов, который мы справедливо ожидали бы увидеть и на собеседовании. Метод Горнера. Данный метод позволяет сэкономить одно умножение на каждой итерации цикла суммирования. Будем считать, что пустой массив нам не передадут: в конце концов, пустой массив коэффициентов означает, что функция, которую мы пытаемся вычислить, не определена. И кстати, заменим-ка int…

double poly(double x, double a[], size_t n)
{
  assert(n != 0);
  double s = a[--n];
  while (n != 0)
    s = s * x + a[--n];
  return s;
}

Движение назад по массиву коэффициентов несколько неприятно, но с этим можно смириться.

Метод Горнера основан на часто встречающейся тернарной операции, называемой умножение-сложение (multiply-add «mad»/«madd», multiply-accumulate «mac»/«macc»). Многие процессоры (CPU и GPU) поддерживают данную операцию на уровне машинных команд (нередко и для целых чисел тоже). Что до плавающей точки, то существует два варианта реализации: с двумя округлениями (последовательное умножение и сложение) и с одним округлением (слитое умножение-сложение, fused multiply-add, «fma»). Второй вариант может быть довольно интересен с точки зрения уменьшения погрешности вычислений и был добавлен в стандарт IEEE-754 2008 года. Процессоры семейства x86 с недавних пор располагают аппаратной реализацией данного действия (наборы команд FMA4, FMA3). Программная реализация возможна, но будет на порядок проигрывать аппаратной с точки зрения быстродействия (которое для машинных команд FMA близко или равно быстродействию обычного умножения). Далее, стандарт языка C принял FMA ещё раньше: ISO C99 math.h включает функции fma, fmaf, fmal, которые могут задействовать соответствующие ресурсы процессора, либо быть основаны на программной реализации. ISO C++11 также предполагает наличие fma в cmath (в MS Visual C++ есть поддержка, начиная с версии 12 (2013)). Аппаратная поддержка fma открывает возможность для эффективных программных реализаций операций деления и взятия квадратного корня с разной степенью допустимой погрешности.

Заметим, что компилятор по умолчанию не имеет права заменять выражения вида (a*b+c) на инструкцию fma или вызов аналогичной функции, так как это влияет на результат (разная погрешность). Поэтому, если мы хотим меньшую погрешность данной операции, то следует явно вызывать функцию fma. Если мы хотели бы меньшую погрешность, но не готовы платить за это резким падением скорости вычислений, то C99 предоставляет макрос FP_FAST_FMA, который должен быть определён, если для вычисления fma используется аппаратная реализация.

Итак, используя fma, получаем вариант

double poly_fma(double x, double a[], size_t n)
{
  assert(n != 0);
  double s = a[--n];
  while (n != 0)
    s = fma(s, x, a[--n]);
  return s;
}

К сожалению, данная реализация далеко не идеальна с точки зрения быстродействия. Современные процессоры обычно включают два и более конвейеров, способных выполнять инструкции, оперирующие значениями в плавающей точке. Однако, чтобы загрузить их работой, выполняемые инструкции должны быть максимально независимы друг от друга. В методе Горнера каждая следующая команда fma использует результат предыдущей (сумму s). Ситуация обостряется при возможности задействовать SIMD-команды, например, Intel SSE или ARM NEON – метод Горнера, увы, сугубо последователен.

Метод Эстрина использует ассоциативность и дистрибутивность операций сложения и умножения, чтобы построить сбалансированное дерево вычисления многочлена. На схеме ниже приведён пример такого дерева для размера восемь (седьмая степень).

Метод Эстрина для многочлена седьмой степени

Это соответствует следующему коду

const auto
  x2 = x * x,
  a10 = fma(a[1], x, a[0]),
  a11 = fma(a[3], x, a[2]),
  a12 = fma(a[5], x, a[4]),
  a13 = fma(a[7], x, a[6]),
  x4 = x2 * x2,
  a20 = fma(a11, x2, a10),
  a21 = fma(a13, x2, a12),
  result = fma(a21, x4, a20);

Параллельное вычисление ветвей дерева позволяет эффективно загрузить конвейеры CPU, либо использовать SIMD-операции. Сложно сказать, насколько эффективна будет, например, реализация данного метода для произвольного массива в рекурсивной форме – накладные расходы могут перевесить. Поэтому при работе с произвольным многочленом представляется разумным явно вычислять лишь нижние уровни дерева, суммируя их результаты тем же методом Горнера, но уже для некоторой степени x (квадрат для нижнего уровня, четвёртая степень для двух нижних уровней и т.д.).

Например, если раскрыть только самый нижний уровень дерева и явно вычислять комбинации вида (a[2*i+1]*x + a[2*i]), то получим примерно следующий код

double poly_estrin2_fma(double x, double a[], size_t n)
{
  assert(n != 0);
  if (n == 1)
    return a[0];

  // один шаг схемы Эстрина + метод Горнера + fma
  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, если был нечётный размер массива
  return n != 2? s: fma(s, x, a[0]);
}

Но по теоретическим прикидкам получается, что для того чтобы загрузить два конвейера, следует раскрыть хотя бы два нижних уровня дерева. Я выразил это в более обобщённом виде с возможностью выбирать способ вычисления умножения-сложения через параметр шаблона.

struct Multiply_add
{
  template <class Float>
  static Float madd(Float a, Float b, Float c)
  {
    return a * b + c;
  }
};

struct Fused_multiply_add
{
  template <class Float>
  static Float madd(Float a, Float b, Float c)
  {
    using std::fma;
    return fma(a, b, c);
  }
};

template <class Madd = Multiply_add>
class Estrin_4 : Madd
{
  using Madd::madd;
public:

  template <class Float, class Float_it>
  static Float compute(Float x, Float_it a, size_t n)
  {
    assert(n != 0);
    Float s, x2 = x * x;
    // начало считаем по разным формулам в зависимости от остатка
    const size_t remainder = n % 4;
    switch (remainder)
    {
    case 1: s = a[n-1]; break;
    case 2: s = madd(a[n-1], x, a[n-2]); break;
    case 3: s = madd(madd(a[n-1], x, a[n-2]), x, a[n-3]); break;
    default:
      s = madd(madd(a[n-1], x, a[n-2]), x2, madd(a[n-3], x, a[n-4]));
      n -= 4;
    }

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

    return s;
  }
};

// пример использования класса
inline double poly_estrin4_fma(double x, double a[], size_t n)
{
  return Estrin_4<Fused_multiply_add>::compute(x, a, n);
}

Beware! Я провёл несколько тестов данного кода с многочленами разного размера и он дал правильный результат, однако я не даю никаких гарантий относительно его корректности и, тем более, реального быстродействия (это я не проверял).

Реклама