使用浮點數最最基本的觀念

pxcwan

转自：冼鏡光 DCview.com達人部落格
http://blog.dcview.com.tw/article.php?a=BTtXMgFpUWA%3D



使用浮點數最最基本的觀念

冼鏡光
2008年5月12日上線

    雖然此地是攝影網站，從本篇起我會偶而貼一些與電子計算機有關的文章，原因是自己懶而且沒時間在電腦網站上另起爐灶，所以攝影同好們請原諒，而玩電腦的朋友（特別是學生與專業程式設計者）們得閒時不妨來逛逛、並且廣為宣傳。我當然知道電腦文章貼在此地點閱率不會高，但目前就只有在這兒張貼的打算。

      本文的引用方式如下：
file:///C:/DOCUME~1/CHING-~1/LOCALS~1/Temp/msohtml1/01/clip_image001.gif
          非學術性請用
          冼鏡光 DCview.com達人部落格 本文的網址
file:///C:/DOCUME~1/CHING-~1/LOCALS~1/Temp/msohtml1/01/clip_image001.gif
          學術性請用
          冼鏡光，使用浮點數最最基本的觀念，DCview達人部落格，本文的網址，讀取時間（本文標題中最新更正時間）




        不知道各位在程式中使用浮點（或實）數的機會有多大，也不知道各位的程式有否因為浮點數問題而不能產生正確結果的困擾，也不知道是否有程式因為沒有用好浮點數而得到錯誤結論而不自知。在科學、工程、統計上幾乎無法避免使用浮點數，所以如何趨吉避凶是每一個常用浮點數的程式員的基本知識。不過妥善運用浮點數的學問很大，學校裡頭教得不多，充其量是在數值方法這門課中學些皮毛，涵蓋面的廣度與深度相當有限，所以打算透過這篇文字提供一些使用浮點數的最基本觀念，而不做深入的理論性探討。雖然本文的程式都用Fortran 90寫作，改用其它語言的效果並無不同；另外，文中很多名詞採用直覺的涵意，而沒有很精確的定義，主要是避免太過艱澀的理論題材。正因為這個出發點，本文所需的預備知識大致上就只有最簡單的程式寫作基礎，比如單精度、倍精度、浮點數的表示方式等等，所以大多數讀者應該都可以看完這篇文字。這篇文字很長，希望您可以耐心看完，相信對您的工作與教學會有些幫助。

        在進入正文之前要聲明一點，以下是我教計算機工程三年級數值方法講義的一部份，未來會出現在本人的著作中，所以未經本人同意請勿抄錄、轉載、或轉貼，敬請尊重智慧財產權。


1. 浮點數的最根本問題

        使用浮點數有兩個最根本的問題：（1）輸入與儲存的值不一定精確、與（2）計算的結果會有誤差。第一點的成因不外乎以下幾項：輸入值本身就不精確、計算機儲存浮點數用的記憶體有限、不同底數之間的轉換本來就不精確、溢位、以及數學上習慣使用的運算定律在浮點數運算時不一定成立。第二，因為使用浮點數進行運算時會有上述第一點所提到的以及其它的誤差，這些誤差經過反覆而且大量運算之後很有可能會蔓延到其它所在，從而產生一發不可收拾的後果。下面各節先就第一點的成因做個介紹，稍後再討論第二點。


1.1  輸入值本身不精確 ***

        因為實數可能有無限位，於是就無法輸入所有位數，當然計算機得到的浮點數就不會精確；譬如說，有理數1/3、1/7、1/13等有無限位數，無理數如2的平方根、3的立方根、圓周率p = 3.1415926……等也有無限位數。所以，我們無法完整地把這些值輸入到計算機中，從而產生誤差。圖一中是一個頂點在(n,0)、(-n,0)與 (0,sqrt(3)*n)的等邊三角形，數學中表示這個三角形的座標不會有任何困難，但因為3的平方根有無限位數，輸入sqrt(3)*n時就不可能準確，所以輸入後儲存在計算機中的三角形會是個非常靠近等邊三角形的等腰三角形，而不是個真正的等邊三角形，於是程式中基於等邊三角形所進行的運算就很可能不精確了。




1.2 儲存浮點數的記憶體有限

        計算機中儲存浮點數都是用小數加乘冪的標準型，以十進位而言就像是這樣： ± 0.xxx…xxx 10±yyy；換言之，實數-123.456會轉換成-0.123456 x 103，而0.00123456會轉換成0.123456 10-2。在計算機中，前面的小數部份±0.xxx…xxx與後面的乘冪±yyy會分開來儲存到記憶體中，就像是圖二那樣。



        在圖中，小數部份與它的符號通常是分開來的，而乘冪則另外儲存，至於乘冪部份的符號要如何儲存則取決於硬體的設計；為了避免太多技術性的討論，我們姑且不去碰IEEE浮點數算術這個標準。因為儲存小數部份與乘冪部份都只有有限個數元（bit），這就產生了兩個問題。

        第一，小數部份不可能完整地儲存一個有無限位的實數（比如，根號2、1/3、圓周率等），所以會有誤差。第二，用來儲存乘冪的數元（含符號部份）比用來儲存小數部份的要少得多，如果一個浮點數的乘冪值比能夠儲存的值來得大，硬體就無法儲存這個乘冪，於是發生溢位（overflow，有時也叫做上溢位）；反之，如果一個浮點數乘冪值比能夠儲存的值來得小，硬體也無法儲存這個乘冪，因而發生下溢位（underflow）。看個例子，如果硬體在儲存乘冪時的最大與最小值分別是+63與-63，於是0.1234 ´ 1068會造成（上）溢位，0.1234 ´ 10-75會造成下溢位。

        在計算過程中如果出現上溢位與下溢位，計算機會停止作業，並且顯示floating-point exception或類似訊息，有些機型可能會明白指出上溢位、下溢位、或用0除。


