15. 浮點數運算：問題與限制

在計算機架構中，浮點數 (floating-point number) 是以基數為 2（二進位）的小數表示。例如說，在十進位小數中 0.625 可被分為 6/10 + 2/100 + 5/1000，同樣的道理，二進位小數 0.101 可被分為 1/2 + 0/4 + 1/8。這兩個小數有相同的數值，而唯一真正的不同在於前者以十進位表示，後者以二進位表示。

不幸的是，大多數十進位小數無法精準地以二進位小數表示。一般的結果為，你輸入的十進位浮點數只能由實際儲存在計算機中的二進位浮點數近似。

在十進位中，這個問題更容易被理解。以分數 1/3 為例，你可以將其近似為十進位小數：

0.3

或者，更好的近似：

0.33

或者，更好的近似：

0.333

依此類推，不論你使用多少位數表示小數，最後的結果都無法精準地表示 1/3，但你還是能越來越精準地表示 1/3。

同樣的道理，不論你願意以多少位數表示二進位小數，十進位小數 0.1 都無法被二進位小數精準地表達。在二進位小數中，1/10 會是一個無限循環小數：

0.0001100110011001100110011001100110011001100110011...

只要你停在任何有限的位數，你就只會得到近似值。而現在大多數的計算機中，浮點數是透過二進位分數近似的，其中分子是從最高有效位元開始，用 53 個位元表示，分母則是以二為底的指數。在 1/10 的例子中，二進位分數為 3602879701896397 / 2 ** 55，而這樣的表示十分地接近，但不完全等同於 1/10 的真正數值。

由於數值顯示的方式，很多使用者並沒有發現數值是個近似值。Python 只會印出一個十進位近似值，其近似了儲存在計算機中的二進位近似值的真正十進位數值。在大多數的計算機中，如果 Python 真的會印出完整的十進位數值，其表示儲存在計算機中的 0.1 的二進位近似值，它將顯示為：

>>> 0.1
0.1000000000000000055511151231257827021181583404541015625

這比一般人感到有用的位數還多，所以 Python 將位數保持在可以接受的範圍，只顯示捨入後的數值：

>>> 1 / 10
0.1

一定要記住，雖然印出的數字看起來是精準的 1/10，但真正儲存的數值是能表示的二進位分數中，最接近精準數值的數。

有趣的是，有許多不同的十進位數，共用同一個最接近的二進位近似分數。例如說：數字 0.1 和 0.10000000000000001 和 0.1000000000000000055511151231257827021181583404541015625，都由 3602879701896397 / 2 ** 55 近似。由於這三個十進位數值共用同一個近似值，任何一個數值都可以被顯示，同時保持 eval(repr(x)) == x。

歷史上，Python 的提示字元 (prompt) 與內建的 repr() 函式會選擇上段說明中有 17 個有效位元的數：0.10000000000000001。從 Python 3.1 版開始，Python（在大部分的系統上）現在能選擇其中最短的數並簡單地顯示為 0.1。

注意，這是二進位浮點數理所當然的特性，並不是 Python 的錯誤 (bug)，更不是你程式碼的錯誤。只要程式語言有支援硬體的浮點數運算，你將會看到同樣的事情出現在其中（雖然有些程式語言預設不會顯示該差異，有些甚至是在所有的輸出模式中都不會顯示。）

為求更優雅的輸出，你可能想要使用字串的格式化 (string formatting) 產生限定的有效位數：

>>> format(math.pi, '.12g')  # give 12 significant digits
'3.14159265359'

>>> format(math.pi, '.2f')   # give 2 digits after the point
'3.14'

>>> repr(math.pi)
'3.141592653589793'

要了解一件很重要的事，在真正意義上，浮點數的表示是一種幻覺：你基本上在捨入真正機器數值所展示的值。

這種幻覺可能會產生下一個幻覺。舉例來說，因為 0.1 不是真正的 1/10，把三個 0.1 的值相加，也不會產生精準的 0.3：

