double - 及其用法

評分：3.1/5 (8 票)

本文介紹 "double" 算術的一些特性，以及你可以和不可以做什麼。
注意：需要 UTF-8 支援數學符號。

0. 動機
=============
在進行浮點運算時，可能會出現一些奇怪的效果，例如

for(double i=0.0; i!=5.5; i+=0.5)
{
        ... 
}

由於舍入誤差，可能永遠不會終止。特別是測試 (不) 等值通常不會得到期望的結果。一個常見的解決方法是使用一個 epsilon

bool is_equal(double d1, double d2) 
{
        if(abs(d1-d2)<epsilon)
                return true;
        return false;
}

然而，這個 "解決方案" 只在少數情況下適用 - 這樣定義的等價性不再是等價關係（失去了傳遞性：a-b<ε 且 b-c<ε 並不能推匯出 a-c<ε - 哎呀）。另一個問題是：epsilon 應該有多大？

計算機應該計算出正確的結果，或者說明它的不確定性。使用 epsilon 通常無法實現這一點。這就是為什麼你應該對本文的其餘部分感興趣 ;-)

1. 數值型別
==================
IEEE 754 double 值由以下幾部分組成：
1 位符號 (+ 或 -)
11 位指數（範圍：0 到 1047，0 和 1047 是 "保留" 的）
52 位有效數字（也稱為 "尾數"）

則表示的數字為

[符號] 2^{指數-偏差[=1023]}*1[有效數字]

其中 1 是一個隱含位，在有效數字位開始之前設定為 1。

對於歸一化數字，有效數字的第一位是 1 且不儲存，因此有效數字擴充套件到 53 位（顯然，前導零可以省略，因此數字以非零開頭，在二進位制系統中碰巧總是 1）。

當且僅當 0 < 指數 < 2047 時，該值才被歸一化。如果指數為零，IEEE 使用 "非歸一化" 數字，則值為 [符號]2^-1022*[有效數字]，但不包含有效數字中的前導 1。請注意，在非歸一化數字範圍內的計算通常非常慢（指數中的 2047 保留給 "NaN (非數字)、±inf 及類似項"）。

如果你出於某種原因需要位訪問，可以在 *大端序* 機器上使用

union ieee_double
{
   double   a;
   struct { unsigned s : 1;
            unsigned e :11;
            unsigned h :20;
            unsigned l :32;
          } b;
};

小端序需要單獨處理，但方式類似。

2. 操作
=============
假設舍入設定為 "舍入到最近" 和 "偶數舍入斷點"（預設）

設 ∘ 表示 {+,-,*,/} 中的一個運算，設 ⊚ 表示相應的機器運算，即 {⊕,⊖,⊛,⊘} 中的一個。

設 x、y 和 z 為 IEEE 754 double 值。
設 a 為實數。
設 fl(a) 為最接近 a 的 IEEE 754 double 值。

IEEE 754 現在保證
(a) fl(x ∘ y) = x ⊚ y
推論
x ⊖ y = 0 ⇔ x = y
x ⊕ y = y ⊕ x
x ⊛ y = y ⊛ x
(b) fl(sqrt(x)) = msqrt(x)，其中 msqrt 是 "機器 sqrt" 操作。
此外，還定義了一些 "非數字" 結果，例如 NaN（非數字）和 +/- 無窮大。

然而，在一般情況下，由於舍入誤差
a ⊚ ( b ⊚ c ) ≠ ( a ⊚ b ) ⊚ c

3. Sterbenz 引理
=================
在進行 double 運算時，你最重要的朋友之一是 Sterbenz 引理，它可以表述為：
如果 y/2 ≤ x ≤ y 則
x ⊖ y = x - y
如果 x ⊖ y 不下溢

或者，用 "普通英語" 表述：如果 y 在 x/2 和 x 之間，那麼機器減法 x ⊖ y 就會給你精確的結果。

4. 誤差界
===============
給定 (1) 和 (2) 中的性質以及由此產生的引理 (3)，可以為使用 +,-,*,/ 和 'sqrt' 的任何公式計算誤差界。然而，這種計算本身很繁瑣且容易出錯，因此通常使用近似值。

利用 Sterbenz 引理，可以證明一些更簡單的多項式求值誤差界。我不會在這裡展示它們，但你可以搜尋例如：
- Fortune-van-Wyk 界
- Mehlhorn-Näher 界
- Mehlhorn-Neumaier 界
來找到它們。（注意：據我所知，所有這些論文都僅證明了它們對於整數運算元（在 doubles 中表示）的正確性。然而，它們對任何 double 輸入都正確。）所有這些公式一個值得注意的美妙之處在於，實際的誤差界可以使用 double 算術來計算。

5. 浮點過濾器
=========================
記住動機示例：等值測試。給定兩個多項式 P 和 Q，現在可以在某些情況下確定 P(x) = Q(x)。

// equality test:
if(abs(P(x1) - Q(x2)) > ErrorBound(max(x1,x2)))
        return false; /* we know that for sure */
else
        // ???

（請注意，如果 P 或 Q 是多元的，max 最好使用 x1 和 x2 的上確界範數...）

如果你在編譯時（或程式啟動時）知道運算元的大小並且可以預先計算它，那麼在這種情況下，你只需要進行一次簡單的測試即可獲得精確性。
問題：'???' 是什麼？有幾種選擇。你可以使用任意精度數字型別的精確算術，你可以使用浮點展開，或者你可以執行惰性自適應求值（請注意，速度從左到右增加，但實現難度（以及想出正確公式的難度）也隨之增加）。
這些方法超出了本文的範圍（如果時間允許，我可能會寫一篇續篇 ;-））。
在某些領域，另一種可能性是使用受控擾動，這意味著輸入被轉換成不會出現等值且可以顯示不等值的情況（不精確，但可能魯棒且足夠，因為我們保證了某種屬性（即，我們不僅僅說 "不相等"，我們還確保它不是）。

6. 結論
=============
雖然計算誤差界需要額外的工作，但執行時效率是一個小問題（在大多數情況下，只需一次比較即可證明結果的正確性）。同時，由於錯誤地假設等值而導致的許多退化情況可以避免，通常可以顯著提高效能。

________________________________________________________________________________
如果您有任何評論、建議或改進，或者發現了拼寫錯誤（或多個），請告知我。

C++

文章

double - 及其用法