1.3 數系轉換並不精確

        我們日常用十進位數，但計算機的浮點數卻常用二進位或十六進位運算與儲存，所以寫在程式中的十進位數得轉換成二進位或十六進位，但是一個簡單的十進位數卻可能轉換成無限位的二進位或十六進位數。

        圖三是把十進位的0.1轉換成二進位（左行）與十六進位（右行）的計算方式。簡單地說，我們用2或16去乘給定的十進位數，取出小數點前方的值（這是結果的一部份），再用2或16去乘剩下來的小數部份等等。在左行我們用2去乘0.1，依次得到0、0、0、1、1、0，接下來的小數部份是0.4，又回到前面出現過的內容，因此0.1（十進位）就轉換成一個二進位的循環小數0.0 0011 0011 0011 0011 ………；同理，十進位的0.1也轉換成一個十六進位的循環小數0.1999999999……。因此，簡單的十進位0.1就無法完整地在計算機中儲存；換句話說，程式中寫了0.1，但儲存到記憶體中的卻不是0.1、而是0.1的近似值，於是算出來的結果就可能不夠精確了。



        如果浮點數有整數與小數部份，這兩部份得分開轉換，再寫成標準型式；譬如十進位的277.31得分成277與0.31兩部份。若要把277.31轉換成十六進位，於是277轉換後得到115（十六進位），而0.31則轉換成0.4 F5C28 F5C28 F5C28 ……，一個無限位的十六進位循環小數；把兩者合併就是

115.4 F5C28 F5C28 F5C28 ……（十六進位）

寫成標準型得到

0.1154 F5C28 F6C28 F5C28 ……… ´ 163

請注意的是，上面的乘冪是用16做底數，存入計算機內就像是圖四這樣：


        不過，計算機不可能儲存有無限位的循環小數，於是上面結果中的後半部就會被切掉。一般而言，在單精確度（single precision）的情況下，32個數元的浮點數儲存方式可以容納六個十六進位的位數，若不計較四捨五入的規律而只用捨去的方式，計算機就只能儲存1154F5這一段，見圖五中的黃框部份：




換句話說，計算機只儲存了0.1154F5 ´ 163這個數。因為0.1154F5 ´ 163等於115.4F5，把115與0.4F5轉換回十進位分別得到277與0.309814453125，所以結果是277. 309814453126，而不是原來的277.31。這表示在程式中使用277.31時，計算機內部卻用的是277. 309814453126，當然計算的結果就會有誤差。


1.4 溢位

        前面提過，在浮點數的標準型中，如果小數部份的長度超過硬體的界限，多出來的位數就會被四捨五入後存入記憶體，這就產生了誤差。但若硬體無法儲存乘冪部份，就會有溢位；如果乘冪的值大於硬體所能容許的值，結果是上溢位（overflow），如果乘冪的值（負數）小於硬體所能容許的值，結果為下溢位（underflow）。計算過程中如果發生了溢位，計算機會顯示floating-point exception的訊息、然後停止作業。

        程式中的四捨五入誤差與溢位是很令人困擾的，因為前者不容易發現、一旦發生了（結果不正確）又很難追蹤；如果後者（溢位）發生了，程式會停止執行，要補救就有點遲了。所以在寫程式時，我們應該留意可能會產生很大或很小的值的所在，稍許改變程式的寫法來避免這個問題。

        讓我們看個極為簡單的例子，計算a平方加b平方後再取平方根（圖六）



大概沒有人不會寫這個程式，不是嗎？但您是否想到過a或b中有一者很大時所帶來的問題？若a的值很大，a的平方更大而可能會溢位；縱使a與b不一定很大，但a平方與b平方的和卻也有可能造成溢位，雖然這個和的平方根不一定會溢位。若要降低溢位的機會，我們可以把a與b的值用圖七的手法調低：



簡單地說，我們把a與b用k > 0去除、再把平方根乘上k就可以得到相同的結果。那一個k值好用呢？比較理想的是k = MAX(|a|,|b|)，也就是|a|與|b|中大的那一個，於是在a/k平方與b/k平方中有一者為1、另一者小於1，當然兩者都不會發生上溢位（但小的那個卻有可能產生下溢位）；它們的和會小於2，開平方後的值在1與2之間，把平方根乘以k時造成溢位的可能性也低些。

