Арифметическое кодирование с использованием
целочисленной арифметики ограниченной точности

1. Вместо вещественных чисел, представлющих
вероятности символов, можно использовать рациональные
числа с ограниченными знаменателями. Вместо вещественного
массива вероятностей символов
    p = [ 0.2, 0.4, 0.4 ]
мы будем использовать целочисленные частоты символов:
    p = [ 2, 4, 4 ]
Обозначим   
    R = 2 + 4 + 4 = 10
Вероятности символов выражаются рациональными
числами prob_i = p_i/R
    prob = [ 2/10, 4/10, 4/10 ]
Установим точность вычислений, ограничив число
двоичных разрядов в записи целых чисел:
    precision = 20
Максимальное число 2**20 = 1024*1024 = 1048576

В предыдущих алгоритмах мы работали с вещественными
числами в диапазоне от 0. до 1. Теперь мы работаем с
целыми числами в диапазоне 0..2^precision - 1
(их можно представлять как числители рациональных
дробей со знаменателем 2^precision, но знаменатели
мы не храним). Числу 1.0 соответствует целое число
    whole = 2^precision
Половинке интервала [0, 1), т.е. числу 1/2 = 0.5,
соответствует целое цисло
    half = whole/2,
четвертинке 1/4 = 0.25 соответствует
    quoter = whole/4
Таким образом, все целые числа в нашем алгоритме,
несмотря на масштабирование, будут оставаться
в диапазоне 0..whole - 1. Перепишем последний алгоритм,
работающий с вещественными числами, но при этом
использующий масштабирование, заменив вещественную
арифметику операциями с целыми числами.

def ariphmEncoder_FinitePrec(p, message):
    '''An arithmetic encoder using finite-precision
    integer arithmetic with the rescaling technique.
    Input: p[i] is the frequency of i-th symbol in
           the alphabet expressed by an integer number.
           Zero symbol terminates a message;
           message is a list of indices of symbols.
    Return the sequence of bits that encodes the message'''

    precision = 20
    whole = 2**precision
    half = whole//2
    quoter = whole//4

    n = len(p)
    c = [0]*(n + 1)
    R = 0
    c[0] = 0
    for i in range(1, n+1):
        c[i] = c[i-1] + p[i - 1]
        R += p[i - 1]
        
    encodedMessage = []

    a = 0; b = whole; s = 0
    for j in message:
        w = b - a
        b = a + (c[j+1]*w)//R
        a = a + (c[j]*w)//R

        while b <= half or half <= a:
            if b <= half:
                # Select the left half and rescale
                a = a*2; b = b*2
                encodedMessage.append(0)
                # Emit 111...1 (s units)
                for i in range(s):
                    encodedMessage.append(1)
                s = 0
            else:
                assert half <= a
                # Select the right half and rescale
                a = (a - half)*2
                b = (b - half)*2
                encodedMessage.append(1)
                # Emit 000...0 (s zeros)
                for i in range(s):
                    encodedMessage.append(0)
                s = 0

        while quoter < a and b < quoter*3:
            s += 1
            # Rescale the segment [1/4, 3/4]
            a = (a - quoter)*2
            b = (b - quoter)*2
        
    assert a <= quoter or quoter <= b
    s += 1
    if a <= quoter:
        # Select [1/4, 1/2)
        # Emit 0111...1 (zero and s units)
        encodedMessage.append(0)
        for i in range(s):
            encodedMessage.append(1)
    else:
        assert 3*quoter <= b
        # Select [1/2, 3/4)
        # Emit 1000...0 (unit and s zeros)
        encodedMessage.append(1)
        for i in range(s):
            encodedMessage.append(0)
        
    return encodedMessage

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

В бесконечной точности при использовании ариметики
вещественных чисел нашему сообщению соответствует
двоичная дробь в диапазоне [0.0, 1.0):
    0.10110101100101000101011011110...1101
Закодированное сообщение состоит из битов дробной
части этого числа:
      10110101100101000101011011110...1101
Но мы при использовании арифметики конечной точности
можем использовать только числа с числом битов,
не большим, чем precision (в нашем примере 20).
Поэтому из этого числа, кодирующего наше сообщение,
мы берем только старшие precision бит:
      10110101100101000101011011110...1101
      ====================            
  Z = 10110101100101000101
При умножении Z на 2, что соответствует сдвигу Z
влево, мы стираем старший бит (в момент умножения
старший бит равен 0) и добавляем очередной 
непрочитанный бит сообщения. При декодировании
мы применяем масштабирование в соответствие с тем,
как оно применялось при кодировании.

