我認為我應該添加一個硬體設計師的視角,因為我設計和建立浮點硬體。了解錯誤的起源可能有助於理解軟體中發生的情況,最終,我希望這能夠解釋浮點錯誤發生的原因,並且似乎隨著時間的推移而累積。
從工程角度來看,大多數浮點運算都會存在一些誤差,因為執行浮點計算的硬體只需要在最後一位上有不到半個單位的誤差。因此,大多數硬體只會停留在精確度上,只需產生小於最後一位的半個單位的誤差,這在浮點除法中尤其成問題。什麼構成了一個單獨的操作取決於單元接受的操作數數量。對於大多數單元來說,是兩個操作數,但有些單元接受3個或更多的操作數。因此,無法保證重複的操作會產生理想的誤差,因為隨著時間的推移,誤差會累積。
大多數處理器遵循IEEE-754標準,但有些使用非規格化或不同的標準。例如,IEEE-754中有一種非規格化模式,允許以犧牲精度為代價表示非常小的浮點數。然而,以下的內容將涵蓋IEEE-754的規格化模式,這是典型的操作模式。
在IEEE-754標準中,硬體設計師可以允許任何誤差/ε值,只要它小於最後一位的半個單位,並且結果僅對於一個操作必須小於最後一位的半個單位。這解釋了為什麼在重複操作時,誤差會累積。對於IEEE-754雙精度,這是第54位,因為有53位用來表示浮點數的數值部分(規格化),也稱為尾數(例如5.3e5中的5.3)。接下來的幾節將更詳細地討論各種浮點運算中硬體錯誤的原因。
浮點除法中誤差的主要原因是用來計算商的除法演算法。大多數電腦系統使用乘法逆元進行除法計算,主要是在Z=X/Y和Z = X * (1/Y)。除法是迭代計算的,即每個週期計算一些商的位數,直到達到所需的精度,對於IEEE-754來說,這是任何誤差小於最後一位的一個單位的情況。 Y的倒數表(1/Y)被稱為慢除法中的商選擇表(QST),商選擇表的位數通常是基數的寬度,或在每次迭代中計算的商的位數加上幾個保護位。對於IEEE-754標準的雙精度(64位元),它將是除法器的基數大小加上幾個保護位元k,其中k>=2。因此,例如,一個典型的每次計算2位商(基數為4)的除法器的商選擇表將是2 2= 4位元(加上一些可選位元)。
Z=X/Y
Z = X * (1/Y)
k>=2
2 2= 4
3.1 除法捨入誤差:倒數的近似
商選擇表中的倒數取決於除法方法:慢除法(如SRT除法)或快除法(如Goldschmidt除法);每個條目根據除法演算法進行修改,以盡可能產生最低的誤差。無論如何,所有倒數都是實際倒數的近似,並引入一定程度的誤差。慢除法和快除法方法都是迭代計算商,即每步計算一定數量的商位數,然後從被除數中減去結果,除數重複這些步驟,直到誤差小於最後一位的半個單位。慢除法方法在每步計算固定數量的商位數,並且通常更便宜,而快除法方法在每步計算可變數量的位數,並且通常更昂貴。除法方法最重要的部分是它們大多依賴對倒數的近似的重複乘法,因此容易出錯。
所有操作中捨入誤差的另一個原因是IEEE-754允許的不同截斷模式。有截斷、向零捨入、四捨五入(預設)、向下捨入和向上捨入。所有方法都會引入一個小於最後一位的單位的誤差,對於單一操作而言。隨著時間和重複操作的進行,截斷也會累積到結果誤差中。這種截斷誤差在指數運算中尤其成問題,其中涉及某種形式的重複乘法。
由於執行浮點計算的硬體只需要在單一操作中產生一個小於最後一位的半個單位的誤差的結果,如果不加以監控,隨著重複操作,誤差將會增長。這就是為什麼在需要有界誤差的計算中,數學家使用方法,例如使用IEEE-754中的四捨五入最後一位的偶數位,因為隨著時間的推移,誤差更有可能相互抵消,以及結合變化的區間算術和IEEE 754舍入模式的變化來預測舍入誤差並修正它們。由於相對誤差較低,與其他舍入模式相比,四捨五入到最近的偶數位(最後一位)是IEEE-754的預設舍入模式。
請注意,預設的捨入模式,四捨五入到最近的偶數位元,保證了一個操作的誤差小於最後一位的半個單位。僅使用截斷、向上舍入和向下
二進位浮點數數學就是這樣的。在大多數程式語言中,它是基於IEEE 754標準。問題的關鍵在於,數字在這種格式中表示為一個整數乘以二的冪;分母不是二的冪的有理數(例如0.1,即1/10)無法精確表示。
0.1
1/10
對於標準的binary64格式中的0.1,其表示可以精確寫成
binary64
0.1000000000000000055511151231257827021181583404541015625
0x1.9999999999999ap-4
相較之下,有理數0.1,即1/10,可以精確寫成
0x1.99999999999999...p-4
...
你的程式中的常數0.2和0.3也將近似它們的真實值。恰好,最接近0.2的double大於有理數0.2,但最接近0.3的double小於有理數0.3。 0.1和0.2的和最終比有理數0.3大,因此與程式碼中的常數不一致。
0.2
0.3
double
關於浮點數算術問題的一個相當全面的處理是每個電腦科學家都應該了解的浮點數算術知識。更容易理解的解釋,請參閱floating-point-gui.de。
普通的十進制(基數10)數字也存在同樣的問題,這就是為什麼像1/3這樣的數字最終變成了0.333333333...
你剛好遇到了一個(3/10)在十進制系統中很容易表示的數字,但在二進制系統中無法表示。反過來也是一樣(在某種程度上):1/16在十進制中是一個不好看的數字(0.0625),但在二進制中它看起來和十進制中的萬分之一一樣整齊(0.0001)* * - 如果我們習慣在日常生活中使用二進制數係統,你甚至會看著那個數字並本能地理解你可以透過不斷地除以2來得到它。
當然,這並不是浮點數在記憶體中儲存的方式(它們使用一種科學計數法的形式)。然而,它確實說明了一個問題,即二進制浮點精度錯誤往往會出現,因為我們通常感興趣的“真實世界”數字往往是十的冪 - 但這僅僅是因為我們日常使用十進制數係統。這也是為什麼我們會說71%而不是「每7個中的5個」(71%是一個近似值,因為5/7無法用任何十進制數精確表示)。
所以不,二進位浮點數並沒有出錯,它們只是像其他任何基數N的數系統一樣不完美 :)
在實務中,這個精確度問題意味著你需要使用舍入函數將浮點數舍入到你感興趣的小數位數,然後再顯示它們。
你還需要用允許一定容差的比較來取代等號測試,這表示:
不要使用if (x == y) { ... }
if (x == y) { ... }
而是用if (abs(x - y) < myToleranceValue) { ... }。
if (abs(x - y) < myToleranceValue) { ... }
其中abs是絕對值函數。 myToleranceValue需要根據你的特定應用選擇,它與你準備允許多少「搖擺餘地」以及你將要比較的最大數字有很大關係(由於精確度損失問題)。要注意你所選的語言中的“epsilon”風格常數。這些常數可以用作容差值,但它們的有效性取決於你正在處理的數字的大小,因為與大數計算可能超過epsilon閾值。
abs
myToleranceValue
硬體設計師的視角
#我認為我應該添加一個硬體設計師的視角,因為我設計和建立浮點硬體。了解錯誤的起源可能有助於理解軟體中發生的情況,最終,我希望這能夠解釋浮點錯誤發生的原因,並且似乎隨著時間的推移而累積。
1. 概述
從工程角度來看,大多數浮點運算都會存在一些誤差,因為執行浮點計算的硬體只需要在最後一位上有不到半個單位的誤差。因此,大多數硬體只會停留在精確度上,只需產生小於最後一位的半個單位的誤差,這在浮點除法中尤其成問題。什麼構成了一個單獨的操作取決於單元接受的操作數數量。對於大多數單元來說,是兩個操作數,但有些單元接受3個或更多的操作數。因此,無法保證重複的操作會產生理想的誤差,因為隨著時間的推移,誤差會累積。
2. 標準
大多數處理器遵循IEEE-754標準,但有些使用非規格化或不同的標準。例如,IEEE-754中有一種非規格化模式,允許以犧牲精度為代價表示非常小的浮點數。然而,以下的內容將涵蓋IEEE-754的規格化模式,這是典型的操作模式。
在IEEE-754標準中,硬體設計師可以允許任何誤差/ε值,只要它小於最後一位的半個單位,並且結果僅對於一個操作必須小於最後一位的半個單位。這解釋了為什麼在重複操作時,誤差會累積。對於IEEE-754雙精度,這是第54位,因為有53位用來表示浮點數的數值部分(規格化),也稱為尾數(例如5.3e5中的5.3)。接下來的幾節將更詳細地討論各種浮點運算中硬體錯誤的原因。
3. 除法的捨入誤差的原因
浮點除法中誤差的主要原因是用來計算商的除法演算法。大多數電腦系統使用乘法逆元進行除法計算,主要是在
Z=X/Y
和Z = X * (1/Y)
。除法是迭代計算的,即每個週期計算一些商的位數,直到達到所需的精度,對於IEEE-754來說,這是任何誤差小於最後一位的一個單位的情況。 Y的倒數表(1/Y)被稱為慢除法中的商選擇表(QST),商選擇表的位數通常是基數的寬度,或在每次迭代中計算的商的位數加上幾個保護位。對於IEEE-754標準的雙精度(64位元),它將是除法器的基數大小加上幾個保護位元k,其中k>=2
。因此,例如,一個典型的每次計算2位商(基數為4)的除法器的商選擇表將是2 2= 4
位元(加上一些可選位元)。3.1 除法捨入誤差:倒數的近似
商選擇表中的倒數取決於除法方法:慢除法(如SRT除法)或快除法(如Goldschmidt除法);每個條目根據除法演算法進行修改,以盡可能產生最低的誤差。無論如何,所有倒數都是實際倒數的近似,並引入一定程度的誤差。慢除法和快除法方法都是迭代計算商,即每步計算一定數量的商位數,然後從被除數中減去結果,除數重複這些步驟,直到誤差小於最後一位的半個單位。慢除法方法在每步計算固定數量的商位數,並且通常更便宜,而快除法方法在每步計算可變數量的位數,並且通常更昂貴。除法方法最重要的部分是它們大多依賴對倒數的近似的重複乘法,因此容易出錯。
4. 其他操作中的捨入誤差:截斷
所有操作中捨入誤差的另一個原因是IEEE-754允許的不同截斷模式。有截斷、向零捨入、四捨五入(預設)、向下捨入和向上捨入。所有方法都會引入一個小於最後一位的單位的誤差,對於單一操作而言。隨著時間和重複操作的進行,截斷也會累積到結果誤差中。這種截斷誤差在指數運算中尤其成問題,其中涉及某種形式的重複乘法。
5. 重複操作
由於執行浮點計算的硬體只需要在單一操作中產生一個小於最後一位的半個單位的誤差的結果,如果不加以監控,隨著重複操作,誤差將會增長。這就是為什麼在需要有界誤差的計算中,數學家使用方法,例如使用IEEE-754中的四捨五入最後一位的偶數位,因為隨著時間的推移,誤差更有可能相互抵消,以及結合變化的區間算術和IEEE 754舍入模式的變化來預測舍入誤差並修正它們。由於相對誤差較低,與其他舍入模式相比,四捨五入到最近的偶數位(最後一位)是IEEE-754的預設舍入模式。
請注意,預設的捨入模式,四捨五入到最近的偶數位元,保證了一個操作的誤差小於最後一位的半個單位。僅使用截斷、向上舍入和向下
二進位浮點數數學就是這樣的。在大多數程式語言中,它是基於IEEE 754標準。問題的關鍵在於,數字在這種格式中表示為一個整數乘以二的冪;分母不是二的冪的有理數(例如
0.1
,即1/10
)無法精確表示。對於標準的
binary64
格式中的0.1
,其表示可以精確寫成0.1000000000000000055511151231257827021181583404541015625
0x1.9999999999999ap-4
相較之下,有理數
0.1
,即1/10
,可以精確寫成0.1
0x1.99999999999999...p-4
,其中...
表示無盡的9的序列。你的程式中的常數
0.2
和0.3
也將近似它們的真實值。恰好,最接近0.2
的double
大於有理數0.2
,但最接近0.3
的double
小於有理數0.3
。0.1
和0.2
的和最終比有理數0.3
大,因此與程式碼中的常數不一致。關於浮點數算術問題的一個相當全面的處理是每個電腦科學家都應該了解的浮點數算術知識。更容易理解的解釋,請參閱floating-point-gui.de。
普通的十進制(基數10)數字也存在同樣的問題,這就是為什麼像1/3這樣的數字最終變成了0.333333333...
你剛好遇到了一個(3/10)在十進制系統中很容易表示的數字,但在二進制系統中無法表示。反過來也是一樣(在某種程度上):1/16在十進制中是一個不好看的數字(0.0625),但在二進制中它看起來和十進制中的萬分之一一樣整齊(0.0001)* * - 如果我們習慣在日常生活中使用二進制數係統,你甚至會看著那個數字並本能地理解你可以透過不斷地除以2來得到它。
當然,這並不是浮點數在記憶體中儲存的方式(它們使用一種科學計數法的形式)。然而,它確實說明了一個問題,即二進制浮點精度錯誤往往會出現,因為我們通常感興趣的“真實世界”數字往往是十的冪 - 但這僅僅是因為我們日常使用十進制數係統。這也是為什麼我們會說71%而不是「每7個中的5個」(71%是一個近似值,因為5/7無法用任何十進制數精確表示)。
所以不,二進位浮點數並沒有出錯,它們只是像其他任何基數N的數系統一樣不完美 :)
在實務中,這個精確度問題意味著你需要使用舍入函數將浮點數舍入到你感興趣的小數位數,然後再顯示它們。
你還需要用允許一定容差的比較來取代等號測試,這表示:
不要使用
#if (x == y) { ... }
而是用
if (abs(x - y) < myToleranceValue) { ... }
。其中
abs
是絕對值函數。myToleranceValue
需要根據你的特定應用選擇,它與你準備允許多少「搖擺餘地」以及你將要比較的最大數字有很大關係(由於精確度損失問題)。要注意你所選的語言中的“epsilon”風格常數。這些常數可以用作容差值,但它們的有效性取決於你正在處理的數字的大小,因為與大數計算可能超過epsilon閾值。