下面圖八是幾何平均數的定義，所有的x值都大於0。很明顯地，若直接引用這個定義計算，當輸入值很大時出現溢位的可能性是相當大的，您可以用上述的手法避免溢位，請試試看。



        硬體只會對浮點數的溢位產生floating-point exception的訊息，整數的溢位通常是不理會的。舉例而言，算階乘的結果增加非常快，在一部32數元的機型上12!大致就是能夠算出沒有誤差的階乘值的極限，13!以上就會發生整數溢位使結果不正確，下面就是個好例子。程式算到12!時得到479001600，最後兩位數是0，但13!時最後兩位卻不是0，這當然是錯的，因為13! = 12! ´ 13，結果的最後兩位應該是0。另外，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


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

        最常見而且最常用的數學運算定律有三個：交換律（commutative law，a+b = b+a與a´b = b´a），結合律（associative law，(a+b)+c = a+(b+c)與(a´b) ´c = a´(b´c)），與分配律（distributive law，(a+b) ´c = a´c + b´c）。在浮點數運算時，這三者中只有交換律成立，結合律與分配律都會出問題。

        圖九是兩個例子，我們使用只有兩位有效數字的運算；圖中左行是說明結合律不成立的情形。首先，0.12´0.34手算是0.0408，因為只有兩位有效數字，0.0408 = 0.408 ´ 10-1就被四捨五入成0.41 ´ 10-1 = 0.041，因而有了四捨五入誤差；最後0.041 ´ 4是0.164，再次四捨五入得到0.16。另一方面，0.34´4 = 1.36，四捨五入得1.4；接著，0.12´1.4 = 0.168，再次四捨五入得到0.17。因為0.16不等於0.17，所以結合律對浮點數運算不成立！\n


圖九右行是分配律。首先，0.12+0.99是1.11，四捨五入後是1.1，所以1.1´5 = 5.5。另一方面，0.12´5是0.6，0.99´5是4.95、四捨五入後得5.0，於是0.6+5.0為5.6。但是5.5不等於5.6，因此分配律也不成立！

        事實上，結合律還可能造成溢位。考慮圖十的式子，左邊括號內的結果是1，所以等號左邊的結果是1´1063；但右邊括號中的結果是1´10126，在一些機型之下乘冪為126時多半會發生上溢位。這個例子再次告訴我們在浮點數運算時結合律不成立，因此寫作運算式時要特別小心。



        讓我們看個實際的計算例子。圖十一中的五道運算式在數學上是完全相同的，它們不過是把第一式用結合律與分配律改寫過而已：



        我們打算寫一個程式，令R = 3.0並且從x0 = 0.5開始計算，把五道式子的結果算出來；也就是說，對每一道式子而言，我們從x0 = 0.5起算出x1，把x1代回原式算出x2，等等。為了節省空間，我們只顯示x100、x200、x300、……、與x1000的結果，這就是圖十二；圖中第一行是i的值，第二到第六行分別是上面第一到第五道式子算出的結果。仔細比較同一列上的數值，相信不難發現五道式子得來的結果有多大差異；所以，這個結果告訴我們結合律與分配律不成立時的影響有多大。




2 再談四捨五入誤差的影響

        前面提過，計算得來的結果會被四捨五入，當然就不精確，如果這個不精確的結果又用來計算其它的值，於是誤差可能會滾雪球般地愈來愈大，最後使程式算出完全不合情理的結果。

        要注意的是，單一四捨五入所造成的誤差並不大，在最好的情況下大約是最後一位的加／減50%。譬如說，在一部有三位有效數字的機型上，如果計算結果是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.1234之間，或者說成在[0.123-0.0005，0.123+0.0005]=[0.1225,0.1235]中、這正好是最後一位的±50%。

        這個細微的捨入誤差很可能會導致不正確的結果。譬如說，我們希望把區間[1,2]分成三等分，並且在x為1、1 1/3、1 2/3、與2時用x的值做進一步的計算。下面是一個Fortran 90的程式，其中a、b、x、與h都是至少有15位有效數字的實數變數（見DOUBLE的定義）。程式看起來很正確、四平八穩誰都寫得出來，找不到什麼漏洞，不是嗎！？

PROGRAM Rounding
IMPLICIT NONE
INTEGER, PARAMETER :: DOUBLE = SELECTED_REAL_KIND(15)
REAL(KIND=DOUBLE) :: a = 1.0_DOUBLE, b = 2.0_DOUBLE
REAL(KIND=DOUBLE) :: x, h
INTEGER              :: n = 3

h = (b - a)/n
x = a
WRITE(*,*) "x = ", x
DO
    IF (x >= b) EXIT
    x = x + h
    WRITE(*,*) "x = ", x
