費氏數列分析

本文改寫自穆信成博士與李佑棠的文章

• Shin-Cheng Mu: 費氏數怎麼算？
• 李佑棠: 關於費式數列的那些事

---

費氏數列

Leonardo Fibonacci 為 12 世紀到 13 世紀的義大利籍數學家，在他的著作《Liber Abaci》(字面上的意思是計算之書) 經提到一個場景：

「若有一隻免子每個月生一隻小免子，一個月後小免子也開始生產。起初只有一隻 免子，一個月後就有兩隻免子，二個月後有三隻免子，三個月後有五隻免子 (小免子投入生殖) 」。

注意新生的小免子需一個月成長期才會投入生產，類似的道理也可以用於植物的生長，這就是 Fibonacci 數列，一般習慣簡稱費氏數列，數列前幾項是:

1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89

用數學表示 Fibonacci 數: (令

$Nat$ 表示自然數的型別)

$$\begin{gathered} fib::Nat\to Nat \\ fib(0)=0 \\ fib(1)=1 \\ fib(n+2)=fib(n+1)+fib(n) \end{gathered}$$

在沒有其他機制（如 memoization）輔助的情況下直接跑上述定義，會得到一個很慢的演算法，因為許多計算重複了。早年的入門程式設計書籍常錯誤地把這當作「遞迴效率欠佳」的例子。

費式數列給一個初學程式的人都能寫得出來，下方是只用 4 行 Python 程式的實作，一方面展現 Python 在大數運算上的實力，一方面展視了它的簡潔，像是 a , b = a + b, a 這種寫法。

a = b = 1
while b < 1000000000000:
    print(b)
    a, b = a + b, a

當然會寫是一回事，深入進去就不簡單：每次遞迴都會做重複的計算，於是計算以指數的方式成長，比如說用 Python 的 timeit 去測一個遞迴的費式數列函式，很快執行時間大幅增長，約到 fib(30) 以上就會跑得很吃力了。

如果定義

$fib2(n)=(fib(n),fib(n+1))$，我們不難推導出如下（也是遞迴的）定義：

$$\begin{gathered} fib2::Nat\to (Nat,Nat) \\ fib2(0)=(0,1) \\ fib2(n+1)=(y,x+y) \\ \text{where }(x,y)=fib2(n) \end{gathered}$$

這個演算法的效率改進不少 —— 當 n 不太大時堪用。計算 fib2(n) 時只需

$O(n)$ 個遞迴呼叫，但這表示這個演算法的時間複雜度是

$O(n)$ 嗎？

費式數列有個所謂的「公式解」，即:

$fib(n)=\frac{1}{\sqrt{5}}(\frac{1+\sqrt{5}}{2}{)}^n-\frac{1}{\sqrt{5}}(\frac{1-\sqrt{5}}{2}{)}^n$

這個「公式解」一般認為由 Jacques P. M. Binet 在 1843 年發現，仍可追溯得更早。習慣上我們仍稱之為 Binet 等式 (Binet’s Formula)。該公式的證明是常見的練習，可參見 Cornell 大學的 CS 280: Induction and Recursion (Joe Halpern, 2005)。許多人可能看到這是一個封閉式，便認為這是最快的算法，甚至有人以為這個式子的右手邊可在常數時間內算出來。但別忘了，

$(\frac{1+\sqrt{5}}{2}{)}^n$ 和

$(\frac{1-\sqrt{5}}{2}{)}^n$ 對電腦 (和人類) 絕非沒有代價。

為什麼不用這個算式算呢？公式解為何不能看待是常數時間呢？

考慮電腦的離散本質，會遇上一系列問題，例如可用 64 位元整數表示最大值的

$fib(93)$ 為例，整數算的解為

$12200160415121876738$。若是 Python 寫的公式解呢？

def fib(n):
    return (math.pow((1 + math.sqrt(5)) / 2, n) -
    math.pow((1 - math.sqrt(5))/2, n)) / math.sqrt(5)
print(int(fib(93)))

會得到 12200160415121913856，顯然這兩個數值不同！何以致此？問題就是浮點數的不精確問題，可參見〈你所不知道的 C 語言：數值系統篇〉理解背景知識。

我們在計算完 sqrt(5) 之後，只能用一個近似的值來表達結果，在 python 內預設是以雙精度浮點數在儲存，它跟真正的

$\sqrt{5}$ 還是有細微的差距，在隨後的 n 次方、除法上，這個細微的誤差都會被慢慢的放大，最終導致這個巨大的誤差。

