浮點算術  作者：ham 于2008-1-17上傳

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


作者：ham

參考資料：The Art of Computer Programming 《計算機程序設計藝術》

[美]DONALD E.KUNTH著

 

在本文中，我們將通過分析奠定浮點數計算基礎的內部機制，來研究對這種數進行運算的基本原理。也許許多讀者對這樣的課題興趣不大，因為他們的計算機內已經含有浮點的指令了（FPU）。但是，事實上，每一個熟練的程序員都應當有浮點運算的基本步驟在執行期間是如何進行的知識。這一課題并不是像大多數人所想像的那樣全然微不足道，它包含多得驚人的有趣信息。

1.浮點計算

A.浮點記法

在科學計算中還有一種計數方法叫定點表述，在這種情況下，程序員知道在正在操作的數中，假定小數點是在什么位置上，這種方法在銀行用了比較多（定點BCD數），為什么要用定點數呢，我們分析浮點數的精確度時將會看到。但是從許多目的的來看，當程序運行時，如果讓小數點靈活“浮動”，并讓每個數帶一個相應的小數點的位置標志，有著很大的方便性。這一想法，在科學計算中已經用了許多年。

在本文中我們大多數情況下考慮使用二進制來表示浮點數，但在有時為了方便理解將會出現十進制表示的浮點數。

現代計算機和計算機程序按照IEEE在1985年制定的標準來處理浮點數，這個標準也為ANSI(the American national standards institute，美國國家標準局)所認可。ANSI/IEEE Std 754-1985稱作IEEE二進制浮點數算術運算標準。它并不像一般標準那樣長，只有18頁，但卻奠定了以方便的方式編碼二進制浮點數的基礎。

IEEE浮點數標準定義了兩個基本格式：單精度格式，需要4個字節（32位）；雙精度格式，需要8個字節（64位）。

除了這兩種表示方法，其實在我們CPU內部集成FPU中處理的都是擴展精度浮點數，這種浮點數用10個字節（80位）表示，這種計法有先天的優點。

浮點數有三部分構成，我們先來看看單精度浮點數：

符號位（1位）
 指數（8位）
 有效數字（23位）
 

接著是雙精度浮點數：

符號位（1位）
 指數（11位）
 有效數字（52位）
 

符號為0表示此數是正數，為1表示是負數。

在這里我們看到單精度浮點數有效數字是23位，事實上它表示的是24位，這個我們過會兒會看到原因，同理雙精度有效數字實際為53位。

單精度指數可以表示的范圍是0～255，在實際計算時需要減去127，才能表示為真實的指數（因為指數采用的是移碼表示，可能是方便判斷特殊情況和簡化FPU的設計），還有指數為0和255，也就是全0與全1時用于特殊用途，這樣在平時指數的范圍就是1～254，減去127就是-126～127了。

雙精度指數可以表示的范圍是0～2047，實際計算時需要減去1023，同樣指數全0與全1時用于特殊用途。

在這種浮點表示結構中它實際表示的意思就是：

單精度真值 = -1^符號位×1.有效數字×2^(指數-127)

雙精度真值 = -1^符號位×1.有效數字×2^(指數-1023)

這里我們看到在小數點前有個1，這個1在浮點結構中是隱含的，因為在規格化的二進制浮點數中的最高位始終是1，所以省略了，這樣實際精度就比“有效數字”多一位了。

由于單精度的浮點數實際有效數字是24位，這樣它的實際精度就是1/(2^24) = 1/16777216，這說明什么呢？這說明它只能完全的表示十進制的7位有效數字，因為最大只有1/10000000的分母小于16777216，8位有效數字時就不能完全表示了。

同樣雙精度的浮點數精度為1/(2^53) = 1/9007199254740992，這樣它能完全表示十進制的15位有效數字。

由上面的分析我們知道單精度浮點的指數范圍是-126～127，那么在規格化的情況下表示正數的范圍就是1.0……0B×2^-126 ～ 1.11……1B(共24個1)×2^127，化成十進制大約是±1.175484×10^-38～±3.028235×10^38，雙精度就是1.0……0B×2^-1022 ～ 1.11……1B（共53個1）×2^1023，化成十進制大約是±2.11507385807201×10^-308～±1.79769313486232×10^308。

如何表示一個確切的數0呢？這是一個特殊的情況：

如果指數的所有二進制位為0，且所有的二進制有效數字為0（單精度23位，雙精度52位，不包括省略的1），則表示數0，當然符號可以為1，這樣就是 -0。

如果指數的所有二進制位為0，且二進制有效數字不為0，那么數仍然有效，只不過是一種非規格化的數，它表示為：

單精度：-1^符號位×0.有效數字×2^(-127)

雙精度：-1^符號位×0.有效數字×2^(-1023)

如果指數的所有二進制位為1，且所有的二進制有效數字為0，則根據符號位來判斷是負無窮大，還是正無窮大。

如果指數的所有二進制位為1，且二進制有效數字不為0，則這個不是一個合法的浮點數，簡寫成NaN。

B.規格化計算

現在使用中的大部分浮點程序幾乎全部處理規格化的數：對程序的輸入假定為規格化的，而且輸出總是規格化的。這樣我們獲得了速度、一致性，以及在計算中對相對誤差給出相對簡單的界的能力。

現在來詳細研究浮點數算數的操作。同時我們可以考慮構造這些操作的子程序，這里假定我們的計算機沒有內部的浮點硬件（FPU）。

浮點算術的匯編語言子程序通常都編寫得非常依賴于機器，它們使用該計算機的許多最原始的特性，我們在這里就只討論如何在Intel 80x86 CPU上編寫浮點運算的子程序。

