留意浮點數運算的陷阱

資料整理: jserv, ranvd

科學、工程與統計等領域往往離不開浮點數的運用，因此，如何趨吉避凶，正是每位經常處理浮點數的程式開發者應具備的基本素養。儘管妥善運用浮點數是門深奧的學問，學校課程中對此著墨不多，頂多在數值方法這門課中略有涉及，其涵蓋的廣度與深度皆有限。本文探討操作浮點數的基本觀念，不涉及理論分析，並使用直覺的詞彙，刻意避開過於嚴謹或艱澀的定義，以降低理解門檻，閱讀本文所需的預備知識僅限於基礎的程式設計概念，例如單精度與倍精度浮點數的表示方式，預期多數讀者皆可順利閱讀。

本文改寫自冼鏡光教授的文章〈使用浮點數最基本的觀念〉，將原文中的 FORTRAN 程式碼改為 C 語言，並對數學式進行重新排版、調整用語和表達，並補充相關資訊。

浮點數的最根本問題

使用浮點數時，面臨二個最根本的問題：

1. 輸入與儲存的數值不一定精確
2. 計算結果會產生誤差

第一點的成因主要包括以下幾項：輸入值本身就可能不夠精確、計算機儲存浮點數所用的記憶體有限、不同底數之間的轉換本質上難以精確、溢位問題，及數學常見的運算定律在浮點數運算時未必成立。至於第二點，因為浮點數運算本身會受到上述因素及其他誤差的影響，這些誤差經過反覆且大量的運算後，往往可能擴散至程式的其他部分，導致難以預料的結果，甚至造成無法收拾的後果。

以下各節將先介紹第一點的成因，稍後再深入探討第二點所衍生的問題。

輸入值本身不精確

因為實數可能具有無限位數，因此無法輸入其所有位數，計算機接收到的浮點數自然也就無法完全精確。例如，像 $\frac{1}{3}$, $\frac{1}{7}$, $\frac{1}{13}$ 等有理數，其小數表示為無限循環小數；而無理數如 $\sqrt{2}$、$\sqrt{3}$, 圓周率 $\pi =3.1415926\dots$ 等，更是具有無限不循環小數。因此，這些數值無法被完整地輸入到計算機中，進而不可避免地產生誤差。

圖一展示一個頂點分別位於 $(n,0)$, $(-n,0)$ 與 $(0,\sqrt{3}n)$ 的等邊三角形。在數學上，這樣的三角形座標表示完全沒有困難，但由於 $\sqrt{3}$ 是無理數、具有無限位數，因此在輸入 $\sqrt{3}$ 時無法做到精確。結果是，儲存在計算機中的三角形將會是個非常接近等邊三角形的等腰三角形，而非真正的等邊三角形。因此，若在程式中基於等邊三角形的假設進行後續運算，就極可能導致不精確的結果。

記憶體對浮點數儲存的限制

在計算機中，浮點數通常以標準型式儲存，類似於十進位的科學記號，例如：$\pm 0.xxx...xxx\times {10}^{\pm yyy}$。換句話說，實數 $-123.456$ 會被轉換為 $-0.123456\times {10}^3$，而 $0.00123456$ 會被表示為 $0.123456\times {10}^{-2}$。

在計算機內部，小數部份 ($\pm 0.xxx...xxx$) 與指數部份 ($\pm yyy$) 是分開儲存的，如圖二所示：

圖中顯示，小數部份與其符號獨立儲存，而指數部份則額外存放。至於指數部份的符號儲存方式，則取決於硬體設計。為了避免過於技術性的討論，此處不深入探討 IEEE 浮點數標準。

由於小數部份與指數部份的儲存空間皆有限 (以有限個位元表示)，因此會產生二個主要問題：

1. 小數部份無法精確表示具有無限位數的實數 (例如 $\sqrt{2}$、$1/3$、圓周率 $\pi$ 等)，因此會產生誤差
2. 儲存指數部份的位數通常比儲存小數部份的少得多，若浮點數的指數超出可儲存範圍，便會發生溢位 (overflow) 或下溢位 (underflow)：
   • 上溢位 (overflow)：若浮點數的指數超過最大可儲存值，則無法表示。例如，若硬體能儲存的指數範圍為 $[-63,+63]$，則 $0.1234\times {10}^{68}$ 會導致上溢位。
   • 下溢位 (underflow)：若浮點數的指數小於最小可儲存值，則無法表示。例如，$0.1234\times {10}^{-75}$ 會導致下溢位。

在計算過程中，若發生上溢位或下溢位，計算機可能會停止運算並顯示 floating-point exception 或類似錯誤訊息。有些系統甚至會具體標示錯誤類型，如上溢位、下溢位，或除以 0 的錯誤。

數系轉換不精確

我們日常用十進位數，但計算機的浮點數運算與儲存通常採用二進位或十六進位。因此，程式中使用的十進位數必須轉換為二進位或十六進位，然而，即使是簡單的十進位數，轉換後可能會變成無限位的二進位或十六進位數。

圖三展示將十進位數 $0.1$ 轉換為二進位 (左方) 與十六進位 (右方) 的計算方式。其過程是：將給定的十進位數乘以 $2$ 或 $16$，取出小數點前的數字（作為結果的一部分），然後對剩餘的小數部分繼續進行相同運算。例如，在左列中，$0.1$ 依次乘以 $2$，產生的結果為 $0,0,0,1,1,0$，接下來的小數部分變為 $0.4$，回到先前的計算狀態，形成循環。因此，十進位數 $0.1$ 在二進位中表示為循環小數 0.0_0011_0011_0011_0011...，而在十六進位中則表示為 0.1_9999_9999_9...。這說明，即使是簡單的十進位數 $0.1$，計算機也無法精確儲存，只能儲存其近似值，這導致計算結果產生誤差。

若浮點數包含整數部分與小數部分，需要分開轉換，然後再組合成標準型式。例如，十進位數 $277.31$ 需拆分為 $277$ 和 $0.31$。

若將 $277.31$ 轉換為十六進位：

• $277$ 轉換為 115 (十六進位)
• $0.31$ 轉換為 0.4_F5C28_F5C28_F5C28 \dots，這是個無限位的十六進位循環小數

合併後，表示為：
$$115.4F5C28F5C28F5C28\dots {}_{16}$$

寫成標準型式則為：
$$0.1154F5C28F5C28F5C28\cdots \times {16}^3$$

值得注意的是，指數部分的底數為 $16$，存入計算機時的方式如圖四所示：

然而，計算機無法儲存無限位的循環小數，因此後半部分將被截斷。通常，在單精確度 (single precision) 情況下，32 位元的浮點數儲存方式可容納最多 6 個十六進位數字。如果僅採用截斷 (忽略四捨五入)，計算機只能儲存 1154F5 這一段，如圖五中黃色框所示：

