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

в 8:15, , рубрики: Алгоритмы, формула Бине, числа

Красивая формула Бине для чисел Фибоначчи содержит иррациональность - квадратный корень из пяти. Это делает ее непригодной для точного вычисления больших чисел Фибоначчи. Это кажется вполне очевидным. Предлагаю способ, как избавиться от зловредного корня и сделать формулу Бине пригодной для точных вычислений.

Интересно? Читайте дальше

Как хорошо известно, числа Фибоначчи – это целочисленная последовательность, первые два члена которой равны единице, а каждый последующий равен сумме двух предыдущих. За 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+bsqrt{5})x=(c+dsqrt{5})

Пусть

x=e+fsqrt{5}

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

(a+bsqrt{5})(e+fsqrt{5})=(c+dsqrt{5})

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

(a+bsqrt{5})(e+fsqrt{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+rsqrt{5}

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

rsqrt{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 сек. На прилагаемом рисунке показаны разультаты тестов:

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

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

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

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

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.

Автор: Борис Файфель

Источник


* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js