Как стать автором
Обновить

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

Алгоритмы
Читая статью об устройстве на работу в 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] приводит матричное тождество следующего вида:
\begin{pmatrix} F_{n+1} & F_n\\ F_n & F_{n-1} \end{pmatrix} = \begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix}^n
Тождество приводится без доказательства, но довольно просто доказывается по индукции.
Матрица справа иногда называется Q-матрицей.
Обозначим:
Q = \begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix}
Из тождества получаем, что F_n = Q^{n-1}\left ( 1,1 \right ). Т.е. для вычисления F_n нам нужно вычислить матрицу Q^{n-1} и взять первый элемент первой строки (нумерация с 1).

Так как вычисление F_n свелось к возведению матрицы в степень, то рассмотрим этот процесс подробнее.
Пусть у нас имеется некоторая матрица M, которую необходимо возвести в степень n\in \mathbb{N}_1. Предположим также, что n является степенью двойки, т.е. n=2^i,i\in \mathbb{N}_0.
M^n можно представить в виде дерева:

Имеется в виду:
M^n=M^{n/2}\cdot M^{n/2}=...=\prod ^{n}_{1}M^1.
Соответственно, для вычисления матрицы M^n нужно вычислить матрицу M^{n/2} и умножить саму на себя. Для вычисления M^{n/2} нужно тоже самое проделать с M^{n/4} и т.д.
Очевидно, что высота дерева равна \log n.
Оценим время вычисления M^n. Матрица M в любой степени имеет постоянный размер. Поэтому умножение двух матриц в любой степени можно выполнить за O\left ( 1 \right ). Всего таких умножений нужно выполнить \log_2n. Следовательно, сложность вычисления M^n равна O\left ( \log n \right ).

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

Теперь же встает вопрос: что делать, если n не является степенью двойки? Любое натуральное число n можно разложить как сумму чисел, являющихся степенью двойки, причем без повторений (мы занимаемся этим каждый раз, когда переводим число из двоичной системы счисления в десятичную). Т.е.
n=\sum_{p\in 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.

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

Оптимизация

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

Код

Перейдем к кодированию. Я реализовал алгоритм на Python по двум причинам:
Получается такое:
class MatrixFibonacci:
    Q = [[1, 1],
         [1, 0]]

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

    def __multiply_matrices(self, M1, M2):
        """Умножение матриц
        (ожидаются матрицы в виде списка список размером 2x2)."""

        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

    def __get_matrix_power(self, M, p):
        """Возведение матрицы в степень (ожидается 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

    def get_number(self, n):
        """Получение n-го числа Фибоначчи
        (в качестве 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]

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=300 (cтоит заметить, что F_{94} уже не влезет в 64 бита без знака). Мне кажется, что можно говорить о непрактичности данного алгоритма, хотя и имеющего хорошую теоретическую скорость.
Полный код с бенчмарком выложен тут. Код для gnuplot для построения графика по данным — тут.

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


  1. Дональд Кнут. Искусство программирования, том 1. Основные алгоритмы. — 3-е изд. — М.: «Вильямс», 2006.
Теги:числа фибоначчи
Хабы: Алгоритмы
Всего голосов 93: ↑77 и ↓16+61
Просмотры61K

Похожие публикации

Лучшие публикации за сутки