換句話說，計算機實際儲存的是 $0.1154F5\times {16}^3$。因為 $0.1154F5\times {16}^3$ 等於 $115.4F5$，將 $115$ 和 $0.4F5$ 轉換回十進位，分別得到 $277$ 和 $0.309814453125$，最終結果為 $277.309814453126$，而非原始值 $277.31$。這表示，在程式中使用 $277.31$ 時，計算機實際上使用的是 $277.309814453126$，這將導致計算誤差。

溢位

前面提過，在浮點數的標準型中，若小數部分的長度超過硬體的界限，超出的位數會被四捨五入後存入記憶體，這將導致誤差。而當硬體無法儲存乘冪部分時，則會發生溢位；若乘冪值超過硬體可容許的範圍，會發生上溢位(overflow)，若乘冪值（負數）小於硬體可容許的範圍，則會發生下溢位(underflow)。

計算過程中若發生溢位，計算機會顯示 floating-point exception 訊息，然後停止執行。程式中的四捨五入誤差與溢位問題相當棘手，因為前者不易察覺，一旦影響結果的正確性，往往難以追蹤；後者則會導致程式直接終止，使得補救變得困難。因此，撰寫程式時應留意可能產生極大或極小值的情境，並適當調整寫法，以避免這類問題。

考慮這個簡單的例子，計算 $a$ 平方加 $b$ 平方後再取平方根：
$$c=\sqrt{a^2+b^2}$$

看似直觀，但若 $a$ 或 $b$ 其中一個數值很大，可能會帶來問題。例如，若 $a$ 的值很大，則 $a^2$ 會更大，甚至可能發生溢位。此外，即使 $a$ 和 $b$ 本身不算大，$a^2+b^2$ 的總和仍可能超出硬體可表示的範圍，導致溢位，儘管最終計算平方根的結果可能不會溢位。

為了降低溢位的風險，可以使用以下方法對 $a$ 和 $b$ 進行縮放處理：
$$c=k\times \sqrt{{\left(\frac{a}{k}\right)}^2+{\left(\frac{b}{k}\right)}^2}$$

簡單來說，先用某個 $k>0$ 來縮小 $a$ 和 $b$，再將平方根的結果乘回 $k$，即可得到與原始計算相同的結果。

那麼，該如何選擇合適的 $k$ 值呢？理想的選擇是 $k=max(|a|,|b|)$，也就是 $|a|$ 和 $|b|$ 中較大者。這樣一來，$a/k$ 和 $b/k$ 之中較大的一項會變為 $1$，較小的一項則小於 $1$。這樣不僅能避免上溢位，還能確保其平方和不超過 $2$，開平方後的結果介於 $1$ 與 $\sqrt{2}$ 之間。最終將平方根乘回 $k$ 時，發生溢位的風險也較低。

下方是幾何平均數的定義，所有 $x$ 值皆為正數。若直接根據此定義計算，當輸入值過大時，發生溢位的可能性相當高。為了避免此問題，可用前述方法進行處理。
$$\sqrt[n]{x_1\times x_2\times x_3\times \cdots \times x_n}$$

硬體僅會在浮點數發生溢位時產生 floating-point exception 訊息，而對於整數溢位通常不會有任何警告。例如，計算階乘的結果會迅速增長，在 32 位元的系統上，$12!$ 是能夠正確計算的上限，而 $13!$ 以上則會發生整數溢位，導致結果錯誤。

考慮以下階乘運算：

• 計算 $12!$ 時，結果為 $479001600$，最後二位數為 $00$
• 計算 $13!$ 時，結果的最後兩位數卻不再是 $00$，顯然不正確，因為 $13!=12!\times 13$，其最後兩位數理應仍為 $00$
• 進一步計算發現 $14!$ 竟然小於 $13!$，而 $19!$ 甚至比 $12!$ 還小，甚至出現負數結果

之所以有上述異常結果，肇因於計算 $13!$ 時發生溢位，使得後續計算的結果皆不正確。然而，計算機並不會因整數溢位而停止運行，甚至在某些應用(如產生整數隨機數) 中，還會刻意利用整數溢位的特性。

 1! =            1
 2! =            2
 3! =            6
 4! =           24
 5! =          120
 6! =          720
 7! =         5040
 8! =        40320
 9! =       362880
10! =      3628800
11! =     39916800
12! =    479001600
13! =   1932053504
14! =   1278945280
15! =   2004310016
16! =   2004189184
17! =   -288522240
18! =   -898433024
19! =    109641728
20! =  -2102132736

數學運算定律對浮點數運算不一定成立

最常見的數學運算定律有三個：

• 交換律 (commutative law)
  $$a+b=b+a,a\times b=b\times a$$
• 結合律 (associative law)
  $$(a+b)+c=a+(b+c),(a\times b)\times c=a\times (b\times c)$$
• 分配律 (distributive law)
  $$(a+b)\times c=a\times c+b\times c$$

然而，在浮點數運算中，這三者中僅有交換律仍然成立，結合律與分配律則可能出現問題。

假設浮點數運算僅保留二位有效數字，圖九左方說明結合律的不成立：

• 計算 $0.12\times 0.34$，手算結果為 $0.0408$，但由於僅保留兩位有效數字，需四捨五入，即 $0.0408=0.408\times {10}^{-1}\approx 0.41\times {10}^{-1}=0.041$ 產生四捨五入誤差
• 接著計算 $0.041\times 4=0.164$，四捨五入後得到 $0.16$

另一方面，若先計算 $0.34\times 4$：

• $0.34\times 4=1.36$，四捨五入後，得到 $1.4$
• 接著計算 $0.12\times 1.4=0.168$，四捨五入後，得到 $0.17 $

由於 $0.16\ne 0.17$，可見結合律在浮點數運算中不成立。

圖九右方展示分配律的問題。首先，計算 $0.12+0.99$，結果為 $1.11$，四捨五入後變為 $1.1$，再乘以 $5$ 得 $1.1\times 5=5.5$

另一方面，分別計算 $0.12\times 5$ 和 $0.99\times 5$：

• $0.12\times 5=0.6$
• $0.99\times 5=4.95$，四捨五入後為 $5.0$
• 相加得 $0.6+5.0=5.6$

由於 $5.5\ne 5.6$，可見分配律在浮點數運算中不成立。

此外，結合律還可能導致溢位。考慮以下算式：
$$(1\times {10}^{-63}\times 1\times {10}^{63})\times 1\times {10}^{63}=1\times {10}^{-63}\times (1\times {10}^{63}\times 1\times {10}^{63})$$

左側括號內的運算結果為 $1$，因此整體結果為 $1\times {10}^{63}$。然而，右側括號內的結果為 $1\times {10}^{126}$，在某些系統上，當指數達到 $126$ 時可能會發生上溢位。