>>> 0.1 + 0.1 + 0.1 == 0.3
False

同時，因為 0.1 不能再更接近精準的 1/10，還有 0.3 不能再更接近精準的 3/10，預先用 round() 函式捨入並不會有幫助：

>>> round(0.1, 1) + round(0.1, 1) + round(0.1, 1) == round(0.3, 1)
False

雖然數字不會再更接近他們的精準數值，但 math.isclose() 函式可以用來比較不精確的值：

>>> math.isclose(0.1 + 0.1 + 0.1, 0.3)
True

或者，可以使用 round() 函式來比較粗略的近似值：

>>> round(math.pi, ndigits=2) == round(22 / 7, ndigits=2)
True

二進位浮點數架構下還有很多這樣的意外。底下的「表示法誤差 (Representation Error)」章節，詳細地解釋了「0.1」的問題。Examples of Floating Point Problems（浮點數問題範例）一文提供了二進位浮點數的作用方式與現實中這類常見問題的摘錄。如果想要其他常見問題的更完整描述，可以參考 The Perils of Floating Point（浮點數的風險）。

如同上文的末段所述，「這個問題並沒有簡單的答案。」不過，也不必對浮點過度擔心！Python 的 float（浮點數）運算中的誤差，是從浮點硬體繼承而來，而在大多數的計算機上，每次運算的誤差範圍不會大於 2**53 分之 1。這對大多數的任務來說已經相當足夠，但你需要記住，它並非十進位運算，且每一次 float 運算都可能會承受新的捨入誤差。

雖然浮點運算確實存在一些問題，但在一般情況下，如果你只是把最終結果的顯示值，以十進位方式捨入至預期的位數，那麼仍會得到你預期的結果。str() 通常就能滿足要求，而若想要更細緻的控制，可參閱格式化文字語法中關於 str.format() method（方法）的格式規範。

對於需要精準十進位表示法的用例，可以試著用 decimal 模組，它可實作適用於會計應用程式與高精確度應用程式的十進位運算。

另一種支援精準運算的格式為 fractions 模組，該模組基於有理數來實作運算（因此可以精確表示像 1/3 這樣的數字）。

如果你是浮點運算的重度使用者，你應該看一下 NumPy 套件，以及由 SciPy 專案提供的許多用於數學和統計學運算的其他套件。請參閱 <https://scipy.org>。

在罕見情況下，當你真的想知道一個 float 的精準值，Python 提供的工具可協助達成。float.as_integer_ratio() method 可將一個 float 的值表示為分數：

>>> x = 3.14159
>>> x.as_integer_ratio()
(3537115888337719, 1125899906842624)

由於該比率是精準的，它可無損地再現該原始值：

>>> x == 3537115888337719 / 1125899906842624
True

float.hex() method 以十六進位（基數為 16）表示 float，一樣可以給出你的電腦所儲存的精準值：

>>> x.hex()
'0x1.921f9f01b866ep+1'

這種精確的十六進位表示法可用於精準地重建 float 值：

>>> x == float.fromhex('0x1.921f9f01b866ep+1')
True

由於該表示法是精準的，因此適用於在不同版本的 Python 之間可靠地傳送數值（獨立於系統平台），並與支援相同格式的其他語言（如 JAVA 和 C99）交換資料。

另一個有用的工具是 sum() 函式，能在計算總和時幫忙減少精確度的損失。當數值被加到運行中的總計值時，它會追蹤「失去的位數 (lost digits)」。這可以明顯改善總體準確度 (overall accuracy)，使得誤差不至於累積到影響最終總計值的程度：

>>> 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 == 1.0
False
>>> sum([0.1] * 10) == 1.0
True

math.fsum() 更進一步將數值被加到運行總計中的總計時追蹤所有「失去的位數」，以便結果僅有一次捨入。這比 sum() 慢，但在不常見的情況下會更準確，在這種情況下，大量輸入大多相互抵消，最終的總和會接近於零：

