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