上述案例說明浮點數運算不符合結合律，因此在撰寫運算式時應特別留意，避免潛在的誤差與溢位問題。

考慮一個實際的計算例子。下列五個運算式在數學上完全相同，它們只是藉由結合律與分配律改寫：
$$\begin{array}{cc}{x}_{i+1}& =(R+1)x_i-R(x_i{x}_i)\\ x_{i+1}& =(R+1)x_i-(R{x}_i)x_i\\ x_{i+1}& =((R+1)-R{x}_i)x_i\\ x_{i+1}& =R{x}_i+(1-R{x}_i)x_i\\ x_{i+1}& =x_i+R(x_i-x_i{x}_i)\end{array}$$

現在，我們撰寫程式，設定 $R=3.0$，並從 $x_0=0.5$ 開始計算，分別使用這 5 個運算式進行迭代。換句話說，對每個運算式，我們從 $x_0=0.5$ 開始計算 $x_1$，再將 $x_1$ 代回原式求得 $x_2$，依此類推。

為了節省空間，僅顯示 $x_{100}$, $x_{200}$, $x_{300}$, …, $x_{1000}$ 的結果。結果如下表所示，第 1 行為 $i$ 的值，第 2 至第 6 行分別對應上方 5 個運算式的計算結果。

仔細比較同一列的數值，不難發現五個運算式所得結果存在顯著差異，這清楚顯示當結合律與分配律不成立時，浮點數運算可能產生多麼嚴重的影響。

   0 0.5000000 0.5000000 0.5000000 0.5000000 0.5000000
 100 0.2907871 1.2975237 0.0494887 0.7379640 0.0497974
 200 0.4477620 1.3313184 0.1603770 0.2955562 0.0096703
 300 0.1246991 0.1883296 0.3310097 0.0102723 0.3746538
 400 1.2000864 0.3369606 1.2219012 1.3320501 0.4422436
 500 1.3012922 1.2063222 1.3024172 1.1021049 0.0557052
 600 0.8624897 0.8175101 0.1605954 1.1790663 0.0637164
 700 1.0700798 0.6316685 0.9783148 0.0618592 0.6890506
 800 1.2988597 0.5392919 1.1552609 1.2842493 0.7373490
 900 1.3333333 0.0690899 1.1314385 1.0307841 0.9684411
1000 1.3308059 0.3915348 1.1785927 0.0624183 0.1420892

浮點數比較

由於浮點數的儲存方式可能導致精度丟失，因此在比較或計算浮點數時，必須考慮精度誤差的累積及其影響。

初學者可能會用以下方式比較浮點數：

bool equal(float a, float b)
{
    return a == b;
}

然而，若變數 a 受到精度丟失的影響，則 a 和 b 可能不再相等。於是，浮點數之間的比較通常要引入額外的誤差容忍值 epsilon，如下：

/* the smallest difference between 1.0 and the next representable float */
#define FLT_EPSILON 1.19209290E-07F
bool equal(float a, float b)
{
    return fabs(a - b) <= FLT_EPSILON;
}

注意：隨著浮點數的數值變大，其精度相對降低，而 FLT_EPSILON 代表的是 $1.0$ 和下個可表示的浮點值之間的差異。因此，當數值大於 1 時，這種比較方式就不夠準確。

為了應對這個問題，可根據浮點數大小動態縮放誤差，對應的程式碼如下：

#include <math.h>
bool equal(float a, float b, float diff)
{
    return fabs(a - b) / fmin(fabs(a), fabs(b)) < diff;
}

但這種方法正確嗎？

考慮一元二次方程式 $x^2+200x-10.0025=0$，其實根為 $0.050003051757812$ 和 $-200.050003051757812$，考慮到單精度浮點數的有效位數，我們期望的根是 $0.05$ 和 $-200.05$。

若用 equal 函式比較 $0.05$ 和 $0.050003051757812$，則計算結果如下：
$$\frac{\text{fabs}(a-b)}{min(\text{fabs}(a),\text{fabs}(b))}$$

計算結果為 $0.0000610622508880615234375$，如果這個值大於 FLT_EPSILON，則比較結果為 false。

另一種比較方法是利用最小精確度單元（units of least precision, ULP），即知道某個浮點數與其相鄰浮點數的差值，再以 fabs(a - b) < ULP 判斷是否相等。

基於 IEEE 754 浮點數標準，可透過位元表示的數值差來計算兩個浮點數之間的 ULP 差異：

bool equal(float a, float b, int ulpsEpsilon)
{
    return ulps(a, b) < ulpsEpsilon;
}

int ulps(float a, float b)
{
    int32_t ia = *(int32_t *) &a;
    int32_t ib = *(int32_t *) &b;
    int32_t distance = ia - ib;
    return distance < 0 ? -distance : distance;
}

回到前述一元二次方程式的根，計算 $0.0500030517578125$ 和 $0.05$ 的 ULP 差異，結果為 819。這是因為 ULP 取決於浮點數的指數大小，當數值接近某些指數邊界時，ULP 的變化可能會影響比較結果。此外，較小數值範圍內的 ULP 間距較小，導致位元變動可能累積成較大的 ULP 差異。然而，在實際計算前，我們並不總是能確定應該選擇多少 ULP 作為誤差範圍。

綜合前述幾種比較方法，究竟何種方式最適合比較浮點數？其實，不存在單一通用的解法，選擇何種方法取決於具體應用情境。

首先，直接比較 return a == b; 不一定會出錯。只要不涉及類型轉換或運算，浮點數的直接比較仍然適用。例如，若將浮點數寫入檔案後再讀取並比較，只要寫入時逐位元儲存，就不會有問題。

其次，return fabs(a - b) <= FLT_EPSILON; 在 a 和 b 都小於 1 時是合理的。

至於 ULP 方法的比較結果相對精確，但其適用性取決於具體場景，特別是在接近指數邊界時，ULP 的變化可能會影響比較結果。此外，浮點數在經過一系列計算後可能產生不可預測的精度損失，因此難以確定適當的 ulpsEpsilon 值。

最後，相對誤差 (relative epsilon) 方法透過計算數值差異與數值本身的比值來進行比較。例如，$0.05$ 與 $0.0500030517578125$ 的 ULP 差異為 819，而其相對誤差為 $0.00006102025508880615234375$，可見相對於 $0.05$ 的誤差佔比極低。

上述方法尚未考慮浮點數的特殊值 (如 NaN 和無窮大等)。完整的比較函式應充分考慮，如下：

#include <math.h>
#include <stdbool.h>
#include <stdint.h>
#include <limits.h>

bool equal(float a, float b, int ulpsEpsilon)
{
    return ulps(a, b) < ulpsEpsilon;
}