>>> arr = [-0.10430216751806065, -266310978.67179024, 143401161448607.16,
...        -143401161400469.7, 266262841.31058735, -0.003244936839808227]
>>> float(sum(map(Fraction, arr)))   # Exact summation with single rounding
8.042173697819788e-13
>>> math.fsum(arr)                   # Single rounding
8.042173697819788e-13
>>> sum(arr)                         # Multiple roundings in extended precision
8.042178034628478e-13
>>> total = 0.0
>>> for x in arr:
...     total += x                   # Multiple roundings in standard precision
...
>>> total                            # Straight addition has no correct digits!
-0.0051575902860057365

15.1. 表示法誤差 (Representation Error)

本節將會詳細解釋「0.1」的例子，並說明你自己要如何對類似案例執行精準的分析。 以下假設你對二進位浮點表示法已有基本程度的熟悉。

Representation error（表示法誤差）是指在真實情況中，有些（實際上是大多數）十進位小數並不能精準地以二進位（基數為 2）小數表示。這是 Python（或 Perl、C、C++、JAVA、Fortran 和其他許多）通常不會顯示你期望的精準十進位數字的主要原因。

為什麼呢？因為 1/10 無法精準地以一個二進位小數來表示。至少自從 2000 年以來，幾乎所有的計算機皆使用 IEEE-754 浮點數運算標準，並且幾乎所有的平台都以 IEEE 754 binary64 標準中的「雙精度 (double precision)」來作為 Python 的 float。IEEE 754 binary64 的值包含 53 位元的精度，所以在輸入時，電腦會努力把 0.1 轉換到最接近的分數，以 J/2**N 的形式表示，此處 J 是一個正好包含 53 位元的整數。可以將：

1 / 10 ~= J / (2**N)

重寫為：

J ~= 2**N / 10

而前面提到 J 有精準的 53 位元（即 >= 2**52 但 < 2**53），所以 N 的最佳數值是 56：

>>> 2**52 <=  2**56 // 10  < 2**53
True

意即，要使 J 正好有 53 位元，則 56 會是 N 的唯一值。而 J 最有可能的數值就是經過捨入後的該商數：

>>> q, r = divmod(2**56, 10)
>>> r
6

由於餘數超過 10 的一半，所以最佳的近似值是透過進位而得：

>>> q+1
7205759403792794

所以，在 IEEE 754 雙精度下，1/10 的最佳近似值是：

7205759403792794 / 2 ** 56

將分子和分母同除以二，會約分為：

3602879701896397 / 2 ** 55

請注意，由於我們有進位，所以這實際上比 1/10 大了一點；如果我們沒有進位，商數將會有點小於 1/10。但在任何情況下都不可能是精準的 1/10！

所以電腦從來沒有「看到」1/10：它看到的是上述的精準分數，也就是它能得到的 IEEE 754 double 最佳近似值：

>>> 0.1 * 2 ** 55
3602879701896397.0

如果將該分數乘以 10**55，則可以看到該值以 55 個十進位數字顯示：

>>> 3602879701896397 * 10 ** 55 // 2 ** 55
1000000000000000055511151231257827021181583404541015625

這表示儲存在電腦中的精準數值等於十進位值 0.1000000000000000055511151231257827021181583404541015625。與其顯示完整的十進位數值，許多語言（包括 Python 的舊版本）選擇將結果捨入至 17 個有效位數：

>>> format(0.1, '.17f')
'0.10000000000000001'

fractions 與 decimal 模組能使這些計算變得容易：

>>> from decimal import Decimal
>>> from fractions import Fraction

>>> Fraction.from_float(0.1)
Fraction(3602879701896397, 36028797018963968)

>>> (0.1).as_integer_ratio()
(3602879701896397, 36028797018963968)

>>> Decimal.from_float(0.1)
Decimal('0.1000000000000000055511151231257827021181583404541015625')

>>> format(Decimal.from_float(0.1), '.17')
'0.10000000000000001'