一、遞歸法

一個看起來很直觀、用起來很恐怖的算法就是遞歸法。根據(jù)Fibonacci的遞推公式，對于輸入的n，直接遞歸地調(diào)用相同的函數(shù)分別求出F(n - 1)和F(n - 2)，二者相加就是結(jié)果。遞歸的終止點就是遞推方程的初值，即n取0或1的時候。

程序（in Python）寫出來那也是相當(dāng)?shù)暮啙嵵庇^（為了跟后面的程序區(qū)分開來，這里取名SlowFibonacci）。

1
2
3
4

def SlowFibonacci(n):
  assert n >= 0, 'invalid n'
  if n < 2: return n  # F(0) = 0, F(1) = 1
  return SlowFibonacci(n - 1) + SlowFibonacci(n - 2)

這個算法的時間復(fù)雜度有著跟Fibonacci類似的遞推方程：T(n) = T(n - 1) + T(n - 2) + O(1)，很容易得到T(n) = O(1.618 ^ n)（1.618就是黃金分割，(1+√5)/2 ）?？臻g復(fù)雜度取決于遞歸的深度，顯然是O(n)。

二、遞推法

雖然只是一字之差，但遞推法的復(fù)雜度要小的多。這個方法就是按照遞推方程，從n = 0和n = 1開始，逐個求出所有小于n的Fibonacci數(shù)，最后就可以算出F(n)。由于每次計算值需要用到前兩個Fibonacci數(shù)，更小的數(shù)就可以丟棄了，可以將空間復(fù)雜度降到最低。算法如下：

1
2
3
4
5
6
7

def NormFibonacci(n):
  assert n >= 0, 'invalid n'
  if n == 0: return 0
  (prev, curr) = (0, 1)  # F(0), F(1)
  for i in xrange(n - 1):
    (prev, curr) = (curr, prev + curr)
  return curr

顯然時間復(fù)雜度是O(n)，空間復(fù)雜度是O(1)。

比較一下遞歸法和遞推法，二者都用了分治的思想——把目標問題拆為若干個小問題，利用小問題的解得到目標問題的解。二者的區(qū)別實際上就是普通分治算法和動態(tài)規(guī)劃的區(qū)別。

三、矩陣法

算Fibonacci數(shù)精確值的最快的方法應(yīng)該就是矩陣法，看過的人都覺得這個方法很好。如果你跟我一樣，曾經(jīng)為記住這個方法中的矩陣而煩惱，那今天就來看看怎么進行推導(dǎo)。其實方法非常簡單，想清楚了也就自然而然地記住了。

我們把Fibonacci數(shù)列中相鄰的兩項：F(n)和F(n - 1)寫成一個2x1的矩陣，然后對其進行變形，看能得到什么：

是不是非常自然呢？把等式最右邊繼續(xù)算下去，最后得到：

因此要求F(n)，只要對這個二階方陣求n - 1次方，最后取結(jié)果方陣第一行第一列的數(shù)字就可以了。

看起來有點兒化簡為繁的感覺，但關(guān)鍵點在于，冪運算是可以二分加速的。設(shè)有一個方陣a，利用分治法求a的n次方，有：

可見復(fù)雜度滿足T(n) = T(n / 2) + O(1)，根據(jù)Master定理可得：T(n) = O(log n)。

在實現(xiàn)的時候，可以用循環(huán)代替遞歸實現(xiàn)這里的二分分治，好處是降低了空間復(fù)雜度（用遞歸的話，空間復(fù)雜度為O(log n)）。下面的Python程序直接利用的numpy庫中的矩陣乘法（當(dāng)然這個庫也實現(xiàn)了矩陣的冪運算，我把它單獨寫出來是為了強調(diào)這里的分治算法）。另外如果不用第三方庫，我也給出了矩陣乘法的簡單實現(xiàn)。

• Using numpy Library

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

from numpy import matrix

def MatrixPower(mat, n):
  assert n > 0, 'invalid n'
  res = None
  temp = mat
  while True:
    if n & 1:
      if res is None: res = temp
      else: res = res * temp
    n >>= 1
    if n == 0: break
    temp = temp * temp
  return res

