Используем быстрое возведение матриц в степень для написания очень быстрого интерпретатора простого языка программирования и...

Источник: habrahabr

SkidanovAlex

Недавно на хабре появилась неплохая статья про вычисление N-ного числа фибоначи за O(log N) арифметических операций. Разумный вопрос, всплывший в комментариях, был: "зачем это может пригодиться на практике". Само по себе вычисление N-ого числа фибоначи может и не очень интересно, однако подход с матрицами, использованный в статье, на практике может применяться для гораздо более широкого круга задач.

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

loop 1000000000
  loop 1000000000
    loop 1000000000
      a += 1
      b += a
    end
  end
end
end

Незамедлительно выведет a = 1000000000000000000000000000, b = 500000000000000000000000000500000000000000000000000000, несмотря на то, что если бы программа выполнялась наивно, интерпретатору необходимо было бы выполнить октиллион операций.

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

Для простоты ограничим наш интерпретатор четырьмя переменными - A, B, C и D. Для представления состояния интерпретатора в заданный момент будем использовать вектор размера пять, первые четыре элемента которого будут содержать значения четырех переменных соответственно, а последний будет на протяжении всей работы интерпретатора равен единице.

(A, B, C, D, 1)

В начале работы интерпретатора будем полагать значения всех переменных равными нулю.

(0, 0, 0, 0, 1)

Допустим, что первая операция в коде программы содержит строку

A += 5

Эффект этой команды заключается в том, что значение переменной A увеличится на пять, в то время как значения остальных трех переменных не изменятся. Это можно претставить в виде следующей матрицы:

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

v[0] * 1 + v[4] * 5

Так как v[0] содержит текущее значение в переменной A, а v[4] всегда равен единице, то

v[0] * 1 + v[4] * 5 = A + 5

Если вектор текущего состояния умножить на эту матрицу, полученный вектор будет соответствовать состоянию, в котором A на пять больше, что и требовалось.
Если матрицу поменять немного, убрав единицу в первом элементе первой строки:

Как и прежде, значения всех элементов кроме первого не изменятся, в то время как первый элемент станет равным v[4] * 5, или просто пяти. Умножение вектора текущего состояния на такую матрицу эквивалентно выполнению команды

A = 5

Посмотрим на такую матрицу:

Единственное отличие ее от единичной матрицы - это второй элемент в четвертой строке, который равен единице. Очевидно, что умножение вектора текущего состояния на эту матрицу не изменит значения в первом и последних трех элементах, в то время как значение второго элемента изменится на

v[1] * 1 + v[3] * 1

Так как v[1] содержит текущее значение переменной B, а v[3] содержит текущее значение переменной D, то умножение вектора состояния на такую матрицу эквивалентно выполнению команды B += D

Аналогично рассуждая можно понять, что умножение вектора состояния на следующую матрицу эквивалентно выполнению команды C *= 7

Перейдем к комбинированию команд. Пусть вектор v задает текущее состояние, матрица Ma соответствует команде A += 5, а матрица Mm соответствует команде A *= 7. Тогда, чтобы получить вектор r для состояния после выполнения этих двух команд, надо сначала умножить v на Ma, а затем на Mm:

r = v * Ma * Mm

Как и ожидается, умножение вектора начального состояния на эти две матрицы приводит в состояние, в котором a=35:

              1 0 0 0 0     7 0 0 0 0
              0 1 0 0 0     0 1 0 0 0
              0 0 1 0 0     0 0 1 0 0
              0 0 0 1 0     0 0 0 1 0
              5 0 0 0 1     0 0 0 0 1

0 0 0 0 1     5 0 0 0 1    35 0 0 0 1

Рассмотрим другой пример - пусть вместо умножения на семь, мы просто хотим прибавить пять к A семь раз.

r = v * Ma * Ma * Ma * Ma * Ma * Ma * Ma

Тут на помощь приходит ассоциативное свойство умножения матриц:

r = v * Ma * Ma * Ma * Ma * Ma * Ma * Ma = 
    v * (Ma * Ma * Ma * Ma * Ma * Ma * Ma) = 
    v * Ma ^ 7

Это пример простого цикла - вместо того, чтобы умножать v на Ma семь раз, достаточно возвести матрицу Ma в седьмую степень, после чего выполнить только одно умножение. Если бы мы хотели выполнить следующие две операции в цикле три раза:

A += 5
B -= C

(Пусть операция B -= C представляется матрицей Mbc), это бы выглядело следующим образом:

r = v * Ma * Mbc * Ma * Mbc * Ma * Mbc = 
    v * ((Ma * Mbc) * (Ma * Mbc) * (Ma * Mbc)) = 
    v * (Ma * Mbc) ^ 3

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

