Fibonacci number

페이스북의 파이썬 코리아 커뮤니티에 어떤 분이 피보나치 수열 100만번째 항을 구하는 중이라는 포스팅을 하셨다. 댓글이 꽤 많이 달리고, 중간에 피보나치 수열의 일반항이나 행렬 계산법 등에 대한 언급도 있었다. 피보나치 수를 구하는건 나도 꽤 좋아하는 문제라서 개발자 면접에서 자주 꺼낸다. 간단하지만 프로그램이 실제로 수행될 때 어떤 일이 일어나는지를 모르면 공간 복잡도나 시간 복잡도에 대해 엉뚱한 이야기를 할 수도 있다. 답이 여러개가 나올 수 있는데다 이론을 넘어 실제로 결과를 보기 위해서는 아주 약간 귀찮은 작업이 필요하기 때문에 면접의 최초 스크리닝 문제로 꽤 괜찮다고 생각한다.


피보나치 수열

일단 피보나치 수에 대한 간단한 설명. 위키피디아에 잘 나와 있지만 …

  1. 처음 두 수는 1, 1 이다.
    • fibo(1) = 1
    • fibo(2) = 1
  2. 세번째 수부터는 앞의 두 수를 더한 값이 된다.
    • fibo(3) = fibo(2) + fibo(1)
    • fibo(4) = fibo(3) + fibo(2)
    • fibo(n) = fibo(n-1) + fibo(n-2) if n>= 3

위키피디아에는 0부터 시작한다고 되어 있는데 n=0 일때다. 사소한거니까 그냥 무시한다 … 위키피디아에 나온 피보나치 수열의 앞쪽 숫자들은

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946...

이렇다. 이후로 0은 무시한다.


반복문 구현

앞의 두 수를 계속 더해가면 된다. 간단하게 만들어 보았다.

In [1]:
def fibo_iterative(n):
    a, b = 1, 1
    for _ in xrange(n-1):
        a, b = b, a+b
    return a

이 코드로 n=1 부터 n=20 까지 피보나치 수를 구해보면 위키피디아에서 확인한 것 과 같은 수열이 나옴을 확인할 수 있다.

Out[2]:
n1234567891011121314151617181920
fibo_iterative(n)11235813213455891442333776109871597258441816765

재귀 구현

피보나치 수열은 반복문으로 구현할 수도 있지만, 점화식으로 쉽게 이해할 수 있기 때문에 재귀(recursion)를 배울 때에 예제로 자주 등장한다. 위에도 언급되었지만 점화식은,

$fibo(1)=1$

$fibo(2)=1$

$fibo(n) = fibo(n-1) + fibo(n-2)$ if $n>2$

이렇게 정리된다. 그대로 코드로 옮겨보자.

In [3]:
def fibo_recursive(n):
    if n == 1 or n == 2:
        return 1
    return fibo_recursive(n-1) + fibo_recursive(n-2)

위의 fibo_iterative 와 같은 결과가 나오는지 확인해보자.

In [4]:
for i in xrange(1, 10):
    print(fibo_iterative(i) == fibo_recursive(i))
True
True
True
True
True
True
True
True
True

일단 두 구현은 같은 값을 내보내는 것 같다.


복잡도

결론부터 얘기하면 n이 커짐에 따라 위 두 구현이 사용하는 메모리 양과 수행시간은,

  1. fibo_iterative: O(1) 메모리, O(n) 시간
  2. fibo_recursive: O(n) 메모리, O(exp(n)) 시간

이다. 재귀 구현이 메모리도, 시간도 더 많이 쓴다. 메모리는 귀찮으니 패스하고, 시간이 어떻게 되나 눈으로 확인해보자. 함수의 수행시간을 측정하기 위한 함수를 만들었다. 아래 elapsed 함수는 계산 결과는 버리고 수행시간을 초 단위로 리턴한다.

In [5]:
def elapsed(f):
    import time
    def _f(*args, **kwargs):
        st = time.time()
        ret = f(*args, **kwargs)
        return time.time() - st
    return _f
In [6]:
fibo_iter_elapsed = elapsed(fibo_iterative)
fibo_recur_elapsed = elapsed(fibo_recursive)

