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

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

${F}_{n}={F}_{n-1}+{F}_{n-2}$

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

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

F_n=(((1+√5)/2)^n-((1-√5)/2))^n)/(√5)

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

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

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

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

x=a+√5*b

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

(a+b√5)+(c+d√5)=((a+c)+(b+d) √5) (a+b√5)(c+d√5)=((ac+5bd)+(ad+bc) √5)

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

1≡(1+√5*0) 0≡(0+ √5*0)

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

(a+b√5)-(c+d√5)=((a-c)+(b-d) √5)

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

$(a+b\sqrt{5})x=(c+d\sqrt{5})$

Пусть

$x=e+f\sqrt{5}$

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

$(a+b\sqrt{5})(e+f\sqrt{5})=(c+d\sqrt{5})$

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

$(a+b\sqrt{5})(e+f\sqrt{5})=(ae+5bf)+(af+be)\sqrt{5}$

Отсюда:

$\left.\begin{matrix} ae+5bf & = & c\\ be+af & = & d \end{matrix}\right\}$

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

$\begin{vmatrix} a & 5b\\ b & a \end{vmatrix} = {a}^{2}-5{b}^{2}$

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

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

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

(a,b)+(c,d)=((a+c),(b+d)) (a,b)-(c,d)=((a-c),(b-d)) (a,b)*(c,d)=((ac+5bd),(ad+bc))

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

$0+r\sqrt{5}$

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

$r\sqrt{5}$

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

A=(1+√5)^n/2^n

и

B=(1-√5)^n/2^n

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

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]). Идея ускорения состоит в том, что для вычисления xn вычисляется цепочка x → x2 → x4 →…→x2k до тех пор, пока 2k<=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.