int32_t ulps(float a, float b)
{
    /* 如果數值完全相等，返回 0 */
    if (a == b) return 0;

    const int32_t max = INT_MAX;

    /* 檢查 NaN 或無窮值 */
    if (isnan(a) || isnan(b)) return max;
    if (isinf(a) || isinf(b)) return max;

    int32_t ia = *(int32_t*)&a;
    int32_t ib = *(int32_t*)&b;

    /* 若符號不同，返回最大值 */
    if ((ia < 0) != (ib < 0)) return max;

    int32_t distance = ia - ib;
    return distance < 0 ? -distance : distance;
}

上方程式碼讓浮點數之間的比較更穩健。

延伸閱讀: Comparing Floating Point Numbers

再談四捨五入誤差的影響

前面提過，計算結果會經過四捨五入，因此不會完全精確。若這樣的結果被反覆用於後續計算，誤差可能會像滾雪球般累積，最終導致程式產生完全不合理的結果。

需要注意的是，單次四捨五入所造成的誤差通常不大，在最佳情況下，大約為最後一位數的 $\pm 50\mathrm{\%}$。舉例而言，假設某計算機保留小數點後三位有效數字：

• 若計算結果為 $0.1234$，則捨去後變為 $0.123$。
• 若結果為 $0.1235$，則進位得到 $0.124$。

反過來看，若數值顯示為 $0.123$，可能是由以下數值捨去而得：

• $0.1230$、$0.1231$、$0.1232$、$0.1233$ 或 $0.1234$。
• 也可能是由 $0.1225$、$0.1226$、$0.1227$、$0.1228$ 或 $0.1229$ 進位而來。

也就是說，$0.123$ 的實際值可能位於 $[0.1225,0.1235]$ 之間，即 $[0.123-0.0005,0.123+0.0005]=[0.1225,0.1235]$。

這正好對應最後一位數的 $\pm 50\mathrm{\%}$ 誤差範圍。

即使是如此細微的捨入誤差，也可能導致錯誤結果。例如，假設希望將區間 $[1,2]$ 平均分成三等分，並在 $x=1,4/3,5/3,2$ 時使用 $x$ 進行進一步計算。

以下是個 C 語言的程式，其中 a, b, x 和 h 均為至少具備 15 位有效數字的浮點變數（見 double 的定義）。這段程式看起來結構正確，設計嚴謹，乍看之下沒有任何明顯漏洞，對吧？

#include <stdio.h>
int main(void)
{
    const double a = 1.0, b = 2.0;
    int n = 3;
    double h = (b - a) / n;
    double x = a;

    printf("x = %lf\n", x);
    while (1) {
        if (x >= b) break;
        x = x + h;
        printf("x = %lf\n", x);
    }
    
    return 0;
}

實際情況並非如此，執行結果如下，出現一個額外的輸出值 $2.33333\dots$ :

x = 1.000000
x = 1.333333
x = 1.666667
x = 2.000000
x = 2.333333

從輸出結果可見，$x$ 由 $1$ 開始，依序經過 $1.333\dots$ (小於 $\frac{4}{3}$), $1.666\dots$ (小於 $\frac{5}{3}$), $1.999\dots$（小於 $2$），最終多執行一次迴圈才滿足 $x\ge 2$ 的條件，導致離開 while 迴圈。這正是捨入誤差所造成的問題。

在小型程式中，這類問題相對容易察覺，但在規模龐大且結構複雜的程式中，找出額外的一次迴圈可能會相當費時。

在多數情況下，藉由以下方式調整 if 條件判斷，可有效避免額外的迴圈執行問題：

if (x >= b - h / 2)
    break;

讓我們再看一個更長且極端的例子：假設有 $n>1$ 個數值 $x_1,x_2,\dots ,x_n$，它們的平均數 $m$ 定義為這些數值的總和除以 $n$，即：
$$m=\frac{x_1+x_2+\cdots +x_n}{n}$$

變異數（variance）則定義如下（式 (1)）：
$$\begin{array}{}\text{(1)}& S^2=\frac{1}{n-1}\sum _{i=1}^n(x_i-m{)}^2\end{array}$$

因此，計算變異數時，首先需要求出平均數 $m$，再計算各個 $(x_i-m{)}^2$ 的總和，最後再除以 $n-1$。許多初學者不喜歡這種方式，因為它需要使用兩個 for 迴圈 —— 第一個用來計算 $x_i$ 的總和，第二個則計算平方和。他們可能會在課本中找到另一種公式（式 (2)），看起來更為「簡潔」，因為它可以在單一 for 迴圈內同時計算 $x_i$ 的總和與平方和，效率上似乎更佳：
$$\begin{array}{}\text{(2)}& S^2=\frac{1}{n-1}[\sum _{i=1}^n{x}_i^2-\frac{1}{n}{\left(\sum _{i=1}^n{x}_i\right)}^2]\end{array}$$

然而，經驗豐富的軟體開發者可能會對這種計算方式有所顧慮。雖然 $x_i$ 的總和本身不一定會導致溢位，但若總和數值較大，計算其平方時可能會產生更大的數值，進而在儲存時造成尾數位元的損失，引發四捨五入誤差。同樣的問題也可能發生在 $x_i$ 的平方和計算過程中。

不相信嗎？我們來探討其衝擊。考慮 $n=3$ 的情境，每個數皆為 7 位數，接近 32 位元儲存範圍的極限，確保數值不會因為四捨五入而發生截斷。我們先手算：

x1 = 9000000
x2 = 9000001
x3 = 9000002

於是 $x_1+x_2+x_3=27000003$

其平均數為 $m=9000001$

若使用式 (1) 計算變異數：
$$\frac{(9000000-9000001{)}^2+(9000001-9000001{)}^2+(9000002-9000001{)}^2}{2}=1$$

若改用式 (2)，則計算如下：
$$\begin{array}{lcl}[1]& x_1+x_2+x_3& =27000003\\ [2]& (x_1+x_2+x_3{)}^2& =729000162000009\\ [3]& \frac{(x_1+x_2+x_3{)}^2}{3}& =243000054000003\\ [4]& x_1^2+x_2^2+x_3^2& =243000054000005\end{array}$$

變異數計算為：
$$\frac{[4]-[3]}{2}=1$$

然而，若僅有 7 位有效數字，[4] 的結果將會四捨五入為：
$$0.2430001\times {10}^{15}=243000100000000$$

同時，[3] 的結果可能也會捨入為相同的值，導致變異數計算結果為 $0$，這顯然是錯誤的。事實上，Microsoft Excel 的早期版本在進行此計算時，確實會產生 $0$ 的結果！

以下是使用二種公式計算變異數的 C 程式，所有變數皆採用約定的單精度 (32 位元) 運算：