일단 fibo_iterative 의 n을 1부터 1000까지 늘려가며 시간을 측정.

In [7]:
n_x = range(1, 1000+1)
elapsed_iter = map(fibo_iter_elapsed, n_x)
plot(n_x, elapsed_iter)
Out[7]:
[<matplotlib.lines.Line2D at 0x7f14c7892c10>]

전체적으로 선형 관계를 따라가지만 종종 튀는 것들이 보인다. 아마도 캐시미스나 gc 등의 영향인 듯 싶다. 그러면 fibo_recursive 는 어떨까? 1부터 1000까지로 했더니 도저히 끝나질 않아서 1부터 30까지만 해 보았다.

In [8]:
elapsed_recur = map(fibo_recur_elapsed, n_x[:30])
plot(n_x[:30], elapsed_recur)
Out[8]:
[<matplotlib.lines.Line2D at 0x7f14c78ac0d0>]

그냥 딱 보기에도 시간이 폭발하고 있다. 두개를 겹쳐 그리면 아래와 같다 … n 이 더 커지면 스케일 문제로 아래쪽의 fibo_iterative의 결과가 보이지 않게 되어 n=15 까지만 그렸다.

In [9]:
plot(n_x[:15], elapsed_iter[:15], label="fibo_iterative")
plot(n_x[:15], elapsed_recur[:15], label="fibo_recursive")
ylabel("elapsed")
legend()
Out[9]:
<matplotlib.legend.Legend at 0x7f14c6e06650>

재귀 최적화

재귀가 오래 걸리는 것은 동일한 함수에 대해 중복 계산이 많기 때문이다. 위키피디아에서 n=5 인 경우에 대한 콜 트리를 가져와 보았다.

n=1 인 경우가 5번이나 계산되고 있다. 이런 중복 계산을 다시 하지 않고 값을 저장해 두었다가 사용할 수 있다. 이를 메모이제이션이라고 한다. 위의 fibo_recursive의 값을 저장했다 사용해 보자. 일반적으로 완벽하게 하려면 매우 귀찮지만 이 경우에는 인자의 첫 번째 값(n)만 알고 있으면 저장해야 하는 값/가져와야 하는 값이 무엇인지 특정할 수 있다. memoizedfibo_recursive에게 기억력을 주는 오즈의 마법사같은 함수다.

In [10]:
def memoized(f):
    cache = {}
    def _f(*args, **kwargs):
        key = args[0]
        cached = cache.get(key)
        if cached is not None:
            return cached
        ret = f(*args, **kwargs)
        cache[key] = ret
        return ret
    return _f

fibo_recursive = memoized(fibo_recursive)

memoized로 변경시킨 fibo_recursivefibo_iterative 와 똑같은 값을 내는지 확인해보자.

In [11]:
iter_result = map(fibo_iterative, range(1, 20+1))
memo_result = map(fibo_recursive, range(1, 20+1))
make_table([
    ["n"] + range(1, 20+1),
    ["fibo_iterative"] + iter_result,
    ["fibo_memoized"] + memo_result
])
Out[11]:
n1234567891011121314151617181920
fibo_iterative11235813213455891442333776109871597258441816765
fibo_memoized11235813213455891442333776109871597258441816765

굿잡!

그럼 수행 시간은 어떨까? 근성있게 fibo_iterative와 비슷한 O(n) 수행시간을 가질거라고 믿고 1~1000 을 모두 돌려서 fibo_iterative와 비교해 보겠다.

In [12]:
n_x = range(1, 1000+1)
fibo_memo_elapsed = elapsed(fibo_recursive)
elapsed_memoized = map(fibo_memo_elapsed, n_x)
elapsed_iter = map(fibo_iter_elapsed, n_x)

plot(n_x, elapsed_memoized, label="memoized")
plot(n_x, elapsed_iter, label="iterative")
legend()
Out[12]:
<matplotlib.legend.Legend at 0x7f14c6d74e10>