END DO
END PROGRAM Rounding

        其實不然，下面是執行的結果，它多了一個輸出值2.33333…！從輸出中很明顯地可以看出，x從1起經過1.333….（小於1 1/3）、1.666….（小於1 2/3）、1.999…..（小於2），以致於多出一次循環才會得到x >= 2的結果從而離開DO，這就是捨入誤差惹來的禍。這個問題在小程式中很容易察覺，但在大而且複雜的程式中要找出多一次循環就會費工夫了。

x = 1.0
x = 1.3333333333333332
x = 1.6666666666666665
x = 1.9999999999999997
x = 2.333333333333333

        在多數情況下，把IF改成下面的寫法可以避免多出一次循環的問題：

IF (x >= b – h/2) EXIT

        附帶一提，在Fortran 77中DO可以用實數，所以上面的廻圈可以寫成

DO x = a, b, (b-a)/n
   WRITE(*,*) “x = “, x
END DO

        Fortran 77編譯程式會從a、b與(b-a)/n的值算出循環的次數，再用整數做計數，但是在Fortran 90/95中己經建議不要使用這項功能，而且很可能在往後的新標準中消失。因此，用C/C++或Java的朋友能否由此學到些教訓？

        讓我們再看一個較長而且極端的例子結束本段。如果有n > 1個值x1、x2、…、xn，它們的平均數m是這n個值的和被n去除的結果（亦即m = (x1+x2+…+xn)/n），它們的變異數（variance）定義如圖十三：



因此算變異數時要先求出平均數m，再求各個xi - m的平方和，最後用n-1去除。很多生手不喜歡這道式子，因為它要用兩個DO，第一個求各個xi的和，第二個求各個xi - m的平方和；他們也許可以在課本中找到另一道式子（圖十四），它很「漂亮」、因為一個DO就可以把xi的和與xi的平方和算出來，效率上也好一些。



        不過老手們應該會對第二式（圖十四）有所疑慮。因為xi的和也許不會造成溢位，但可能會很大，於是算和的平方時值就更大，從而在存入記憶體時會損失不少在尾部的位數造成四捨五入的誤差；這個問題也會在算xi的平方和時出現。不相信嗎？我們來看看為什麼。

        我們就用三個（n=3）數說明，而且每一個數都有七位、大約是32個數元所能儲存而不致於產生四捨五入的值；我們先用手算：

x1 = 9000000
x2 = 9000001
x3 = 9000002

很明顯地，我們有

x1 + x2 + x3 = 27000003
平均數    = 9000001

如果用第一式（圖十三）算變異數，我們有

變異數={(9000000–9000001)2 + (9000001-9000001)2 + (9000002-9000001)2}/2 = 1

如果用第二式，那就是：

[1] x1+ x2+ x3 = 27000003
[2] (x1+ x2+ x3)2 = 729000162000009
[3] (x1+ x2+ x3)2/3 = 243000054000003
[4] x12 + x22 + x32 = 81000000000000 + 81000018000001 + 81000036000004 = 243000054000005
變異數 = ([4] – [3])/2 = 1

        但若只有七位有效數字可用，[4]的結果會被四捨五入成七位後得到0.2430001 ´ 1015 = 243000100000000，而[3]的則可能也是同一個值，於是變異數為0，這當然是錯的。事實上，Microsoft Excel的早期版本算出來的結果就是0！

        下面是一個用兩道式子算變異數的Fortran 90程式，所有REAL變數都用約定的單精度32數元運算：\n
PROGRAM Variance
IMPLICIT NONE
INTEGER, PARAMETER :: SIZE = 10
REAL, DIMENSION(1:SIZE) :: x
REAL                    :: sumx, sumx2
REAL                    :: mean, var1, var2
INTEGER                 :: n, i
CHARACTER(LEN=*), PARAMETER :: FMT = "(1X,A,F30.10)"

READ(*,*) n, (x(i), i=1, n)
WRITE(*,*) "Number of data items --> ", n
WRITE(*,*)
DO i = 1, n
    WRITE(*,"(1X,A,I2,A,F25.10)") "x(", i, ") = ", x(i)
END DO

sumx = 0.0
sumx2 = 0.0
DO i = 1, n
    sumx = sumx + x(i)
    sumx2 = sumx2 + x(i)**2
END DO
WRITE(*,*)
WRITE(*,FMT) "Sum of x --------------> ", sumx
WRITE(*,FMT) "(Sum of x)**2 ---------> ", sumx**2
WRITE(*,FMT) "[(Sum of x)**2]/n -----> ", sumx**2/n
WRITE(*,FMT) "Sum of x**2 -----------> ", sumx2

WRITE(*,*)
mean = sumx/n
WRITE(*,FMT) "Mean ------------------> ", mean

var1 = (sumx2 - sumx**2/n)/(n-1)
WRITE(*,FMT) "One-Pass Variance -----> ", var1

var2 = 0.0
DO i = 1, n
    var2 = var2+ (x(i) - mean)**2
END DO
var2 = var2/(n-1)
WRITE(*,FMT) "Two-Pass Variance -----> ", var2

