Я участвовал в онлайн-группе чтения книги Thirty-three Miniatures: Mathematical and Algorithmic Applications of Linear Algebra математика Иржи Матушека. Это самая нетрадиционная книга о математике, с которой мне приходилось сталкиваться. Первые две главы посвящены способам быстрого нахождения чисел Фибоначчи. Традиционный, или итеративный метод нахождения чисел Фибоначчи (основанный на хранении промежуточных значений в памяти), который мы изучали на курсах программирования, линеен по времени. Но в книге представлена методика их вычисления приблизительно с логарифмической временной сложностью. Возможно, кто-то из вас знает эту методику, но для меня она была новой, и я решил, что ею стоит поделиться.

▍ Используем стратегию линейной алгебры для нахождения чисел Фибоначчи


Традиционный рекурсивный алгоритм нахождения чисел Фибоначчи с хранением значений в памяти выглядит так:

table = {}
def fib(n):
    if n in table:
        return table[n]
    
    if n == 0 or n == 1:
        result = n
    else:
        result = fib(n - 1) + fib(n - 2)
    
    table[n] = result
    return result

Это алгоритм линейного времени, потому что для вычисления n-ного числа Фибоначчи нам нужно вычислить все предыдущие n - 1 числа в последовательности. При этом возникает привычный вопрос: можно ли добиться лучшего?

Давайте разберёмся. Одна из проблем вышеприведённого алгоритма заключается в том, что для вычисления произвольного значения в последовательности необходимо вычислить все предыдущие значения. Следовательно, нам нужно найти более обобщённое решение, способное вычислять n-ное значение напрямую.

Одно из наблюдений: мы знаем значения базового случая (0 и 1), так что если мы сможем выразить вычисление n-ного числа Фибоначчи, исходя из базовых случаев, то достигнем своей цели. Давайте разовьём это наблюдение.

Базовые случаи последовательности:

F(0) = 0
F(1) = 1

Вычисление F(2) при помощи этих двух случаев выполняется тривиально:

F(2) = F(1) + F(0)

Аналогично, мы можем также выразить F(1), исходя из F(1) и F(0):

F(1) = 1 * F(1) + 0 * F(0)

Можно удобно выразить вычисление F(2) и F(1) в единой операции линейной алгебры:

$\begin{align} \begin{bmatrix} F(2) \\ F(1) \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} F(1) \\ F(0) \end{bmatrix} \\ \\ Или, \begin{bmatrix} F(2) \\ F(1) \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} \end{align}$


Давайте попробуем сделать то же самое для F(3):

F(3) = F(2) + F(1)
F(2) = F(1) + F(0)

Представив это в виде матричной операции, мы получаем следующее:

$\begin{bmatrix} F(3) \\ F(2) \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} F(2) \\ F(1) \end{bmatrix}$


Однако мы знаем значение вектора [F(2), F(1)] из предыдущего этапа; подстановка этого значения даёт нам следующее:

$\begin{bmatrix} F(3) \\ F(2) \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix}$


Давайте упростим:

$Обозначим\ матрицу\ \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} как\ M$


Это даёт нам:

$\begin{bmatrix} F(3) \\ F(2) \end{bmatrix} = M^2 \begin{bmatrix} 1 \\ 0 \end{bmatrix}$


Похоже, мы создаём паттерн для вычисления чисел относительно матрицы M и базовых случаев F(1) и F(0). Посмотрим, сохраняется ли этот паттерн для F(4)

F(4) = F(3) + F(2)
F(3) = F(2) + F(1)

Снова представив это в виде матричных операций, получим:

$\begin{align} \begin{bmatrix} F(4) \\ F(3) \end{bmatrix} &= \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} F(3) \\ F(2) \end{bmatrix} \\ \\ Или, \begin{bmatrix} F(4) \\ F(3) \end{bmatrix} &= M \begin{bmatrix} F(3) \\ F(2) \end{bmatrix} \end{align}$


Подставив значение вектора [F(3), F(2)] из предыдущего вычисления, получаем следующий результат:

$\begin{bmatrix} F(4) \\ F(3) \end{bmatrix} = M . M^2 \begin{bmatrix} 1 \\ 0 \end{bmatrix}$


Что преобразуется в следующее:

$\begin{bmatrix} F(4) \\ F(3) \end{bmatrix} = M^3 \begin{bmatrix} 1 \\ 0 \end{bmatrix}$


В общем случае мы можем увидеть, что этот паттерн повторяется и обобщается до следующего вида:

$\begin{bmatrix} F(n+1) \\ F(n) \end{bmatrix} = M^n \begin{bmatrix} 1 \\ 0 \end{bmatrix}$


Вы также можете доказать для себя, что это истинно, при помощи математической индукции.

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

▍ Быстрое вычисление степеней


Наивный алгоритм вычисления n-ной степени требует n умножений, а следовательно, линеен по времени. Но для этой операции существуют более быстрые алгоритмы. Например, если n — это степень двойки, то можно вычислить M^n многократным возведением в квадрат M, что позволяет сделать это за время log2(n).

Во втором томе «Искусства программирования» Дональд Кнут показывает более общий алгоритм вычисления x^n. В этом алгоритме используется свойство, что любое число можно разбить на суммы степеней двойки. Например, 5 можно записать как 2^2 + 2^0, а 9 можно записать как 2^3 + 2^0. Следовательно, можно записать x^n следующим образом:

$x^n = x^{k_1} * x^{k_2} * x^{k_3}...x^{k_m}$


Где k1, k2, …, km — степени двойки, сумма которых равна n.

Хотя эта идея проста, сам алгоритм для выполнения этого чуть сложнее, потому что не существует способа подобного разбиения напрямую. Вместо этого нужно многократно делить число на два и смотреть на остаток от результата. В третьей редакции «Искусства программирования», томе втором (страница 462) есть листинг следующего алгоритма для вычисления x^n:

Алгоритм для вычисления x^n. Основан на алгоритме, опубликованном на странице 462 «Искусства программирования», том второй, третья редакция

Если вы предпочитаете изучать код вместо абстрактного псевдокода, то метод linalg.matrix_power в numpy реализует именно этот алгоритм (показанный на изображении ниже).

Реализация функции linalg.matrix_power в Numpy

Какова сложность этого алгоритма? Цикл выполняется не более log2(n) шагов, потому что на каждом шаге мы уменьшаем n вдвое. И на каждом шаге цикла мы умножаем Z на себя, что приводит к log2(n) умножениям. Кроме того, когда мы получаем остаток 1, то также умножаем Z на Y, увеличивая количество умножений. Формально, общее количество умножений, выполняемое этим алгоритмом, будет таким:

$\lfloor \log_2(n) \rfloor + v(n)$


Где v(n) — количество битов со значением 1 в двоичной форме n.

Кроме того, это сложность нашего алгоритма нахождения чисел Фибоначчи, потому что единственная выполняемая им операция, зависящая отn — это вычисление M^n.

Получается, мы получили более эффективный алгоритм. Но прежде чем праздновать победу, нужно запустить алгоритмы и сравнить их скорость. Однако перед этим давайте рассмотрим ещё один альтернативный способ нахождения чисел Фибоначчи.

▍ Аналитическое решение для чисел Фибоначчи


Существует формула аналитического решения для прямого вычисления n-ного числа Фибоначчи, в которой используется золотое сечение:

$F(n) = \frac{1}{\sqrt{5}} \left(\left(\frac{1 + \sqrt{5}}{2}\right)^n - \left(\frac{1 - \sqrt{5}}{2}\right)^n\right)$


Вероятно, это гораздо более известная формула, выведенная в «Искусстве программирования», том первый (страница 79). В книге Thirty-three Miniatures представлено альтернативное доказательство этой формы с использованием линейной алгебры. Из соображений времени и объёма статьи я не буду повторять здесь эти доказательства, если интересно, можете изучить их в книгах. Я просто хотел показать методику, чтобы у нас был более полный бенчмарк. Давайте рассмотрим его.

▍ Сравнение разных алгоритмов нахождения чисел Фибоначчи


Итак, у нас есть три алгоритма для нахождения чисел, настало время проверить их и сравнить скорость. Производительность этих трёх методик показана ниже.

Время выполнения трёх алгоритмов нахождения чисел Фибоначчи для растущих значений n

Здесь:

  • fib(n) — традиционный алгоритм с линейным временем;
  • fib_fast(n) — ускоренный алгоритм, изученный нами во втором разделе; вычисляет M^n.
  • fib_closed_form(n) — формула аналитического решения на основании золотого сечения.

Это отличный пример различия скоростей алгоритмов с линейным и логарифмическим временем. По сравнению со скоростью алгоритма с линейным временем другие алгоритмы выглядят как плоские линии.

