Monthly Archives: Ноябрь 2014

Xor-Δ-кодирование

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

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

  • a xor 0 ≡ a;
  • a xor a ≡ 0.

Откуда легко выводится тождество (a xor b) xor b ≡ a.

Таким образом, если мы определим множество функций { xora | a ∈ Z, xora(x) ≡ a xor x }, то единица этого множества — тождественное преобразование xor0, и каждый элемент этого множества является обратным к самому себе относительно операции композиции функций (в частности, множество таких функций является абелевой группой относительно операции композиции).

Благодаря своим свойствам, операция xor очень популярна в области криптографических алгоритмов.

Итак, я написал следующий C++-код, заменяющий операции вычитания и сложения в дельта-кодировании/декодировании на операцию поразрядного исключающего или.

void xor_delta_encode(char *begin, size_t n)
{
  char prev = begin[0], t;
  for (size_t i = 1; i < n; ++i)
  {
    t = begin[i];
    begin[i] ^= prev;
    prev = t;
  }
}

void xor_delta_decode(char *begin, size_t n)
{
  char prev = begin[0];
  for (size_t i = 1; i < n; ++i)
    prev = begin[i] ^= prev;
}

Тестировочный код целиком: Bitbucket.

Если запустить этот код, то можно увидеть довольно забавное свойство (которое и сподвигло меня написать этот пост).

Повторное xor-дельта-кодирование приводит к повторению префиксов исходной последовательности. А именно (это можно показать «на бумажке», сокращая элементы благодаря свойствам операции xor): после 2n повторений имеем совпадение с исходной последовательностью префикса длины до 2n. Через 2N повторений операции кодирования, где N = ceil(log2(длина последовательности)), имеем восстановление исходной последовательности, «круг замкнулся». Повторение обратной операции (xor-дельта-декодирования) перебирает те же 2N возможных последовательностей, но уже в обратном порядке.

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

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

Итак, пусть дан массив чисел, представляющий набор коэффициентов многочлена по целым степеням 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! Я провёл несколько тестов данного кода с многочленами разного размера и он дал правильный результат, однако я не даю никаких гарантий относительно его корректности и, тем более, реального быстродействия (это я не проверял).