Некоторое время назад я заметил, что большое количество работников индустрии информационных технологий имеют крайне смутное представление о том, что такое информация, что одновременно смешно и печально. Кто-то по этому вопросу ударяется в неконструктивную философию, кто-то вспоминает полторы формулы из теории информации, не особо вникая в их смысл. Между тем, мне кажется, есть довольно хороший путь к пониманию, и этот путь пролегает через сжатие данных. Когда речь заходит о сжатии, большинство вспоминают про RLE, алгоритмы семейства LZ и коды Хаффмана (на них основана zlib, а что еще нужно?). Но это все равно что сказать "Программирование? Да-да, знаю. Фортран и Кобол".
Поскольку мне наивно кажется, что в сжатии данных я кое-что понимаю, я решил сделать образовательную серию постов, в которой мы, пользуясь школьной математикой и несколькими банальными идеями и наблюдениями, выведем некоторые важные факты теории информации и напишем на чистом функциональном языке компрессор, который будет сжимать лучше, чем zip, rar, bzip2 и 7z в их дефолтных режимах. Потом оглядимся вокруг, в том числе на физику, и в результате, надеюсь, что-то поймем об информации.

Часть 1. Изобретаем арифметик

Пусть у нас есть сообщение, которое мы хотим записать или передать, а потом прочитать. Например, "HELLO". Эти пять букв в кодировке ASCII или UTF-8 представлены пятью байтами, которые в шестнадцатеричной системе выглядят как 48 45 4C 4C 4F. Банальное наблюдение номер один: если мы запишем эти числа подряд, а перед ними поставим "0.", то получим вещественное число между 0 и 1:
0.48454C4C4F
Мы можем так записать любое сколь угодно длинное сообщение, и получить вещественное число между 0 и 1. Поскольку в вещественном числе бесконечно много цифр (мы просто привыкли опускать бесконечную последовательность нулей в конце), чтобы знать какие из нулей входят в сообщение, а какие нет, вместе с числом надо также передать количество значащих цифр, т.е. длину исходного сообщения. Так любое сообщение (хоть все содержимое интернета) можно закодировать двумя числами, а потом восстановить. Причем, если нам известна длина сообщения, то само сообщение можно представить бесконечным количеством чисел - незначимые цифры могут быть любыми. Наше "HELLO" может быть представлено любым числом от 0.48454C4C4F0000... до 0.48454C4C4FFFFF... Фактически, сообщению соответствует интервал [0.48454C4C4F, 0.48454C4C50). Кодируя сообщение символ за символом, мы фактически занимаемся тем, что постепенно сужаем и уточняем интервал: изначально у нас есть интервал [0, 1), первая шестнадцатеричная цифра сужает его до [0.4, 0.5), вторая до [0.48, 0.49), третья до [0.484, 0.485) и т.д. Да простят меня за столь подробное изложение столь очевидных вещей. Вот так процесс выглядит:

"A way to skateboard" (Watch on YouTube)

(Если флэш-ролик не показывается, то вот GIF 4.5 MB или AVI 2.6 MB)
Цветом фона в ролике показано двоичное разбиение вещественной оси, видно, как получается двоичная запись числа. На каждую из шестнадцатеричных цифр уходит по 4 бита.

А теперь давайте чуть-чуть модифицируем процесс. Будем точно так же последовательно уточнять интервал, выбирая один из 16 подинтервалов, соответствующий очередной цифре, но подинтервалы сделаем разного размера. Раньше они все были размером 1/16 от текущего интервала, а теперь назначим цифре 4 интервал размером 1/2, цифре С - 1/8, цифрам 5, 8 и F - 1/16, а между остальными цифрами разделим оставшееся пространство примерно поровну. И вот что получится:

"Richmond" (Watch on YouTube)

(GIF 3.5 MB или AVI 1.6 MB)
После кодирования 10 шестнадцатеричных цифр у нас получился интервал размером 1/2 * 1/16 * 1/2 * 1/16 * 1/2 * 1/8 * 1/2 * 1/8 * 1/2 * 1/16 = 1/8388608 = 1/(2^23). Мы видим, что на цифры с интервалом 1/16 у нас по-прежнему уходит по 4 бита, но на цифру С с интервалом 1/8 тратится только 3 бита, а на четверку с ее интервалом 1/2 всего один бит. Чтобы закодировать исходное сообщение, нам помимо его длины достаточно 23 битов числа из получившегося интервала. Сравните это с 40 битами в предыдущем случае. Поздравляю, мы изобрели арифметическое сжатие! Оно же арифметическое кодирование, оно же просто арифметик.

