15. 浮點(diǎn)算術(shù):爭(zhēng)議和限制?

Floating-point numbers are represented in computer hardware as base 2 (binary) fractions. For example, the decimal fraction 0.125 has value 1/10 + 2/100 + 5/1000, and in the same way the binary fraction 0.001 has value 0/2 + 0/4 + 1/8. These two fractions have identical values, the only real difference being that the first is written in base 10 fractional notation, and the second in base 2.

不幸的是,大多數(shù)的十進(jìn)制小數(shù)都不能精確地表示為二進(jìn)制小數(shù)。這導(dǎo)致在大多數(shù)情況下,你輸入的十進(jìn)制浮點(diǎn)數(shù)都只能近似地以二進(jìn)制浮點(diǎn)數(shù)形式儲(chǔ)存在計(jì)算機(jī)中。

用十進(jìn)制來(lái)理解這個(gè)問(wèn)題顯得更加容易一些。考慮分?jǐn)?shù) 1/3 。我們可以得到它在十進(jìn)制下的一個(gè)近似值

0.3

或者,更近似的,:

0.33

或者,更近似的,:

0.333

以此類推。結(jié)果是無(wú)論你寫(xiě)下多少的數(shù)字,它都永遠(yuǎn)不會(huì)等于 1/3 ,只是更加更加地接近 1/3 。

同樣的道理,無(wú)論你使用多少位以 2 為基數(shù)的數(shù)碼,十進(jìn)制的 0.1 都無(wú)法精確地表示為一個(gè)以 2 為基數(shù)的小數(shù)。 在以 2 為基數(shù)的情況下, 1/10 是一個(gè)無(wú)限循環(huán)小數(shù)

0.0001100110011001100110011001100110011001100110011...

在任何一個(gè)位置停下,你都只能得到一個(gè)近似值。因此,在今天的大部分架構(gòu)上,浮點(diǎn)數(shù)都只能近似地使用二進(jìn)制小數(shù)表示,對(duì)應(yīng)分?jǐn)?shù)的分子使用每 8 字節(jié)的前 53 位表示,分母則表示為 2 的冪次。在 1/10 這個(gè)例子中,相應(yīng)的二進(jìn)制分?jǐn)?shù)是 3602879701896397 / 2 ** 55 ,它很接近 1/10 ,但并不是 1/10 。

大部分用戶都不會(huì)意識(shí)到這個(gè)差異的存在,因?yàn)?Python 只會(huì)打印計(jì)算機(jī)中存儲(chǔ)的二進(jìn)制值的十進(jìn)制近似值。在大部分計(jì)算機(jī)中,如果 Python 想把 0.1 的二進(jìn)制對(duì)應(yīng)的精確十進(jìn)制打印出來(lái),將會(huì)變成這樣

>>>
>>> 0.1
0.1000000000000000055511151231257827021181583404541015625

這比大多數(shù)人認(rèn)為有用的數(shù)字更多,因此Python通過(guò)顯示舍入值來(lái)保持可管理的位數(shù)

>>>
>>> 1 / 10
0.1

牢記,即使輸出的結(jié)果看起來(lái)好像就是 1/10 的精確值,實(shí)際儲(chǔ)存的值只是最接近 1/10 的計(jì)算機(jī)可表示的二進(jìn)制分?jǐn)?shù)。

有趣的是,有許多不同的十進(jìn)制數(shù)共享相同的最接近的近似二進(jìn)制小數(shù)。例如, 0.1 、 0.10000000000000001 、 0.1000000000000000055511151231257827021181583404541015625 全都近似于 3602879701896397 / 2 ** 55 。由于所有這些十進(jìn)制值都具有相同的近似值,因此可以顯示其中任何一個(gè),同時(shí)仍然保留不變的 eval(repr(x)) == x 。

在歷史上,Python 提示符和內(nèi)置的 repr() 函數(shù)會(huì)選擇具有 17 位有效數(shù)字的來(lái)顯示,即 0.10000000000000001。 從 Python 3.1 開(kāi)始,Python(在大多數(shù)系統(tǒng)上)現(xiàn)在能夠選擇這些表示中最短的并簡(jiǎn)單地顯示 0.1

請(qǐng)注意這種情況是二進(jìn)制浮點(diǎn)數(shù)的本質(zhì)特性:它不是 Python 的錯(cuò)誤,也不是你代碼中的錯(cuò)誤。 你會(huì)在所有支持你的硬件中的浮點(diǎn)運(yùn)算的語(yǔ)言中發(fā)現(xiàn)同樣的情況(雖然某些語(yǔ)言在默認(rèn)狀態(tài)或所有輸出模塊下都不會(huì) 顯示 這種差異)。