def FastFibonacci(n):
  assert n >= 0, 'invalid n'
  if n < 2: return n  # F(0) = 0, F(1) = 1
  mat = matrix([[1, 1], [1, 0]], dtype=object)
  mat = MatrixPower(mat, n - 1)
  return mat[0, 0]

• Without numpy Library

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

def DotProduct(x, y):
  n = len(x)
  assert len(y) == n, 'x and y must have the same length'
  s = 0
  for i in xrange(n):
    s += x[i] * y[i]
  return s

def MatrixMultiply(x, y):
  # x is a m*a matrix, y is a a*n matrix.
  # x * y is a m*n matrix.
  m = len(x)
  n = len(y[0])
  a = len(x[0])
  assert len(y) == a

  # transpose y
  y = [[y[i][j] for i in xrange(a)] for j in xrange(n)]

  res = [[DotProduct(x[j], y[i]) for i in xrange(n)] for j in xrange(m)]
  return res

def MatrixPower(mat, n):
  assert n > 0, 'invalid n'
  res = None
  temp = mat
  while True:
    if n & 1:
      if res is None: res = temp
      else: res = MatrixMultiply(res, temp)
    n >>= 1
    if n == 0: break
    temp = MatrixMultiply(temp, temp)
  return res

def FastFibonacci(n):
  assert n >= 0, 'invalid n'
  if n < 2: return n  # F(0) = 0, F(1) = 1
  mat = [[1, 1], [1, 0]]
  mat = MatrixPower(mat, n - 1)
  return mat[0][0]

二階方陣相乘一次可以看成是常數(shù)時間（雖然這個常數(shù)會比較大），因此整個算法的時間復(fù)雜度是O(log n)，空間復(fù)雜度是O(1)。

四、運行時間大比拼

至此，我們得到的時間復(fù)雜度分別是O(1.618 ^ n)、O(n)和O(log n)的算法，讓我們來直觀地比較比較它們。

用Python的timeit模塊對以上三個算法的運行時間進行了測量，記錄了每個算法對于不同的n的每千次運算所消耗的時間（單位是秒），部分數(shù)據(jù)記錄在fibonacci_data。利用Mathematica可以很方便地對這些數(shù)據(jù)進行擬合，對于較小的n，用三個復(fù)雜度表達式分別去擬合，得到的效果都非常好。尤其值得注意的是，對于第一個算法，我用a * b ^ n去擬合，結(jié)果得到b等于1.61816，這與黃金分割數(shù)的正確值相差無幾。

• 遞歸法擬合結(jié)果：0.000501741 * 1.61816 ^ n，RSquare = 0.999993。
• 遞推法擬合結(jié)果：0.000788421 + 0.000115831 * n，RSquare = 0.999464。
• 矩陣法擬合結(jié)果：-0.0114923 + 0.0253609 log(n)，RSquare = 0.986576。

下圖是n <= 35時，三種算法的千次運行耗時比較。其中紅色為O(1.618 ^ n)的遞歸法；藍色為O(n)的遞推法；綠色為O(log n)的矩陣法。散點為實際測量到的運行時間，實線為擬合方程的曲線。

三種算法的運行時間比較

當(dāng)n > 10的時候，指數(shù)時間就已經(jīng)超出畫面范圍了。另外在這張圖里，身為對數(shù)時間復(fù)雜度的矩陣法似乎沒有任何優(yōu)勢，其耗時遠遠高于線性時間復(fù)雜度的遞推法。這是因為n還不夠大，體現(xiàn)不出log(n)的優(yōu)勢。在考慮更大的n之前，先來看看指數(shù)時間復(fù)雜度會增大到什么程度。

三種算法的運行時間比較（對數(shù)坐標軸）

五、大整數(shù)情況下的復(fù)雜度

Python內(nèi)置了大整數(shù)支持，因此上面的程序都可以直接接受任意大的n。當(dāng)整數(shù)在32位或64位以內(nèi)時，加法和乘法都是常數(shù)時間，但大整數(shù)情況下，這個時間就不能忽略了。

