費波拉契數列的演算法分析

Python · 發表 2018-10-06 11:37:00

摘要：　　版權申明：本文為博主窗戶(Colin Cai)原創，歡迎轉帖。如要轉貼，必須註明原文網址　　http://www.cnblogs.com/Colin-Cai/p/9717119.html 　　作者：窗戶　　QQ/微信：6679072 　　E-mail：[email protected]...

　　版權申明：本文為博主窗戶(Colin Cai)原創，歡迎轉帖。如要轉貼，必須註明原文網址

　　http://www.cnblogs.com/Colin-Cai/p/9717119.html 

　　作者：窗戶

　　QQ/微信：6679072

　　E-mail：[email protected]

看過我其他一些文章的人，可能想象不出我會寫一篇關於費波拉契數列的文章。因為可能會感覺1,1,2,3…這樣一個數列能講出什麼高深的名堂?嗯，本篇文章的確是關於費氏數列，但我的目的還是為了說一些應該有95%以上程式員不明白的東西。如果能夠跟著我弄明白文中分析的手法，其好處是不言而喻的。請聽我細細道來。

費波拉契數列

什麼叫費波拉契數列呢？

數學家費波拉契在自己的著作中用兔子繁殖模型引入了這樣一個數列：1，1，2，3，5，8，13…

這個數列的第1項和第2項都為1，以後的項都是前面兩項之和。

寫成遞推公式如下：

f(n) = 1 n=1,2

f(n) = f(n-1)+f(n-2) n>2

通項公式

我們都知道對於費波拉契數列，隨著n趨向於無窮大，f(n+1)/f(n)是存在極限的，我們在此略去在不知道通項公式的情況下對上述比值存在極限的數學證明( 數學裡的一切並不是想當然 )。

假設有極限，我們知道

而f(n+2)=f(n+1)+f(n)，所以

從而

又因為這個極限大於0，解方程得

費波拉契數列不僅有遞推公式，也有通項公式。

既然上述極限存在，我們完全可以猜測它是多個等差數列的和，然後用待定係數法算出來。

通項公式如下：

樹遞迴

現在我們就開始本節的重點， 如何計算費波拉契數列的第n項 。

既然已經知道費波拉契數列的遞推公式，那麼很容易就給出一個遞迴函式的版本，因為涉及到大數，我們可以採用Python來描述，本文後續主要採用Python：

def f(n):
if n<3:
return 1
return f(n-1)+f(n-2)

JavaScript版本當然如下（只可惜不是自身支援大數，需要實現，否則我更願意用js描述）：

function f(n) {
if(n<3)
return 1;
return f(n-1)+f(n-2);
}

如果用Scheme或Common Lisp來表示除了字首表示式古怪一點，看起來也差別不大：

(define (f n)
 (if (< n 2) 1
(+ (f (- n 1)) (f (- n 2)))
 )
)

上面是Scheme的描述；Common Lisp只要把開頭的define改成defun。

我們用Python來計算一下數列的前40項：

print(map(f, range(1,41)))

執行好慢啊，我的機器上運行了1分多鐘才出來了結果

[1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811, 514229, 832040, 1346269, 2178309, 3524578, 5702887, 9227465, 14930352, 24157817, 39088169, 63245986, 102334155]

運算量如此之大，這是怎麼回事呢？

我們來看看機器是怎麼計算這段遞迴的，我們就以f(5)的計算為例子：

展開就是這麼一個樣子，樹中的每個節點都在計算過程中出現，樹的規模是指數級(f(6)比f(5)多了6個節點)，也就是運算時間是 指數級 。

這個太誇張了，難怪這麼慢！

帶快取的遞迴

我們發現上面計算f(5)的遞迴計算的樹裡，f(3)是重複展開計算了。從而推斷，之所以樹遞迴計算規模這麼大，原因就在於出現了大量的重複計算。

於是我們自然就想到了，如何避免計算過的值再次遞迴展開重複計算呢？換句話說，怎麼讓計算機知道誰計算過了誰沒有計算過？

其實我們只需要引入一段快取，就可以解決這個問題，程式每當要計算一個值時，先到快取裡查一下就可以避免重複展開了。

我們用Python生成一個數組來做這段cache，沒有計算過的是None。程式碼如下：