Рассмотренных примеров достаточно, чтобы начать работать над интерпретатором простого языка, поддерживающего присваивание, сложение, вычитание, умножение (только на константу) и циклы. Для этого мы научимся представлять любую такую программу в виде матрицы размера N+1 на N+1, где N - это количество переменных, которыми программа оперирует, после чего будем просто умножать вектор с начальным состоянием на эту матрицу.

Правила представления программы в виде матрицы очень просты:
1. Каждая отдельная команда представляется в виде матрицы, отличающейся от единичной одним элементом (или двумя для операции присваивания). Примеры таких матриц рассмотрены выше в этой статье.
2. Несколько подряд идущих команд представляются в виде матрицы, равной произведению матричного представления каждой отдельной команды.
3. Цикл представляется в виде матрицы, представляющей тело цикла, возведенной в степень количества итераций цикла.

Если у нас есть функция identity, возвращающая единичную матрицу:

def identity():
    return [[1 if i == j else 0 for j in range(REGS + 1)] for i in range(REGS + 1)]

То фукнция, строящая матрицу для команды r1 += r2 (где r1 и r2 - переменные) может выглядеть так:

def addreg(r1, r2):
    ret = identity()
    ret[r2][r1] = 1
    return ret

А для команды r += val (r - переменная, val - константа) вот так:

def addval(r, val):
    ret = identity()
    ret[REGS][r] = val
    return ret

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

Интерпретатор без циклов теперь пишется очень просто - пусть матрица mat соответствует уже прочитанному коду. В начале она равна единичной матрице, потому что пустая программа не меняет состояния. Затем мы считываем команды по одной, разбиваем их на три элемента (левый операнд, оператор, правый операнд), и в зависимости от оператора домножаем матрицу, соответствующую всей программе, на матрицу, соответствующую текущей команде:

def doit():
    mat = identity()
    while True:
        line = sys.stdin.readline().lower()
        tokens = line.split()
        if tokens[0] == 'loop':
            # тут будет код для циклов
        elif tokens[0] == 'end':
            return mat
        else:
            r1 = reg_names.index(tokens[0])
            try:
                r2 = reg_names.index(tokens[2])
            except:
                r2 = -1
            if tokens[1] == '+=':
                if r2 == -1: cur = addval(r1, long(tokens[2]))
                else: cur = addreg(r1, r2)
            elif tokens[1] == '-=':
            ....
        mat = matmul(mat, cur)

Осталось дело за малым - добавить поддержку циклов. Цикл возводит матрицу тела цикла в степень количества итераций цикла. Возведение в степень, как известно, требует только O(log N) операций, где N - это степень, в которую матрица возводится. Алгоритм возведения в степень очень прост:
1. Если степень равна нулю, вернуть единичную матрицу.
2. Если степень четная, пусть 2N, то можно рекурсивно вычислить M^N, а затем вернуть квадрат получившейся матрицы.
3. Если степень нечетная, пусть 2N+1, то достаточно рекурсивно посчитать значение M^2N, и вернуть полученную матрицу, умноженную на M.

Так как каждые две итерации степень сокращается в двое, сложность такого алгоритма логарифмическая.

def matpow(m, p):
    if p == 0: return identity()
    elif p % 2 == 0:
        tmp = matpow(m, p / 2)
        return matmul(tmp, tmp)
    else: return matmul(m, matpow(m, p - 1))

В интерпретаторе теперь осталось добавить одну строку:

        ...
        if tokens[0] == 'loop':
            cur = matpow(doit(), long(tokens[1]))
        ...

И интерпретатор готов.

Пример интерпретатора доступен на гитхабе. Весь код занимает меньше 100 строк.

Для теста скорости можно вернуться к уже упомянутым числам фибоначи. Например, такой код:

A = 1
B = 1
loop 100
  C = A
  C += B
  A = B
  B = C
end
end

Вычислит 101-ое и 102-ое числа фибоначи:

A = 573147844013817084101, B = 927372692193078999176

Замена 100 на 1000000 вычислит миллион первое и миллион второе числа за четыре секунды. Выполнение такой программы в лоб заняло бы гораздо больше, потому что программе приходится оперировать многотысячезначными числами. Если написать код, которому не приходится оперировать большими числами, например код для вычисления суммы арифметической прогрессии, приведенный в начале статьи, то количество итераций может уходить за рамки разумного, но код будет выполняться за доли секунды

loop 1000000000000000000000000000000000000000000000
  loop 1000000000000000000000000000000000000000000000
    loop 1000000000000000000000000000000000000000000000
      a += 1
      b += a
    end
  end
end
end

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