使用 Decimal 提高浮點運算精確度

以費氏數列公式解為例

第四篇提到費氏數列的公式解會遇到浮點數問題,只能精準算到 122 位。因為浮點數是利用二進位的小數來做儲存,對於非二進位的數字會有誤差。而且其儲存位數有上限,對於無理數的運算,超過一定的大小就會出現精確度問題。

所以這篇就要來簡單使用 Python 內建的 Decimal module 來拉高浮點數運算的精確度。

以費氏數列公式解為例

首先回顧一下公式:

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

可以注意到裡面有無理數,也就是黃金比例的部分。對無理數做 n 次方運算,當 n 放大後很明顯誤差也會樂乘越大。為了提高精確度,我們可以提高儲存的位元。

Decimal module

Python 提供了一個內建 module:Decimal,是一個十進位運算模型,可以避免一些浮點數誤差問題:

In [1]: 0.1 * 3 == 0.3
Out[1]: False

In [2]: from decimal import Decimal
    ...: 

In [3]: Decimal('0.1') * 3 == Decimal('0.3')
Out[3]: True

注意我們是使用 Decimal('0.1') 而不是 Decimal(0.1),因為後者會把不精確的 float 0.1 丟到 Decimal,還是變得不精確:

In [4]: Decimal(0.1)
Out[4]: Decimal('0.1000000000000000055511151231257827021181583404541015625')

套用到費氏數列公式解

原本的 code:

def fib_binet(n):
    golden_ratio = (1 + 5 ** 0.5) / 2
    golden_ratio_alt = (1 - 5 ** 0.5) / 2
    return round((golden_ratio ** n - golden_ratio_alt ** n) / 5 ** 0.5)

我們把數字用 Decimal 包起來:

def accurate_fib_binet(n):
    sqrt5 = Decimal(5) ** Decimal('0.5')
    golden_ratio = (1 + sqrt5) / Decimal(2)
    golden_ratio_alt = (1 - sqrt5) / Decimal(2)
    return round((golden_ratio ** Decimal(n) - golden_ratio_alt ** Decimal(n)) / sqrt5)

並拿矩陣解法fib_fast_double_iter 來比較:

x = 0
while True:
    x += 1
    diff = decimal_fib_binet(x) - fib_fast_double_iter(x)
    if diff != 0:
        print(x, diff)
        break
$ python binet_de.py
123 -1

嗯?跟原本一樣?

主要是因為我們的問題最主要是儲存位數不夠,畢竟我們原本的數值就是無理數,用十進位也無法精確表示。

然而 Decimal module 提供了我們提升精確度的方法:

以空間換取精確度

from decimal import getcontext

getcontext().prec = 10000  # default is 28

把上面這段 code 加進去之後,原本的 while loop 就會跑很久…,也就代表數字到很大的時候,還是可以算準。我自己測試約可以到四萬附近。

不過精確度越大,就會算越久,可以從 MAX_PREC 這個常數看最大可以設到多少:

In [1]: from decimal import MAX_PREC
   ...: MAX_PREC
Out[1]: 999999999999999999

不要真的設到這麼大,不然會卡住 …

但實際上還是要看應用情境吧,對我來說如果不知道進來的數字會多大,不如直接使用公式解以外的算法比較保險。

Reference