Напомню схему Горнера: как вычислить число по
записи его цифр (число в системе счисления
с основанием, допустим, 10 -- это многочлен по
степеням десятки, коэффициенты которого -- это
цифры числа, записанные по убыванию степеней):
    1024 = ((1*10 + 0)*10 + 2)*10 + 4
Вычисляя целое число Z, соответствующее дробной части
закодированного сообщения y в двоичной записи, при
добавлении справа очередного бита y[i] мы умножаем Z
на 2 и прибавляем этот бит:
    z = z*2 + y[i]  

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

def ariphmDecoder_FinitePrec(p, encodedMessage):
    '''An arithmetic dencoder using finite-precision
    integer arithmetic with the rescaling technique.
    Input: p[i] is the frequency of i-th symbol in
           the alphabet expressed as an integer number.
           Zero symbol terminates a message;
           message is a list of indices of symbols.
    Return the sequence of bits that encodes the message'''

    precision = 20
    whole = 2**precision
    half = whole//2
    quoter = whole//4

    n = len(p)
    c = [0]*(n + 1)
    R = 0
    c[0] = 0
    for i in range(1, n+1):
        c[i] = c[i-1] + p[i - 1]
        R += p[i - 1]

    # print("p =", p, "R =", R)
    # print("c =", c)
        
    m = len(encodedMessage)
    z = 0
    i = 0
    while i < precision:
        z = z*2
        if i < m and encodedMessage[i] == 1:
            z += 1
        i += 1

    # print("z =", z, "binary =", bin(z))

    decodedMessage = []

    a = 0; b = whole
    while True:
        w = b - a
        for j in range(n):
            b0 = a + (c[j+1]*w)//R
            a0 = a + (c[j]*w)//R
            if a0 <= z < b0:
                # print("a0 =", a0, "b0 =", b0, "j =", j)
                decodedMessage.append(j)
                if j == 0:
                    return decodedMessage
                break

        a = a0; b = b0

        while b <= half or half <= a:
            if b <= half:
                # Select the left half and rescale
                a *= 2; b *= 2
                z *= 2
                if i < m and encodedMessage[i] == 1:
                    z += 1
                i += 1
            else:
                assert half <= a
                # Select the right half and rescale
                a = (a - half)*2
                b = (b - half)*2
                z = (z - half)*2
                if i < m and encodedMessage[i] == 1:
                    z += 1
                i += 1

        while quoter < a and b < quoter*3:
            # Rescale the segment [1/4, 3/4]
            a = (a - quoter)*2
            b = (b - quoter)*2
            z = (z - quoter)*2
            if i < m and encodedMessage[i] == 1:
                z += 1
            i += 1

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

>>> p = [2, 4, 4]
>>> x = [2, 1, 2, 2, 1, 2, 1, 1, 0]
>>> y = ariphmEncoder_FinitePrec(p, x)
>>> y
[1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1]
>>> xx = ariphmDecoder_FinitePrec(p, y)
>>> xx
[2, 1, 2, 2, 1, 2, 1, 1, 0]
>>> x
[2, 1, 2, 2, 1, 2, 1, 1, 0]
>>> x == xx
True

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


Почему арифметическе кодирование лучше, чем
коды Хаффмана?

Замечания

1. В кодах Хаффмана длина кода для каждой буквы --
это целое число, в частности, код любого символа
имеет длину не меньше единицы (энтропия 1-го символа
округлется до целого числа в сторону увеличения).
Насколько я помню, средняя длина закодированного
сообщения y = E(x) при использовании кодов Хаффмана
    len(y) <= H(p)*n + n,
где H(p) -- энтропия распределения вероятностей p,
n = len(x) -- длина исходного сообщения x.
При арифметическом кодировании средняя длина
закодированного сообщения длины n оценивается как
    len(y) <= H(p)*n + 2.

2. Надо иметь в виду, что мы рассматривали модель
исходных данных i.i.d. -- 
Indepedent Identical Distribution (каждая следующая
буква сообщения -- это случайная величина, не зависящая
от предыдущих букв, и распределенная так же, как и все
предыдущие).

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

Реальные языки не соответствуют модели i.i.d.,
там вероятность следующей буквы зависит от предыдущих
(в частном случае, от одной предыдущей). Но
арифметическое кодирование легко может учитывать
условные вероятности! Нужно просто делить отрезок
после предыдущей буквы или нескольких букв
в соответствии не с глобальными, а с условными
вероятностями (используем не массив вероятностей букв,
а, например, матрицу условных вероятностей букв при
учете одной предыдущей буквы).

3. Легко строится адаптивная модель арифметческого
кодирования, когда вероятности учитывают частоты
символов в уже закодированной части (намного удобнее,
чем, к примеру, адаптивные коды Хаффмана, которые требуют
нетривиальной процедуры перестройки кодирующего
дерева, как в алгоритме Виттера).