James L Holloway 在 1988 年發表的碩士論文〈Algorithms for computing Fibonacci numbers quickly 〉回顧好幾個計算 Fibonacci 數的方法，包括會重複的遞迴、反覆的加法、Binet's Formula（即大家所說的「公式解」）、好幾種基於矩陣乘法（因乘法的遞移律，可以做 divide and conquer）的做法，生成函數等等，有理論分析，也有執行程式驗證。

一些有趣的結果如下。

用

$F(n)$ 表示第

$n$ 個 Fib 數。理論分析方面，我們得先知道

$F(n)$ 大約需要多少個 bit 來表示。若

$N(n)$ 是

$F(n)$ 需要的 bit 數，大約可推算出

$N(n)=n\times 0.69424$，也就是說，

$F(n)$ 和

$N(n)$ 是線性關係，和我們對於沒調整過的

$F(n)$ 運算以指數速度增長的認知符合。

一般我們會說用反覆的加法（也就是用兩個變數，存最近的兩個 Fibonacci 數的手法，即 memoization）算

$F(n)$ 的演算法的時間複雜度為

$O(n)$。但這當然是假設加法運算只需

$O(1)$。如果考慮大數字的 bit 數，以這個演算法算

$F(n)$ 需要的 bit 運算數目大約是

$N(n^2/2-n/2)$ (注意: 這邊說的都不是 big

$O$）

Binet's Formula (公式解) 呢？

$\sqrt{5}$ 可以用牛頓法來算，假設每個

$n$ bit 除法需要

$n^2$ 個 bit 運算，但每一輪只需要留下前面

$Nn+\log n$ 個 bit. 整體上說來需要的 bit 運算數中最大宗的大約是

$\log n\times (Nn+\log n{)}^2$, 再加上一些較不顯著的運算 —— 其實沒有明顯地比反覆加法快。

Vorobev 等式

另有一系列演算法只需

$O(\log n)$ 次遞迴呼叫便可算出

$fib(n)$。要理解它們，可由這個問題出發

$fib(n+k)$ 與

$fib(n)$ 和

$fib(k)$ (及它們附近的數)有關聯嗎？

如果可找出些關聯，即可由

$fib(n)$ 算出

$fib(n+n)$，有望設計出

$O(\log n)$ 次遞迴呼叫的演算法。

當

$n>=1$ 時，存在以下特性 (以下將 fib 簡寫為

$f$ ):

$$f(n+k)=f(n-1)\times f(k)+f(n)\times f(k+1)$$

如果有這個式子，首先令

$n:=n+1,k:=n$ ，可以得到

$f(2n+1)=(f(n+1){)}^2+(f(n){)}^2$

再令

$n:=n+1,k:=n+1$，可以得到

$f(2n+2)=2\times f(n)\times f(n+1)+(f(n+1){)}^2$

接著根據費式數列的特性

$f(2n)=f(2n+2)-f(2n+1)$ ，則可以由上面兩式得到

$f(2n)=2\times f(n+1)\times f(n)-(f(n){)}^2$

那麼我們便可以由

$f(n)$ 和

$f(n+1)$ 算出

$f(2n)$ 和

$f(2n+1)$，不難由此做出一個只需

$O(\log n)$ 次遞迴呼叫的演算法（但這不表示執行時間會是

$O(\log n)$! Holloway 的碩論顯示差遠了）。

解決運算精準度的問題

上述因浮點數運算累積的誤差，當然有現成的解決方案，可運用 GMP 或 MPFR 這這兩個 GNU 函式庫是所謂的「無限」位數的整數跟「無限」精確度的浮點數，當然他們不是真的無限，只是完全壓榨記憶體來記錄儘可能多的位數以求精確，理論上記憶體撐上去就能把精確度逼上去。究竟這個函式庫有多麼的強大呢？我們可做簡單試驗，例如用 C 程式計算黃金比例：

mpfr_t phi;
unsigned long int precision, x=5;
uint64_t digits = DIGITS;

precision = ceill((digits+1)*logl(10)/logl(2));
mpfr_init2(phi, precision);
mpfr_sqrt_ui(phi, x, MPFR_RNDN);
mpfr_add_ui(phi, phi, 1.0, MPFR_RNDN);
mpfr_div_ui(phi, phi, 2.0, MPFR_RNDN);
mpfr_printf("%.10000Rf\n", phi);
mpfr_clear(phi);

唯一要注意的是 MPFR 內部用的 precision 是以 2 進位為底，所以我們在十進位需要的精確度，要先換算為 2 進位的數值，再來就能直接算出

$\varphi$，試著算過 50000 位數再對個網路上找的答案，數字是完全一樣。

這個函式庫算得非常快，一百萬位的

$\varphi$ 也是閃一下就出來了，一億位在 64 bit Linux, 2 GHz AMD Ryzen 5 需時 37 秒，相比

$e$ 和

$\pi$ 這類超越數，

$\varphi$ 只需要

$\sqrt{5}$ 容易許多。

如果我們要用 MPFR 這個函式庫，利用公式解來算

$fib(93)$，要怎麼做呢？

$fib(93)$ 到底有多少位數呢？我們可用

$2^n$ 作為

$F(n)$ 的上界，最後所需的位數至少就是 ceill(n*log(2))，相對應的我們運算中的浮點數精確度的要求，

$2^n$ 這個上界有點粗糙但有用，頂多會浪費點記憶體，最後除出來的小數點後面多幾個零而已，如果能套用更精確的上界當然更好。

MPFR 的函式庫設計精良，呼叫上非常直覺，這段程式碼其實就是寫公式解，應該滿好懂的，程式碼在此。有了這機制，至少可確保 fib(10000) 的運算是正確的。

但公式解真的有比較快嗎？

答案是否定的，我們同樣可以用 Fibonacci 快速解搭配 GMP 函式庫來計算，因為都是整數的運算可以做到非常快，測試程式碼在此：

同樣是計算 fib(100,000,000)：

formulafib.c: 57.39s user 2.04s system 97% cpu 1:01.09 total
fastfib.c: 4.70s user 0.20s system 75% cpu 6.524 total

$O(\log n)$ 的 Fibonacci 快速解遠比

$O(1)$ 的公式解來得快。

問題就在於，到了所謂的大數區域，本來我們假定

$O(1)$ 的加法、乘法都不再是常數時間，而是與數字的長度 k (位元數)有關。而上面我們有提到，基本上可以用

$2^n$ 作為費式數列的上界，也因此費式數列的數字長度 k ~= n，加法、乘法複雜度就會視實作方式上升到

$O(n)$ 跟

$O(n^2)$ 或

$O(n\log n)$ 左右。

在 Fibonacci 快速解，我們需要做

$O(\log n)$ 次的 iteration，每次三個乘法兩個加減；公式解雖然沒有 iteration，但需要計算兩次次方運算，也等於是

$O(\log n)$ 次的乘法跟加法，然後還有除法，我們運算的又不是整數而是浮點數，這又需要更多的成本，一來一往之間就抵消了公式解直接算出答案的優勢。

在通常的應用上以及現今電腦的實作，我們仍可假設整數的加減乘都能在近乎常數時間內結束，這樣我們才能討論資料結構與演算法的複雜度。費氏數列的問題在於，在數字仍小，就不用考慮運算複雜度，公式解和

$O(\log n)$ 的 Fibonacci 快速解看不出差異，等到

$n$ 終於大到看得出

$O(\log n)$ 跟

$O(1)$ 的差異時，已把運算複雜度納入考量。

結語

理論上我們當然可以假設有個計算模型，無論有多少位的數字，無論浮點數有多少精確度要求，四則運算與次方都能在常數時間內結束，這時公式解就能來到

$O(1)$，但這樣的狀況一如停機問題所假設的萬能機器，不存在真實世界中。

運用 GMP 和 MPFR 這樣的函式庫，插滿記憶體甚至把硬碟當記憶體來用、把記憶體當 cache 用，消耗幾週光陰跟一堆電力，我們可把無理數算到小數點下一億位、十億位，甚至更多，取決於當下的硬體水準，這是前人們精心為我們建的巨塔，可是數字還是無窮無盡，站在巨塔上反而才看得出我們跟無限有多麼遙遠，誠然人腦可以透過思考一窺數學之奧妙，但不代表我們能超脫數學的嚴格限制浮空而起。

本文透過兩個實作，讓大家體會一下所謂公式解並不一定是

$O(1)$，背後一定有對應的成本。

參考資料

• James L Holloway 和指導老師 Paul Cull 合寫一篇期刊論文 (1989 年):
  Computing fibonacci numbers quickly。
• Dijkstra 的 EWD654 "In honor of Fibonacci" (1978) 推導出了一個
  $O(\log n)$ 次遞迴呼叫的演算法。
• Ron Knott (University of Surrey) 的網頁 有很詳盡的資訊，值得一讀
• 你不知道的費氏數列