費氏數怎麼算?

2019 年初,臉書上有位仁兄與所有人對計算費氏數的複雜度展開激烈辯論,並提出了應該把圓周率與自然對數等等重要常數存在月球(?!)上以便查詢的想法。於是我也跟風查了些關於費氏數的有趣資料。


原圖:sciencefreak @ pixabay

大家都知道如下的 Fibonacci 數定義(令 Nat 表示自然數的型別):

‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍fib :: Nat -> Nat
   fib 0     = 0
‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍fib 1     = 1
‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍fib (n+2) = fib (n+1) + fib n

在沒有其他機制(如 memoization)輔助的情況下直接跑上述定義,會得到一個很慢的演算法,因為許多計算重複了。早年的入門程式教科書常錯誤地把這當作「遞迴效率欠佳」的例子。如果定義 fib2 n = (fib n, fib (n+1)),我們不難推導出如下的(也是遞迴的)定義:

‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍fib2 :: Nat -> (Nat, Nat)
   fib2 0     = (0, 1)
‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍fib2 (n+1) = (y, x+y)
      where (x,y) = fib2 n 

這個演算法的效率改進不少 — 當 n 不太大時是很堪用的。計算 fib2 n 時只需 O(n) 個遞迴呼叫,但這表示這個演算法的時間複雜度是 O(n) 嗎?

如果還想更快呢?許多人會說,費氏數不是有個公式解嗎?

‍‍‍‍‍‍ ‍‍  fib n = (((1+√5)/2)^n - ((1-√5)/2)^n) /√5

這個「公式解」一般認為由 Jacques P. M. Binet在 1843 年發現,但事實上還可追溯得更早。習慣上我們仍稱之為 Binet 等式 (Binet’s Formula). 該公式的證明是常見的練習,我推薦大家看看這篇詳細的教材(CS 280 @ Cornell, by Joe Halpern, 2005)。許多人可能看到這是一個封閉式,便認為這是最快的算法,甚至有人以為這個式子的右手邊可在常數時間內算出來。但別忘了, ((1+√5)/2)^n((1-√5)/2)^n 可不總是一彈手指就算得出來的!

事實上,當數字一大,我們便不能假設加法、乘法等等基本運算只花常數時間。「怎麼算大數乘法」這個問題甚至到最近還有新進展 — 推薦大家一讀郭君逸教授最近在泛科學的文章

Vorobev 等式與 O(log n) 次遞迴呼叫的演算法

另有一系列演算法只需 O(log n) 次遞迴呼叫便可算出 fib n. 要理解它們,可由這個問題出發:fib (n+k)fib nfib k (以及它們附近的數)有關係嗎?如果可以找出些關聯,也許我們便可以由 fib n 算出 fib (n+n), 有希望設計出 O(log n) 次遞迴呼叫的演算法了。

確實,對 n >= 1,我們有這樣的性質:

‍‍‍‍‍‍‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍fib (n+k) = fib (n-1) * fib k + fib n * fib (k+1) .      -- (Vor)

這個性質目前可追溯到 Nikolai. N. Vorobev,姑且便稱之為「Vorobev 等式」。我們待會兒會試著證明,目前先看看如何用它。由 Vorobev 等式我們可導出幾種只需 O(log n) 次遞迴呼叫的演算法。例如,將 (Vor) 中的 n, k 代換為 n+1, n,我們可得到:

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍   fib (2n+1) = (fib (n+1))^2 + (fib n)^2                   -- (1)

將 (Vor) 中的 n, k 代換為 n+1, n+1 可得到:

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍   fib (2n+2) = 2 * fib n * fib (n+1) + (fib (n+1))^2       -- (2)

把 (2) 減去 (1), 我們得到

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍fib 2n = 2 * fib n * fib (n+1) - (fib n)^2               -- (3)

式子 (1) 和 (3) 之中,等號左邊的參數一個是奇數,一個是偶數,兩者的等號右邊都只需要 fib nfib (n+1)。若定義 fib2v n = (fib n, fib (n+1)),我們可得到如下的程式:

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍  fib2v :: Nat -> (Nat, Nat)
  fib2v 0 = (0, 1)
  fib2v n | n `mod` 2 == 0 = (c,d)
          | otherwise      = (d, c + d)
     where (a, b) = fib2v (div n 2)
           c      = 2 * a * b - a * a
           d      = a * a + b * b