END PROGRAM Variance

        如果我們用n=3，9000000、9000001、與9000002做輸入，就會得到下面的結果；其中One-Pass是第二式（只用一個DO），而Two-Pass是第一式（用兩個DO）。程式算出來的資料的和為27000002，與實際的結果差1，但求平均數時由於四捨五入的緣故，結果還是正確的。然而因為各數的和稍有誤差，它的平方（有15位）只有前6位是完全正確的，剩下的9位並不可靠，所以除以3之後的結果也只有七位是可靠的。同理，平方和243000055693312也只有前八位可靠，所以，後者減前者的結果中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

        我們或許會想到用倍精度來增加位數、從而解決這個問題，不過我們仍然可以找到三個足夠大的數使倍精度的程式發生同樣的困擾。所以，根本的做法是採用一個好的算法，而不是肓目地提昇精確度；以這個算變異數的題目為例，就是用兩個DO的寫法。

        上面也提過算平均數時得箕各資料的和，這個和也許會很大、使計算機無法容納所有位數，從而產生四捨五入誤差。但平均數會在最大與最小的值之間，所以在能夠相當精確地儲存最大與最小值的前提下，也應該可以很精確地儲存平均數。事實上，求平均數與變異數都有所謂的「更新」公式，它把目前算出的平均數（或變異數）與下一個資料結合而算出新的平均數（或變異數），於是就避免了求各項和時產生很大的值，從而降低誤差。


3 失去精確位數

        另一個讓人頭痛的浮點數問題是使用減法時可能會失去精確位數。下面的Fortran程式中x與y的值都有16位，在一部有15位精確度的機型中最多只有最後兩位不很精確，它們的差（手算）是0.000000000004111，於是兩個原來有十五位精確位數的值相減後只留下四位，因此精確位數最多也不過四位；換言之，兩個值幾乎相等的浮點數經過減法運算後會失去很多精確位數。

PROGRAM Cancellation
IMPLICIT NONE
INTEGER, PARAMETER :: DOUBLE = SELECTED_REAL_KIND(15)
REAL(KIND=DOUBLE) :: x = 3.141592653589793_DOUBLE
REAL(KIND=DOUBLE) :: y = 3.141592653585682_DOUBLE
REAL(KIND=DOUBLE) :: z
CHARACTER(LEN=*), PARAMETER :: FMT = "(1X, A, E23.17)"

z = x - y
WRITE(*,*) "x     = ", x
WRITE(*,*) "y     = ", y
WRITE(*,*) "x = y = ", z
END PROGRAM Cancellation

        下面是程式的輸出，兩者的差仍然有15位，但這15位中只有最前面三位411是可信的，雖然從第5位的9經過四捨五入後會得到4111的結果，但是在程式執行時我們不一定能夠了解在何處四捨五入，所以程式的計算結果雖然有足夠的15位（倍精度），但可靠的卻只有三到四位；這個失去精確位數的成因，就是兩個幾乎相等的值經過減法運算的結果。另外，在前面討論計算變異數時，我們也談到過失去精確位數的現象，得來的結果全然不可信。

x     = 3.141592653589793
y     = 3.141592653585682
x - y = 4.11093381558203E-12

把兩個幾乎相等的值相減常常不是程式的本意，而是從其它沒有處理好的計算所衍生出來的結果，正因為如此，失去精確位數的成因往往不容易發現。下面就是一個簡單的例子，求一元二次方程式ax2 + bx + c = 0的根（見圖十五）：



程式十分容易，相信在學初級程式設計時大概都寫過，下面是個Fortran 90的例子，所有實數變數都用單精度（約7位精確位數）：

PROGRAM Degree2_Eqn
IMPLICIT NONE
REAL :: a, b, c
REAL :: d, r1, r2
CHARACTER(LEN=*), PARAMETER :: FMT = "(1X,A,E15.7)"

READ(*,*) a, b, c
WRITE(*,FMT) "a      = ", a
WRITE(*,FMT) "b      = ", b
WRITE(*,FMT) "c      = ", c

d = SQRT(b*b - 4.0*a*c)
r1 = (-b + d)/(2.0*a)
r2 = (-b - d)/(2.0*a)
WRITE(*,*)
WRITE(*,FMT) "d      = ", d
WRITE(*,FMT) "Root 1 = ", r1
WRITE(*,FMT) "Root 2 = ", r2
END PROGRAM Degree2_Eqn

        這個程式正確嗎？是的。這個程式寫得好嗎？以初級程式設計而言，很好；就專業角度而言，不好！這個不好並不是因為沒有考慮到a = 0與b2 – 4ac 等等的狀況，而是它很容易發生失去精確度的問題（因為分子中有減法）。避免b2 – 4ac中的減法複雜些，所以姑且不論，剩下來在根號前面的±就是問題所在，我們看看為什麼。

        如果b的值很大而a與c都很小，於是b2就更大，當然b2 – 4ac的結果差不多就是b2，所以b2 – 4ac的平方根就與b非常靠近，經過減法後分子很可能是0。譬如說，若a = c = 1而b = 20000，於是b2 = 400000000，b2 – 4ac = 399999996 = 0.399999996 ´ 109。如果用單精度（7位）計算，這個數會被四捨五入成0.4000000 ´ 109，因此b2 – 4ac的平方根為20000，於是一個根是(-20000+20000)/2 = 0，另一個根則是(-20000-20000)/2 = -20000！這當然不可能是正確的結果，因為ax2 + bx + c = 0兩個根的和是-b/a、兩個根的積是c/a；然而輸入的方程式為x2 + 20000x + 1 = 0，不可能有為0的根！下面就是程式的執行結果。