想要更美觀的輸出,你可能會(huì)希望使用字符串格式化來(lái)產(chǎn)生限定長(zhǎng)度的有效位數(shù):

>>>
>>> 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'

必須重點(diǎn)了解的是,這在實(shí)際上只是一個(gè)假象:你只是將真正的機(jī)器碼值進(jìn)行了舍入操作再 顯示 而已。

一個(gè)假象還可能導(dǎo)致另一個(gè)假象。 例如,由于這個(gè) 0.1 并非真正的 1/10,將三個(gè) 0.1 的值相加也不一定能恰好得到 0.3:

>>>
>>> .1 + .1 + .1 == .3
False

而且,由于這個(gè) 0.1 無(wú)法精確表示 1/10 的值而這個(gè) 0.3 也無(wú)法精確表示 3/10 的值,使用 round() 函數(shù)進(jìn)行預(yù)先舍入也是沒(méi)用的:

>>>
>>> round(.1, 1) + round(.1, 1) + round(.1, 1) == round(.3, 1)
False

雖然這些小數(shù)無(wú)法精確表示其所要代表的實(shí)際值,round() 函數(shù)還是可以用來(lái)“事后舍入”,使得實(shí)際的結(jié)果值可以做相互比較:

>>>
>>> round(.1 + .1 + .1, 10) == round(.3, 10)
True

Binary floating-point arithmetic holds many surprises like this. The problem with "0.1" is explained in precise detail below, in the "Representation Error" section. See The Perils of Floating Point for a more complete account of other common surprises.

正如那篇文章的結(jié)尾所言,“對(duì)此問(wèn)題并無(wú)簡(jiǎn)單的答案。” 但是也不必過(guò)于擔(dān)心浮點(diǎn)數(shù)的問(wèn)題! Python 浮點(diǎn)運(yùn)算中的錯(cuò)誤是從浮點(diǎn)運(yùn)算硬件繼承而來(lái),而在大多數(shù)機(jī)器上每次浮點(diǎn)運(yùn)算得到的 2**53 數(shù)碼位都會(huì)被作為 1 個(gè)整體來(lái)處理。 這對(duì)大多數(shù)任務(wù)來(lái)說(shuō)都已足夠,但你確實(shí)需要記住它并非十進(jìn)制算術(shù),且每次浮點(diǎn)運(yùn)算都可能會(huì)導(dǎo)致新的舍入錯(cuò)誤。

雖然病態(tài)的情況確實(shí)存在,但對(duì)于大多數(shù)正常的浮點(diǎn)運(yùn)算使用來(lái)說(shuō),你只需簡(jiǎn)單地將最終顯示的結(jié)果舍入為你期望的十進(jìn)制數(shù)值即可得到你期望的結(jié)果。 str() 通常已足夠,對(duì)于更精度的控制可參看 格式字符串語(yǔ)法str.format() 方法的格式描述符。

對(duì)于需要精確十進(jìn)制表示的使用場(chǎng)景,請(qǐng)嘗試使用 decimal 模塊,該模塊實(shí)現(xiàn)了適合會(huì)計(jì)應(yīng)用和高精度應(yīng)用的十進(jìn)制運(yùn)算。

另一種形式的精確運(yùn)算由 fractions 模塊提供支持,該模塊實(shí)現(xiàn)了基于有理數(shù)的算術(shù)運(yùn)算(因此可以精確表示像 1/3 這樣的數(shù)值)。

如果你是浮點(diǎn)運(yùn)算的重度用戶則你應(yīng)當(dāng)了解一下 NumPy 包以及由 SciPy 項(xiàng)目所提供的許多其他數(shù)字和統(tǒng)計(jì)運(yùn)算包。 參見(jiàn) <https://scipy.org>。

Python 也提供了一些工具,可以在你真的 想要 知道一個(gè)浮點(diǎn)數(shù)精確值的少數(shù)情況下提供幫助。 例如 float.as_integer_ratio() 方法會(huì)將浮點(diǎn)數(shù)表示為一個(gè)分?jǐn)?shù):

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

由于這是一個(gè)精確的比值,它可以被用來(lái)無(wú)損地重建原始值:

>>>
>>> x == 3537115888337719 / 1125899906842624
True