#include <stdio.h>
#define SIZE 10
int main(void)
{
    int n, i;
    float x[SIZE];
    float sumx = 0.0f, sumx2 = 0.0f;
    float mean, var1, var2 = 0.0f;

    scanf("%d", &n);
    if (n > SIZE) {
        printf("Exceed size limit (%d)\n", SIZE);
        return 1;
    }
    for (i = 0; i < n; i++)
        scanf("%f", &x[i]);
    
    printf("Number of data items --> %d\n\n", n);
    for (i = 0; i < n; i++)
        printf(" x(%2d) = %25.10f\n", i+1, x[i]);
    
    for(i = 0; i < n; i++) {
        sumx  += x[i];
        sumx2 += x[i] * x[i];
    }
    printf("\n");
    printf(" Sum of x --------------> %30.10f\n", sumx);
    printf(" (Sum of x)**2 ---------> %30.10f\n", sumx * sumx);
    printf(" [(Sum of x)**2]/n -----> %30.10f\n", (sumx * sumx) / n);
    printf(" Sum of x**2 -----------> %30.10f\n", sumx2);
    
    mean = sumx / n;
    printf("\n");
    printf(" Mean ------------------> %30.10f\n", mean);
    
    if (n > 1) {
        var1 = (sumx2 - (sumx * sumx) / n) / (n - 1);
        printf(" One-Pass Variance -----> %30.10f\n", var1);
    } else {
        printf(" One-Pass Variance -----> Undefined (n must be > 1)\n");
    }
    
    if (n > 1) {
        for(i = 0; i < n; i++)
            var2 += (x[i] - mean) * (x[i] - mean);
        var2 = var2 / (n - 1);
        printf(" Two-Pass Variance -----> %30.10f\n", var2);
    } else {
        printf(" Two-Pass Variance -----> Undefined (n must be > 1)\n");
    }
    
    return 0;
}

若設 $n=3$，並以 $9000000$、$9000001$、$9000002$ 作為輸入，則程式的計算結果如下。其中，One-Pass 代表式 (2)（僅用一個 for 迴圈），而 Two-Pass 則為式 (1)（使用二個 for 迴圈）。

程式計算出的數值總和為 $27000002$，與實際結果相差 $1$。由於四捨五入的影響，平均數的計算仍然正確。然而，總和略有誤差，導致其平方 (共有 15 位數) 中僅前 6 位是完全正確的，其餘 9 位不可靠。因此，將其除以 $3$ 後，結果的可靠位數僅有 7 位。同理，平方和 $243000055693312$ 亦僅有前 8 位可靠。因此，當後者減去前者時，所得的 $16777216$ 幾乎沒有任何可靠的位數 (儘管數值看似很大)，這導致變異數的計算結果異常（$8388608$），顯然不可能正確。

Number of data items --> 3

 x( 1) =        9000000.0000000000
 x( 2) =        9000001.0000000000
 x( 3) =        9000002.0000000000

 Sum of x -------------->            27000002.0000000000
 (Sum of x)**2 --------->     729000099971072.0000000000
 [(Sum of x)**2]/n ----->     243000038916096.0000000000
 Sum of x**2 ----------->     243000055693312.0000000000

 Mean ------------------>             9000001.0000000000
 One-Pass Variance ----->             8388608.0000000000
 Two-Pass Variance ----->                   1.0000000000

或許會想到透過倍精度運算 (double precision) 來提高數值精確度，以解決此問題。然而，即使使用倍精度，也能找到三個足夠大的數，使計算仍然遭遇相同困境。因此，根本的解決方案應是選擇適當的演算法，而非單純提升精度。在變異數的計算問題中，最佳做法是採用上述 Two-Pass 方法，即要用到二個 for 迴圈的計算方式。

此外，前述計算平均數的過程需累加所有數值，這可能導致總和過大，超出計算機可容納範圍，從而產生四捨五入誤差。然而，平均數始終介於最大與最小值之間，因此，在能夠精確儲存最大與最小值的前提下，平均數亦該能被精確儲存。

事實上，計算平均數與變異數皆可透過「更新公式」來處理，該方法透過將目前計算出的平均數 (或變異數) 與新輸入的數值結合，以逐步更新結果，從而避免累積總和過大，並降低計算誤差。

失去精確位數

另一個令人頭痛的浮點數問題是，減法運算可能導致精度損失。當兩個數值極為接近時，相減後的結果可能會失去大量有效數字。

以下 C 程式中，x 和 y 均為 16 位數。在具有 15 位精確度的系統上，這兩個數最多只有最後兩位可能不精確。它們的真實差值（手算結果）為 $ x - y = 0.000000000004111$。因此，兩個原本具有 15 位有效數字的數值相減後，結果僅剩下 4 位有效數字，也就是說，當浮點數相減時，如果兩個數值非常接近，計算結果可能會嚴重喪失精度。

#include <stdio.h>
int main(void)
{
    double x = 3.141592653589793;
    double y = 3.141592653585682;
    double z = x - y;

    printf("x     = %23.17E\n", x);
    printf("y     = %23.17E\n", y);
    printf("x - y = %23.17E\n", z);

    return 0;
}

程式執行輸出如下，計算結果仍然保留 15 位，但其中僅小數點後的前三位數 (411) 可信：

x     = 3.14159265358979310E+00
y     = 3.14159265358568200E+00
x - y = 4.11100000000000000E-12

雖然 z 仍保留 15 位有效數字，但其中只有最前 3 到 4 位可靠，因為後續位數受到了四捨五入的影響。然而，在程式執行過程中，我們無法確定在哪一位發生四捨五入，因此結果雖然表面上具有 15 位 (符合倍精度的規格)，但實際上只有 3 到 4 位是可靠的。這種精度損失的主因，正是當兩個幾乎相等的浮點數相減時，大部分有效數字相互抵消，僅剩下極少數可用位數。前面討論計算變異數時，也曾提及這類精度喪失問題，當計算結果依賴這類減法運算時，可能會導致完全不可信的結果。

當兩個幾乎相等的數值相減時，這通常並非程式設計的本意，而是來自其他計算過程中未妥善處理的結果。正因如此，精度損失的根本原因往往不易察覺。以下是個典型的例子：求解一元二次方程式 $a{x}^2+bx+c=0$ 的根（見以下公式）：
$$x=\frac{-b\pm \sqrt{b^2-4ac}}{2a}$$

這段程式邏輯簡單，相信多數人在學習初級程式設計時都曾撰寫過。以下為 C 語言的實作範例，其中所有浮點變數皆採用單精度 (約 7 位有效數字)：

#include <stdio.h>
#include <math.h>