a      =   0.1000000E+01
b      =   0.2000000E+05
c      =   0.1000000E+01

d      =   0.2000000E+05
Root 1 =   0.0000000E+00
Root 2 = -0.2000000E+05

        如何去掉分子的減號呢？請記位一點，兩個根的積是c/a，所以從根的式子算出一者，令為r1（假設不是0，如果是0就容易了），另一個根就是(c/a)/r1。我們可以這樣想：如果b ，-b就是正數，所以-b+SQRT(b2-4ac)是兩個正數的和；反之，若b > 0，則-b小於0，-b-SQRT(b2-4ac)是兩個負數的和。這樣，減法運算就沒有了。下面的程式片段是根據這個原則稍許修改而成：

d = SQRT(b*b - 4.0*a*c)
IF (b > 0) THEN
    r1 = -b - d
ELSE
    r1 = -b + d
END IF
r1 = r1/(2*a)
r2 = (c/a)/r1

        以下是計算的結果，兩個根是-20000與-0.00005。但是，這兩個根的和是-20000.00005，而不是-20000，結果仍然不十分準確；不過就浮點數運算而言，在單精度（7位）運算下，20000.00005會被四捨五入成20000，所以結果是正確的！雖然有稍許誤差，這總是要比完全錯誤的結果好太多了。

a      =   0.1000000E+01
b      =   0.2000000E+05
c      =   0.1000000E+01

d      =   0.2000000E+05
Root 1 = -0.2000000E+05
Root 2 = -0.5000000E-04

        去掉減法不是件容易的工作，而且還不一定做得到，不過有一個技巧是常常用到的，這就是圖十六的式子。如果式子中出現a-b，我們可以乘上(a+b)/(a+b)得到(a2-b2)/(a+b)，若a2-b2的部份可以化簡到沒有減號，那就大功告成了。但是這個技巧不是萬靈丹，如果不成功就得用其它技巧。要注意的是，只有當a與b的值幾乎相等而會造成損失精確位數時，才需要做這道手續；如果已知a與b的值相差不小而不致於發生損失精確位數時，做這樣的轉換就沒多大意義了。



        讓我們看看如何算圖十七的式子，這道式子只有在x很大時才會出問題。當x很大時，因為四捨五入的原因，x+1的值可能非常靠近x，於是x+1的平方根就與x的平方根幾乎相等甚或相等，有造成用0除的可能性。縱使不是這兩種情況，因為兩個幾乎相等的值相減後會損失精確位數，得來的結果也不會十分準確。



        如果我們用圖十八的手法，圖十七就可以轉化成x+1的平方根與x的平方根的和，當然減法就消失了。不過要注意的是，不是每一道式都可以用這個手法處理，有時還得用其它技巧。



        我們不妨寫個用單精度（7位）的程式來証實上面的分析，這就是圖十九。這張表的第一行是十個從小到大的x值，它們至少都有七位。第二行是x+1的平方根減x的平方根，從表中可以看出，當x為7777777、9999999、與123456789時兩者的差為零，當然算倒數時就會發生用0除的floating-point exception；另外，因為失去了精確位數，從3333333起算出來的差如果不是零就是同一值，當然都不正確。第三行是第二行的倒數，正是圖十七中式子的值；由於第二行的值不準確，第三行的也不準確。第四行是採用圖十八的算法，不妨用手上的計算機驗算，看看是否準確得多！？最後，因為x很大時x+1與x很靠近，當然它們的平方根也很靠近，於是第四行的結果（在x很大時）就差不多是x的平方根的兩倍，這是第五行的結果，從x = 4444444起兩行的結果就幾乎一樣了。從這張表我們可以看出失去精確位數的影響、以及轉換式子去掉減號的重要性。



        下面是幾道練習題，有興趣的朋友不妨試一試：

(1)        (1-cos(x))/sin(x)， x非常靠近0
(2)        (1-cos(x))/x^2， x非常靠近0（提示：考慮半角公式）
(3)        sin(x)-cos(x)， x非常靠近PI/4（提示：考慮半角公式）


