# 費氏數列分析 本文改寫自[穆信成博士](http://www.iis.sinica.edu.tw/~scm/)與[李佑棠](https://yodalee.blogspot.com)的文章 * [Shin-Cheng Mu: 費氏數怎麼算?](https://scm.iis.sinica.edu.tw/ncs/2019/06/how-to-compute-fibonacci/) * [李佑棠: 關於費式數列的那些事](https://yodalee.blogspot.com/2019/02/fibonacci.html) --- ## 費氏數列 [Leonardo Fibonacci](https://en.wikipedia.org/wiki/Fibonacci) 為 12 世紀到 13 世紀的義大利籍數學家,在他的著作《Liber Abaci》(字面上的意思是計算之書) 經提到一個場景: > 「若有一隻免子每個月生一隻小免子,一個月後小免子也開始生產。起初只有一隻 免子,一個月後就有兩隻免子,二個月後有三隻免子,三個月後有五隻免子 (小免子投入生殖) 」。 注意新生的小免子需一個月成長期才會投入生產,類似的道理也可以用於植物的生長,這就是 Fibonacci 數列,一般習慣簡稱費氏數列,數列前幾項是: > 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89 用數學表示 Fibonacci 數: (令 $Nat$ 表示自然數的型別) $$ fib :: Nat \to Nat \\ fib(0) = 0 \\ fib(1) = 1 \\ fib(n+2) = fib(n+1) + fib(n) $$ 在沒有其他機制(如 memoization)輔助的情況下直接跑上述定義,會得到一個很慢的演算法,因為許多計算重複了。早年的入門程式設計書籍常錯誤地把這當作「遞迴效率欠佳」的例子。 費式數列給一個初學程式的人都能寫得出來,下方是只用 4 行 Python 程式的實作,一方面展現 Python 在大數運算上的實力,一方面展視了它的簡潔,像是 `a , b = a + b, a` 這種寫法。 ```python a = b = 1 while b < 1000000000000: print(b) a, b = a + b, a ``` 當然會寫是一回事,深入進去就不簡單:每次遞迴都會做重複的計算,於是計算以指數的方式成長,比如說用 Python 的 [timeit](https://docs.python.org/3/library/timeit.html) 去測一個遞迴的費式數列函式,很快執行時間大幅增長,約到 fib(30) 以上就會跑得很吃力了。 如果定義 $fib2(n) = (fib(n), fib(n+1))$,我們不難推導出如下(也是遞迴的)定義: $$ fib2 :: Nat \to (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) = \frac{1}{\sqrt{5}}(\frac{1+\sqrt{5}}{2})^n-\frac{1}{\sqrt{5}}(\frac{1-\sqrt{5}}{2})^n$ 這個「公式解」一般認為由 [Jacques P. M. Binet](https://en.wikipedia.org/wiki/Jacques_Philippe_Marie_Binet) 在 1843 年發現,仍可[追溯得更早](http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/fibFormula.html)。習慣上我們仍稱之為 Binet 等式 (Binet’s Formula)。該公式的證明是常見的練習,可參見 Cornell 大學的 [CS 280: Induction and Recursion](http://www.cs.cornell.edu/courses/cs280/2005fa/induction.pdf) (Joe Halpern, 2005)。許多人可能看到這是一個封閉式,便認為這是最快的算法,甚至有人以為這個式子的右手邊可在常數時間內算出來。但別忘了,$(\frac{1+\sqrt{5}}{2})^n$ 和 $(\frac{1-\sqrt{5}}{2})^n$ 對電腦 (和人類) 絕非沒有代價。 為什麼不用這個算式算呢?公式解為何不能看待是常數時間呢? 考慮電腦的離散本質,會遇上一系列問題,例如可用 64 位元整數表示最大值的 $fib(93)$ 為例,整數算的解為 $12200160415121876738$。若是 Python 寫的公式解呢? ```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 語言:數值系統篇](https://hackmd.io/@sysprog/c-numerics) 理解背景知識。 我們在計算完 `sqrt(5)` 之後,只能用一個近似的值來表達結果,在 python 內預設是以雙精度浮點數在儲存,它跟真正的 $\sqrt{5}$ 還是有細微的差距,在隨後的 n 次方、除法上,這個細微的誤差都會被慢慢的放大,最終導致這個巨大的誤差。 James L Holloway 在 1988 年發表的碩士論文 [Algorithms for computing Fibonacci numbers quickly](https://ir.library.oregonstate.edu/downloads/t435gg51w) 回顧好幾個算 Fibonacci 數的方法,包括會重複的遞迴、反覆的加法、[Binet's Formula](https://proofwiki.org/wiki/Euler-Binet_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 運算,但每一輪只需要留下前面 $N\ \ n + log n$ 個 bit. 整體上說來需要的 bit 運算數中最大宗的大約是 $log n \times (N n + 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, k := n+1$ ,可以得到 $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](https://gmplib.org/) 或 [MPFR](https://www.mpfr.org/) 這這兩個 GNU 函式庫是所謂的「無限」位數的整數跟「無限」精確度的浮點數,當然他們不是真的無限,只是完全壓榨記憶體來記錄儘可能多的位數以求精確,理論上記憶體撐上去就能把精確度逼上去。究竟這個函式庫有多麼的強大呢?我們可做簡單試驗,例如用 C 程式計算黃金比例: ```cpp 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](https://www.mpfr.org/) 內部用的 precision 是以 2 進位為底,所以我們在十進位需要的精確度,要先換算為 2 進位的數值,再來就能直接算出 $\phi$,試著算過 50000 位數再對個[網路上找的答案](https://www2.cs.arizona.edu/icon/oddsends/phi.htm),數字是完全一樣。 這個函式庫算得非常快,一百萬位的 $\phi$ 也是閃一下就出來了,一億位在 64 bit Linux, 2 GHz AMD Ryzen 5 需時 37 秒,相比 $e$ 和 $\pi$ 這類超越數,$\phi$ 只需要 $\sqrt{5}$ 容易許多。 如果我們要用 [MPFR](https://www.mpfr.org/) 這個函式庫,利用公式解來算 $fib(93)$,要怎麼做呢? $fib(93)$ 到底有多少位數呢?我們[可用 $2^n$ 作為 $F(n)$ 的上界](https://math.stackexchange.com/questions/2971350/show-that-log-fib-n-is-thetan),最後所需的位數至少就是 `ceill(n*log(2))`,相對應的我們運算中的浮點數精確度的要求,$2^n$ 這個上界有點粗糙但有用,頂多會浪費點記憶體,最後除出來的小數點後面多幾個零而已,如果能套用更精確的上界當然更好。 [MPFR](https://www.mpfr.org/) 的函式庫設計精良,呼叫上非常直覺,這段程式碼其實就是寫公式解,應該滿好懂的,[程式碼在此](https://gist.github.com/yodalee/4e221b081be4b367e9c7ef328ada7db5)。有了這機制,至少可確保 fib(10000) 的運算是正確的。 但公式解真的有比較快嗎? 答案是否定的,我們同樣可以用 Fibonacci 快速解搭配 [GMP](https://gmplib.org/) 函式庫來計算,因為都是整數的運算可以做到非常快,[測試程式碼在此](https://gist.github.com/yodalee/4e221b081be4b367e9c7ef328ada7db5): 同樣是計算 fib(100,000,000): ```shell 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)$,但這樣的狀況一如[停機問題](https://en.wikipedia.org/wiki/Halting_problem)所假設的萬能機器,不存在真實世界中。 運用 [GMP](https://gmplib.org/) 和 [MPFR](https://www.mpfr.org/) 這樣的函式庫,插滿記憶體甚至把硬碟當記憶體來用、把記憶體當 cache 用,消耗幾週光陰跟一堆電力,我們可把無理數算到小數點下一億位、十億位,甚至更多,取決於當下的硬體水準,這是前人們精心為我們建的巨塔,可是數字還是無窮無盡,站在巨塔上反而才看得出我們跟無限有多麼遙遠,誠然人腦可以透過思考一窺數學之奧妙,但不代表我們能超脫數學的嚴格限制浮空而起。 本文透過兩個實作,讓大家體會一下所謂公式解並不一定是 $O(1)$,背後一定有對應的成本。 ## 參考資料 * James L Holloway 和指導老師 Paul Cull 合寫一篇期刊論文 (1989 年): [Computing fibonacci numbers quickly](https://doi.org/10.1016/0020-0190(89)90015-X)。 * Dijkstra 的 [EWD654 "In honor of Fibonacci"](http://www.cs.utexas.edu/users/EWD/ewd06xx/EWD654.PDF) (1978) 推導出了一個 $O(\log{n})$ 次遞迴呼叫的演算法。 * Ron Knott (University of Surrey) 的[網頁](http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/fibFormula.html) 有很詳盡的資訊,值得一讀。