Также здесь стоит отметить производительность алгоритмов fib_fast и fib_closed_form. По сути, оба эти алгоритма вычисляют n-ную степень выражения, а значит, их сложность относительно n одинакова. Однако алгоритм на основании золотого сечения вычисляет n-ную степень скалярного значения, что гораздо быстрее, чем вычисление n-ной степени матрицы 2X2, отсюда и различия в их реальной скорости времени выполнения.

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

▍ Подведём итог


На этом мы заканчиваем экскурсию по алгоритмам нахождения чисел Фибоначчи. Когда я впервые увидел формулу линейной алгебры для нахождения чисел Фибоначчи, она показалась мне магией. Автор книги Thirty-three Miniatures Иржи Матушек сообщает, что источник этой методики неизвестен, поэтому мы можем никогда не узнать, как он был открыт. Я уверен, что экспериментируя с подобными последовательностями, вы обнаружите паттерны, способные повести вас по пути открытий и, возможно, способы обнаружения таких алгоритмов могут быть различными. Я показал вам лишь один из способов их нахождения, но могут существовать и способы получше. Если у кого-то возникнут более качественные рассуждения, то поделитесь ими в комментариях.

▍ Благодарности


Я бы хотел поблагодарить профессора Нилдхару Мисру за проведение сессии по этой теме и за рассуждения о задействованном математическом аппарате. Если вам интересны обсуждения тем из CS и математики, то можете подписаться на неё в X/Twitter.

Также я хотел бы поблагодарить Алехандро Пиада Морфиса за вычитку и отзывы по статье. Он пишет прекрасные статьи о computer science в своём блоге Mostly Harmless Ideas.

▍ Дополнительная литература и ссылки



Узнавайте о новых акциях и промокодах первыми из нашего Telegram-канала ????

Комментарии (15)


  1. deafcafe
    13.11.2023 14:32
    +5

    Карацубу забыли :)

    Fast Fibonacci algorithms (nayuki.io)


    1. mcferden
      13.11.2023 14:32

      Python и так использует Карацубу для больших чисел.


      1. Alexandroppolus
        13.11.2023 14:32

        в JS упоролись чуть больше и выбирают алгоритм в зависимости от длины


  1. Fedorkov
    13.11.2023 14:32

    (удалено)


  1. CrazyOpossum
    13.11.2023 14:32
    +3

    Как написан fib_closed_form? Возводим double в степень, вычитаем, делим на double, округляем? Тогда вопрос - начиная с какого n даблы начнуть давать ошибки?


    1. Pshir
      13.11.2023 14:32

      На самом деле, вторую скобку можно просто выкинуть, потому что там в n-ную степень возводится число по модулю меньше единицы. А если оценить погрешность, то получится, что ошибки начнутся в районе n=70.


    1. yizraor
      13.11.2023 14:32

      Красота этой формулы в том, что корни всегда сокращаются, а результат оказывается целым числом.
      Поэтому есть возможность рассчитать её без 'double' (к сожалению, ценой скатывания алгоритма в линейную сложность).


      1. wataru
        13.11.2023 14:32

        Можно и без этого. Правда по сути решение скатиться все в ту же линейную алгебру. Надо просто вести вычисления в кольце алгебраического расширения целых чисел через sqrt(5). По-человечески это значит, что числа у нас будут в форме a+b\sqrt5. Вот уже число - это вектор длины 2. Домножение на такое число - это умножение на матрицу 2x2 определенного вида. Вот и опять можно логарифмически возводить в степень матрицу.


  1. Brotherofken
    13.11.2023 14:32
    +5

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

    В аналитической формуле числа, возводимые в степень n, не из воздуха берутся.

    Т.е. ещё чуть-чуть текста добавить и было бы гораздо более полно.


    1. slonpts
      13.11.2023 14:32

      Допишете?



      1. erik-mv
        13.11.2023 14:32

        Если автор это сделал, то пришел бы к 3 случаю. Чтобы понять это, достаточно взглянуть на собственные значения матрицы https://matrixcalc.org/ru/vectors.html#eigenvectors({{1,1},{1,0}})


  1. bigcrush
    13.11.2023 14:32

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

    M = P \cdot D \cdot P^{-1}


  1. Gudd-Head
    13.11.2023 14:32
    +1

    не существует способа подобного разбиения напрямую. Вместо этого нужно многократно делить число на два и смотреть на остаток от результата

    Эээ... Целочисленное число в двоичном виде уже есть сумма двоек разной степени


    1. ksotar
      13.11.2023 14:32

      Да. От этого умножение/деление на 2 и выполняется сверхбыстрыми в сравнении с остальным операциями сдвига.