(495) 925-0049, ITShop интернет-магазин 229-0436, Учебный Центр 925-0049
  Главная страница Карта сайта Контакты
Поиск
Вход
Регистрация
Рассылки сайта
 
 
 
 
 

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

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

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

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

Ссылки по теме


 Распечатать »
 Правила публикации »
  Обсудить материал в конференции Embarcadero »
Обсудить материал в конференции Microsoft »
Обсудить материал в конференции Дизайн, графика, обработка изображений »
Написать редактору 
 Рекомендовать » Дата публикации: 07.08.2012 
 

Магазин программного обеспечения   WWW.ITSHOP.RU
Microsoft Office 365 Персональный 32-bit/x64. 1 ПК/MAC + 1 Планшет + 1 Телефон. Все языки. Подписка на 1 год.
Microsoft Office 365 Бизнес. Подписка на 1 рабочее место на 1 год
Delphi Professional Named User
Microsoft Office для дома и учебы 2019 (лицензия ESD)
Microsoft Office 365 для Дома 32-bit/x64. 5 ПК/Mac + 5 Планшетов + 5 Телефонов. Подписка на 1 год.
 
Другие предложения...
 
Курсы обучения   WWW.ITSHOP.RU
 
Другие предложения...
 
Магазин сертификационных экзаменов   WWW.ITSHOP.RU
 
Другие предложения...
 
3D Принтеры | 3D Печать   WWW.ITSHOP.RU
 
Другие предложения...
 
Новости по теме
 
Рассылки Subscribe.ru
Информационные технологии: CASE, RAD, ERP, OLAP
Безопасность компьютерных сетей и защита информации
Новости ITShop.ru - ПО, книги, документация, курсы обучения
Программирование на Microsoft Access
CASE-технологии
Один день системного администратора
Все о PHP и даже больше
 
Статьи по теме
 
Новинки каталога Download
 
Исходники
 
Документация
 
Обсуждения в форумах
Пишу программы на заказ профессионально (3248)
Пишу программы на заказ на языках Pascal (численные методы, списки, деревья, прерывания) под...
 
Написание программ для микроконтроллеров AVR, PIC, ARM, STM32 (25)
Напишу любую программу на любом искусственном языке. Профессиональный программист. Основная...
 
Разработка устройств на микроконтроллерах (39)
Профессиональный программист. Основная специализация: МИКРОКОНТРОЛЛЕРЫ, АССЕМБЛЕР для любых...
 
Пишу программы на заказ для студентов (254)
Пишу для студентов на с, с++, паскаль в средах ms visual studio, qt, builder, borland c, delphi....
 
Ремонт холодильников (1)
Если нужно починить срочно холодильник в Харькове то советую обращаться в эту...
 
 
 



    
rambler's top100 Rambler's Top100