先來看一下Fibonacci數(shù)的二進制位數(shù)。我們知道Fibonacci數(shù)的通項公式是：

當(dāng)n充分大（其實都不需要很大）的時候，第二項就可以忽略不計了。把第一項對2取對數(shù)，就可以得到Fibonacci數(shù)的二進制位數(shù)的近似表達式，大概是log21.618×n0.5log25=log21.618×n1.161=O(n) 。由此可以算出，F(xiàn)(47)是32位無符號整數(shù)可以表達的最大的Fibonacci數(shù)，F(xiàn)(93)是64位無符號整數(shù)可以表達的最大的Fibonacci數(shù)。上面圖中的n在36以內(nèi)，不需要動用大整數(shù)運算，復(fù)雜度也比較符合之前的結(jié)論。但對于更大的n，之前的復(fù)雜度就不再適用了。

指數(shù)復(fù)雜度的算法就不管了，還不等用到大整數(shù)，它就已經(jīng)慢到不行了。

來看看O(n)時間復(fù)雜度的遞推法。每次遞推的時候都要計算兩個Fibonacci數(shù)之和，第i次運算時，這兩個Fibonacci數(shù)分別有O(i)個二進制位，完成加法需要O(i)的時間。因此總的時間大約是：

可見對于很大的n，遞推法的時間復(fù)雜度實際上是O(n ^ 2)的，空間復(fù)雜度是O(n)用來存儲Fibonacci數(shù)的各個二進制位。

再看矩陣法，注意到矩陣運算中有乘法，兩個長度為n的大整數(shù)相乘，傳統(tǒng)算法是O(n ^ 2)時間復(fù)雜度，較好的Karatsuba算法是O(n ^ (log 3 / log 2))時間，更快的快速傅立葉變換法是O(n log n)時間。Python 2.5中使用的是Karatsuba算法（Python 3里面似乎是快速傅立葉變換法）（參見Python源碼中的算法分析 之 大整數(shù)乘法）。以Karatsuba算法為例，矩陣法的時間復(fù)雜度遞推方程為：T(n)=T(n/2)+O(nlog23) ，應(yīng)用Master定理求得T(n)=O(nlog23) 。因此對于很大的n，矩陣法的時間復(fù)雜度為O(n ^ 1.585)，空間復(fù)雜度O(n)。

利用Mathematica對大n情況下這兩種算法每千次運行時間進行擬合，分別得到：

• 遞推法大整數(shù)擬合結(jié)果：0.0131216 + 0.000102101 * n + 2.44765 * 10 ^ -7 * n ^ 2，RSquare = 0.999482。
• 矩陣法大整數(shù)擬合結(jié)果：0.171487 + 9.74496 * 10 ^ -7 * n ^ 1.51827，RSquare = 0.998395。

看一下n在4000以內(nèi)時，兩種復(fù)雜度的對比情況：

遞推法（藍色）與矩陣法（綠色）運行時間比較（大整數(shù)）

從圖中可以看出，遞推法的增長速度也是很快的，當(dāng)n增大到60多的時候，它的運行時間就超過矩陣法了。矩陣法的增長速度非常慢，看起來像是線性的，讓我們把n調(diào)的更大來看一下。

矩陣法的運行時間（更大的n）

六、更快的算法？

試了試Mathematica中的Fibonacci函數(shù)，發(fā)現(xiàn)其運算速度相當(dāng)驚人，估計時間復(fù)雜度在O(n log n)上下，而且對于相同的n，運算速度遠遠高于我的矩陣法?？上疫€不了解它的算法，只是在幫助文檔里看到：

Fibonacci[n] uses an iterative method based on the binary digit sequence of n.

來看看它到底有多快：

矩陣法（綠色）與Mathematica Fibonacci函數(shù)（橙色）運行時間比較

好吧，這個問題留待以后慢慢研究。

最后相關(guān)的Mathematica命令文件放在這里：fibonacci_timecost