4 幾個「無心」但後果嚴重的浮點數過失

        在計算機的歷史中有不少因為浮點數的錯誤而產生嚴重後果的例子，以下舉幾個供大家參考。

        第一，為了避免被敵軍偵測到位置，美國愛國者（Patriot）反飛彈系統在設計時就打算在每一個點只停留並且運轉若干小時後就移往它處。來襲飛彈的速度是用一個浮點數表示，但愛國者飛彈內部的時鐘卻是一個整數，1相當於1/10秒；雖然這個時鐘的解析度不高（許多計算機中1代表1/100秒），但用整數值計時是常見的。這個整數時間會乘上0.1成為一個24個數元的二進位浮點數後再使用，但是0.1的十進位在轉換成二進位時是個循環小數（我們在前面提過），所以轉換後的誤差大約是9.5 ´ 10-8。

        來襲飛彈位置的誤差與來襲飛彈的速度和系統運轉時間的乘積成正比。如果一個愛國者飛彈系統運轉100小時，來襲飛彈的速度是每秒1676公尺，位置的誤差就有573公尺。在第一次波斯灣戰爭時，1991年2月15日伊拉克向在沙烏地阿拉伯達蘭市（Dhahran）的美軍射出一枚飛毛腿（Scud）飛彈，但由於上述的計算誤差使愛國者飛彈沒有成功攔截到來襲飛彈，造成28名美軍陣亡。

        第二，溫哥華股市在每筆交易後都會更新它的股票指數，這個指數有三位小數，在更新之後把第四位之後的位數捨去而留下三位。經過22個月後，股票指數從起初的1000跌到524.881，但是正確算出來的指數卻是1098.881（上漲）。這是什麼時代，還會發生這種簡單的錯誤？很抱歉，這是近代的事，1982年！

        第三，1996年歐州太空機構（European Space Agency）在法屬幾內亞（French Guiana）的Kourou地方發射Ariane 5無人火箭，40秒後火箭失去控制必須引爆摧毀。事後追查原因，居然是一位程式員把一個代表某種控制訊號的64個數元的浮點數存入一個有16個數元的帶符號整數（最大值為32767），所以訊號全都錯亂，使火箭在空中失去控制。專業程式員會犯這種錯，當初是如何測試系統的，是不是覺得很奇怪！？

        第四，1997年9月，美國海軍USS Yorktown號的水兵在操作時不小心輸入了0，於是造成用0除的問題；不但如此，錯誤的結果不斷蔓延，終於把艦上的推進系統關掉（停俥），使軍艦在海面上靜止2小時45分鐘，還好這不是在戰時。想想一個程式不撿查除數是否為0、不檢查輸入是否正確，這是否有點不可思議！？

        話又說回來，我們在寫程式時是否已經儘力避免這些細微但明顯的錯誤？有沒有因為類似的錯誤而使公司災情慘重，或是使公司甚或政府做出錯誤的決策？我常在想一個問題，計算機科系因為開了很多本科的課，而不再講授數值方法（因為本科的學生用不著，真的嗎？我不知道），甚至於連本文的基本觀念都不提，所以本科系學生在浮點數運算的基本知識上很可能是一片空白。非本科的學生大概也學不到這些內容，要用的話就找一套程式庫（IMSL、NAG等）或Matlab等，這當然解決了難的部份，但是在用這些程式庫之前，自己總得先整理資料，如果在過程中出現誤差或錯誤，垃圾進垃圾出就會發生。然而日常程式寫作中總是會用到浮點數，為了避免產生與浮點數有關的錯誤，我們還是得要了解運用浮點數的基本原則的，不論您是本科還是非本科的學生。


5 綜合討論

        前面把幾個使用浮點數時容易出錯的地方分開來用例子做了說明，不過浮點數的問題並不單獨存在，而是與在出問題處所使用的資料有關；換言之，在很多情況下，出問題之前所算出的值可能才是關鍵所在。就用減法為例，出問題所在的減法可能沒有問題，真正問題也許是計算該減法兩個運算元的值的所在。另外，如果減法運算元的值是準確的（譬如輸入值），這個減法也非做不可，所以一些使用大量減法的程式也不會出什麼大差錯（因為它們用原始的輸入值）。有時候，不論如何轉換也無法去掉減法，甚至於去掉了減法後得來的結果仍然不精確，成因很可能是所處理題目本身就很不容易算得精確（ill-conditioned），因此更換演算法才是正途。還有，有些減法是無害的；譬如說，若a與b幾乎相等，但a – b的值卻比c小很多，於是(a-b)+c中的減法就不會有什麼影響（請問為什麼？）。

        四捨五入的誤差又如何呢？幸運的是，在大多數情況下捨去與進位的誤差會相互抵消，而且由捨入誤差造成的困擾很可能只源自若干點、而不一定是由大量反覆運算而得，於是找對了源頭就有辦法找出對策。圖二十是一個十分簡單的例子，您也許在微積分中學到過；Euler証明了這個極限值是e = 2.71828.....，我們打算寫一個程式讓n愈來愈大，從而算出極限值。



        程式十分容易寫如下，我們用n = 10, 100, 1000, 10000, … 等nn個數值代入上式看看是否會得到近似2.71828…的結果。請注意的是，這個程式片段中的四捨五入只侷限在算E的式子而不會蔓延。

n = 10
EE = exp(1.0)
DO i = 1, nn
    E = (1.0 + 1.0/n)**n
    write(*,*) n, E, abs(E-EE)
    n = n*10