由于浮點計算的不精確性，我們使用

＋、－、×、÷

來區別于精確的+、-、*、/。

現在我們來討論單精度浮點加法，這是我們討論的第一個算法（而且是最最困難的一個）。

算法A：規定兩個單精度浮點數u和v，進行u＋v運算，這里只考慮一般情況，不考慮溢出和特殊狀態。

A1.分離u和v的指數部分與有效數字，有效數字要補上舍去的1。

A2.把指數大的有效數字的放在a中，把指數小的有效數字放到b中，然后對b進行邏輯右移，移動的位數等于a的指數減去b的指數。

A3.檢測u與v的符號位是否相同，相同就跳到A7

A4.計算a - b的值，如果結果為負就把結果取反。

A5.對結果規格化，也就是把有效數字的最高位的1移動到第23位（從0開始），然后去掉第23位作為隱含值，同時調整指數，最后結果的符號位由指數大的那個數的符號位來確定，當a-b的值為負數時，就將結果符號取反。

A6.跳到A9。

A7.計算a + b的值。

A8.對結果規格化，也就是把有效數字的最高位的1移動到第23位（從0開始），然后去掉第23位作為隱含值，同時調整指數，符號位就根據u或v的符號卻定，如果負，那么結果就是負，如果正，那么結構就是正。

A9. 合并指數、有效數字、符號

下面是我用匯編語言描述的子程序，如發現不當的地方指正一下：

這里的floatU就是u，floatV就是v

eax就是a,ebx就是b

返回值在eax中

 

_fadd           proc   floatU,floatV

        mov   eax,floatU

        mov   ecx,eax

        mov   edi,eax;用來保存u的符號位

        mov   esi,eax

        mov   ebx,floatV

        mov   edx,floatV

        shl   ecx,1;去除u的符號位

        shl   edx,1;去除v的符號位

        shr   ecx,24;得到u指數

        shr   edx,24;得到v指數

        push  ecx;我們需要這個來設置結果的指數,

                   ;指數由指數大的數提供,我們先保存u的指數

        and   eax,7fffffh;得到u的23位有效數字

        and   ebx,7fffffh;得到v的23位有效數字

        or    eax,800000h;補上u隱含的1,此時eax就是a

        or    ebx,800000h;補上v隱含的1,此時ebx就是b

        sub   ecx,edx;此時的ecx就是兩個數的指數差,用于對b邏輯右移

        jae   @f;當u的指數小于v的指數時,我們必須交換eax,ebx

                  ;讓指數大的有效數字保存在a中,然后專心對b邏輯右移

                  ;并且需要對ecx取反,方便b移位

        xor   ecx,-1

        xchg  eax,ebx;交換a與b，讓指數大的有效數字保存在a中

        inc   ecx

        mov   edi,floatV;此時指數大的數就是v了,我們需要它的符號位

        mov   [esp],edx;現在我們替換掉剛剛保存的u的指數,保存v的指數

        @@:

        xor   esi,floatV;查看u和v的符號位是否相同,相同就直接加法

                          ;不同就做減法

        test  esi,esi

        jns   lp

        ;-------------------------

        ;進行減法計算

        shr   edi,31;得到結果的符號位

        pop   edx;得到結果指數的原始值

        shr   ebx,cl

        sbb   eax,ebx

        mov   ebx,800000h

        jnc   @f

        xor   eax,-1

        xor   edi,1;a-b為負數時對結果的符號位取反

        inc   eax

        jmp   @f

        ;===================

        lp1:

        ;對結果規格化

        shl   eax,1

        dec   edx

        jz    lp2;這個主要用來判斷結果是否為0,

                   ;如果指數為0,那么有效數字也不需要規格化了

        @@:

        test  eax,ebx

        jz    lp1

        ;===================

        xor   eax,ebx;去掉隱含的1

        lp2:

        shl   edi,31

        shl   edx,23

        or    eax,edx

        or    eax,edi

        ret

        ;--------------------------------

        ;進行加法計算

        lp:

        shr   ebx,cl

        adc   eax,ebx

        pop   ecx;得到結果指數的原始值

        test  eax,1000000h

        jz    @f

        ;對結果規格化,因為假如第24位為0,那么第23位肯定為1

        shr   eax,1

        inc   ecx

        @@:

        mov   edi,floatV

        and   eax,7fffffh;去掉隱含的1

        shl   ecx,24

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

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

        or    eax,ecx

        ret

_fadd           endp

 

減法其實就很簡單了，比如u－v，就直接把v的符號位取反，然后調用加法子程序。

_fsub   proc    floatU,floatV

        mov     eax,floatV

        mov     ecx,80000000h

        xor     eax,ecx;把v的符號位取反

        push    eax

        push    floatU

        call    _fadd

        ret

_fsub   endp

現在我們再開始分析單精度浮點數的乘法運算，乘法運算很簡單，在這里我們不考慮溢出以及特殊狀態。

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

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

B2.把u和v的指數部分相加，然后減去127，因為u和v的指數都加了127，所以要減去一個127，吧結果存入s。

B3.計算a*b。

B4.規格化有效數字：因為a*b結果最多會變成47位或48位，所以對48位先邏輯右移23位，然后檢測第24位（從0位開始）是否為1，不是1，那么第23位肯定是1了，如果是1，那么把結果再邏輯右移一位，并把指數s加1，最后去掉第23位的隱含1。

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

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

 

現在來看看我寫的用匯編語言處理浮點乘法的子程序:

_fmul   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

        or    ebx,ecx

        mul   ebx

        shl   esi,1

        shl   edi,1

        shr   esi,24

        shr   edi,24

        lea   esi,[esi+edi-127]

        ;sub  esi,23