N-е число Фибоначчи за O(log N)

в 10:38, , рубрики: Алгоритмы, числа фибоначчи, метки:

Читая статью об устройстве на работу в ABBYY, встретил в ней упоминание задачи:

быстро – за O( log N ) арифметических операций над числами – найти N-е число Фибоначчи

Я задумался над ней и понял, что сходу в голову приходят только решения, работающие за время O(N). Однако позже решение было найдено.

Для удобства я перейду от обозначения N к n. Также я буду использовать обозначения множеств: mathbb{N}_0=left { 0, 1, 2, 3, ... right } — целые неотрицательные числа, mathbb{N}_1=left { 1, 2, 3, ... right } — целые положительные числа. Дело в том, что в разных математических традициях множество натуральных чисел может включать или не включать 0. Поэтому сейчас в международных математических текстах предпочитают указывать это явно.

Итак, решение

Кнут [1, с. 112] приводит матричное тождество следующего вида:
image
Тождество приводится без доказательства, но довольно просто доказывается по индукции.
Матрица слева иногда называется Q-матрицей.
Обозначим:
image
Из тождества получаем, что F_n=Q^{n-1}left (% 1,1 right ). Т.е. для вычисления F_n нам нужно вычислить матрицу Q^{n-1} и взять первый элемент первой строки (нумерация с 1).

Так как вычисление F_n свелось к возведению матрицы в степень, то рассмотрим этот процесс подробнее.
Пусть у нас имеется некоторая матрица M, которую необходимо возвести в степень nin mathbb{N}_1. Предположим также, что n является степенью двойки, т.е. n=2^i,iin mathbb{N}_0.
image можно представить в виде дерева:
N е число Фибоначчи за O(log N)
Имеется в виду:
image.
Соответственно, для вычисления матрицы M^n нужно вычислить матрицу image и умножить саму на себя. Для вычисления image нужно тоже самое проделать с image и т.д.
Очевидно, что высота дерева равна log n.
Оценим время вычисления M^n. Матрица M в любой степени имеет постоянный размер. Поэтому умножение двух матриц в любой степени можно выполнить за Oleft ( 1 right ). Всего таких умножений нужно выполнить log_2n. Следовательно, сложность вычисления M^n равна Oleft ( logn right ).

А если n — не степень двойки?

Теперь же встает вопрос: что делать, если n не является степенью двойки? Любое натуральное число n можно разложить как сумму чисел, являющихся степенью двойки, причем без повторений (мы занимаемся этим каждый раз, когда переводим число из двоичной системы счисления в десятичную). Т.е.
n=sum_{pin P}2^p,
где P_n subset mathbb{N}_0 — множество степеней, через которые выражается конкретное n. Если вспомнить, что n — это степень матрицы, то мы получим:
M^n=prod_{p in P_n}M^{2^p}.
Хотя в общем произведение матриц не коммутативно, т.е. порядок операндов при умножении важен, но для т.н. перестановочных матриц коммутативность соблюдается. Матрица M^i является перестановочной для M^j, i,j in mathbb{N}_0. Следовательно, нам не придется иметь в виду порядок операндов при умножении, что несколько облегчает дело.

Итак, алгоритм вычисления M^n можно представить в виде следующих шагов:

  1. Разложить n на сумму степеней двойки в виде множества P_n.
  2. Вычислить все элементы множества S=left {M^p|p in P_n right }.
  3. Вычислить M^n=prod_{s in S}s.

Оценим время работы данного алгоритма.
Первый шаг выполняется за время Oleft ( left lfloor log_2n rfloor+1 right )=Oleft ( log n right ), где lfloor log_2 n rfloor + 1 — количество двоичных разрядов в n.
Второй шаг выполняется за Oleft ( left left ( lfloor log_2 n rfloor + 1 right ) cdot log_2 n right )=Oleft ( left log^2 n right ), т.к. нам нужно выполнить не более lfloor log_2 n rfloor + 1 возведение матрицы в степень.
Третий шаг выполняется за Oleft ( left lfloor log_2 n rfloor right )=Oleft ( left log n right ), т.к. нам нужно выполнить умножение матриц не более lfloor log_2 n rfloor раз.

Оптимизация

