Leonardo Fibonacci 為 12 世紀到 13 世紀的義大利籍數學家,在他的著作《Liber Abaci》(字面上的意思是計算之書) 經提到一個場景:
「若有一隻免子每個月生一隻小免子,一個月後小免子也開始生產。起初只有一隻 免子,一個月後就有兩隻免子,二個月後有三隻免子,三個月後有五隻免子 (小免子投入生殖) 」。
注意新生的小免子需一個月成長期才會投入生產,類似的道理也可以用於植物的生長,這就是 Fibonacci 數列,一般習慣簡稱費氏數列,數列前幾項是:
1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89
用數學表示 Fibonacci 數: (令 表示自然數的型別)
在沒有其他機制(如 memoization)輔助的情況下直接跑上述定義,會得到一個很慢的演算法,因為許多計算重複了。早年的入門程式設計書籍常錯誤地把這當作「遞迴效率欠佳」的例子。
費式數列給一個初學程式的人都能寫得出來,下方是只用 4 行 Python 程式的實作,一方面展現 Python 在大數運算上的實力,一方面展視了它的簡潔,像是 a , b = a + b, a
這種寫法。
當然會寫是一回事,深入進去就不簡單:每次遞迴都會做重複的計算,於是計算以指數的方式成長,比如說用 Python 的 timeit 去測一個遞迴的費式數列函式,很快執行時間大幅增長,約到 fib(30) 以上就會跑得很吃力了。
如果定義 ,我們不難推導出如下(也是遞迴的)定義:
這個演算法的效率改進不少 —— 當 n 不太大時堪用。計算 fib2(n) 時只需 個遞迴呼叫,但這表示這個演算法的時間複雜度是 嗎?
費式數列有個所謂的「公式解」,即:
這個「公式解」一般認為由 Jacques P. M. Binet 在 1843 年發現,仍可追溯得更早。習慣上我們仍稱之為 Binet 等式 (Binet’s Formula)。該公式的證明是常見的練習,可參見 Cornell 大學的 CS 280: Induction and Recursion (Joe Halpern, 2005)。許多人可能看到這是一個封閉式,便認為這是最快的算法,甚至有人以為這個式子的右手邊可在常數時間內算出來。但別忘了, 和 對電腦 (和人類) 絕非沒有代價。
為什麼不用這個算式算呢?公式解為何不能看待是常數時間呢?
考慮電腦的離散本質,會遇上一系列問題,例如可用 64 位元整數表示最大值的 為例,整數算的解為 。若是 Python 寫的公式解呢?
會得到 12200160415121913856
,顯然這兩個數值不同!何以致此?問題就是浮點數的不精確問題,可參見〈你所不知道的 C 語言:數值系統篇〉理解背景知識。
我們在計算完 sqrt(5)
之後,只能用一個近似的值來表達結果,在 python 內預設是以雙精度浮點數在儲存,它跟真正的 還是有細微的差距,在隨後的 n 次方、除法上,這個細微的誤差都會被慢慢的放大,最終導致這個巨大的誤差。
James L Holloway 在 1988 年發表的碩士論文〈Algorithms for computing Fibonacci numbers quickly 〉回顧好幾個計算 Fibonacci 數的方法,包括會重複的遞迴、反覆的加法、Binet's Formula(即大家所說的「公式解」)、好幾種基於矩陣乘法(因乘法的遞移律,可以做 divide and conquer)的做法,生成函數等等,有理論分析,也有執行程式驗證。
一些有趣的結果如下。
用 表示第 個 Fib 數。理論分析方面,我們得先知道 大約需要多少個 bit 來表示。若 是 需要的 bit 數,大約可推算出 ,也就是說, 和 是線性關係,和我們對於沒調整過的 運算以指數速度增長的認知符合。
一般我們會說用反覆的加法(也就是用兩個變數,存最近的兩個 Fibonacci 數的手法,即 memoization)算 的演算法的時間複雜度為 。但這當然是假設加法運算只需 。如果考慮大數字的 bit 數,以這個演算法算 需要的 bit 運算數目大約是 (注意: 這邊說的都不是 big )
Binet's Formula (公式解) 呢? 可以用牛頓法來算,假設每個 bit 除法需要 個 bit 運算,但每一輪只需要留下前面 個 bit. 整體上說來需要的 bit 運算數中最大宗的大約是 , 再加上一些較不顯著的運算 —— 其實沒有明顯地比反覆加法快。
另有一系列演算法只需 次遞迴呼叫便可算出 。要理解它們,可由這個問題出發
與 和 (及它們附近的數)有關聯嗎?
如果可找出些關聯,即可由 算出 ,有望設計出 次遞迴呼叫的演算法。
當 時,存在以下特性 (以下將 fib 簡寫為 ):
如果有這個式子,首先令 ,可以得到
再令 ,可以得到
接著根據費式數列的特性 ,則可以由上面兩式得到
那麼我們便可以由 和 算出 和 ,不難由此做出一個只需 次遞迴呼叫的演算法(但這不表示執行時間會是 ! Holloway 的碩論顯示差遠了)。
上述因浮點數運算累積的誤差,當然有現成的解決方案,可運用 GMP 或 MPFR 這這兩個 GNU 函式庫是所謂的「無限」位數的整數跟「無限」精確度的浮點數,當然他們不是真的無限,只是完全壓榨記憶體來記錄儘可能多的位數以求精確,理論上記憶體撐上去就能把精確度逼上去。究竟這個函式庫有多麼的強大呢?我們可做簡單試驗,例如用 C 程式計算黃金比例:
唯一要注意的是 MPFR 內部用的 precision 是以 2 進位為底,所以我們在十進位需要的精確度,要先換算為 2 進位的數值,再來就能直接算出 ,試著算過 50000 位數再對個網路上找的答案,數字是完全一樣。
這個函式庫算得非常快,一百萬位的 也是閃一下就出來了,一億位在 64 bit Linux, 2 GHz AMD Ryzen 5 需時 37 秒,相比 和 這類超越數, 只需要 容易許多。
如果我們要用 MPFR 這個函式庫,利用公式解來算 ,要怎麼做呢?
到底有多少位數呢?我們可用 作為 的上界,最後所需的位數至少就是 ceill(n*log(2))
,相對應的我們運算中的浮點數精確度的要求, 這個上界有點粗糙但有用,頂多會浪費點記憶體,最後除出來的小數點後面多幾個零而已,如果能套用更精確的上界當然更好。
MPFR 的函式庫設計精良,呼叫上非常直覺,這段程式碼其實就是寫公式解,應該滿好懂的,程式碼在此。有了這機制,至少可確保 fib(10000) 的運算是正確的。
但公式解真的有比較快嗎?
答案是否定的,我們同樣可以用 Fibonacci 快速解搭配 GMP 函式庫來計算,因為都是整數的運算可以做到非常快,測試程式碼在此:
同樣是計算 fib(100,000,000):
的 Fibonacci 快速解遠比 的公式解來得快。
問題就在於,到了所謂的大數區域,本來我們假定 的加法、乘法都不再是常數時間,而是與數字的長度 k (位元數)有關。而上面我們有提到,基本上可以用 作為費式數列的上界,也因此費式數列的數字長度 k ~= n,加法、乘法複雜度就會視實作方式上升到 跟 或 左右。
在 Fibonacci 快速解,我們需要做 次的 iteration,每次三個乘法兩個加減;公式解雖然沒有 iteration,但需要計算兩次次方運算,也等於是 次的乘法跟加法,然後還有除法,我們運算的又不是整數而是浮點數,這又需要更多的成本,一來一往之間就抵消了公式解直接算出答案的優勢。
在通常的應用上以及現今電腦的實作,我們仍可假設整數的加減乘都能在近乎常數時間內結束,這樣我們才能討論資料結構與演算法的複雜度。費氏數列的問題在於,在數字仍小,就不用考慮運算複雜度,公式解和 的 Fibonacci 快速解看不出差異,等到 終於大到看得出 跟 的差異時,已把運算複雜度納入考量。
理論上我們當然可以假設有個計算模型,無論有多少位的數字,無論浮點數有多少精確度要求,四則運算與次方都能在常數時間內結束,這時公式解就能來到 ,但這樣的狀況一如停機問題所假設的萬能機器,不存在真實世界中。
運用 GMP 和 MPFR 這樣的函式庫,插滿記憶體甚至把硬碟當記憶體來用、把記憶體當 cache 用,消耗幾週光陰跟一堆電力,我們可把無理數算到小數點下一億位、十億位,甚至更多,取決於當下的硬體水準,這是前人們精心為我們建的巨塔,可是數字還是無窮無盡,站在巨塔上反而才看得出我們跟無限有多麼遙遠,誠然人腦可以透過思考一窺數學之奧妙,但不代表我們能超脫數學的嚴格限制浮空而起。
本文透過兩個實作,讓大家體會一下所謂公式解並不一定是 ,背後一定有對應的成本。