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

Недавно на хабре появилась неплохая статья про вычисление 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


На практике этот подход может применяться, например, в оптимизирующих компиляторах, которые могут таким образом сворачивать циклы с большим количеством итераций, оперирующие на небольшом количестве переменных.
Метки:
Поделиться публикацией
Похожие публикации
Комментарии 55
  • +9
    Это просто замечательно. Никогда не подозревал ранее, что благодаря матрицам можно вытворять такое. Спасибо за полезные примеры.
    • +2
      Действительно, хорошее направление для развития оптимизирующих компиляторов.
      Напомнило Задачу о порядке перемножения матриц
      • +3
        У вас нотация довольно необычная. Обычно в математике используются векторы-столбцы и умножают матрицы на векторы а не наоборот. Потратил минуту чтобы понять, почему же матрица выглядит именно так.

        А вообще довольно интересно.
        • 0
          Вообще-то в математике есть много чего, в том числе и row-vector наравне с column-vector (обычно их определяют как матрицы 1xN и Nx1).
          • +3
            Более того, столбцы можно считать обычными векторами, а строки — ковекторами (функционалами на векторах), тогда произведение строки на столбец будет соответствовать применению функционала к вектору, а столбца на строку — взятию тензорного произведения вектора на функционал (и получится линейный оператор).
            Таки тензоры рулят:)
            • 0
              Можно, только без векторного пространства и преобразований координат в нём они ничем не отличаются от матриц 1xN и Nx1 )
              • 0
                Ну это один из способов записи (при фиксированном базисе). А тензорное произведение — это обобщение произведения матриц, которое, в свою очередь, обобщение, например, применения функционала к вектору.
                А транспонирование матрицы — это свертка с \delta_i^j.
                В общем, мне лично нравятся такие теории, которые выявляют общее в на первый взгляд существенно различных вещах.
                Но это были так, мысли вслух)
                • 0
                  Да, в такие моменты задумываешся, почему не ходил на пары с векторного и тензорного анализа… интересно блин, а понять ничего не могу… )
                  • +2
                    Если интересно — почитайте «Алгебру» Винберга, там хорошо эти темы изложены.
        • –4
          а теперь, то же самое для

          loop 100000000
             psw = md5(psw + "salt")
          end


          :D
          • 0
            Гениально…

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

            p.s. а теперь добавьте сюда параллельные вычисления ведь одновременные вычисления с разными переменными неплохо вписываются в одну операцию с матрицами…
            • +2
              Я думаю, что способа достигнуть катарсиса нет. Мы можем свести к перемножению матриц только линейные операторы в пространстве линейных форм над <1, a, b, c, d...>. А, скажем, умножение переменной на переменную уже дает не линейную, а квадратичную форму.
            • +1
              Спасибо за статью. Есть вопрос: можно ли реализовать(имеется ввиду через матричные преобразования) битовые операции: NOT, OR, AND (все операнды — переменные)?
              • +1
                если вместо сложения взять XOR, вместо умножения — AND, алгебраические тождества будут выполняться (получим вычисления в поле вычетов по модулю 2). NOT можно получить формулой NOT x = 1 XOR x

                другой любопытный пример — пусть сложение OR, а умножение AND.
                если взять матрицу A смежности графа (A[i,j]=1 если между i и j есть ребро), то A^2 — матрица достижимости пути за 2 шага, A^3 — за 3 и т.д… A^n — матрица связности (An[i,j]=1 если есть путь от i к j, n — число вершин графа, любой кратчайший путь от i к j не длинее n, поэтому степень выше n рассматривать избыточно). возведение в степень n для поиска матрицы связности тоже можно делать этим быстрым алгоритмом
                • 0
                  Вот только матрицы перемножаются дольше, чем O(n^2), так что этот алгоритм будет работать дольше, чем O(n^2 log n), в то время как банальный поиск в ширину даст сложность O(n+m) для списка ребер (где m — количество ребер), или, для матрицы смежности, O(n^2).
                  • 0
                    но поиск в ширину/глубину — последовательный алгоритм (стек, очередь). распараллелить можно попробовать, но теории нет, чтобы сделать это формально, не задумываясь
                    • 0
                      Зато у него константа куда лучше, чем у быстрого перемножения матриц. И пригодные для практического использования алгоритмы умножения вроде бы имеют сложность O(n^{2+epsilon}) для, как правило, достаточно большого (>0.1) эпсилон.
                      Так что, поскольку это всё есть смысл обсуждать только для очень больших графов (десятки тысяч и больше вершин), то распараллеливание едва ли спасет алгоритм с худшей ассимптотикой и константой.
                      К тому же поиск пусть не в ширину, а в глубину вроде бы вполне красиво параллелится.
                      • 0
                        и всё-таки, у этого подхода есть практическое применение, опережающее обход в ширину.
                        это когда нужно найти матрицу расстояний для каждой пары i, j

                        только уже не булевские операции надо использовать, а сложение обычное, а умножение — min
                        • +1
                          Ну да, для матрицы расстояний это получается обобщением алгоритма Флойда (который получается, если вместо быстрого перемножения матриц использовать обычное). Но для матрицы достижимости такой метод всё-таки неинтересен.
                          • 0
                            Кстати, а алгоритмы быстрого перемножения матриц не завязаны на, например, дистрибутивность умножения (я этого не помню, а смотреть и думать сейчас лень)?
                            • 0
                              в системах операций OR/AND, XOR/AND дистрибутивность выполняется
                              • 0
                                Я про минимум (для нахождения матрицы расстояний).
                    • 0
                      Флойдом можно найти матрицу достижимости за куб — на практике это должно быть всегда быстрее, чем возведение матрицы смежности в степень.
                      • 0
                        Нет, конечно. Быстрое умножение матриц имеет не такую уж большую константу.
                        • 0
                          Даже если при больших N быстрое возведение в степень с быстрым умножением будут быстрее, чем флойд, я бы все равно остановился на флойде, хотя бы из-за его простоты.
                          Вариант с поисками в ширину тоже приятный — у него тоже кубическая сложность на плотных графах (а на разряженных еще меньше), но зато его можно легко запустить на нескольких процессорах сразу. Аналогично, если надо найти не матрицу достижимости, а все кратчайшие пути, на современном железе, наверное, N дийкстр будут очень эффективны, потому что их можно тоже легко распараллелить, сохраняя все ту же кубическую сложность.
                          • 0
                            Угу. А пузырек вместо быстрой сортировки из-за его простоты Вы не выберете? Как и перемножение матриц «в лоб» вместо быстрых алгоритмов?)
                            Простота реализации важна в некоторых случаях, но далеко не во всех. Вы, думаю, не обрадуетесь видеокартам в 5 раз дороже и в 10 раз медленнее из-за того, что их разработчики, как и Вы, предпочитают медленные, но простые алгоритмы.
                            А все стандартные алгоритмы быстрого перемножения матриц вроде бы прекрасно параллелятся.
                            • 0
                              Я думаю причиной этого небольшого разногласия является мое устаревшее представление о быстром умножении. Лучшее, что я знаю, работает за степень примерно 2.7. На графе с 10000 вершин разница с кубом будет в 15 раз, но быстрое возведение в степень потребует еще логарифм, который от 10000 даст 14, и смысл потеряется, потому что у быстрого умножения костанта выше, чем у флойда.

                              Если есть алгоритм быстрого умножения заметно быстрее, чем за (N^2.7), то тогда мои размышления о флойде, вероятно, заметно отличались бы.
                              • 0
                                Из применимых на практике вроде бы быстрее всего прямые суммы Шенхаге, которые вроде бы работают примерно за n^2.52. Что дает выигрыш в 63 раза по сравнению с лобовым.
                                С учетом логарифма — выигрыш в 6 раз на 10000 и в 15 на 100. Что, скорее всего, с лихвой перекрывает константу.
                                • 0
                                  Вообще я, быть может, высказался немного резко, приношу извинения. Просто мне очень уж не нравятся фразы вида «даже если по-другому будет лучше, всё равно сделаю как проще».
                                  • 0
                                    Асимптотически лучший алгоритм перемножения матриц имеет экспоненту 2.373 (т.е. сложность n2.373)
                                    • +1
                                      Вот только константа и затраты памяти там, вроде бы, совершенно ужасные.
                                      • 0
                                        Не «вроде бы», а точно :) Даже у Штрассена константа раз в 10-100 (зависит от реализации) больше, чем у алгоритма по формуле.
                          • +1
                            Есть задачи, которые решаются возведением матрицы смежности графа в степень.
                            Например: codeforces.ru/contest/147/problem/B.
                            • 0
                              Если конкретнее, то задача звучит так: кратчайшие пути между всеми всеми вершинами, состоящие ровно из N рёбер :)
                          • 0
                            Зависит от способа записи:)
                            NOT X при любом разумном способе записывается как 2^k — 1 — X.
                            А вот OR и AND, скорее всего, ни при какой разумной записи выразить будет нельзя. Разве что хранить вектор значений всех возможных функций наших переменных:)
                            • 0
                              я имел ввиду матрицы не из чисел, а из bool
                              • 0
                                Ну, в Z_2 AND вообще нелинейная операция и выразить ее линейно нельзя)
                                • 0
                                  не понял утверждения.
                                  умножение в действительных числах тоже не выражается через сложение.
                                  • 0
                                    Вопрос был «можно ли реализовать(имеется ввиду через матричные преобразования) битовые операции».
                                    Насколько я понимаю, имелся в виду тот же метод, что и в статье, т.е. вектор из переменных (и, быть может, побочных значений) мы линейным преобразованием переделываем в вектор функций от этих переменных (и, быть может, побочных значений).
                                    Утверждение было: преобразование (x, y) -> (x AND y, z(x, y)) нелинейно над Z_2^2 ни для какой функции z.
                                    • 0
                                      да, теперь понял. для данного примера нужно искать систему операций, где AND выполняла бы роль сложения
                                  • 0
                                    Z_2 операция AND и есть умножение
                            • 0
                              Произведение любых матриц ассоциативно, не только квадратных, разве нет?
                            • +5
                              Оригинальный подход — за него однозначный плюс :) Хоть и малоприменим на практике, так как позволяет брать только уж совсем специфичные случаи.
                              • 0
                                Простая идея, простая реализация, молодец! Подобный DSL очень бы пригодился на всяких codeforces/topcoder :)
                                • +1
                                  Отличная статья, приятно видеть такие на хабре. Она почти вернула мне веру в то, что математика в программировании таки небесполезна :)
                                  • +1
                                    Всегда знал, что матрицами можно «сворачивать» вычисления, использовал это и в расчетах четырехполюсников, и когда с 3D разбирался.

                                    Но никогда не понимал, ПОЧЕМУ это работает? ПОЧЕМУ блин так происходит, А?

                                    То есть, я пользуюсь аппаратом, но не понимаю, почему он работает.

                                    Я уже считаю, что это какое-то непознанное свойство материи или там пространства-времени, которое так отображается в математике нашего мира.
                                    • 0
                                      «Почему» — очень сложный вопрос. То, что он работает, можно доказать. Является ли это свойством реальности или лишь свойством нашей теории — сказать, скорее всего, невозможно.
                                      Причем описание матрицами линейных преобразований — это еще очень простой аппарат. Вот то, что работают дифференциальные формы, тензорный анализ на гильбертовых пространствах, обобщенные функции и т.д. — гораздо удивительнее.
                                      Пользуясь случаем, в очередной раз прорекламирую статью «Непостижимая эффективность математики в естественных науках».
                                    • +1
                                      Это действительно замечательно и я всерьёз намерен использовать это для ускорения вычислений. Жаль, правда, что такой язык не обладает полнотой, т.к. может работать только с циклами, для которых точно известно количество итераций, что даёт очень узкий класс решаемых задач (это даже не примитивно-рекурсивные, для которых известно максимальное количество итераций — ещё уже). И поэтому, наверное, языком программирования это вообще называться не может. Или я что-то проглядел?
                                      • +2
                                        Алгоритм кстати на столько прост, что его можно автоматически внедрять в код, подставляя его в код вместо циклов… ведь во время выполнения больше случаев, когда количество итераций известно…

                                        p.s. я считаю идею нужно развивать, ведь математические выражения можно упрощать не только через возведения в степень (правда алгоритм получается не совсем тривиальным)… основная польза от этого — нереально крутые возможности по распараллеливанию кода (матричные вычисления этому очень способствуют).
                                      • 0
                                        Любобытно.
                                        Возникла идея: а что если расширить цикл и на нецелые числа и вычислять степени матриц с помощью разложения в ряд Тейлора? Тогда возможно будет вычислять корни.
                                        Например вот так: loop 0.5. Не уверен, что это корректная операция, но попробовать стоит.
                                        • 0
                                          Вычислять корни не получится — умножать-то можно только на константу.
                                          умножение (только на константу)
                                          • –2
                                            Вы внимательно прочитали, что я написал?
                                            Я говорил про возведение матриц в дробную степень и умножение переменных на константу тут совсем не при чем.
                                            • +1
                                              Ну расскажите мне, как умножив на какую-либо константу половину раза, Вы сможете вычислить корень из переменной (корни из констант, очевидно, не нужны — это те же константы).

                                              Предположим, что мы как-то придумали возведение в степень. Скомбинировали тысячу различных извращений (включая ряды Тейлора, Маклорена, Лорана — кого пожелаете) и получили его (а точнее — её, чудесно-волшебную матрицу, возводящую элементы любого вектора в некоторую степень, допустим, 3). Т.е. мы получили для вектора v (пусть он будет таким: [ a b 1 ]) матрицу M, такую, что v * M = [ a³ … … ]. Тогда имеем для произвольного λ:
                                              λ[ a b 1 ] * M = [ λa λb λ ] * M = [ (λa)³ … … ]
                                              с другой стороны
                                              λ[ a b 1 ] * M = λ ([ a b 1 ] * M) = λ [ a³ … … ] = [ λa³ … … ]
                                              Таким образом, λa³ = (λa)³, a = 0 (возводить нули в степень скучно) => λ = λ³, но у нас-то λ любое!

                                              Очевидно, что если бы такая матрица M существовала, но она должна была бы зависеть от текущих переменных, т.е. являться функцией вектора M(v). А это ломает всё, что мы здесь придумали — цепочку вычислений
                                              v = v * M1v
                                              v = v * M2v

                                              v = v * Mnv
                                              (для упрощения скобки у вызовов функции опущены) Мы уже не сможем «свернуть» так же просто, как в статье. Разве что

                                              v = v * M1v * M2(v * M1v) * M3(v * M1v * M2(v * M1v)) * … * Mn(v * M1v * M2(v * M1v) * … * Mn-1(…))
                                              однако вычислять это в таком виде — очевидно, феерическая глупость.
                                        • +3
                                          Привет — два года спустя пользователь hx0 написал подробную статьи про свой модуль оптимизации байт-кода Python, основанный на алгоритме из твоей статьи. Оставляю этот комментарий, чтобы поместить линк на более новую статью: «Автоматическая оптимизация алгоритмов с помощью быстрого возведения матриц в степень».

                                          Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.