費波拉契數列的算法分析
版權申明：本文為博主窗戶(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(n2)。於是，我們也就可以理解雖然1000000並不是大到離譜，但叠代計算時間卻如此之慢了。
最終算法
我們回頭去看看費波拉契數列的通項公式，是可以由兩個等比數列合成。
關於求整數次冪顯然有快速的算法，乘法的次數為對數級，這個我在之前好幾篇博文裏都有說到過，可以認為這個是基本算法。
an是n個a相乘，平凡的算法下我們要計算n-1個乘法。
但我們註意到如下的計算：
a×a=a2
a2×a2=a4
a4×a4=a8
a8×a8=a16
a16×a16=a32
a32×a32=a64
…
以上短短的幾個乘法就得到了a2、a4、a8、a16、a32、a64 …這些指數都是2的整數次冪。
而我們所要算的冪的指數顯然可以表示為二進制，從而表示為1、2、4、8、16、32、64…這些2的整數次冪的一部分之和。
例如我們要算a57，
57用二進制表示為111001，也就是
57=1+8+16+32
所以a57=a×a8×a16×a32
以上的分析表明，快速的求冪算法是存在的。當然，我們甚至沒有必要把指數先展開成二進制方式，而是選擇直接叠代就可以做出來。
其實，對於
b×an 當然，這裏b不為0，n也為非負整數。
如果n=0，那麽
b×an = b
此外如果n是偶數，那麽
b×an = b×(a×a)n/2
剩下的情況，n是奇數，那麽
b×an = (a×b)×(a×a)(n-1)/2
以上三條規則不斷叠代1×an自然可以得到最後的解。
以上的分析寫成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)
T2: (a,b)->(a+b,a+2b)
T3: (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)時間復雜度問題倒是真的。
最後那個對數級的快速算法，我當初是先知道它的存在，然後硬是不看答案甚至不看提示自己推出來，推的過程乃至代碼描述和本文有很大的不同。但在寫這篇文章的時候，推著推著，覺得文中這樣的推理過程更加合理，從而更能讓人接受，於是就按文所寫了。
費波拉契數列的算法分析