;把兩指數相加再減去127得到的就是結果指數的原始值

        ;本來需要再減去23的，因為此時的小數點停留在第22與23位之間，

        ;但是由于后面的邏輯右移指數又需要加上23，所以就省略了

        shrd  eax,edx,23

        ;add  esi,23

        test  eax,1000000h;測試第24位是否為1，原因見上文

        jz    @f

        shr   eax,1

        inc   esi

        @@:

        xor   eax,ecx;去掉隱含的1

        shl   esi,24

        mov   ebx,floatU

        mov   ecx,floatV

        xor   ebx,ecx;設置結果的符號位

        shl   ebx,1;取出符號位放入CF標志位中

        rcr   esi,1;從CF中得到符號位

        or    eax,esi

        ret

_fmul   endp

 

最后我們來看如何寫浮點除法的子程序：

算法C：規定兩個單精度浮點數u和v，進行u÷v運算。

C1.分離u和v的指數部分與有效數字，有效數字要補上舍去的1，然后把u的有效數字放入a中，v的有效數字放入b中。

C2.把a左移24位，這是為了與b相除時，可以得到24或25位有效數字。

C3.把u的指數和v的指數相減，然后加上127這個偏移值，結果存入s。

C4.計算a/b。

C5.規格化浮點數，因為a左移24位，因為結果的低24位都是小數部分，而正常結果只有23位有效數字，最高位是隱含位，s要減去1，檢測a/b的結果的第24位（從0開始數）是否為1，不是1的話，那第23位就是1了，如果是，那么a/b的結果就邏輯右移1位，同時把s也加1。

C6.結果的符號位可以有u和v的符號位“異或”處理得到。

C7.合并符號位，指數，有效數字。

最后來看看我寫的用匯編語言處理浮點除法的子程序：

_fdiv   proc    floatU,floatV

        mov   ecx,800000h

        mov   eax,floatU

        mov   esi,floatU

        mov   ebx,floatV

        mov   edi,floatV

        and   eax,7fffffh

        and   ebx,7fffffh

        or    eax,ecx

        mov   edx,eax

        shl   eax,24;把u的有效數字左移24位,這樣可以得到25位

        or    ebx,ecx

        shr   edx,8;或24位的結果,因為結果的低24位都是小數部分,

        ;而正常結果只有23位小數,最高位是隱含位,所以在下面的

        ;結果指數要減去1

        div   ebx

        shl   esi,1

        shl   edi,1

        shr   esi,24

        shr   edi,24

        sub   esi,edi

        add   esi,126;下面兩條指令的合成

        ;add  esi,127

        ;因為u和v的指數都加了127,當兩個相減后抵消了,所以要加127

        ;dec  esi;指數減去1(原因見上注釋)

        test  eax,1000000h

        jz    @f

        shr   eax,1

        inc   esi

        @@:

        xor   eax,ecx;去掉隱含的1

        shl   esi,24

        mov   ebx,floatU

        mov   ecx,floatV

        xor   ebx,ecx;設置結果的符號位

        shl   ebx,1;取出符號位放入CF標志位中

        rcr   esi,1;從CF中得到符號位

        or    eax,esi

        ret

_fdiv   endp

 

看懂了單精度的算數運算程序，那么寫雙精度的也就沒問題了，思路是完全一樣的。

 

2.浮點算術的精確度

就本質來說，浮點算術是不精確的，而且程序員們很容易就會濫用它，從而使計算的答案幾乎全尤“噪聲”所組成。數值分析的主要問題之一，是確定某個數值方法的結果的精度如何。許多認真的數學家曾試圖對于浮點操作的序列給出嚴格的分析，結果發現這個任務實在大得驚人，只好做些含糊其詞的論證來聊以自慰。

當然，對于誤差分析技術的徹底考察，超出了本文的范圍，但我們將在本文中研究浮點技術的誤差的某些特征。我們的目標是來發現以什么樣一種方法實施浮點算術，使得誤差傳播的合理分析盡可能簡便。

浮點運算特征的一個粗糙的（但相當有用的）方式，可以以“有效位”或相對誤差的概念為基礎。粗略的講，浮點乘法和除法的運算并不使相對誤差擴大許多；但近乎相等的量的加法（當u為正數，v為負數時）則可以大大增加其相對誤差。所以我們決定從加法還有減法中去尋找主要的精度損失，而不是乘法和除法。

還有一點，這一點似乎不太合常人的想法，但的確需要適當的理解，因為“不好”的加法和減法是以完全的精度被執行的。

現在我們以X=x(1＋E)來表示在一臺計算機內的一個確切的實數x，則E=(X-x)/x稱為該近似值的相對誤差。這里的E是單詞ERROR（誤差）的首寫字母

浮點加法可能一個不可靠的結果之一，是破壞了結合律：

(u＋v)＋w≠u＋(v＋w)              (1)

例如：

我們為方便而使用十進制來表示，假如有效位數有8位。下面的E是指數的意思表示為×10^*。