어라? 왜 상수시간 O(1) 처럼 보이지? 비밀은 테스트를 위해서 n=1 부터 n=1000 까지 계속 함수 호출을 하는데, 메모이제이션을 위한 캐시가 유지되기 때문에 이전 함수 호출의 결과를 다음 함수 호출에서 사용하기 때문이다. 피보나치 수를 리턴하는 함수는 입력이 같으면 언제나 같은 값을 내보내는 참조 투명성이 있기 때문에, 캐시가 함수 호출 사이에 공유되어도 상관 없다. 매 함수 호출마다 캐시를 리셋하면 fibo_iterative와 비슷한 모양이 나올 것으로 예상된다. 해보자.

아름답진 않지만 전역 공간에 cache 를 두고 동작하는 코드를 만들어 보겠다. 그리고 이건 재귀로 동작하기 때문에 recursion depth 를 높혀둘 필요가 있으므로 좀 높혀 두겠다.

In [13]:
import sys
sys.setrecursionlimit(8192)
In [14]:
n_x = range(1, 200+1)
memo_cache = {}

@elapsed
def fibo_memo(n):
    ret = memo_cache.get(n)
    if ret is not None:
        return ret
    if n == 1 or n == 2:
        return 1
    memo_cache[n] = fibo_memo(n-1) + fibo_memo(n-2)
    return memo_cache[n]

elapsed_memo_reset = []
for i in n_x:
    memo_cache = {}  # 캐시 리셋
    elapsed_memo_reset.append(fibo_memo(i))
    
elapsed_iter = map(fibo_iter_elapsed, n_x)
plot(n_x, elapsed_memo_reset, color="green", label="memo_reset")
plot(n_x, elapsed_iter, color="red", label="iterative")
legend()
Out[14]:
<matplotlib.legend.Legend at 0x7f14c6c6bed0>

기울기는 많이( …) 차이가 나지만 잘 보면 둘다 선형이다. memo_reset 이 같은 O(n) 임에도 성능이 훨씬 안나오는 이유는 아마도 캐시 참조와 함수 호출 오버헤드 때문일 것 같다. 하지만 위에도 말했듯 참조 투명성이 있으니 굳이 캐시를 지울 필요는 없다.


행렬 표현

위의 방법은 반복문이건, 메모이제이션이건 n번째 수를 처음 구할 때에는 O(n) 만큼의 스텝이 필요하다. 원 포스팅의 페이스북 댓글에도 언급이 있고, SICP 에도 나오는 예인데 n 번째 피보나치 수를 선형보다 빨리 구하는 방법이 있다. 피보나치 수의 일반항을 식으로 써낼 수 있는데, 이건 일단 귀찮으니 스킵. (SICP보세요) 행렬 계산으로 바꾸는 방식을 보자. 일단 fibo_iteration 의 식을 보면,

a, b = b, a+b

다시 잘 보면, 매 반복마다

$$a’ = 0 \times a + 1 \times b$$ $$b’ = 1 \times a + 1 \times b$$

이렇게 새 값을 만들고 있다. 이는 원래의 a, b 값에 2x2 행렬

[ [0, 1]
  [1, 1] ]

을 곱하는 연산이다. 즉, 다음 피보나치 수를 구하는 연산은 위 행렬을 한번 곱하는 것과 같다. 위 행렬을 M 이라 하면 n번째 피보나치 수를 구하기 위해서는 $M^{n}$ 을 구하면 된다. 이제 문제는 행렬의 곱하기를 빨리 하는 것으로 바뀌었다. 일단 이 방법이 맞는지 검증해 보자.

행렬을 latex 플러그인으로 썼더니 nbconvert 가 깨먹는다!!!!
Out[33]:
n12345678910
M^n[[0,1] [1,1]][[1,1] [1,2]][[1,2] [2,3]][[2,3] [3,5]][[3,5] [5,8]][[5,,8] [8,13]][[8,13] [13,21]][[13,21] [21,34]][[21,34] [34,55]][[34,55] [55,89]]

잘 보면 $M^{n}$ 의 [0, 1] 엘리먼트가 n번째 피보나치 수임을 알 수 있다.

In [16]:
table = [["n"], ["(M^n)[0,1]"]]
for i in xrange(1, 20+1):
    table[0].append(i)
    table[1].append((M**i)[0,1])