到底哪種算法比較快?

這麼多演算法中,到底哪個比較快? 2019 年初李祐棠(Lee You Tang)親自做了許多實驗,相當有趣。其實,在 1988 年,已有位 James L Holloway 寫了一整篇碩士論文談費氏數的計算。Holloway 回顧了好幾個算費氏數的方法,包括有許多重複的遞迴法(相當於本篇的第一個 fib)、反覆加法(相當於 fib2)、Binet 等式、好幾種基於矩陣乘法(因乘法的遞移律,可以做 divide and conquer)的做法,生成函數… 有做理論分析,也有跑程式。

一些有趣的結果如下。

理論分析方面,我們得先知道 fib n 大約需要多少個 bit 來表示。令 N n 為存放 fib n 所需要的 bit 數,Holloway 推算出的近似值是 N n = n * 0.69424fib nN n 是線性關係,和我們對於 fib n 以指數速度增長的直覺符合。

演算法 fib2 只需要 O(n) 個遞迴呼叫。但這不表示計算時間也是 O(n). 以這個演算法算 fib n 需要的 bit 運算數目大約是 N (n^2/2 - n/2). 注意:此處說的已經不是 big O,而是實質的上限。

Binet 等式呢?√5可用牛頓法來算,假設每個 n bit 除法需要 n^2 個 bit 運算,但每一輪只需要留下前面 N n + log n 個 bit. 整體上說來需要的 bit 運算數中最大宗的大約是 log n * (N n + log n)^2, 再加上一些較不顯著的運算 — 其實並沒有很明顯地比 fib2 快。

幾個矩陣演算法需要的 bit 運算數目大都是 (N n)^2 乘以一個常數,常數大小各有不同。

基於 Vorobev 等式的 fib2v 表現得相當不錯,需要大約 1/2 * (N n)^2 個 bit 運算。

實測時的效率呢?Holloway 嘗試每個演算法,如果花時間超過五分鐘就不做了。一個測試中,Binet等式解在 log n = 7 的時候就超出了五分鐘。fib2 倒是可以撐到 log n = 15. 另一個測試(只考慮 n 為二的乘冪的較簡單演算法)中,Binet 等式解在 log n = 13 時超過時限。此時幾個矩陣演算法都只花了大約 3 – 5 秒的時間,fib2v 大約兩秒。

Vorobev 等式怎麼證?

最後我們來看看 Vorobev 等式怎麼證。我們在 n 上作歸納。當 n := 12 的 case 很容易成立。假設 n:= n+1 的 case 也成立(以下我們把fib簡寫為f):

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍f (n+1+k) = f n * f k + f (n+1) * f(k+1)      -- (Vor')

我們來證明 n:=n+2 的 case:

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ ‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍f (n+2+k)
  = { f 之定義 }
  ‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍f (n+k) + f (n+k+1)
  = { (Vor) & (Vor') }
  ‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍f (n-1) * f k + f n * f (k+1) +
  ‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍f n * f k + f (n+1) * f(k+1)
  = { f (n+1) = f n + f (n-1) }
  ‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍f (n+1) * f k + f n * f (k+1) + f (n+1) * f (k+1)
  = { f (n+2) = f (n+1) + f n } 
  ‍‍‍‍‍‍ ‍‍ ‍‍‍‍‍‍ ‍‍f (n+1) * f k + f (n+2) * f (k+1) .


如此證出了 (Vor).

參考資料

Holloway 的碩士論文可在這兒找到。

此外他和指導老師 Paul Cull 也有出一篇 journal 版:

這位 James Holloway 近況不明,只知並沒有待在學術界。不知是否去了月球。

Dijkstra 的 EWD654 “In honor of Fibonacci” (1978) 也推導出了一個 O(log n) 次遞迴呼叫的演算法。

Ron Knott (University of Surrey) 的這個網頁有很詳盡的資訊,值得一讀。

2019 年初 李祐棠親自做了許多實驗,考慮到了浮點數精確度、挑選適合的浮點運算函式庫等等實際問題,也很值得一讀。

This entry was posted in 計算算計 and tagged , . Bookmark the permalink. Post a comment or leave a trackback: Trackback URL.

Post a Comment

Your email is never published nor shared. Required fields are marked *

You may use these HTML tags and attributes <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

*
*