(1.1111113E7＋(-1.1111111E7)＋7.5111111E0

= 2.0000000E0＋7.5111111E0

= 9.5111111E0

而：

1.1111113E7＋((-1.1111111E7)＋7.5111111E0)

= 1.1111113E7＋(-1.1111103E7)

= 1.0000000E1

發現了嗎？兩個結果相差將近0.5，所以程序員們必須特別小心，不得暗中假定結合律的正確性。

在開頭我們提到銀行通常使用定點BCD數，就是因為這個原因，因為有時浮點可以精確到厘，而有時卻不能精確到元，通常銀行做得最多的就是加法運算，浮點運算如此的誤差，顯然是不能接受的！

還好，交換律是成立的：

u＋v = v＋u                       (2)

可別小看這一定律，她對程序設計和程序分析來說，在概念上不失為一種有價值的財富。

我們應該尋找一些為浮點數的＋、－、×、÷所滿足的重要定律；而且，我們應該把浮點程序設計成為滿足盡可能多的數學定律。顯然，如果有更多的公理成立，則就更容易寫出好的程序來。

現在我們來看看其它的一些基本定律，它們對于規格化的浮點運算是正確的：

u－v = u＋(-v)                    (3)

-(u＋v) = -u＋(-v)                (4)

當 v = -u 時 u＋v = 0             (5)

u＋0 = u                          (6)

進一步得到：

u－v = -(v－u)                    (7)

這些定律我們可以從算法A中導出。

另外：

如果u≤v 則 u＋w≤v＋w            (8)

仔細分析算法A的設計原理可以說明這一點。

浮點運算應滿足：

u＋v = round(u+v)

u－v = round(u-v)

u×v = round(u*v)

u÷v = round(u/v)                 (9)

round(x)代表x的最好的浮點數近似。

顯然  round(-x) = -round(x)      (10)

x≤y  蘊含  round(x)≤round(y)   (11)

根據這些基本性質。我們還可以寫出：

u×v = v×u

(-u)×v = -(u×v)

1×v = v

當u = 0或v = 0時 u×v = 0

(-u)÷v = u÷(-v) = -(u÷v)

0÷v = 0

u÷1 = u

u÷u = 1

如果u≤v和w＞0，則 u×w≤v×w 且 u÷w≤v÷w

當v≥u＞0時，w÷u≥w÷v

當u＋v = u+v，則(u＋v)－v = u；且如果u×v = u*v ≠ 0，則(u×v)÷v = u

當然，還有若干熟悉的代數規則仍然不存在。

浮點乘法的結合律不是嚴格正確的。而且稍微更糟的是＋和－之間的分配率不成立。設u = 2.0000000E4，v = -6.0000000E0，w = 6.0000003E0 則

(u×v)＋(u×w) = -1.2000000E5＋1.2000001E5 = 1.0000000E-2

u×(v＋w) = 2.0000000E4×3.0000000E-7 = 6.0000000E-3

所以

u×(v＋w)≠(u×v)＋(u×w)       (12)

當b是浮點數的進制時，我們有b×(v＋w) = (b×v)＋(b×w)，因為

round(bx) = b round(x)          (13)

嚴格來講，我們這里的恒等式和不等式隱含地假定不出現指數的上溢和下溢。當|x|太大或太小是函數round(x)是無定義的，而且像(13)這樣的等式僅當兩邊有定義時才成立。

下面是一個在現代教科書上求標準差的公式：

    (14)

經驗不多的程序員經常發現他們用此公式取得是一個負數的平方根！用浮點算術平均值和標準差的好的多的方式是使用遞推式：

，           (15)

，       (16)

其中2≤k≤n，。

很明顯，在這個公式中決不能是負的，下面的例子將使我們看到這個問題的嚴重性，以至于我們不得不采用(15)和(16)的公式。

現在我們假設有n=1000000，且對于所有的k，我們有=1.1111111E0，那么在8位精度的十進制下我們來計算1000000個1.1111111E0的和，我們先使用最原始的公式：

當：



因此：

這里我們取時使用最接近的近似：1.2345679E0

當n=1000000時：



乘以n后得到1.2247821E12。

1.2090991E6平方后得到1.2301008E12

結果發現我們取的是：

(1.2247821E12－1.2301008E12)÷(n(n-1)) = -5.3187053E-3的平方根！

在這種情況下使用(15)和(16)的公式將不會有這個問題。

 

雖然傳統的代數法則不一定正確，但也不是偏離了太離譜。當在b進制系統中，x假如是真值，我們用e來表示指數(不是很好，因為容易與自然底數e混淆)，那么2^(e-1)≤x＜2^e。

設round(x) = x + ρ(x)，其中≤，p表示有效位數

且有round(x) = x(1+δ(x))

那么與x無關的相對誤差的界δ(x)：

≤≤＜

由此式我們得到：

u＋v = (u+v)(1+δ(u+v))

的相對誤差。

我們現在使用此式來估計乘法結合律的誤差。

在一般情況下(u×v)×w≠u×(v×w)，但可以肯定比加法的結合律好得多。

不考慮上溢和下溢對于

δ1,δ2,δ3,δ4,

我們有：

(u×v)×w = ((uv)(1+δ1))×w = uvw(1+δ1)(1+δ2)

u×(v×w) = u×((vw)(1+δ3)) = uvw(1+δ3)(1+δ4)

因此：



其中，

那

因此乘法的結合律的相對誤差大約在以內。

程序員在做浮點運算時應該避免測試兩個被計算的值是否相等，雖然這極不可能出現。

好，關于浮點算術，這期的內容就寫到這里，下期還有更精彩的，待續！



--------------------------------------------------------------------------------

歡迎訪問AoGo匯編小站:http://www.aogosoft.com/