Можно ли улучшить время работы данного алгоритма? Да, можно. Заметим, что во втором шаге в общем-то не оговаривается метод вычисления матриц, входящих во множество S. Нам известно о них то, что для каждой из них p является степенью двойки. Если вернуться к рисунку с деревом, это означает, что они все лежат на тех или иных ярусах дерева и для вычисления больших необходимо вычисление меньших. Если мы применим методику мемоизации и будем сохранять вычисленные матрицы на всех ярусах дерева, то сможем сократить время работы второго шага до Oleft ( log n right ) и время выполнения всего алгоритма также до Oleft ( log n right ). Именно это и требуется нам по условию задачи.

Код

Перейдем к кодированию. Я реализовал алгоритм на Python по двум причинам:

Получается такое:

class MatrixFibonacci:
    Q = [[1, 1],
         [1, 0]]

    def __init__(self):
        self.__memo = {}

    # умножение матриц
    # ожидаются матрицы в виде списка список размером 2x2
    def __multiply_matrices(self, M1, M2):
        a11 = M1[0][0]*M2[0][0] + M1[0][1]*M2[1][0]
        a12 = M1[0][0]*M2[0][1] + M1[0][1]*M2[1][1]
        a21 = M1[1][0]*M2[0][0] + M1[1][1]*M2[1][0]
        a22 = M1[1][0]*M2[0][1] + M1[1][1]*M2[1][1]
        r = [[a11, a12], [a21, a22]]
        return r

    # возведение матрицы в степень
    # ожидается p равная степени двойки
    def __get_matrix_power(self, M, p):
        if p == 1:
            return M
        if p in self.__memo:
            return self.__memo[p]
        K = self.__get_matrix_power(M, int(p/2))
        R = self.__multiply_matrices(K, K)
        self.__memo[p] = R
        return R

    # получение n-го числа Фибоначчи
    # в качестве n ожидается неотрицательное целое число
    def get_number(self, n):
        if n == 0:
            return 0
        if n == 1:
            return 1
        # разложение переданной степени на степени, равные степени двойки
        # т.е. 62 = 2^5 + 2^4 + 2^3 + 2^2 + 2^0 = 32 + 16 + 8 + 4 + 1
        powers = [ int(pow(2, b)) for (b, d) in enumerate(reversed(bin(n-1)[2:])) if d == '1' ]
        # тоже самое, но менее pythonic: http://pastebin.com/h8cKDkHX

        matrices = [ self.__get_matrix_power(MatrixFibonacci.Q, p) for p in powers ]
        while len(matrices) > 1:
            M1 = matrices.pop()
            M2 = matrices.pop()
            R = self.__multiply_matrices(M1, M2)
            matrices.append(R)
        return matrices[0][0][0]<source lang="ruby">

mfib = MatrixFibonacci()
for n in range(0, 128):
    num = mfib.get_number(n)
    print(num)

Бенчмарк

Я хотел сравнить производительность моего алгоритма с классическим итерационным алгоритмом, основанным на последовательном вычислении F_1,F_2,...,F_n с сохранением предыдущих двух чисел. Я реализовал его так:

def get_number(self, n):
    if n == 0:
        return 0
    a = 0
    b = 1
    c = 1
    for i in range(n-1):
        c = a + b
        a = b
        b = c
    return c

Применялась следующая методика тестирования производительности. Число a принимает значения a=10,20,...,1000. Для каждого a создаются объекты классов MatrixFibonacci и IterationFibonacci (класс, вычисляющий числа Фибоначчи итерационно). Далее 10 000 раз генерируется случайное целое n in left [ a-5,a+5 right ). Каждый из объектов вычисляет F_n с замером времени. На выходе теста получается множество строк вида: n <tab> T1 <tab> T2. По данным строкам строится график.
N е число Фибоначчи за O(log N)
Как видно из графика, матричный алгоритм начинает уверенно обгонять итерационный начиная где-то с n=300 (cтоит заметить, что F_{94} уже не влезет в 64 бита без знака). Мне кажется, что можно говорить о непрактичности данного алгоритма, хотя и имеющего хорошую теоретическую скорость.
Полный код с бенчмарком выложен тут. Код для gnuplot для построения графика по данным — тут.

P.S. Хотел попросить прощения за то, что TeX-формулы в виде картинок не всегда адекватно выравниваются внутри строчки. И похоже, что на Хабре нельзя сделать выравнивание руками.


  1. Дональд Кнут. Искусство программирования, том 1. Основные алгоритмы. — 3-е изд. — М.: «Вильямс», 2006.

Автор: Ivanhoe

* - обязательные к заполнению поля