費氏數怎麼算？

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 n 和 fib 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 n 和 fib (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.69424 — fib n 和 N 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 := 1 和 2 的 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 的碩士論文可在這兒找到。

• James L Holloway. Algorithms for Computing Fibonaci Numbers Quickly. Master Thesis, Oregon State University, 1988.

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

• Paul Cull, James L. Holloway. Computing fibonacci numbers quickly. Information Processing Letters, 32(3), pp 143-149. 1989.

這位 James Holloway 近況不明，只知並沒有待在學術界。不知是否去了月球？無論你做的是什麼題目，只要懷抱著熱情，多年後總會被另一些懷抱同樣熱情的人發掘 — 這樣的想法令我覺得很溫馨與安心。

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

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

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