Python: 자크 비네의 피보나치 수열 방정식

Jacques Philippe Marie Binet

 

1. 소개

자크 필리프 마리 비네의 n번째 피보나치 수열 생성 방정식은 아래와 같다.

선형 동차 점화식을 해결하는 과정에서 특성 방정식을 풀고 그 해를 사용해 일반 해를 구하면 이 식이 도출되는 듯하다.

루트를 사용하기 때문에 72번째부터는 값이 다르다.

그래서 O(logn)이면서 큰 수를 구할 때는 행렬 거듭제곱을 사용할 수 있겠

 

(n0 = 0, n1 = 1)

@wikipedia

 

2. 파이썬

from math import sqrt


def fibo_binet(n):
    sqrt5 = sqrt(5)
    return int((((1 + sqrt5) ** n - (1 - sqrt5) ** n) / (2 ** n * sqrt5)))

 

꼼수로 250번 째 까지 작동하게 하려면 precision을 바꿔야 한다.

from decimal import Decimal, getcontext

# the larger n, the higher precision required
getcontext().prec = 75

def fibo_binet(n):
    sqrt5 = Decimal(5).sqrt()
    phi = (1 + sqrt5) / 2
    psi = (1 - sqrt5) / 2
    return int((phi ** n - psi ** n) / sqrt5)

# Example
print(fibo_binet(250))

 

피보나치 최적화 버전 vs 비넷 버전

비넷 버전이 조금 더 빠르다.

from math import sqrt
from timeit import default_timer


def fibo_python(n):
    counter = 0
    first, second = 0, 1
    while counter < n:
        first, second = second, first + second
        counter += 1
    return first


def fibo_binet(n):
    sqrt5 = sqrt(5)
    return int((((1 + sqrt5) ** n - (1 - sqrt5) ** n) / (2 ** n * sqrt5)))

start = default_timer()
for i in range(2, 71):
    fibo_python(i)  # 최적화 버전
print(f"{default_timer() - start:.5f}s")

start = default_timer()
for i in range(2, 71):
    fibo_binet(71)  # 비넷 버전
print(f"{default_timer() - start:.5f}s")


""" 결과

0.00015s
0.00010s

"""

 

3. 번외

피보나치 수열을 계산하는 가장 빠른 수학적 방법은 행렬 거듭제곱을 사용하는 것이다.

이 방법은 \(O(\log n)\)의 시간 복잡도를 가진다.

아래의 행렬을 사용한다

\[
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\]

이 행렬을 \(n-1\)번 거듭제곱하면, 그 결과 행렬의 왼쪽 상단 원소가 \(F(n)\)이 되는데

수학적으로 표현하면, 행렬 \( A \)를

\[
A = \begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\]

라고 할 때, \( A^{(n-1)} \)의 왼쪽 상단 원소가 \( F(n) \)이다.

from math import sqrt

def fibo_binet(n):
    sqrt5 = sqrt(5)
    return int((((1 + sqrt5) ** n - (1 - sqrt5) ** n) / (2 ** n * sqrt5)))

def mat_mult(a, b):
    return [
        [a[0][0]*b[0][0] + a[0][1]*b[1][0], a[0][0]*b[0][1] + a[0][1]*b[1][1]],
        [a[1][0]*b[0][0] + a[1][1]*b[1][0], a[1][0]*b[0][1] + a[1][1]*b[1][1]]
    ]

def mat_pow(mat, n):
    if n == 1:
        return mat
    if n % 2 == 0:
        half_pow = mat_pow(mat, n // 2)
        return mat_mult(half_pow, half_pow)
    else:
        return mat_mult(mat, mat_pow(mat, n - 1))

def fast_fibo(n):
    if n == 0:
        return 0
    if n == 1:
        return 1
    
    initial_mat = [[1, 1], [1, 0]]
    result_mat = mat_pow(initial_mat, n - 1)
    return result_mat[0][0]

# Test the function
print(fast_fibo(10), fast_fibo(20), fast_fibo(100))
print(fibo_binet(10), fibo_binet(20), fibo_binet(30))