#http://cnblogs.com/colin-cai
def f2(n,cache):
if cache[n]!=None:
return cache[n]
ret = f2(n-1, cache) + f2(n-2, cache)
cache[n] = ret
return ret

def f(n):
if n<3:
return 1
cache = [0,1,1]+[None]*(n-1)
return f2(n, cache)

當然，上面快取剛產生的時候，第0個元為0其實沒有用到，算是浪費了，其實去掉這個0再做一點點簡單的修改就可以了，這個留給讀者吧。

我們還是運算一下數列的前40個元素：

print(map(f, range(1,41)))

結果出的非常之快。

迭代

試想一下，如果讓我們在黑板上寫出費波拉契數列的前40項，我們會怎麼做？

我們會先把數列的第1項和第2項都寫成1；

然後從第3項開始，每一項都用前兩項相加，算出一項往黑板上寫一項，一項一項往前推進，直到寫完。

每一項的產生在是相互關聯的，而我們之前用Python裡的map函式生成數列的前40項，過程中每次呼叫f都是孤立的。

原來，如果我們的目的是生成費波拉契數列的前n項，剛才寫黑板的演算法就已經非常棒。用Python描述一下如下：

#http://cnblogs.com/colin-cai
def list_f(n):
if n<3:
return [1]*n
a = [1,1]+[None]*(n-2)
for i in range(2,n):
a[i]=a[i-1]+a[i-2]
return a

這個for迴圈就是從第3項起不斷的算一項寫一項的過程。

我們在黑板上寫各項的過程中，每次寫新項都只關心它前面兩項，而在此之前的項根本不關心。我們要生成新的一項，只要每次保持住前面的兩項即可。

我們每次都只需要記得當前是第幾項，已寫數列的最後兩項是什麼，然後算出新項，然後再只需要記得當前的項數以及最後兩項……如此不斷反覆，直到寫完。

這給了我們一個啟示，我們其實是在不斷的推進狀態：

當前項數n	第n-1項	第n項
2	1	1
3	1	2
4	2	3
5	3	5
…	…	…

我們要計算費波拉契數列的第n項，就不斷的如此狀態轉換，直到當前項數達到我們的要求則停止計算。

於是就有了以下的Python迭代版本：

#http://cnblogs.com/colin-cai
def f(n):
if n<3:
return 1
stat = {"n":2, "f(n-1)":1, "f(n)":1};
while stat["n"]<n:
stat["f(n-1)"],stat["f(n)"] = stat["f(n)"],stat["f(n-1)"]+stat["f(n)"]
stat["n"] += 1
return stat["f(n)"]

這個看似有著不錯的效率了。

上述迭代的效率

我們試圖用上述迭代版本的f計算費波拉契的第1000000項，結果我的計算機花了半分鐘以上。

版本中使用了字典，可能效率低一些。於是改用陣列來代替字典如下：

#http://cnblogs.com/colin-cai
def f(n):
if n<3:
return 1
stat = [2,1,1]
while stat[0]<n:
stat[1],stat[2] = stat[2],stat[1]+stat[2]
stat[0] += 1
return stat[2]

甚至於，有人可能覺得完全基於狀態推演的迭代看起來不那麼親切，那麼咱們換個寫法：

#http://cnblogs.com/colin-cai
def f(n):
if n<3:
return 1
x,y = 1,1
for m in range(2,n):
x,y = y,x+y
return y

結果發現折騰來折騰去並沒有明顯的好轉。

這個迭代的演算法稍有演算法功底的人就可以看的出來這是一個線性時間複雜度的演算法。但是，我們之所以說是線性時間複雜度是建立在裡面的所有加法、比較、賦值操作都是常數級時間。

然而，此處認為加法是常數級時間未必合適。在迭代的過程中，我們的stat[1]、stat[2]的長度應該是在上升的，準確的說，長度應該是線性的，也就是O(n)。從而，加法也是O(n)級的。我們基於此，可以認為整個迭代是O(n ² )。於是，我們也就可以理解雖然1000000並不是大到離譜，但迭代計算時間卻如此之慢了。

最終演算法

我們回頭去看看費波拉契數列的通項公式，是可以由兩個等比數列合成。

