一、遞歸法

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

程序（in Python）寫出來那也是相當(dāng)?shù)暮?jiǎn)潔直觀（為了跟后面的程序區(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)

這個(gè)算法的時(shí)間復(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ù)雜度要小的多。這個(gè)方法就是按照遞推方程，從n = 0和n = 1開始，逐個(gè)求出所有小于n的Fibonacci數(shù)，最后就可以算出F(n)。由于每次計(jì)算值需要用到前兩個(gè)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

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

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

三、矩陣法

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

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

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

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

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

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

在實(shí)現(xiàn)的時(shí)候，可以用循環(huán)代替遞歸實(shí)現(xiàn)這里的二分分治，好處是降低了空間復(fù)雜度（用遞歸的話，空間復(fù)雜度為O(log n)）。下面的Python程序直接利用的numpy庫中的矩陣乘法（當(dāng)然這個(gè)庫也實(shí)現(xiàn)了矩陣的冪運(yùn)算，我把它單獨(dú)寫出來是為了強(qiáng)調(diào)這里的分治算法）。另外如果不用第三方庫，我也給出了矩陣乘法的簡(jiǎn)單實(shí)現(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í)間（雖然這個(gè)常數(shù)會(huì)比較大），因此整個(gè)算法的時(shí)間復(fù)雜度是O(log n)，空間復(fù)雜度是O(1)。

四、運(yùn)行時(shí)間大比拼

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

用Python的timeit模塊對(duì)以上三個(gè)算法的運(yùn)行時(shí)間進(jìn)行了測(cè)量，記錄了每個(gè)算法對(duì)于不同的n的每千次運(yùn)算所消耗的時(shí)間（單位是秒），部分?jǐn)?shù)據(jù)記錄在fibonacci_data。利用Mathematica可以很方便地對(duì)這些數(shù)據(jù)進(jìn)行擬合，對(duì)于較小的n，用三個(gè)復(fù)雜度表達(dá)式分別去擬合，得到的效果都非常好。尤其值得注意的是，對(duì)于第一個(gè)算法，我用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時(shí)，三種算法的千次運(yùn)行耗時(shí)比較。其中紅色為O(1.618 ^ n)的遞歸法；藍(lán)色為O(n)的遞推法；綠色為O(log n)的矩陣法。散點(diǎn)為實(shí)際測(cè)量到的運(yùn)行時(shí)間，實(shí)線為擬合方程的曲線。

三種算法的運(yùn)行時(shí)間比較

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

三種算法的運(yùn)行時(shí)間比較（對(duì)數(shù)坐標(biāo)軸）

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

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

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

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

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

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

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

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

利用Mathematica對(duì)大n情況下這兩種算法每千次運(yùn)行時(shí)間進(jì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)時(shí)，兩種復(fù)雜度的對(duì)比情況：

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

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

矩陣法的運(yùn)行時(shí)間（更大的n）

六、更快的算法？

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

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

來看看它到底有多快：

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

好吧，這個(gè)問題留待以后慢慢研究。

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