Try   HackMD

費氏數列分析

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


費氏數列

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

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

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

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

用數學表示 Fibonacci 數: (令

Nat 表示自然數的型別)
fib::NatNatfib(0)=0fib(1)=1fib(n+2)=fib(n+1)+fib(n)

在沒有其他機制(如 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)),我們不難推導出如下(也是遞迴的)定義:
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)=15(1+52)n15(152)n

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

(1+52)n
(152)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 內預設是以雙精度浮點數在儲存,它跟真正的

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×0.69424
,也就是說,
F(n)
N(n)
是線性關係,和我們對於沒調整過的
F(n)
運算以指數速度增長的認知符合。

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

F(n) 的演算法的時間複雜度為
O(n)
。但這當然是假設加法運算只需
O(1)
。如果考慮大數字的 bit 數,以這個演算法算
F(n)
需要的 bit 運算數目大約是
N(n2/2n/2)
(注意: 這邊說的都不是 big
O

Binet's Formula (公式解) 呢?

5 可以用牛頓法來算,假設每個
n
bit 除法需要
n2
個 bit 運算,但每一輪只需要留下前面
N  n+logn
個 bit. 整體上說來需要的 bit 運算數中最大宗的大約是
logn×(Nn+logn)2
, 再加上一些較不顯著的運算 —— 其實沒有明顯地比反覆加法快。

Vorobev 等式

另有一系列演算法只需

O(logn) 次遞迴呼叫便可算出
fib(n)
。要理解它們,可由這個問題出發

fib(n+k)
fib(n)
fib(k)
(及它們附近的數)有關聯嗎?

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

fib(n) 算出
fib(n+n)
,有望設計出
O(logn)
次遞迴呼叫的演算法。

n>=1 時,存在以下特性 (以下將 fib 簡寫為
f
):

f(n+k)=f(n1)×f(k)+f(n)×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×f(n)×f(n+1)+(f(n+1))2

接著根據費式數列的特性

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

f(2n)=2×f(n+1)×f(n)(f(n))2

那麼我們便可以由

f(n)
f(n+1)
算出
f(2n)
f(2n+1)
,不難由此做出一個只需
O(logn)
次遞迴呼叫的演算法(但這不表示執行時間會是
O(logn)
! Holloway 的碩論顯示差遠了)。

解決運算精準度的問題

上述因浮點數運算累積的誤差,當然有現成的解決方案,可運用 GMPMPFR 這這兩個 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 進位的數值,再來就能直接算出

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

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

ϕ 也是閃一下就出來了,一億位在 64 bit Linux, 2 GHz AMD Ryzen 5 需時 37 秒,相比
e
π
這類超越數,
ϕ
只需要
5
容易許多。

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

fib(93),要怎麼做呢?
fib(93)
到底有多少位數呢?我們可用
2n
作為
F(n)
的上界
,最後所需的位數至少就是 ceill(n*log(2)),相對應的我們運算中的浮點數精確度的要求,
2n
這個上界有點粗糙但有用,頂多會浪費點記憶體,最後除出來的小數點後面多幾個零而已,如果能套用更精確的上界當然更好。

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(logn) 的 Fibonacci 快速解遠比
O(1)
的公式解來得快。

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

O(1) 的加法、乘法都不再是常數時間,而是與數字的長度 k (位元數)有關。而上面我們有提到,基本上可以用
2n
作為費式數列的上界,也因此費式數列的數字長度 k ~= n,加法、乘法複雜度就會視實作方式上升到
O(n)
O(n2)
O(nlogn)
左右。

在 Fibonacci 快速解,我們需要做

O(logn) 次的 iteration,每次三個乘法兩個加減;公式解雖然沒有 iteration,但需要計算兩次次方運算,也等於是
O(logn)
次的乘法跟加法,然後還有除法,我們運算的又不是整數而是浮點數,這又需要更多的成本,一來一往之間就抵消了公式解直接算出答案的優勢。

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

O(logn) 的 Fibonacci 快速解看不出差異,等到
n
終於大到看得出
O(logn)
O(1)
的差異時,已把運算複雜度納入考量。

結語

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

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

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

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

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

參考資料