Пирамида

Иван Колесников в комментарии к прошлом посту подсказал мне использовать пирамиду для подсчета частот символов в реализации арифметического сжатия. Его идея оказалось очень удачной и позволила существенно ускорить процесс как кодирования, так и декодирования данных.

Дело в том, что в арифметическом сжатии при кодирования каждого символа нужно выполнить работу не только по его преобразованию в сжатый вид, но перед этим получить вероятность его появления в тексте, а также, в случае использования адаптивной статистической модели, еще и обновить саму модель.
В общем-то, обе операции не очень сложные, но вызываются очень часто, поэтому составляют изрядную долю времени работы алгоритма.
В принципе, можно использовать не адаптивную, а постоянную статистическую модель, но этот вариант в общем случае представляется мне менее удобным, так как в этом варианте данные сжимаются за два прохода. На первом проходе нужно собрать статистику для построения модели, а на втором – собственно, сжать сами данные.
Конечно, если у нас уже есть готовая статистическая модель данных, то мы можем сразу приступить к сжатию и оно, естественно, будет быстрее. Но это не избавляет нас от второго недостатка – необходимости в этом случае вместе со сжатым файлом хранить также и модель исходных данных, что, конечно, несколько снижает коэффициент сжатия. 

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

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

При арифметическом сжатии нам нужно, во-первых, получить две частоты с накоплением: кодируемого символа и предыдущего, а во вторых, обновить модель, то есть увеличить, как минимум элементы второго массива на единицу. Если поддерживается два массива, то также увеличить на единицу частоту кодируемого символа.
Сложность первого процесса крайне мала и равна O(1), а вот со вторым дело обстоит немного хуже – там сложность O(n), где n – размер алфавита.
Обычно в такой реализации процессы разнесены по времени: первый выполняется до кодирования символа, а второй – после. В принципе, их не так уж сложно объединить, но особой необходимости в этом нет, а реализация станет сложней и чуть-чуть медленее.

Теперь попробуем ускорить второй процесс (первый процесс с его O(1) ускорять уже некуда) с помощью пирамиды. Вот так она будет выглядеть для показанного выше примера.

Собственно, пирамида – это идеально сбалансированное дерево постоянной структуры, обычно частично упорядоченное. В нашем случае значение в каждом узле, кроме листовых, является суммой двух его сыновей. А листовые узлы хранят исходную таблицу частот (без накопления).
Постоянство структуры дает одно преимущество – пирамиду можно хранить без дополнительных расходов на описание структуры, что в двоичном дереве обычно реализуется с помощью двух указателей на сыновей. В пирамиде же ее просто можно сохранить в массиве, не тратя дополнительной памяти на указатели.
В этом примере голубые числа вверху слева от вершин – индексы в массиве. Переход от родителя (определяемого, например, индексом i) к сыновьям очень прост и эффективен:

  Left := i*2+1;
  Right := i*2+2;

В случае пирамиды операция поиска частоты с накоплением имеет сложность O(log₂n), такую же, как и изменение модели. Если выполнять оба процесса последовательно, то получим O(2log₂n), что лучше, чем при использовании обычного массива O(n+1).
Но в данном случае имеет смысл объединить оба процесса в один, и тогда общая сложность будет в два раза меньше: O(log₂n).

В результате реализации этих процессов при сжатии данных получаем компактный и красивый код:

function TArithmeticCoder.GetLimit(Symbol: byte; var Prev: cardinal): cardinal;
var
  i : cardinal;
begin
  i := Symbol + 255;
  Prev := 0;
  Result := Limits[i];
  Limit255 := Limits[0];
  while i <> 0 do
  begin
    if (i and 1) = 0 then // правый узел
      Prev := Prev + Limits[i-1]; // добавляем левый узел
    inc(Limits[i]);
    i := (i-1) shr 1; // переходим на уровень вверх
  end;
  inc(Limits[i]);
  Result := Prev+Result;
end;

Функция возвращает частоту с накоплением переданного ей в параметре символа, а частоту предыдущего возвращает в переменной Prev.
Так как при кодировании код символа известен, то проще идти в обратном направлении, от листьев к корню.
Limit255 – общее количество закодированных символов, кроме последнего. Необходимость этой переменной неоднозначна, в принципе в корне дерева (Limits[0]), лежит нужное нам значение, правда, в связи с тем, что мы при обходе дерева сразу же пересчитали модель, увеличенное на единицу. Так как к корню мы обращается на каждой итерации сжатия, то он гарантированно лежит в кэше, правда, из него надо вычесть единицу. В общем, не понятно.

Обратная функция, используемая в процессе декодирования:

