費氏數列 O(LogN) 的解法

fcamel

Published in

fcamel的程式開發心得

5 min readOct 6, 2017

在網路上看到可以用矩陣算出第 N 個費氏數列，時間複雜度是 O(LogN)。於是自己推想了一下。在使用矩陣的提示下，想出解法不會太難。

題目定義

F(n) = F(n-1) + F(n-2), n > 1
F(0) = F(1) = 1
求 F(n) 的值

O(LogN) 的解

假設我們能用矩陣表示出 F(n)，即:

v1 = [F(1)]
     [F(0)]m 是 2x2 矩陣m * v1 = [F(2)]
         [F(1)]
m^n * v1 = [F(n+1)]
           [F(n)]

那麼，只要我們能快速求出 m^n，就能快速求出 F(n)。

另一方面，對任何數字或矩陣來說，可以用 divide-and-conquer 的方式用 O(LogN)的時間求出它的 n 次方。以矩陣 m 為例，作法如下:

G(n) = G(n/2) * G(n/2), n > 1 且 n 是偶數
G(n) = G(n/2) * G(n/2) * G(1), n > 1 且 n 是奇數
G(1) = m

求 G(n) 的虛擬碼如下:

M G(n) {
  if (n <= 1)
    return m;  M t = G(n / 2);
  t = t * t;  if (n % 2)
    t *= m;
  return t;
}

費氏數列的 O(LogN) 解

於是，只要我們能找出表示 F(n) 的矩陣 m, 套入 O(LogN) 的求矩陣 n 次方的解，自然就有 O(LogN) 的費氏數列解。

稍微推敲一下即可得知:

m = [1 1]
    [1 0]v1 = [F(1)]
     [F(0)]m * v1 = [F(1) + F(0)] = [F(2)]
         [F(1)       ]   [F(1)]

推廣公式

更進一步來看，對於 H(n) 來說，只要 H(n) 能用 H(n-1), H(n-2), …, H(n-k) 的線性組合表示，就可以將 H(n) 轉成用矩陣表示的式子，從而得到 O(LogN) 的解。

太多代數看起來太抽象了，這裡以 k = 3 為例:

H(n) = a*H(n-1) + b*H(n-2) + c*(H-3)
H(0), H(1), H(2) 為某個常數

於是我們可以這樣表示 H(n):

v1 = [H(2)]
     [H(1)]
     [H(0)]m = [a b c]
    [1 0 0]
    [0 1 0]m * v1 = [a*H(2) + b*H(1) + c*H(0)] = [H(3)]
         [H(2)                    ]   [H(2)]
         [H(1)                    ]   [H(1)]

一樣，用 divide-and-conquer 的方式求出 m ^ (n — 2)，再乘上 v1，就能得到 H(n)，時間複雜度是 O(LogN)。

時間複雜度裡常數的影響

雖說矩陣乘法的版本是 O(LogN)，但矩陣乘法的計算量比一次加法多。推測 N 不大的時候，反而是 O(N) 的迴圈解比較快。我實作了幾個版本作比較，程式接收參數 n, k，會計算 F(n) k 次，然後輸出結果和花費的秒數。

$ clang++ -std=c++11 fib.cpp -o fibO0 -O0
$ clang++ -std=c++11 fib.cpp -o fibO2 -O2
$ ./fibO0 70 1000000
fib           (loop)   308061521170129 0.161764
fib by matrix (linear) 308061521170129 0.298222
fib by matrix (log(N)) 308061521170129 0.087025
$ ./fibO2 70 1000000
fib           (loop)   308061521170129 0.011287
fib by matrix (linear) 308061521170129 0.010232
fib by matrix (log(N)) 308061521170129 0.027675

fib (loop) 是一般的迴圈版。
fib by matrix (linear) 是用矩陣的迴圈版。
fib by matrix (log(N)) 是上面說的 O(LogN) 求矩陣 N 次方的版本。
fibO0 是用 O0 編譯，也就是 compiler 沒最佳化的版本。
fiO2 是用 O2 編譯，也就是 compiler 有最佳化的版本。

雖然矩陣乘法的版本用 O0 編譯的版本計算 F(70) 的時候比較快，但改用 O2 編的時候，反而比較慢了。觀察組譯後程式得知，fib (loop) O0 的版本在迴圈計算的過程裡，反覆使用 stack 上的區域變數，而 fib (loop) O2 的版本都用暫存器。這可以解釋為什麼 O2 快上十倍。至於為何矩陣乘法 log(N) 的版本較慢，但 linear 版卻和 loop 版一樣快，有機會再來研究看看。

$ ./fibO0 1000 1000000
fib           (loop)   9079565065540428013 2.37288
fib by matrix (linear) 9079565065540428013 4.1127
fib by matrix (log(N)) 9079565065540428013 0.141836
$ ./fibO2 1000 1000000
fib           (loop)   9079565065540428013 0.248233
fib by matrix (linear) 9079565065540428013 0.246838
fib by matrix (log(N)) 9079565065540428013 0.04743

F(1000) 已超出 long long 的最大值，所以計算結果不對。不過可藉此看出 O(LogN) 的效率差別，即使是 O2 編的版本，也是快上不少。