float.hex() 方法會(huì)以十六進(jìn)制(以 16 為基數(shù))來(lái)表示浮點(diǎn)數(shù),同樣能給出保存在你的計(jì)算機(jī)中的精確值:

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

這種精確的十六進(jìn)制表示法可被用來(lái)精確地重建浮點(diǎn)值:

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

由于這種表示法是精確的,它適用于跨越不同版本(平臺(tái)無(wú)關(guān))的 Python 移植數(shù)值,以及與支持相同格式的其他語(yǔ)言(例如 Java 和 C99)交換數(shù)據(jù).

另一個(gè)有用的工具是 math.fsum() 函數(shù),它有助于減少求和過(guò)程中的精度損失。 它會(huì)在數(shù)值被添加到總計(jì)值的時(shí)候跟蹤“丟失的位”。 這可以很好地保持總計(jì)值的精確度, 使得錯(cuò)誤不會(huì)積累到能影響結(jié)果總數(shù)的程度:

>>>
>>> sum([0.1] * 10) == 1.0
False
>>> math.fsum([0.1] * 10) == 1.0
True

15.1. 表示性錯(cuò)誤?

本小節(jié)將詳細(xì)解釋 "0.1" 的例子,并說(shuō)明你可以怎樣親自對(duì)此類情況進(jìn)行精確分析。 假定前提是已基本熟悉二進(jìn)制浮點(diǎn)表示法。

表示性錯(cuò)誤 是指某些(其實(shí)是大多數(shù))十進(jìn)制小數(shù)無(wú)法以二進(jìn)制(以 2 為基數(shù)的計(jì)數(shù)制)精確表示這一事實(shí)造成的錯(cuò)誤。 這就是為什么 Python(或者 Perl、C、C++、Java、Fortran 以及許多其他語(yǔ)言)經(jīng)常不會(huì)顯示你所期待的精確十進(jìn)制數(shù)值的主要原因。

為什么會(huì)這樣? 1/10 是無(wú)法用二進(jìn)制小數(shù)精確表示的。 目前(2000年11月)幾乎所有使用 IEEE-754 浮點(diǎn)運(yùn)算標(biāo)準(zhǔn)的機(jī)器以及幾乎所有系統(tǒng)平臺(tái)都會(huì)將 Python 浮點(diǎn)數(shù)映射為 IEEE-754 “雙精度類型”。 754 雙精度類型包含 53 位精度,因此在輸入時(shí),計(jì)算會(huì)盡量將 0.1 轉(zhuǎn)換為以 J/2**N 形式所能表示的最接近分?jǐn)?shù),其中 J 為恰好包含 53 個(gè)二進(jìn)制位的整數(shù)。 重新將

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

寫(xiě)為

J ~= 2**N / 10

并且由于 J 恰好有 53 位 (即 >= 2**52< 2**53),N 的最佳值為 56:

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

也就是說(shuō),56 是唯一的 N 值能令 J 恰好有 53 位。 這樣 J 的最佳可能值就是經(jīng)過(guò)舍入的商:

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

由于余數(shù)超過(guò) 10 的一半,最佳近似值可通過(guò)四舍五入獲得:

>>>
>>> q+1
7205759403792794

這樣在 754 雙精度下 1/10 的最佳近似值為:

7205759403792794 / 2 ** 56

分子和分母都除以二則結(jié)果小數(shù)為:

3602879701896397 / 2 ** 55

請(qǐng)注意由于我們做了向上舍入,這個(gè)結(jié)果實(shí)際上略大于 1/10;如果我們沒(méi)有向上舍入,則商將會(huì)略小于 1/10。 但無(wú)論如何它都不會(huì)是 精確的 1/10!

因此計(jì)算永遠(yuǎn)不會(huì)“看到”1/10:它實(shí)際看到的就是上面所給出的小數(shù),它所能達(dá)到的最佳 754 雙精度近似值:

>>>
>>> 0.1 * 2 ** 55
3602879701896397.0

如果我們將該小數(shù)乘以 10**55,我們可以看到該值輸出為 55 位的十進(jìn)制數(shù):

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

這意味著存儲(chǔ)在計(jì)算機(jī)中的確切數(shù)值等于十進(jìn)制數(shù)值 0.1000000000000000055511151231257827021181583404541015625。 許多語(yǔ)言(包括較舊版本的 Python)都不會(huì)顯示這個(gè)完整的十進(jìn)制數(shù)值,而是將結(jié)果舍入為 17 位有效數(shù)字:

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

fractionsdecimal 模塊可令進(jìn)行此類計(jì)算更加容易:

>>>
>>> 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'