END DO

        下面是計算結果，當n到達10000時的近似值為2.718597、與真正e的值差3.1518936E-4，但自此之後算出來的值與真正的e值差距愈來愈大，最後結果居然是1，很明顯的錯誤！問題何在？正是在算不會蔓延的E時的捨入誤差。前面提過，0.1的十進位轉換成二進位或十六進位時會是循環小數，當然1/10n也是如此，於是存入記憶體就會有捨入誤差。不過根本問題並不在此，而是在加上1的結果！如果計算機有七位精確位數，1/107 = 0.0000001，加上1之後為1.0000001，寫成標準型是0.10000001 ´ 101，四捨五入到七位後得到0.1 ´ 101，這就是最終結果為1的原因。所以，捨入誤差並不一定得蔓延才會產生問題，關鍵在找出產生問題的所在！以這個程式而言，解決之道就是改寫這道式子來降低捨入誤差（這是個頗有挑戰性的問題，數學好的朋友不妨一試）。

10    2.593743        0.12453866
100    2.7048113       0.013470411
1000    2.7170507       1.2309551E-3
      10000    2.718597        3.1518936E-4
     100000    2.7219622       3.6804676E-3
    1000000    2.5952267       0.12305498
   10000000    3.2939677       0.575686
  100000000    1.0             1.7182817
1000000000    1.0             1.7182817

        另一個看起來相當「明顯」的做法是改用倍精度來增加精確位數，這是可行的；事實上，當n = 100000000時，倍精度運算得到2.7182817983473577，與真正e的誤差是3.011168780986395E-8！但是，用倍精度的計算只得到8位精確位數而不是15位，是不是有點奇怪呢？用這個方法可能得到15位精確位數嗎？

        遺憾的是，提高精確度並不是萬靈丹，無法解決所有捨入誤差所帶來的問題。您不妨找個計算機、輸入一個正數，再按開方鍵若干次得到一個數，接著按同樣多次的平方鍵，看看會有什麼結果。一般而言，結果會與輸入值相差不大，不過若按開方鍵次數足夠多，您會發現結果會愈來愈靠近1，於是再平方回來時與原本的輸入就有很大的差距；當然，一旦開方的結果變成1，不論如何平方都不可能得到原來的值。這個現象在任何計算機、使用任何精確度都成立，所差的就是何時變成1而已。

        我們也可以寫一個程式驗証這個現象；下面是程式的（部份）輸出，這個程式求2的平方根若干次，再把結果求平方相同次數。首先，我們開方十次再平方十次，於是第十次開方得到1.0006771088：

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

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

    0   1.0006771088
       ………
    5   1.0218964815
       ………
   10   1.9999581575

如果把開方次數增加到25次，下面是結果，到第23次結果就是1，當然平方回來就不可能得到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

   改用倍精度又如何？因為單精度（七位）差不多在開方23到25次就得到1，有十五位的倍精度大約用50（或多幾）次開方也會得到1。所以，就這個問題而言，增加精確度只是延後問題出現的時間而已，當然，如果在問題發生之前就己經得到可用的結果，那就提到精確度吧！


6 結論

        前面用不少例子說明使用浮點數的最基本課題，但我們沒有做嚴謹的定義、也沒有用嚴格的論証，只用到最基本寫程式的知識，這樣做的用意是希望一般程式員（不論是入門或資深）甚至於學生都看得懂。以下是幾點最粗淺的建議：

第一，儘可能避免使用兩個值幾乎相同的運算元的減法，除非這兩個運算元的值是精確的、或是減法的結果對往後的運算影響不大。\n第二，留神溢位或近似溢位的狀況。一般而言，如果算出來的中間結果比最終結果（與輸入資料）大很多或小很多，就是可能會產生溢位的徵候，在這個情況下不妨改用其它式子或把值做比例放大或縮小，以避免溢位。
第三，留神捨入誤差，特別是把兩個相差極大的值相加或相減時，小的那個值可能完全沒有貢獻，從而得到不正確的結果。
第四，不要隨便抄課本的演算法或程式在緊要的應用中使用。課本中的演算法與程式的重點是在解釋觀念，除非作者很明白地指出該演算法或程式已經考量到數值計算的穩定度（stable）等課題，或者是該演算法或程式沒有穩定度的問題；如果不是如此，當程式的輸入比較「怪異」時，該演算法或程式的輸出很可能不正確。

        最後，希望這篇長文對您在未來使用浮點數上頭有點助益。



本文的引用方式如下：
file:///C:/DOCUME~1/CHING-~1/LOCALS~1/Temp/msohtml1/01/clip_image001.gif
          非學術性請用
          冼鏡光 DCview.com達人部落格 本文的網址
file:///C:/DOCUME~1/CHING-~1/LOCALS~1/Temp/msohtml1/01/clip_image001.gif
          學術性請用
          冼鏡光，使用浮點數最最基本的觀念，DCview達人部落格，本文的網址，讀取時間（本文標題中最新更正時間）


本文来自ChinaUnix博客，如果查看原文请点：http://blog.chinaunix.net/u2/63543/showart_697260.html