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

Недавно на хабре появилась неплохая статья про вычисление 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 увеличится на пять, в то время как значения остальных трех переменных не изменятся. Это можно претставить в виде следующей матрицы:

1 0 0 0 0

0 1 0 0 0

0 0 1 0 0

0 0 0 1 0

5 0 0 0 1

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

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

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

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

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

Если матрицу поменять немного, убрав единицу в первом элементе первой строки:

0 0 0 0 0

0 1 0 0 0

0 0 1 0 0

0 0 0 1 0

5 0 0 0 1

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

A = 5

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

1 0 0 0 0

0 1 0 0 0

0 0 1 0 0

0 1 0 1 0

0 0 0 0 1

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

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

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

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

1 0 0 0 0

0 1 0 0 0

0 0 7 0 0

0 0 0 1 0

0 0 0 0 1

Перейдем к комбинированию команд. Пусть вектор 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

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