int main(void)
{
    float a, b, c;
    float d, r1, r2;

    scanf("%f %f %f", &a, &b, &c);

    printf(" a      = %15.7E\n", a);
    printf(" b      = %15.7E\n", b);
    printf(" c      = %15.7E\n", c);

    /* 計算判別式與根 */
    d = sqrt(b * b - 4.0f * a * c);
    r1 = (-b + d) / (2.0f * a);
    r2 = (-b - d) / (2.0f * a);

    printf("\n");
    printf(" d      = %15.7E\n", d);
    printf(" Root 1 = %15.7E\n", r1);
    printf(" Root 2 = %15.7E\n", r2);

    return 0;
}

這段程式正確嗎？是的。那麼，這段程式寫得好嗎？從初級程式設計的角度來看，該寫法不差，但從專業角度而言，仍然存在潛在問題。問題不是沒有考慮 $a=0$ 或 $b^2-4ac<0$ 等特殊情況，而是由於可能發生嚴重的精度損失，也就是分子中涉及減法運算。要避免 $b^2-4ac$ 計算中的減法相對困難，因此暫且不論，關鍵問題出現在平方根前的 $\pm$ 符號。

假設 $b$ 很大，而 $a$ 和 $c$ 都很小，則 $b^2$ 會更大，因此 $b^2-4ac$ 的結果幾乎等於 $b^2$，導致其平方根與 $b$ 非常接近，進而在減法運算時，會產生顯著的精度損失。

舉例而言，假設 $a=c=1$，$b=20000$，則
$$b^2=400000000,b^2-4ac=399999996=0.399999996\times {10}^9$$

若使用單精度運算 (7 位有效數字)，此數值會被四捨五入為 $0.4000000\times {10}^9$

因此，$\sqrt{b^2-4ac}$ 近似為 $20000$，導致計算出的兩個根為：
$$x_1=\frac{-20000+20000}{2}=0,x_2=\frac{-20000-20000}{2}=-20000$$

這顯然是不正確的結果！

根據二次方程式性質，其根的和應為 $-b/a$，根的積應為 $c/a$。然而，給定方程式 $x^2+20000x+1=0$ 的根不可能為 $0$，但程式卻計算出錯誤的結果。

以下是程式的實際執行結果：

 a      =   1.0000000E+00
 b      =   2.0000000E+04
 c      =   1.0000000E+00

 d      =   2.0000000E+04
 Root 1 =   0.0000000E+00
 Root 2 =  -2.0000000E+04

如何去除分子中的減法運算？要注意，二個根的積為 $c/a$，因此當已知其中一個根 (記為 $r_1$，假設不為 0，若根為 0 則問題會更簡單)，另一個根可表示為 $\frac{c/a}{r_1}$。

接著來思考：若 $b<0$，則 $-b$ 為正數，因此 $-b+\sqrt{b^2-4ac}$ 是兩個正數的和；反之，若 $b>0$，則 $-b$ 為負數，使得 $-b-\sqrt{b^2-4ac}$ 是兩個負數相加。如此一來，計算過程中不再涉及減法運算，從而避免精度損失。

以下程式片段採取該策略，進行修正：

    d = sqrt(b * b - 4.0f * a * c);
    if (b > 0)
        r1 = -b - d;
    else
        r1 = -b + d;

    r1 = r1 / (2.0f * a);
    r2 = (c/a) / r1;

對應的計算結果指出，二個根分別為 $-20000$ 與 $-0.00004999\dots$。然而，這二個根的和為 $-20000.00004999\dots$，與理論值 $-20000$ 略有偏差，結果仍稱不上完全準確。不過，就浮點數運算而言，在單精度 (7 位有效數字) 的條件下，$20000.00004999\dots$ 會被四捨五入為 $20000$，因此最終結果仍然是正確的。

雖然仍存在些許誤差，但相比於先前導致嚴重錯誤的計算方式，這樣的結果顯然要可靠得多。

 a      =   1.0000000E+00
 b      =   2.0000000E+04
 c      =   1.0000000E+00

 d      =   2.0000000E+04
 Root 1 =  -2.0000000E+04
 Root 2 =  -4.9999999E-05

去除減法運算並非易事，且非所有情況都能適用。不過，有個常見的技巧是利用以下式 (3) 進行轉換。當數式中出現 $a-b$ 時，可藉由乘上 $\frac{a+b}{a+b}$，將其轉換為 $\frac{a^2-b^2}{a+b}$。若 $a^2-b^2$ 能進一步化簡，使減法消失，於是問題可迎刃而解。然而，這種技巧並非萬能，若無法成功，則需考慮其他方法。

需要注意的是，只有當 $a$ 與 $b$ 的數值接近，可能導致精度損失時，才有必要進行這種轉換；若已知 $a$ 和 $b$ 相差足夠大，不至於發生精度損失，則該變換意義不大。
$$\begin{array}{}\text{(3)}& (a+b)(a-b)=a^2-b^2\end{array}$$

我們來看如何處理以下算式 (式 (4))：當 $x$ 很大時，它才會出現問題。由於四捨五入的影響，當 $x$ 很大時，$x+1$ 與 $x$ 的數值極為接近，因此 $\sqrt{x+1}$ 和 $\sqrt{x}$ 幾乎相等，甚至可能相等，這可能導致除數為零，進而引發「除以零」的錯誤。即便未發生除零錯誤，兩個幾乎相等的數相減仍會導致精度損失，使結果不夠準確。
$$\begin{array}{}\text{(4)}& \frac{1}{\sqrt{x+1}-\sqrt{x}}\end{array}$$

利用上述技巧，式 (4) 可轉換為 $\sqrt{x+1}+\sqrt{x}$，使得減法運算消失，也就是式 (5)。
$$\begin{array}{}\text{(5)}& \frac{1}{\sqrt{x+1}-\sqrt{x}}=\frac{1}{\sqrt{x+1}-\sqrt{x}}\times \frac{\sqrt{x+1}+\sqrt{x}}{\sqrt{x+1}+\sqrt{x}}=\sqrt{x+1}+\sqrt{x}\end{array}$$

我們可撰寫單精度 (7 位有效數字) 的程式來驗證上述分析，結果如下表。

$x$	$\sqrt{x+1}-\sqrt{x}$	$\frac{1}{\sqrt{x+1}-\sqrt{x}}$	$\sqrt{x+1}+\sqrt{x}$	$2\sqrt{x}$
1111111	$4.8828125\times {10}^{-4}$	2048	2108.1855	2108.185
2222222	$3.6621093\times {10}^{-4}$	2730.6667	2981.4249	2981.4238
3333333	$2.4414062\times {10}^{-4}$	4096	3651.484	3651.4836
4444444	$2.4414062\times {10}^{-4}$	4096	4126.37	4216.37
5555555	$2.4414062\times {10}^{-4}$	4096	4714.045	4714.045
6666666	$2.4414062\times {10}^{-4}$	4096	5163.9775	5163.9775
7777777	9	FPE	5577.7333	5577.7333
8888888	$2.4414062\times {10}^{-4}$	4096	5962.8476	5962.8476
9999999	0	FPE	6324.555	6324.555
123456789	0	FPE	22222.223	22222.223