關於求整數次冪顯然有快速的演算法，乘法的次數為對數級，這個我在之前好幾篇博文裡都有說到過，可以認為這個是基本演算法。

a ⁿ 是n個a相乘，平凡的演算法下我們要計算n-1個乘法。

但我們注意到如下的計算：

a×a=a ²

a ² ×a ² =a ⁴

a ⁴ ×a ⁴ =a ⁸

a ⁸ ×a ⁸ =a ¹⁶

a ¹⁶ ×a ¹⁶ =a ³²

a ³² ×a ³² =a ⁶⁴

…

以上短短的幾個乘法就得到了 a ² 、 a ⁴ 、 a ⁸ 、 a ¹⁶ 、 a ³² 、 a ⁶⁴ …這些指數都是 2的整數次冪 。

而我們所要算的冪的指數顯然可以表示為二進位制，從而表示為 1 、 2 、 4 、 8 、 16 、 32 、 64 …這些2的整數次冪的一部分之和。

例如我們要算a ⁵⁷ ，

57用二進位制表示為111001，也就是

57=1+8+16+32

所以a ⁵⁷ =a×a ⁸ ×a ¹⁶ ×a ³²

以上的分析表明，快速的求冪演算法是存在的。當然，我們甚至沒有必要把指數先展開成二進位制方式，而是選擇直接迭代就可以做出來。

其實，對於

b×a ⁿ 當然，這裡b不為0，n也為非負整數。

如果 n=0 ，那麼

b×a ⁿ = b

此外如果 n是偶數，那麼

b×a ⁿ = b×(a×a) ^n/2

剩下的情況， n是奇數，那麼

b×a ⁿ = (a×b)×(a×a) ^(n-1)/2

以上三條規則不斷迭代1×a ⁿ 自然可以得到最後的解。

以上的分析寫成Python如下：

#http://cnblogs.com/colin-cai
def exp(a,n):
def exp2(b,a,n):
if n==0:
return b;
elif n%2==0:
return exp2(b,a*a,n/2)
return exp2(a*b,a*a,(n-1)/2)
return exp2(1,a,n)

既然乘方存在快速演算法，那麼我們就有理由懷疑費波拉契數列求項也存在快速演算法。可惜，通項公式裡那兩個底數都是無理數，近似的算算沒大問題，但是不能稱之為真正解決的演算法。而把冪展開來來推演演算法實在太麻煩( 悄悄地說，這個真的可行 )，在這裡不作講解。

我們把之前迭代中每一次向後推一項，狀態的轉換稱之為 T變換 ，也就是

T: (a,b)->(b,a+b)

這是一個函式，

狀態我們就用一個元組(tuple)來表示，現在我們用Python改寫一下最後出現的那個迭代版本，結果如下：

#http://cnblogs.com/colin-cai
def f(n):
if n<3:
return 1
T = lambda (a,b):(b,a+b)
r = (1,1)
for m in range(2,n):
r = T(r)
return r[1]

順帶說一句，Python就不能學學js那樣，用箭頭來表示lambda?多簡潔啊。

我們可以構造一個運算元 mul_f 來計算兩個函式的積，然後通過 mul_f 運算元來構造 exp_f 運算元來計算函式的冪。

那麼我們求數列第2項之後的第n項就相當於T變換的n-2次冪作用於(1,1)。

於是我們可以根據這個改寫一下程式碼:

#http://cnblogs.com/colin-cai
def f(n):
def mul_f(f1,f2):
return lambda x:f1(f2(x))
def exp_f(func, n):
if n==0 :
return lambda x:x;
return mul_f(func, exp_f(func, n-1))
if n<3:
return 1
return exp_f(lambda (a,b):(b,a+b),n-2)((1,1))[1]

帶有運算元、lambda的程式碼對於不是很習慣於函數語言程式設計的可能看起來有一點難懂，不過這裡只要理解了函式的冪即可。

我們利用mul_f函式乘積運算元還可以模仿之前快速 求冪演算法 ，程式碼如下：

#http://cnblogs.com/colin-cai
def f(n):
def mul_f(f1,f2):
return lambda x:f1(f2(x))
def f2(T1,T2,n):
if n==0:
return T1
elif n%2==0:
return f2(T1,mul_f(T2,T2),n/2)
return f2(mul_f(T1,T2),mul_f(T2,T2),(n-1)/2)