make_table(table)
Out[16]:
n1234567891011121314151617181920
(M^n)[0,1]11235813213455891442333776109871597258441816765

야호!

행렬의 n제곱을 구하기 위해서는 실제로 n번을 곱할 필요가 없다. 쉬운 예를 들면 $M^{16}$을 구하기 위해서는,

$M$ -> $M^2$ -> $M^4$ -> $M^8$ -> $M^{16}$

이렇게 4번의 곱셈만 하면 된다. 이를 SICP 에서는 SuccessiveSquaring이라고 설명하고 있다. 대략적으로 아주 큰 n 이 있다고 하면 O(logn) 만큼의 행렬 곱셈이 필요하다. 한번 해보자. 위키피디아 페이지에 있는 수도코드를 그대로 파이썬 코드로 옮겨 보았다.

… 사실 numpy의 ** 연산자에서 잘 해주고 있긴 할텐데, 뒤에 아마(?)필요할테니 해 두겠다.

In [17]:
def squaring(matrix, n):
    if n == 0:
        return np.matrix([[1, 0], [0, 1]])
    if n == 1:
        return matrix
    half = n/2
    squared = matrix*matrix
    if n%2 == 0:  # even
        return squaring(squared, half)
    else:  # odd
        return matrix*squaring(squared, half)
    
table = [["n"], ["M**n"], ["squaring(M,n)"]]
for i in xrange(1, 20+1):
    table[0].append(i)
    table[1].append((M**i)[0,1])
    table[2].append(squaring(M, i)[0,1])
make_table(table)
Out[17]:
n1234567891011121314151617181920
M**n11235813213455891442333776109871597258441816765
squaring(M,n)11235813213455891442333776109871597258441816765

야호!

과연 기대대로 빠를까? fibo_matrix(=squaring사용)과 fibo_matrix2(=numpy ** 사용)를 만들어서 fibo_iterative와 비교해 보자.

In [18]:
def fibo_matrix(n):
    return squaring(M, n)[0,1]

def fibo_matrix2(n):
    return (M**n)[0,1]

n_x = range(1, 1000+1)
fibo_matrix_elapsed = elapsed(fibo_matrix)
fibo_matrix2_elapsed = elapsed(fibo_matrix2)
elapsed_matrix = map(fibo_matrix_elapsed, n_x)
elapsed_matrix2 = map(fibo_matrix2_elapsed, n_x)
elapsed_iter = map(fibo_iter_elapsed, n_x)
plot(n_x, elapsed_matrix, label="squaring")
plot(n_x, elapsed_iter, label="iterative")
plot(n_x, elapsed_matrix2, label="np **")
legend()
Out[18]:
<matplotlib.legend.Legend at 0x7f14c6b99d10>

행렬 계산은 확실히 로그같은 모양을 보여주고 있다. 처음엔 느리다가 대충 n=500 근처에서 반복문을 사용한 계산보다 빨라지고 있다.

그런데, 이 값이 과연 큰 수에 대해서도 올바를까? 사용한 numpy 라이브러리는 C 타입을 사용하기 때문에 어느 정도 큰 수가 나오면 오버플로우가 일어날 가능성이 있다. 한번 값을 보자. 위의 값 확인은 n=20 까지만 했는데 1000까지 모두 같은 값이 나오나 확인해보자.

In [19]:
for i in range(1, 1000):
    a = fibo_iterative(i)
    b = fibo_matrix(i)
    if a != b:
        print("fibo_iterative and fibo_matrix generate different results at n=%d: %s %s" % (i, a, b))
        break
else:
    print("Looks good")
fibo_iterative and fibo_matrix generate different results at n=93: 12200160415121876738 -6246583658587674878

n=93 에서 행렬 계산의 값이 음수가 나왔다! ㅠㅠ

어디서부터 잘못되었을까?

In [20]:
f92 = fibo_iterative(92)
f93 = fibo_iterative(93)
int64_max = 0x7fffffffffffffff
print f92 < int64_max
print f93 < int64_max
True
False

