Формула Бине без плавающей точки

Как хорошо известно, числа Фибоначчи – это целочисленная последовательность, первые два члена которой равны единице, а каждый последующий равен сумме двух предыдущих. За 500 лет, прошедших с момента ввода этой последовательности в математический обиход, она основательно изучена. Открыто много интереснейших формул с участием чисел Фибоначчи. Но одной из “непреходящих” учебных задач является вычисление чисел Фибоначчи. Для этого придумано много способов: от прямой рекурсии, основанной на формуле:

до матричного метода, описанного, например, в книге Д.Кнута [1].  Большая часть этих подходов (кроме матричного метода Кнута) основаны на рекуррентных свойствах последовательности Фибоначчи и позволяют вычислить величину Fn в лучшем случае за  время O(n). Матричный метод Кнута (использующий матричное возведение в степень) позволяет вычислить число Фибоначчи за логарифмическое время [2].

Особняком в этом ряду алгоритмов располагается формула Бине (известная еще Муавру) имеющая вид:

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

Сказанное означает, что вычисления не будут точными; в них вносится погрешность ограничения. Мне однажды попалась на глаза публикация [3] в которой использовалась формула Бине для вычисления очень большого числа Фибоначчи, но реализация предполагала использование плавающей арифметики сверхвысокой разрядности (с тем, чтобы нужное число полностью уместилось в мантиссу).

Мы пойдем совсем другим путем!

Для начала, рассмотрим множество чисел вида:

при целых a и b. Легко убедиться в том, что это множество алгебраически замкнуто относительно операций обычного сложения и умножения:

Более того, умножение и сложение будут коммутативными, что тоже легко проверяется непосредственно. Кроме того, ноль и единица хорошо "вписываются" в рассматриваемое множество:

Вполне натурально реализуется и вычитание подобных чисел:

Можно определить и деление (разумеется, лишь в случае, когда делитель отличен от нуля). Результат деления можно определить как корень уравнения:

Пусть

Тогда предыдущее равество эквивалентно следующему:

Раскрывая произведение в левой части последнего равенства, получим систему линейных уравнений для определения неизвестных e и f:

Отсюда:

Главный определитель этой системы равен:

Поскольку a и b здесь целые числа, то значение определителя всегда отлично от нуля, а значит, система имеет единственное решение и деление определяется корректно. Впрочем, мы увлеклись. Деление нам не понадобится.

Мы пришли к тому, что рассматриваемое множество с операциям сложения и умножения образует кольцо [4].

А теперь - самое главное! Зачем нам нужен корень из пяти? Никто не мешает реализовать арифметику на множестве пар (a,b), в которой сложение, вычитание и умножение будет описываться формулами:

Таким образом, можно “благополучно забыть” про корень из пяти  и реализовать прямое вычисление по формуле Бине. Для того, чтобы значение оказалось целым (и корень из пяти сократился), нужно, чтобы числитель дроби формулы Бине представлял собой число вида:

которое мы отождествим с "обычным" числом

Деление этого "обычного" иррационального числа на корень из пяти и даст нам искомый целый результат. Естественно, что в действительности делить не требуется, достаточно вычислить (используя описанную выше арифметику пар) два бинома:

и

а потом произвести вычитание. Не представляет труда реализовать этот подход на любом языке программирования. Мы сделаем это на Питоне (привлекает неограниченная разрядность целых в этом языке).

def prod_pairs(a,b):
    return (a[0]*b[0]+5*a[1]*b[1],a[0]*b[1]+a[1]*b[0])
def sub_pairs(a,b):
    return (a[0]-b[0],a[1]-b[1])
def pow_pair(a,n):
    c=a
    for _ in range(n-1):
        c=prod_pairs(c,a)
    return c
def fib_bine(n):
    x1=pow_pair((1,1),n)
    x2=pow_pair((1,-1),n)
    z=sub_pairs(x1,x2)
    return z[1]//(2**n)

Комментарии излишни - все очень просто. Сразу же возникает вопрос, а можно ли ускорить этот код? Очевидно, что "узким местом" здесь являтся возведение пары в целую степень. Для ускорения этой операции имеется стандартный прием - быстрое возведение в степень (этим же приемом воспользовался и автор [2]). Идея ускорения состоит в том, что для вычисления xⁿ вычисляется цепочка x -> x² -> x⁴ ->…->x^2k до тех пор, пока 2^k<=n, а затем аналогичным образом вычисляется x^(n-2k).

Теперь реализуем быстрое возведения пары в целую степень:

def pow_pair(a,n):
if (n==1):
    return a
c=copy(a)
k=1
while k*2<=n:
    if k<=n:
       c=prod_pairs(c,c)
       k=k*2
    p=n-k
    if p>=1:
       tmp=pow_pair(a,p)
       return prod_pairs(tmp,c)
    else:
       return c

Использование этого приема позволяет вычислять числа Фибоначчи за время, близкое к логарифмическому по формуле Бине и без использования арифметики с плавающей точкой. Для сравнения производительности предлагаемого метода с методом, основанным на простой рекурсии, написан следующий простейший код:

def fib_ite(n):
    c,p=0,1
    for _ in range(n):
        c,p=c+p,c
    return c

И что же? Несмотря на очевидную простоту кода fib_ite, функция fib_bine показывает значительно лучшие результаты. Так, на компьютере автора четырехсоттысячное число Фибоначчи по описываемому алгоритму вычисляется примерно за 2 сек, а прямыми итерациями – за 27 сек. На прилагаемом рисунке показаны разультаты тестов:

По горизонтальной оси откладывается номер рассчитываемого числа Фибоначчи, по вертикальной – время в секундах.

Получается, что формула Бине вполне пригодна для быстрых и точных вычислений чисел Фибоначчи.

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

1. Д.Кнут Искусство программирования на ЭВМ, т.1, Основные алгоритмы. – М: Вильямс, - 2017. - 720 C.

2. N-е число Фибоначчи за O(log N) https://habr.com/ru/post/148336/

3. Расчет миллионного числа Фибоначчи https://habr.com/ru/company/skillfactory/blog/555914/

4. С.Ленг Алгебра. М.: Наука, - 1965. - 431 C.