表中第一欄展現 10 個從小到大的 $x$ 值，它們至少具有 7 位有效數字。第二欄顯示 $\sqrt{x+1}-\sqrt{x}$ 的計算結果，從表中可見，當 $x$ 為 $7777777$, $9999999$ 或 $123456789$ 時，二者的差為 $0$，導致計算倒數時發生「除以零」的 floating-point exception。

此外，由於精度損失，從 $x=3333333$ (第 3 行) 開始，計算出的差值若非零，則皆為相同數值，顯然不正確。第三欄展現第二欄的倒數，即式 (4) 的計算結果。由於第二欄的數值不準確，第三欄的結果自然也不可靠。

第四欄則採用式 (5) 進行計算，讀者不妨利用計算機驗證，觀察其結果是否更準確。值得注意的是，當 $x$ 很大時，$x+1$ 與 $x$ 極為接近，二者的平方根亦接近，因此第四欄的結果（當 $x$ 很大時）與 $2\sqrt{x}$ 十分接近。這一點可從第五欄的結果觀察到，當 $x\ge 4444444$ (第 4 行) 時，二個欄位的數值幾乎相同。

透過這張表，可直觀地看到精度損失的影響，及利用數式轉換來消除減法運算的重要性。然而，需要強調的是，並非所有運算都能用此方法處理，有時仍需借助其他技巧。

以下是練習題，留給讀者：

1. 計算 $\frac{1-\cos x}{\sin x}$，其中 $x$ 非常接近 $0$。
2. 計算 $\frac{1-\cos x}{x^2}$，其中 $x$ 非常接近 $0$。（提示：考慮半角公式）
3. 計算 $\sin x-\cos x$，其中 $x$ 非常接近 $\frac{\pi }{4}$。（提示：考慮半角公式）

案例探討

前文已藉由若干個範例，說明浮點數運算中常見的錯誤情境。然而，浮點數的問題並非孤立存在，而是與計算過程中的數值密切相關。換言之，在許多情況下，錯誤並非直接源自浮點數的運算，而可能來自於該運算之前的計算結果。例如，以減法為例，導致問題的減法本身可能沒有錯，真正的問題可能出現在減法運算。

此外，若減法運算元的值本身是準確的 (例如來自輸入的數值)，且減法運算無法避免，那麼即便程式包含大量減法運算，仍不會導致嚴重誤差，因為運算仍基於可靠的原始數值。有時，即使經過數式轉換仍無法消弭減法，甚至去除減法後結果仍然不精確，這可能是因為問題本身的性質難以計算精確 (ill-conditioned)。在這類情況下，更換演算法方為上策。此外，某些減法是無害的，例如，若 $a$ 與 $b$ 幾乎相等，但 $a-b$ 的值遠小於 $c$，則 $(a-b)+c$ 中的減法影響甚微。（請思考其中的數學原理）

至於四捨五入誤差，大多數情況下，捨去與進位的誤差會相互抵消。此外，由捨入誤差引發的問題往往僅發生於特定運算，而非來自大量反覆計算，因此只要能追溯誤差的源頭，通常就能找到解決方案。

以下公式是個簡單的例子，數學家 Leonhard Euler 證明此極限的值為 $e=2.71828\dots$，我們藉由電腦程式，讓 $n$ 逐漸增大，觀察是否能計算出接近 $e$ 的值。
$$\underset{n\to \mathrm{\infty }}{lim}{(1+\frac{1}{n})}^n$$

程式碼相當直觀，如下所示。

#include <stdio.h>
#include <math.h>

int main(void)
{
    int nn = 9;
    int i;
    int n = 10;
    float EE = expf(1.0);
    float E;
    for(i = 0; i < nn; i++){
        E = powf(1.0 + 1.0/n, n);
        printf("%10d %.10f %.10f\n", n, E, fabsf(E - EE));
        n *= 10;
    }

    return 0;
}

我們將 $n=10,100,1000,10000,\dots$ 代入公式，檢視結果是否趨近於 $2.71828\dots$。請注意，程式中的四捨五入影響僅限於 E 的計算，誤差不會進一步擴散。

以下是計算結果：

        10 2.5937430859 0.1245386600
       100 2.7048113346 0.0134704113
      1000 2.7170507908 0.0012309551
     10000 2.7185969353 0.0003151894
    100000 2.7219622135 0.0036804676
   1000000 2.5952267647 0.1230549812
  10000000 3.2939677238 0.5756859779
 100000000 1.0000000000 1.7182817459
1000000000 1.0000000000 1.7182817459

當 $n$ 達到 10000 時，計算出的近似值為 $2.718597$，與真實值 $e$ 相差 $3.1518936\times {10}^{-4}$。然而，自此之後，計算結果與 $e$ 的差距不減反增，最終結果竟然收斂至 1，這顯然是錯誤的！問題究竟出在哪裡？

根本原因在於捨入誤差，尤其是在計算 E 時產生的誤差。前文提及，$0.1$ 在十進位轉換為二進位或十六進位時，會形成循環小數，而 $\frac{1}{{10}^n}$ 亦同理，導致儲存至記憶體時必然產生捨入誤差。然而，真正的問題並非來自這些捨入誤差本身，而是來自加法運算！

假設計算機具有 7 位有效數字，則 $1/{10}^7=0.0000001$，加上 1 後得到 $1.0000001$，以標準型式表示為 $0.10000001\times {10}^1$。然而，當捨入至 7 位數時，結果變為 $0.1\times {10}^1$，即回到 1，這正是最終計算結果錯誤的根本原因。

這個例子清楚說明，捨入誤差不一定需要累積才能導致問題，關鍵在於識別誤差產生的關鍵環節。在此程式中，解決方案在於重新改寫公式，以降低捨入誤差對計算結果的影響。這是個頗具挑戰性的數學問題，以下藉由 Spigot 演算法，避免傳統方法累加後再取小數而造成的精度損失，而 Spigot 演算法則每次直接產生下位數字，程式碼如下：

#include <stdio.h>
#define N 300
int main()
{
    int arr[N + 9] = {0}, idx = 0;
    char digits[N + 1];

    for (int n = 1; n <= N; n++) arr[n] = 1;

    for (int m = N; m > 0; m--) {
        int x = 0;
        for (int n = m; n > 1; n--) {
            x += arr[n] * 10;

            /* Calculate arr[n] = x % n */
            int mod = x;
            while (mod >= n) mod -= n;
            arr[n] = mod;

            /* Calculate x /= n */
            int div = 0;
            for (int sum = x; sum >= n; sum -= n)
                div++;
            x = div;
        }
        digits[idx++] = x + '0';
    }

    digits[idx] = '\0';
    printf("e = 2.%s\n", digits);
    return 0;
}

這段程式碼的功能是輸出 $e$ 的前 300 位有效數字，讀者可嘗試將巨集 N 替換為其他數值進行測試。該程式不僅完全避免浮點數的精度問題，還省略除法運算，藉由整數運算與進位處理來逐位展開數學常數，有效確保數值準確性。

另一種看似直觀的方法是，改用倍精度來提高計算的精確度，例如，當 $n=100000000$ 時，倍精度計算得到 $2.7182817983473577$，與真實值 $e$ 的誤差僅為 $3.011168780986395\times {10}^{-8}$。然而，倍精度理應提供約 15 位有效數字，但該方法僅獲得 8 位，這是否有些奇怪？還是說，即便改用倍精度計算，仍無法真正達到 15 位的精確度？

遺憾的是，提高精度並非萬靈丹，無法解決所有由捨入誤差引發的問題。讀者不妨嘗試這實驗：拿一台計算機，輸入一個正數，反覆按開平方鍵多次，然後再按相同次數的平方鍵，觀察最終結果。一般而言，結果與原始輸入值相差不大，但若開平方的次數足夠多，結果將逐漸趨近於 1，進而在平方回復時與原始值產生顯著誤差。一旦計算過程中開平方的結果變為 1，無論再平方多少次，最終都無法回到原始輸入值。這種現象無論在何種計算機、何種精度下皆會發生，差異僅在於何時趨近於 1 而已。

我們可透過撰寫程式來驗證本現象。以下為程式的部分輸出結果，此程式對 $2$ 反覆進行平方根運算若干次，然後再進行相同次數的平方運算。首先，當開平方 10 次後再平方 10 次，第 10 次的開平方結果為 $1.0006771088$：

    0   2.0000000000
       ………     
    5   1.0218970776
       ………
   10   1.0006771088

再平方 10 次，結果是 $1.9999581575$，誤差是 0.00004 左右，還好。

    0   1.0006771088
       ……
    5   1.0218964815
       ……
   10   1.9999581575

若把開方次數增加到 25 次，下面是結果，到第 23 次結果就是 $1.0$，當然平方回來就不可能得到 $2$。

    0   2.0000000000
       ……
   10   1.0006771088
       ……
   15   1.0000211000
       ……
   20   1.0000005960
   21   1.0000002384
   22   1.0000001192
   23   1.0000000000
   24   1.0000000000
   25   1.0000000000

於是，改用倍精度是否能解決問題？由於單精度 (約 7 位精度) 在開平方約 23 至 25 次後便會收斂至 $1$，而倍精度 (約 15 位精度) 則大約在 50 次 (或稍多) 開平方後也會趨近於 1。因此，在這個問題中，增加精度只是延後問題的發生，不能真正解決問題。

不過，若在精確度耗盡前，已獲得可用的結果，提高精度仍然是可行的選擇。

真實世界中因浮點數誤差導致的問題

過往因浮點數錯誤、導致嚴重後果的事件，在人類短暫的數位化過程中屢見不鮮，以下選出三個值得參考的案例。

1991 年波斯灣戰爭期間，美國愛國者 (Patriot) 反飛彈系統因浮點運算誤差未能攔截來襲的飛毛腿 (Scud) 飛彈，導致 28 名美軍陣亡。該系統的內部時鐘以整數計時，每個單位代表 0.1 秒，而來襲飛彈的速度則使用浮點數表示。當整數時間乘以 0.1 轉換為 24 位元浮點數時，因 $0.1$ 在二進位中為無限循環小數，轉換時產生約 $9.5\times {10}^{-8}$ 的誤差。此誤差與飛彈速度及運轉時間成正比，當系統運行 100 小時後，來襲飛彈的位置誤差累積至 573 公尺，導致攔截失敗。

延伸閱讀: Patriot Missile Software Problem

1996 年，歐洲太空總署 (European Space Agency) 在法屬圭亞那 Kourou 發射 Ariane 5 火箭，然而發射後 40 秒火箭失去控制，被迫引爆摧毀，火箭爆炸解體令 2 名法國士兵當場死亡。事後調查發現，程式將一個 64 位元浮點數轉換為 16 位元有號整數 (最大值 32767)，導致整數溢位，使火箭控制系統錯亂，進而失控。

延伸閱讀: Ariane flight V88

1997 年 9 月，美國海軍 USS Yorktown 號驅逐艦因水兵誤輸入 0，導致除零錯誤。錯誤結果迅速擴散，使艦上的推進系統關閉，導致軍艦在海上靜止長達 2 小時 45 分鐘。雖然事件未發生於戰時，但這類錯誤顯示基礎程式設計缺陷可能帶來嚴重後果。未檢查除數是否為 0，甚至未驗證輸入數值的正確性，這些都是可避免的錯誤。

延伸閱讀: When Smart Ships Divide By Zer0 — Stranding the USS Yorktown

更多案例可見〈軟體缺失導致的危害〉，遺憾之餘，令我們不禁反思：撰寫程式時，是否已盡力避免這些細微但影響深遠的錯誤？是否曾因類似問題導致企業遭受嚴重損失，甚至影響政府決策？

浮點數計算無可避免地影響日常程式開發。無論是否為電腦科學專業，都應理解浮點數的基本原則，以避免重蹈歷史錯誤，確保程式計算結果的可靠。

結語

本文藉由多個範例，闡述使用浮點數時常見的問題，但我們未進行嚴格的定義或數學論證，而僅依賴最基本的程式設計知識，這是便於程式開發者的理解。

以下是幾點基本建議：

1. 儘量避免對二個幾乎相等的數值進行減法運算，除非這兩個運算元的值是完全精確的，或減法的結果對後續運算影響不大。
2. 注意溢位或接近溢位的狀況。一般而言，若計算過程中的中間結果遠大於或遠小於最終結果 (或輸入數值)，則可能存在溢位風險。在這種情況下，可嘗試改寫公式，或對數值進行適當的縮放，以降低溢位風險。
3. 留意捨入誤差，特別是在處理數值差距極大的加法或減法運算時，較小的數值可能會完全失去影響力，致使錯誤的計算結果。
4. 不要在關鍵應用場域中，直接套用教材中的演算法或程式。教學用的範例主要用於解釋概念，除非作者明確指出該演算法已考量數值計算的穩定性議題，否則當輸入數值較為特殊時，演算法的輸出可能是不可靠的。

這些原則有助於降低浮點數誤差對計算結果的影響，確保電腦程式的強健與可靠。