function TArithmeticCoder.GetLimitIndex(Limit: cardinal; var Prev, Current : cardinal): byte;
var
  i : cardinal;
  Left, Right : cardinal;
begin
  i := 0;
  Limit255 := Limits[i];
  Prev := 0;
  while (i < 255) do
  begin
    Left := (i shl 1)+1;
    Right := (i shl 1)+2;
    inc(Limits[i]);
    if Limit < Prev+Limits[Left] then
      i := Left
    else
    begin
      Prev := Prev + Limits[Left];
      i := Right;
    end;
  end;
  Current := Prev + Limits[i];
  inc(Limits[i]);
  Result := i-255;
end;

Принимает в качестве параметра вероятность символа, возвращает код символа, а также вероятность этого и предыдущего символа в соответствии с моделью, ну и по ходу дела меняя саму модуль.

Тем не менее, этот код был написан последним, чисто для сравнения и не является рабочим. Дело в том, что я очень "жадный". И если посмотреть на объем памяти, занимаемый пирамидой, то видно что он почти в два раза больше, чем первоначальный массив.
Возможный механизм того, как можно сэкономить память, описал Иван, когда описывал свою идею в комментарии. Так как родитель и его два непосредственных потомка связаны простым соотношением

Root := Left + Right;

то очевидно, что необязательно хранить все три значения. Достаточно хранить либо правую, либо левую ветвь, особой разницы нет. Я выбрал левую, потому что для расчета частот с накоплением используется именно она. Но для обхода пирамиды все равно нужно рассчитывать значение правого потомка, так что подходы абсолютно одинаковы. Пример такого дерева приведен на рисунке ниже.

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

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

Теперь пирамида помещается ровно в ту же память, что и массив в исходном алгоритме с обычной таблицей и у меня есть вариант ее обхода:

function TArithmeticCoder.GetLimit(Symbol: byte; var Prev: cardinal): cardinal;
var
  i, Base, Size : integer;
  Root, Left, Right : cardinal;
begin
  i := 0; // индекс корня
  Root := Limits[i];
  Limit255 := Root;
  inc(Limits[0]);
  Base := 128;
  Size := 128;
  Prev := 0;
  while (Size <> 0) do
  begin
    Left := Limits[i+1];
    Right := Root - Left;
    if Symbol < Base then // идем влево
    begin
      inc(i);
      Size := Size shr 1;
      Root := Left;
      inc(Limits[i]);
      Base := Base - Size;
    end
    else // идем вправо
    begin
      i := i + Size;
      Size := Size shr 1;
      Prev := Prev + Left;
      Root := Right;
      Base := Base + Size;
    end;
  end;
  if Symbol < Base then
    Result := Prev + Left
  else
    Result := Prev + Right;
end;

И вариант для усеченной пирамиды, используемый в процессе декодирования:

function TArithmeticCoder.GetLimitIndex(Limit: cardinal; var Prev, Current : cardinal): byte;
var
  i, Base, Size : integer;
  Root, Left, Right : cardinal;
  LM : cardinal;
begin
  i := 0; // индекс корня
  Root := Limits[i];
  Limit255 := Root;
  inc(Limits[0]);
  Base := 128;
  Size := 128;
  Prev := 0;
  while {(Level <> 256)} (Size <> 0) do
  begin
    Left := Limits[i+1];
    Right := Root - Left;
    if Limit < Prev + Left then // идем влево
    begin
      inc(i);
      Size := Size shr 1;
      Root := Left;
      inc(Limits[i]);
      Base := Base - Size;
      LM := 1;
    end
    else // идем вправо
    begin
      i := i + Size;
      Size := Size shr 1;
      Prev := Prev + Left;
      Root := Right;
      Base := Base + Size;
      LM := 0;
    end;
  end;
  Result := Base - LM;
  if LM = 1 then
    Current := Prev + Left
  else
    Current := Prev + Right;
end;

Как видно из кода, функции сильно разбухли и стали менее понятны. Тем не менее, на x64 они работают немного быстрее, чем вариант с полной пирамидой.
Этому помогает, во-первых, меньший размер массива, а значит, и меньшее количество промахов в кэше L1. Ну а во-вторых, большое количество регистров в этой архитектуре, что позволяет все локальные переменный держать в регистрах. Не уверен, что такой код будет более быстр на x86 с её куцыми 8 регистрами общего назначения.

Стоит ли овчинка выделки? В смысле, надо ли использовать более сложный код, если выигрыш не так велик? Я думаю что да, ведь основная цель – реализация варианта с 16-битным алфавитом, в котором размер пирамиды увеличится в 256 раз, что существенно  увеличит вероятность промахов в кэше. Но вариант с 16-битном алфавите на полной пирамиде я не тестировал, так что привести численное сравнение не могу.