결과가 signed int64 의 최대값을 넘어가는 순간 음수로 바뀐걸로 보아 integer overflow 문제가 맞다.

이를 극복하기 위해서 그냥 파이썬의 정수계산의 힘을 빌려보겠다. 순수한 파이썬 코드로 2x2 행렬과 곱셈을 구현한다.


행렬 곱셈 직접 하기

위의 numpy 행렬과 구분하기 위해 내가 만든 피보나치 계산 행렬은 MM 이라고 하겠다. 2x2행렬 곱하기가 헷갈리면 여기를 참조하면 된다.

In [21]:
class Matrix22:
    def __init__(self, a, b, c, d):
        """ [[a, b]
             [c, d]] """
        self.a = a
        self.b = b
        self.c = c
        self.d = d
        
    def __mul__(self, that):
        a = self.a * that.a + self.b * that.c
        b = self.a * that.b + self.b * that.d
        c = self.c * that.a + self.d * that.c
        d = self.c * that.b + self.d * that.d
        return Matrix22(a, b, c, d)
        
    def __repr__(self):
        return "[[%s %s]\n [%s %s]]" % (self.a, self.b, self.c, self.d)

MM = Matrix22(0, 1, 1, 1)
In [22]:
def fibo_my_matrix(n):
    return squaring(MM, n).b

이제 fibo_iterativefibo_my_matrix 의 값을 비교해 보자.

In [23]:
for i in range(1, 1000):
    a = fibo_iterative(i)
    b = fibo_my_matrix(i)
    if a != b:
        print("fibo_iterative and fibo_my_matrix generate different results at n=%d: %s %s" % (i, a, b))
        break
else:
    print("Looks good")
Looks good

n=93 에서의 오버플로는 해결되었다. 이제 수행시간을 확인한다.

In [24]:
n_x = range(1, 1000+1)

fibo_my_matrix_elapsed = elapsed(fibo_my_matrix)

elapsed_matrix = map(fibo_matrix_elapsed, n_x)
elapsed_my_matrix = map(fibo_my_matrix_elapsed, n_x)
elapsed_iter = map(fibo_iter_elapsed, n_x)

plot(n_x, elapsed_matrix, label="np matrix")
plot(n_x, elapsed_my_matrix, label="my matrix")
plot(n_x, elapsed_iter, label="iteration")
legend()
Out[24]:
<matplotlib.legend.Legend at 0x7f14c6af87d0>

이제 로그 시간에 n번째 피보나치 수를 가져오는 함수를 손에 넣었다. 계산 가능한 수의 크기는 (아마도) 파이썬의 정수 계산 한계에 바운딩될텐데, 이건 메모리 크기와 관련있다고 알고 있다.


N을 키워본다

자, 이제 나도 100만번째, 1000만번째, 그리고 궁금했던 1억번째 피보나치 수를 구해보자.

In [25]:
a = fibo_iterative(100*10000) 
b = fibo_my_matrix(100*10000)
assert(a==b)

print "100만번째 피보나치 수: %s ... ... ... %s" % (str(a)[:32], str(a)[-32:])
print "Elapsed(iter): %.3f" % fibo_iter_elapsed(100*10000)
print "Elapsed(my_matrix): %.3f" % fibo_my_matrix_elapsed(100*10000)
100만번째 피보나치 수: 19532821287077577316320149475962 ... ... ... 01719893411568996526838242546875
Elapsed(iter): 9.942
Elapsed(my_matrix): 0.304

100만번째 피보나치 수를 구하는데 fibo_iterativefibo_my_matrix 의 수행 시간 차이는 약 30배. 수가 커지면 이 차이는 더 벌어진다.

ipython notebook 에서 실행하긴 여의치 않아서 + 너무 길어서 천만번째와 1억번째 수는 터미널에서…

1억번째 피보나치 수는 계산을 마치고 elapsed time은 떳는데, 아직 디스크에 쓰기 위한 직렬화(Serialization)을 하느라 CPU를 태우고 있는 것 같다. 큰 숫자를 문자열로 바꾸는게 생각보다 느린 작업인가보다. 현재 링크는 깨져 있고, 답이 나오면 업데이트 될 예정… :’(

Comments