Что будет, если размер подинтервала для какого-то символа будет больше 1/2? На него уйдет меньше одного бита, что для кодов Хаффмана уже немыслимо. Например, если первый символ сообщения закодируется интервалом размером 2/3, а второй - интервалом размером 3/4, то получится интервал размером 2/3 * 3/4 = 1/2, на который нужен один бит. В общем случае для интервала S на него уходит log2(1/S) бит. Так, log2(1 / (2/3)) + log2(1 / (3/4)) = 0.585 + 0.415 = 1 = log2 (1/ (1/2)).

Как можно реализовать это все в программе, не прибегая к манипуляциям вещественными числами неограниченной длины (понятно, что float'ов тут не хватит)? Заметим, что если в процессе сужения интервала его верхняя граница оказалась меньше 1/2, то первый бит лежащего в конечном интервале числа будет нулем. И наоборот: если нижняя граница поднялась до 1/2, то она уже не опустится, и первый бит будет единицей. Этот бит можно вывести в файл, а рабочий интервал увеличить вдвое. Так размер рабочего интервала всегда будет между 1/2 и 1, и огромная точность для манипуляций уже не потребуется. Можно вообще обойтись одними целыми числами, работая как бы с дробями в fixed point: числа от 0 до 0.FFFFFF представить целыми числами от 0 до FFFFFF. И поскольку выводить результат побитово не слишком удобно и эффективно, можно получать его побайтово: выводить старший байт и увеличивать рабочий интервал, когда он будет размером менее 1/256 и старший байт обеих границ окажется одинаковым. Такая целочисленная побайтовая версия арифметика называется range coder. В его реализации есть пара тонких моментов, например, что делать, если интервал уже сильно сузился, но его верхняя и нижняя границы еще недостаточно похожи (например, 0.4FF8 и 0.5011). В такой ситуации можно видеть, что первый байт будет либо 4F, либо 50, поэтому можно запомнить значение 4F в переменной-кэше, не выводя пока в файл, интервал все-таки увеличить, а потом, когда границы таки сойдутся и судьба того байта определится, вывести сперва кэш (плюс один если надо), а потом очередной байт.

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

Напишем range coder, который при кодировании будет принимать интервалы (два рациональных числа a/b и c/b в виде трех целых чисел a,b,c) и записывать результат (биты длинного вещественного числа из конечного интервала) в буфер, а при декодировании будет уметь по данным из буфера и переданному ему размеру текущего интервала говорить, где в этом интервале лежит заветное число, что позволит на каждом шаге выбирать нужный подинтервал и тем самым извлекать закодированный символ. Состояние кодера можно описать такой структурой:

:: *RCState = { buf   :: !*Buffer
              , low   :: !Int //нижняя граница рабочего интервала
              , range :: !Int //размер рабочего интервала
              , val   :: !Int //текущее число (используется при декодировании)
              , ffnum :: !Int
              , cache :: !Int // переменные для "сложных" ситуаций
              , carry :: !Int
              }

newEncoder :: Int -> RCState
newEncoder sz = { buf = newBuffer sz, low = 0,
  range = 0xFFFFFF, val = 0, ffnum = 0, cache = 0, carry = 0 }

Восклицательный знак означает строгое значение (не откладывать его вычисление лениво), звездочка - атрибут уникальности, дающий право менять содержимое по месту, не копируя. См. мой недавний пост про Clean и уникальные типы.
Функция encode кодирует один интервал (один символ). Нижняя граница интервала задается числом cumFreq/totFreq, размер интервала - freq/totFreq. Почему они так называются, будет ясно дальше - это обычно частоты символов: cumulative frequency, frequency, total frequency. Нижняя граница rc.low рабочего интервала растет на cumFreq/totFreq от его размера,а сам размер интервала rc.range уменьшается до freq/totFreq от его старого значения.

encode :: !*RCState !.Int !.Int !.Int -> *RCState
encode rc cumFreq freq totFreq = renorm rc1 where
    quant = rc.range / totFreq
    low1 = (rc.low + cumFreq * quant) bitand 0xFFFFFF
    rc1  = { rc & low = low1
           , range = quant * freq
           , carry = if (low1 < rc.low) (rc.carry+1) rc.carry
           }

renorm rc
  | rc.range << TOP = renorm (shiftLow rc)
                    = rc

shiftLow rc = if (rc.low < Thres || rc.carry > 0) rc1 rc2 where
  range1 = rc.range << 8
  low1 = (rc.low bitand 0xFFFF) << 8
  rc1 = { rc & range = range1, low = low1, ffnum = 0, cache = rc.low >> 16,
          buf = buf1, carry = 0 }
  buf1 = writeBytes [toChar (rc.cache + rc.carry) : repeatn rc.ffnum (toChar (rc.carry-1))] rc.buf
  rc2 = { rc & low = low1, range = range1, ffnum = rc.ffnum+1 }

Полученный интервал подвергается перенормировке: пока он меньше порогового значения TOP, он увеличивается в 256 раз, а старший байт записывается в буфер или кэш. Между записью в кэш и моментом, когда границы интервала достаточно сойдутся, чтобы определиться со значением невыведенного в буфер байта, может потребоваться увеличивать интервал несколько раз. Например, при рабочем интервале [0.4FFFFFF0, 0.50000002). Тогда после записанного в кэш значения могут идти несколько байтов FF или 00. Их количество хранится в переменной ffnum, а выбор между FF и 00 (и между 4F и 50 в этом примере) - в переменной carry.

Когда все символы закодированы, остается вывести в буфер данные о текущем интервале, и можно вернуть его содержимое в качестве результата сжатия.

finishEncode :: *RCState -> .{#Char}
finishEncode rc = bufferContents rc1.buf where
  rc1 = iter 4 shiftLow rc

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

newDecoder :: {#Char} -> *RCState
newDecoder str = iter 4 loadByte rc where
  rc = { buf = (buf,0), low = 0, range = 0xFFFFFF,
         val = 0, ffnum = 0, cache = 0, carry = 0 }
  buf = {c \\ c <-: str}

loadByte rc
  # (c, buf1) = readByte rc.buf
  = { rc & buf = buf1, val = ((rc.val bitand 0xFFFF) << 8) + toInt c }

При декодировании нам нужно, имея длинное вещественное число, понять, в какие интервалы оно попадает. Для этого для каждого символа нужно сначала получить точку в текущем интервале, по которой мы определим текущий подинтервал (его номер будет декодированным символом). Для этого передадим в декодер знаменатель totFreq и получим числитель freq.

getFreq :: !*RCState !.Int -> *(Int,*RCState)
getFreq rc totFreq
  #! rc = { rc & range = rc.range / totFreq }
  #! freq = rc.val / rc.range
  = (freq, rc)

Потом, уже зная подинтервал, нужно сузить рабочий интервал до него:

decode :: !*RCState !.Int !.Int -> *RCState
decode rc cumFreq freq = load rc1 where
  rc1 = { rc & val = rc.val - cumFreq * rc.range, range = rc.range * freq }
  load rc
    | rc.range < TOP = load (loadByte {rc & range = rc.range << 8})
                     = rc

Если при сужении размер рабочего интервала оказался меньше TOP, то его надо увеличить, а текущее значение rc.val сдвинуть влево на то же число бит и младший байт пополнить из буфера.

Итак, мы умеем закодировать сообщение последовательностью интервалов. Но какие интервалы выбрать? Мы видим, что чем больше интервал для символа, тем меньше на него тратится битов. Предположим, что разбиение на интервалы фиксировано, т.е. не меняется в процессе кодирования. Рассмотрим для начала вариант с алфавитом из двух символов: 0 и 1. Пусть в сообщении N символов: N0 нулей и N1 единиц, N0+N1=N. Пусть ноль кодируется интервалом размера S, а единица - интервалом размера 1-S. Тогда нули займут N0*log2(1/S) бит, а единицы - N1*log2(1/(1-S)). Суммарный размер закодированного сообщения:

Он минимален, когда аргумент логарифма минимален, т.е. когда обратная к нему величина максимальна. Приравняем ее производную нулю и найдем S:

Получается, что размер закодированного сообщения из нулей и единиц будет минимален, когда размер интервала для нуля S = N0/N, а для единицы - N1/N. Т.е. когда интервалы совпадают с вероятностями символов. Если у нас в алфавите не два символа, а много, то можно считать, что N1 - количество всех символов, отличных от некоторого выбранного, а S-1 - их совместный интервал. Тогда рассуждения, формулы и результат получаются те же самые: размер минимален, когда интервалы соответствуют вероятностям. Если Pi = Ni / N - вероятность i-го символа, то среднее количество бит на символ в оптимально закодированном сообщении получается