Python 費氏數列解法(四):公式解與推導

費氏數列有 O(1) 解法嗎?

第二篇介紹了費氏數列的矩陣解法,不過費氏數列其實是可以直接用公式算出第 n 項的值的,這邊就來介紹並推導一下公式解,順便幫自己複習一下數學 XD

公式解:Binet’s Formula

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

是的,$\frac{1 + \sqrt{5}}{2}$,傳說中的黃金比例。

但奇怪,費氏數列不是都是整數嗎?為什麼會跑出一堆根號呢?不信邪的話自己寫程式去跑,就會發現解剛好都是整數:

# binet.py
def fib_binet(n):
    golden_ratio = (1 + 5 ** 0.5) / 2
    golden_ratio_alt = (1 - 5 ** 0.5) / 2
    return (golden_ratio ** n - golden_ratio_alt ** n) / 5 ** 0.5

    # 一般來說,還是會轉換型態成 int:
    # return round((golden_ratio ** n - golden_ratio_alt ** n) / 5 ** 0.5)

if __name__ == '__main__':
    for i in range(30):
        print(fib_binet(i))
$ python binet.py
0.0
1.0
1.0
2.0
3.0
5.0
8.0
13.0
21.0
...

公式推導

這邊需要一些基本的線性代數知識,主要是對角化的部分。

對於從量子化學領域跳過來的我,解 eigenvalues 就像回母校一樣充滿了熟悉感 XD

下面就只附上大概的流程,計算細節就都省略。

讓我們從第二篇的矩陣表示法開始:

$$
\begin{bmatrix}
F(n+1)\\
F(n)
\end{bmatrix} = \begin{bmatrix}
1 & 1\\
1 & 0
\end{bmatrix}^{n} \begin{bmatrix}
1\\
0
\end{bmatrix}
$$

中間的矩陣我們令其為 A:

$$
A = \begin{bmatrix}
1 & 1\\
1 & 0
\end{bmatrix}
$$

因為要算 A 的 n 次方,可以將其對角化,對角矩陣的 n 次方即是其每個值的 n 次方:

$$
A^n = (Q\lambda Q^{-1})^n = Q\lambda^nQ^{-1}
$$

而 A 的兩個 eigenvalues 為:

$$
\frac{1 \pm \sqrt{5}}{2}
$$

對應的 eigenvectors 分別為($+$對$+$,$-$對$-$):

$$
t \begin{bmatrix}
1 \\
\frac{- 1 \pm \sqrt{5}}{2}
\end{bmatrix}
$$

所以原本的 $A^n$ 就可以利用對角化展開:

$$
A^n = \begin{bmatrix}
1 & 1\\
1 & 0
\end{bmatrix}^n = \begin{bmatrix}
1 & 1\\
\frac{- 1 + \sqrt{5}}{2} & \frac{- 1 - \sqrt{5}}{2}
\end{bmatrix} \begin{bmatrix}
\frac{1 + \sqrt{5}}{2} & 0\\
0 & \frac{1 - \sqrt{5}}{2}
\end{bmatrix}^n \begin{bmatrix}
1 & 1\\
\frac{- 1 + \sqrt{5}}{2} & \frac{- 1 - \sqrt{5}}{2}
\end{bmatrix}^{-1}
$$

帶回原本的式子:

$$
\begin{bmatrix}
F(n+1)\\
F(n)
\end{bmatrix} = \begin{bmatrix}
1 & 1\\
\frac{- 1 + \sqrt{5}}{2} & \frac{- 1 - \sqrt{5}}{2}
\end{bmatrix} \begin{bmatrix}
\frac{1 + \sqrt{5}}{2} & 0\\
0 & \frac{1 - \sqrt{5}}{2}
\end{bmatrix}^n \begin{bmatrix}
1 & 1\\
\frac{- 1 + \sqrt{5}}{2} & \frac{- 1 - \sqrt{5}}{2}
\end{bmatrix}^{-1} \begin{bmatrix}
1\\
0
\end{bmatrix}
$$

剩下就是高中數學的程度了,乘開簡化之後就可以得到公式解了。

公式解的 issues

浮點數精確度

可以注意到,公式解多了很多浮點數運算。

浮點數就是二進位的科學記號,所以當數值大的時候就會有精確度問題出現。因為我們最後會過一個 int,數值小的時候這個誤差其實是不會影響的,但數值夠大就可以看到誤差。

# ./binet.py
# ...
x = 0
while True:
    x += 1
    diff = fib_binet(x) - fib_fast_double_iter(x)
    if diff != 0:
        print(x, diff)
        break
$ python binet.py
123 -1

實際寫個 code 去跑就可以發現,在第 124 項之後就會出現誤差。

其實用上一篇的fast doubling來計算,不但沒有慢上多少,也不用擔心誤差問題。

常數時間?

不是程式碼沒有迴圈就是常數時間,這個公式裡面有黃金比例次方:

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

如果是二的次方,可以靠 shift 來達到常數時間運算,但這是個浮點數次方,以前篇提到的快速冪 fast doubling 來計算也是要 $O(\lg n)$ 的複雜度,即便這個增長很小,在理論上也不是常數時間。

有趣的小故事

其實會來研究 Python 與費氏數列的,與幾年前一串 Python Taiwan 的討論串有關,可惜最熱鬧的那則已經被刪掉了。

故事大致上就是有人在底下說了費氏數列公式解是 $O(1)$ 的言論,甚至提出只要把 $\pi$ 和 $e$ 等無理數都存在月球上的話,就能直接 $O(1)$ 解決很多問題,嗯…

其實 $O(1)$ 在這種時候已經變成迷思了,不然我大可以寫個 function,不管輸入是多少,我都故弄玄虛,算到一個超級大的值,然後再把相對應個輸出丟出去就好,就可以宣稱我是 $O(1)$ 解。

def fib_fake_o1(n):
    max_n = 300
    if n >= max_n:
        raise Exception(f"n must be < {max_n}")

    f1, f2 = 1, 1
    res = f2
    for i in range(2, max_n):
        f1, f2 = f2, f1 + f2
        if i == n:
            res = f2
    return res

Big O 複雜度分析是把常數忽略的,但實務上那個常數也是重要的。

小結

綜規以上來看,其實公式解也沒這麼厲害,而且也不是真的很省事(誤差問題)。

不過實務上通常不會需要這麼大的 fibonacci number 吧…,程式是死的,人是活的,還是要看情況決定哪種解法好。

大部分時候甚至寫個 O(n) 解法就夠用了。(實際上我面試那天在寫出 O(n)解法之後也沒有再被追問了,畢竟也只是眾多面試題中的其中一個小問題。)

參考資料: