不熟悉浮點數的人最容易犯的錯誤之一,就是直接用 == 或 != 比較兩個浮點數。以最常見的 IEEE 754 浮點數來說,下面這樣的判斷式竟然不成立:

if (0.1 + 0.2 == 0.3)

原因在於以二進位表示的浮點數並沒有辦法精確儲存 0.1、0.2、0.3 這些十進位實數,只能以最接近的浮點數表示,和原本的數值有微小的誤差。三個各自帶有誤差的數字要碰巧讓整個等式成立,實在是相當困難的一件事。基於同樣的理由,在採用 IEEE 754 的環境下,以下程式片段陷入無窮迴圈也就沒什麼好奇怪的了:

double x = 0;
while (x != 5.0) {
    x += 0.1;
}

不同的計算環境可能採用不同的浮點數系統,但進制轉換、儲存、運算的誤差必然存在。使用 <、> 不會有問題,但遇到 == 或 != 就要變得疑神疑鬼,浮點數不會輕易讓等號成立。

對於一般應用我們並不指望兩數完完全全相等,只要兩數差在可接受範圍內就好了。最簡單的方法是計算絕對誤差:

if (fabs(a – b) <= 0.001)

所謂的「可接受範圍」通常與應用領域有關,有個偷懶選擇是直接利用DBL_EPSILON、FLT_EPSILON 做為誤差範圍的門檻,但這並不總是合理的做法,對於離 0 太遠的數值範圍甚至是很詭異的選擇,以下就用實驗來說明。

我手邊的編譯器 FLT_EPSILON 大約是 1.19e-7,精確數字不重要,讀者只要對這個數量級(0.0000001)有概念就行了。假如我要比較 50000.0f 和 50000.0039f 兩數,並且以 FLT_EPSILON 做為絕對誤差門檻,計算絕對誤差的程式碼如下:

float a = 50000.0f;
float b = 50000.0039f;

if (fabs(a - b) <= FLT_EPSILON) {
    puts("close enough");
}

問題來了,0.0039 比 FLT_EPSILON 大上 3 萬多倍,顯然沒有資格落入足夠接近的範圍。但在 float 表示法中,50000.0039f已經是所有大於 50000.0f的數字中,和 50000.0f 最接近的數字了。

要知道離 0 越遠,接鄰兩浮點數之間的差就越大。所有可正常表示的 float 當中,只有在 [-1, 1] 區間內,接鄰兩數差才會小於 FLT_EPSILON;若在 [1, 2] 區間,接鄰兩數差恰好為 FLT_EPSILON,負數亦然。隨著數字增加,FLT_EPSILON 很快就會小到不像話,到了 16.0f 時,接鄰兩數差將達到 FLT_EPSILON 的 16 倍;到了16777216.0f 時,接鄰兩浮點數差已經大於 1.0。

這就是為什麼在 [-1, 1] 區間外使用 DBL_EPSILON、FLT_EPSILON 並不合理,就如同死守公分做為誤差單位,對於精密機械而言大得誇張,但對於計算交通路程來說卻又小得可笑,還不如實際考慮需求而設定誤差門檻。

在數值範圍變化大的場合,更好的作法是隨著尺度調整誤差門檻,也就是運用相對誤差。其中一個非常簡單的想法概念碼如下:

int AlmostEqual(float a, float b, float tolerance)
{
    float absA = fabsf(a);
    float absB = fabsf(b);
    return fabsf(a - b) / fmaxf(absA, absB) <= tolerance;
}

為了避免除以零的狀況,移項之後得到:

return fabsf(a - b) <= tolerance * fmaxf(absA, absB);

上面這種比較法只需要誤差相對於 a、b 其中一個足夠小即成立,更嚴格的比較則要求誤差同時對兩者都夠小。由於乘法有機會造成 underflow,因此也有人比較喜歡除法,畢竟運算是否造成 underflow 並不像除以零那麼容易事先檢查。至於 tolerance 該如何選擇呢?這個問題沒有固定答案,通常會選擇 FLT_EPSILON 的整數倍。

浮點數的不連續性衍生出另一種實現相對誤差想法,就是以兩數中較大者為比較基準,測試另一個數是否在基準的前後 N 個浮點數之內。C99 內建的 nextafter() 函數可以用來列舉接鄰的浮點數,然而效能未必令人滿意。進一步的加速巧門必須直接操作底層表示法,而且得手工處理可能出現的浮點數特殊值,細節在此不多說明。

以上方法各自有適用情境,這裡還必須提醒一件事,帶有誤差的比較都會失去遞移性,這是和等號不同的地方。不論採用哪一種方法,最關鍵的考量還是計算範圍與精度兩方面的需求。以上面的例子來說,若計算範圍可達 50000,而且精度需達小數點以下三位,則 float 的有效位數根本不可能滿足需求,任何比較方法都無濟於事。

接下來介紹一個不良示範,似乎有人認為可以使用取巧的方法來比較兩個浮點數陣列:

int ArrayEqual(double* m1, double* m2, size_t count)
{
    return !memcmp(m1, m2, sizeof(double) * count);
}

除非函數目的確實是要比較記憶體內容,否則這可以說是非常糟糕的主意,甚至比直接用 == 還糟糕。因為許多浮點數系統允許特殊數值存在多種不同的表示法,例如 0 和 NaN。就拿 0 來當例子吧,IEEE 754 有 +0 和 -0 兩種表示法,這兩種 0 會被 == 視為相等,計算特性也幾乎相同,但內部的表示法就是不一樣。有興趣的話可以做個實驗:

double a = 0.0;
double b = -1.0 * 0.0;

printf("%d \n", a == b);
printf("%d \n", memcmp(&a, &b, sizeof a));

在大部分的計算環境底下輸出應該會是:

1
-1

經由不同計算途徑得到的 0 很可能使用不同的內部表示法,儘管從數值的觀點來看都代表0。


關於浮點數比較的基本觀念: http://www.boost.org/doc/libs/1_34_1/libs/test/doc/components/test_tools/floating_point_comparison.html

這篇比較偏向實作: http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/

不良示範取材自:http://edisonx.pixnet.net/blog/post/88826345

novus 發表在 痞客邦 PIXNET 留言(1) 人氣()


留言列表 (1)

發表留言
  • edisonx
  • [-1,1] 裡才取用 DBL_EPS / FLT_EPS,這點還真漏考慮。你沒翻出來,我都忘了以前發生過那事。