Back to Basic: 談浮點數的比較 - novus log - 痞客邦

不熟悉浮點數的人最容易犯的錯誤之一，就是直接用== 或!= 比較兩個浮點數。

以最常見的IEEE 754 浮點數來說，下面這樣的判斷式竟然不成立： if (0.1 + ... novuslog 跳到主文 . 部落格全站分類：生活綜合 相簿 部落格 留言 名片 Jun08Sat201301:40 BacktoBasic:談浮點數的比較 不熟悉浮點數的人最容易犯的錯誤之一，就是直接用==或!=比較兩個浮點數。

以最常見的IEEE754浮點數來說，下面這樣的判斷式竟然不成立： if(0.1+0.2==0.3) 原因在於以二進位表示的浮點數並沒有辦法精確儲存0.1、0.2、0.3這些十進位實數，只能以最接近的浮點數表示，和原本的數值有微小的誤差。

三個各自帶有誤差的數字要碰巧讓整個等式成立，實在是相當困難的一件事。

基於同樣的理由，在採用IEEE754的環境下，以下程式片段陷入無窮迴圈也就沒什麼好奇怪的了： doublex=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做為絕對誤差門檻，計算絕對誤差的程式碼如下： floata=50000.0f; floatb=50000.0039f; if(fabs(a-b)<=FLT_EPSILON){ puts("closeenough"); } 問題來了，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並不合理，就如同死守公分做為誤差單位，對於精密機械而言大得誇張，但對於計算交通路程來說卻又小得可笑，還不如實際考慮需求而設定誤差門檻。

在數值範圍變化大的場合，更好的作法是隨著尺度調整誤差門檻，也就是運用相對誤差。

其中一個非常簡單的想法概念碼如下： intAlmostEqual(floata,floatb,floattolerance) { floatabsA=fabsf(a); floatabsB=fabsf(b); returnfabsf(a-b)/fmaxf(absA,absB)<=tolerance; } 為了避免除以零的狀況，移項之後得到： returnfabsf(a-b)<=tolerance*fmaxf(absA,absB); 上面這種比較法只需要誤差相對於a、b其中一個足夠小即成立，更嚴格的比較則要求誤差同時對兩者都夠小。

由於乘法有機會造成underflow，因此也有人比較喜歡除法，畢竟運算是否造成underflow並不像除以零那麼容易事先檢查。

至於tolerance該如何選擇呢？這個問題沒有固定答案，通常會選擇FLT_EPSILON的整數倍。

浮點數的不連續性衍生出另一種實現相對誤差想法，就是以兩數中較大者為比較基準，測試另一個數是否在基準的前後N個浮點數之內。

C99內建的nextafter()函數可以用來列舉接鄰的浮點數，然而效能未必令人滿意。

進一步的加速巧門必須直接操作底層表示法，而且得手工處理可能出現的浮點數特殊值，細節在此不多說明。

以上方法各自有適用情境，這裡還必須提醒一件事，帶有誤差的比較都會失去遞移性，這是和等號不同的地方。

不論採用哪一種方法，最關鍵的考量還是計算範圍與精度兩方面的需求。

以上面的例子來說，若計算範圍可達50000，而且精度需達小數點以下三位，則float的有效位數根本不可能滿足需求，任何比較方法都無濟於事。

接下來介紹一個不良示範，似乎有人認為可以使用取巧的方法來比較兩個浮點數陣列： intArrayEqual(double*m1,double*m2,size_tcount) { return!memcmp(m1,m2,sizeof(double)*count); } 除非函數目的確實是要比較記憶體內容，否則這可以說是非常糟糕的主意，甚至比直接用==還糟糕。

因為許多浮點數系統允許特殊數值存在多種不同的表示法，例如0和NaN。

就拿0來當例子吧，IEEE754有+0和-0兩種表示法，這兩種0會被==視為相等，計算特性也幾乎相同，但內部的表示法就是不一樣。

有興趣的話可以做個實驗： doublea=0.0; doubleb=-1.0*0.0; printf("%d\n",a==b); printf("%d\n",memcmp(&a,&b,sizeofa)); 在大部分的計算環境底下輸出應該會是： 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 novuslog novus發表在痞客邦留言(1)人氣() E-mail轉寄 全站分類：不設分類個人分類：C&C++此分類上一篇：volatile是不夠的 此分類下一篇：在編譯期選擇演算法 上一篇：volatile是不夠的 下一篇：一張流程圖教你解讀陰謀論 ▲top 留言列表 發表留言 Category Computer(6) 軟體(8)Python(4)C&C++(51)Linux(4)心得(11)圖形(10) 剪貼簿(2)生態環保(4)Funstuff(39)超自然揭秘(12)Science&Math(29)雜文(69)繪畫與工藝(7)版務(6)未分類文章(1) Recent Comment Archive Archive 2018十一月(1) 2018十月(1) 2017八月(3) 2017七月(3) 2017六月(1) 2017五月(1) 2017一月(1) 2016八月(1) 2016七月(1) 2016三月(1) 2016二月(1) 2015十一月(1) 2015十月(1) 2015九月(3) 2015八月(1) 2015七月(1) 2015五月(2) 2015四月(1) 2015三月(1) 2015二月(3) 2014十月(1) 2014九月(2) 2014八月(5) 2014七月(2) 2014六月(3) 2014四月(3) 2014三月(1) 2014二月(1) 2014一月(2) 2013十二月(1) 2013十一月(3) 2013八月(3) 2013七月(3) 2013六月(3) 2013五月(3) 2013四月(1) 2013三月(2) 2013二月(3) 2013一月(4) 2012十二月(3) 2012十一月(1) 2012十月(1) 2012八月(1) 2012七月(1) 2012六月(3) 2012五月(4) 2012四月(1) 2012三月(2) 2012二月(4) 2012一月(1) 2011十二月(2) 2011十一月(2) 2011十月(2) 2011九月(1) 2011八月(5) 2011七月(7) 2011六月(3) 2011五月(4) 2011四月(2) 2011二月(4) 2011一月(4) 2010十一月(6) 2010十月(5) 2010九月(3) 2010八月(6) 2010七月(8) 2010六月(8) 2010五月(11) 2010四月(4) 2010三月(2) 2010二月(2) 2010一月(1) 2009十二月(7) 2009十一月(5) 2009十月(4) 2009九月(2) 2009八月(1) 2009七月(2) 2009六月(3) 2009五月(3) 2009四月(1) 2009三月(2) 2009二月(2) 2009一月(3) 2008十二月(5) 2008十一月(3) 2008十月(1) 2008九月(5) 2008八月(10) 2008七月(5) 2008六月(1) 2008五月(1) 2008三月(1) 2008一月(1) 所有文章列表 Search RSS Links Visited 本日人氣： 累積人氣： 回到頁首 回到主文 免費註冊 客服中心 痞客邦首頁 ©2003-2022PIXNET 關閉視窗