if n<3:
return 1
return f2(lambda (a,b):(a,b), lambda (a,b):(b,a+b), n-2)((1,1))[1]

函式在這裡運算，程式碼有點難懂，再繼續忍耐忍耐吧，馬上就要OK了。

上面的程式碼裡面有一堆函式的運算，但是最終運算還是要落實到實際機器上，從而不僅不減少運算量，甚至還多了一點函式的運算。

說白了，那個mul_f實在太“虛”了。如果我們可以讓這裡的函式乘積變成真正的數字計算，那麼我們上述的思路就成立了。

我們嘗試著去手動計算T變換的冪：

T: (a,b)->(b,a+b)

T ² : (a,b)->(a+b,a+2b)

T ³ : (a,b)->(a+2b,2a+3b)

…

T的任何次冪，都是把(a,b)轉換為a、b的兩個線性組合組成的元組。

於是我們受到了啟發：從現在開始，我們完全可以用一個四元組 (m,n,p,q) 來表示T的任意次冪，也就是上述任何的mul_f的結果。(m,n,p,q)代表

(a,b)->(m*a+n*b,p*a+q*b)

函式的乘法mul_f我們需要仔細想想，寫出來似乎長了一點。其他地方簡單的改寫一下，最後函式求值函式採用lambda的編寫就很簡單了。最終程式碼如下：

#http://cnblogs.com/colin-cai
def f(n):
mul_f = lambda (m1,n1,p1,q1),(m2,n2,p2,q2) \
: (m2*m1+n2*p1, m2*n1+n2*q1, p2*m1+q2*p1, p2*n1+q2*q1)
def f2(T1,T2,n):
if n==0:
return T1
elif n%2==0:
return f2(T1,mul_f(T2,T2),n/2)
return f2(mul_f(T1,T2),mul_f(T2,T2),(n-1)/2)
if n<3:
return 1
return (lambda (m,n,p,q),(a,b) : p*a+q*b)\
(f2((1,0,0,1), (0,1,1,1), n-2),(1,1))

運算一下f(1000000)，在我的電腦不到5秒就算了出來。不過看來，計算真未必是Python的長項，它還是做粘合劑比較合適。

當然，最後還是再配上我喜歡的Scheme的實現版本：

(define (f n)
 (define mul
(lambda (a b)
(let ((s1 (list (car a) (caddr a)))
(s2 (list (cadr a) (cadddr a)))
(s3 (list (car b) (cadr b)))
(s4 (cddr b))
(c (lambda (x y) (apply + (map * x y)))))
(list
(c s1 s3)
(c s2 s3)
(c s1 s4)
(c s2 s4)))))
 (define (f2 T1 T2 n)
(cond
((zero? n) T1)
((even? n) (f2 T1 (mul T2 T2) (/ n 2)))
(else (f2 (mul T1 T2) (mul T2 T2) (/ (- n 1) 2)))))
(if (< n 3) 1
((lambda (T init) (apply + (map * (cddr T) init)))
(f2 '(1 0 0 1) '(0 1 1 1) (- n 2))
'(1 1))))

後記

本文我決定使用Python作為問題描述的主要語言，而不以我以往描述問題習慣使用的C、Scheme、bc，誰叫Python流行呢。

C語言作為“高階彙編”，描述演算法確實很吃力，本文中涉及到大數運算，C語言描述並不那麼直接；當然，C++也一樣不能直接描述大數運算，雖然C++有著很多的抽象能力，比如模板，但它身上有“高階彙編”的根，這點無論怎麼包裝都是瞞不住的。我是越來越不傾向於用這樣底層的語言來描述演算法。不過，從底層更加容易理解我這裡提到的大數加法的O(n)時間複雜度問題倒是真的。

最後那個對數級的快速演算法，我當初是先知道它的存在，然後硬是不看答案甚至不看提示自己推出來，推的過程乃至程式碼描述和本文有很大的不同。但在寫這篇文章的時候，推著推著，覺得文中這樣的推理過程更加合理，從而更能讓人接受，於是就按文所寫了。

費波拉契數列的演算法分析

通項公式

帶快取的遞迴

迭代

上述迭代的效率

後記

您可能也會喜歡…