N-e число обобщённых Фибоначчи за O(log N)

В курсовой работе потребовалось написать алгоритм с логарифмической сложностью, который будет находить N-е число из последовательности Фибоначчи.

Алгоритм

Нашёл несколько статей по этой теме, во всех них рассматривалась классический ряд чисел Фибоначчи. Для него можно применять данную формулу:

Но у меня в работе использовались обобщённые ряды, в которых первые два числа — это ноль и некоторый параметр. Спустя часы поисков наткнулся на интересный комментарий, в котором автор приводит формулу циклического нахождения любых линейно-рекуррентных последовательностей (в т.ч. ряда чисел Фибоначчи).

Здесь Q это матрица 2×2, элементы которой можно запросто вычислить.

Подставив любые несколько чисел Фибоначчи, узнаём, что матрица Q = [[0,1], [1,1]].

Итоговую формулу матрицы, в которой содержится N-е число обобщённого ряда Фибоначчи, можно записать так так:

Fn — искомое число Фибоначчи,
C — ключ,
n — порядковый номер числа

Очевидно, что для эффективности всего алгоритма необходимо наиболее быстрым алгоритмом возвести в степень n матрицу Q. С этим мне помогла справиться данная статья. В ней объясняется, что для возведения матрицы в степень n можно разбить n на степени двойки и затем возводить матрицы только в эти степени. Такой подход даёт сложность O(log N).

Далее остаётся только умножить на матрицу [[C, C], [C, 0]] и получить элемент с индексом [0,1].

Реализация на Python

class FiboMatrix:     key = 0     matrix_cur = [[0,0], [0,0]]     matrix_one = [[0,1], [1,1]]      def __init__(self, key):         self.key = key         self.matrix_cur = [[key, key], [key, 0]]         self.PowsBuffer = {}         # Словарь для хранения возведённых в степень матриц      def MultiplyMatrix(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]         return [[a11, a12], [a21, a22]]      def PowMatrix(self, M: list, p: int):         """Возведение матрицы в степень.         Ожидаются матрица 2x2 в виде списка и степень.         """          if p == 1:             return M         if p in self.PowsBuffer:             return self.PowsBuffer[p]          K = self.PowMatrix(M, int(p / 2))         R = self.MultiplyMatrix(K, K)         self.PowsBuffer[p] = R         return R      def GetNum(self, n):         if n == 0:             return 0         if n == 1:             return 1                  # Разложение переданного номера номера числа n на степени двойки         powers = []         p = 0         for digit in reversed(bin(n)[2:]):             if digit == '1':                 powers.append(int(pow(2, p)))             p += 1          matrices = [self.PowMatrix(FiboMatrix.matrix_one, p)                     for p in powers]          # Возведение матриц в необходимые степени          while len(matrices) > 1:             M1 = matrices.pop()             M2 = matrices.pop()             R = self.MultiplyMatrix(M1, M2)             # Перемножение полученных матриц             matrices.append(R)          self.matrix_cur = self.MultiplyMatrix(self.matrix_cur, matrices[0])         # Умножение матрицы с F1 и F2 на матрицу, полученную возведением в степень         return self.matrix_cur[0][1]

Для сравнения эффективности был написан простейший аналог этого алгоритма с использованием цикла:

def Fib(num, key):     fib_1 = 0     fib_2 = key      for dec in range(num):         fib_1, fib_2 = fib_2, fib_1+fib_2      return fib_2

Бенчмарки подтверждают наши ожидания: классический алгоритм уже на 10000-м числе затрачивает в несколько раз больше времени.