十進制浮點數的表示方法

使用十進制浮點數,可以避免二進制浮點數與我們習慣的十進制數之間的表示誤差.這個在金融領域是非常重要的.但是計算機基本都只能對二進制浮點數進行計算,也就是IEEE754格式表示的浮點數.很多程序都會自己模擬十進制浮點數的計算.為了統一,IEEE754做了擴展,包括了十進制的浮點數.

IEEE 754-2008里面規定了十進制浮點數的一些規范.不過里面沒有說具體的二進制表示方法.只是規定了32位,64位,128位的十進制浮點數的表示范圍和有效位數. 因為具體一個浮點數的二進制里面每個位表示啥,都是每個機器自己決定的.不需要跟外界一致.只是在傳輸的時候要保證數據的精度和范圍一致就行了.下表來自wikipedia,列出了每種浮點數的有效位數,指數的范圍.

實際的系統中,十進制浮點數有兩種表示方法,分別是Densely Packed Decimal(密集十進制數)和Binary Integer Decimal(二進制整數表示的十進制數).

DPD表示方便轉換成十進制的浮點數字符串,但是需要專門的計算單元來做計算,軟件模擬比較麻煩.

而BID表示更直觀,轉換到二進制會比較容易.很方便用二進制的整數運算單元來計算.

所以Power6上有了硬件的十進制浮點計算單元,就用DPD表示.而在x86 x64 cpu上沒有十進制計算單元, 各種軟件實現的十進制浮點庫默認大都用BID方式表示.比如Intel就實現了一個開源的c 語言的十進制浮點數庫。http://software.intel.com/en-us/articles/intel-decimal-floating-point-math-library/

十進制浮點的意義，在于更符合人們的習慣，比如下面的例子

1. #include < stdio . h >
   

2.   int main ( )
   

3.   {
   

4.      double a = 7 . ;
   

5.      double b = 0 . 00007 ;
   

6.      printf ( "%d/n" , a = = b * 100000 ) ;
   

7.   }

正確的輸出應該是1，但是實際的輸出結果是0，在做相等比較的時候，還不得不考慮一下這個誤差了。而某些時候誤差會在計算過程中累計，變成比較明顯的錯誤了。

如果用intel的十進制浮點庫賴做這個計算，結果就會不同了。intel這個庫明顯還在試驗階段，用起來比較麻煩。

1. int main ( )
   

2.   {
   

3.      Decimal64 a , b , c ;
   

4.      _IDEC_round my_rnd_mode = _IDEC_dflround ;
   

5.      _IDEC_flags my_fpsf = _IDEC_allflagsclear ;
   

6.      a = bid64_from_int32 ( 7 ) ;
   

7.     b = bid64_from_string ( "0.00007" , my_rnd_mode , & my_fpsf ) ;
   

8.      c = bid64_mul ( b , bid64_from_int32 ( 100000 ) , my_rnd_mode , & my_fpsf ) ;
   

9.     printf ( "%d/n" , bid64_quiet_equal ( a , c , & my_fpsf ) ) ;
   

10.      return 0 ;
    

11.   }

使用和double位數相同的Decimal64，結果就是1了。這里顯然不是精度的問題，而是十進制浮點數能絲毫不變的表示十進制的小數。

我們可以看到這里使用的是BID的表示方法。函數名前面都帶個bid前綴。

接下來，我們來具體看看BID的表示方法，我們可以把剛才程序中的a和c按照十六進制輸出

    printf("%llx/n%llx/n",a,c);

結果是

31c0000000000007
31200000000aae60

可見，兩個相等的十進制浮點數的BID表示不一定是相同的。也就是說，一個數有多種表示方法。

a的表示里，最低位的16進制數就是7，而c的表示里，最低的5位15進制數aae60，其實就是十進制的700000?？磥磉@后面的就是有效數字部分了。查一下BID的表示方法，還是比較復雜的，有6種情況。最高位是符號位，這里當然是0.符號位后面的兩位是00，01，或10時，64位BID每個位的意義是這樣的，s后面的2位和之后的8位是指數部分，之后53位T和t都是有效數字部分

s 00eeeeeeee TTTtttttttttttttttttttt tttttttttttttttttttttttttttttt
s 01eeeeeeee TTTtttttttttttttttttttt tttttttttttttttttttttttttttttt
s 10eeeeeeee TTTtttttttttttttttttttt tttttttttttttttttttttttttttttt
而如果符號位后面的兩位是11，那么每一位的意義是
s 11 00eeeeeeee (100)Ttttttttttttttttttttt tttttttttttttttttttttttttttttt
s 11 01eeeeeeee (100)Ttttttttttttttttttttt tttttttttttttttttttttttttttttt
s 11 10eeeeeeee (100)Ttttttttttttttttttttt tttttttttttttttttttttttttttttt
這時，有效數字前面就加了隱含的100.
這個BID表示的數的值就是 (-1)^S *T*10^(E-398) ,其中T 是實際的有效數字（就是說如果有隱含的100需要加上后計算），E是指數，T，E都是2進制表示的
還是回到我們的例子
a的二進制數
0 0110001110 00000 00000000 00000000 00000000 00000000 00000000 00000111
指數部分就是0110001110，也就是398，所以a就是 7*10^(398-398) ,也就是7

而c的二進制是
0 0110001001 00000 00000000 00000000 00000000 00001010 10101110 01100000
指數部分是 0110001001，也就是393， 所以c的值是 700000*10^(393-398)， 還是7.
這就能看明白為啥同樣是7，二進制表示卻不同。這也是十進制浮點和二進制浮點一個不同之處，十進制浮點沒有規定一定要是哪一種表示。這也給相等比較帶來了一點麻煩。

power6 里面內置了十進制浮點計算單元，而power6上面的編譯器也就支持了內置的十進制浮點類型。前面已經說了，power上面的十進制浮點才用的是DPD表示方法。還是看個程序吧。下面這個程序在一個使用Power6的P520機器上，操作系統是AIX5.3 ML6， 用xlc 10.2編譯。_Decimal64就是64位的十進制浮點。

1. int main ( int argc , char * * argv )
   

2. {
   

3. long i , count ;
   

4. double dfund , dinterest ;
   

5. _Decimal64 Dfund , Dinterest ; / * 定義十進制浮點類型的變量 * /
   

6.  long long value ;
   

7.  union trans {
   

8.   _Decimal64 dv ;
   

9.    int av [ 2 ] ;
   

10.   } transTemp ;
    

11.  dfund = atof ( argv [ 1 ] ) ;
    

12.  dinterest = atof ( argv [ 2 ] ) ;
    

13.  Dfund = atodecimal ( argv [ 1 ] ) ;
    

14.  Dinterest = atodecimal ( argv [ 2 ] ) ;
    

15.  count = atoi ( argv [ 3 ] ) ;
    

16. / * 下面把_Decimal64 類型的Dinterest轉換成兩個int，然后按照十六進制格式顯示 * /
    

17.  transTemp . dv = Dinterest ;
    

18.  printf ( "value=%#x,%#x/n" , transTemp . av [ ] , transTemp . av [ 1 ] ) ;
    

19.  printf ( "double  fund=%20.10f interest=%40.30f/n" , dfund , dinterest ) ;
    

20.  printf ( "Decimal fund=%20.10Df interest=%40.30Df/n" , Dfund , Dinterest ) ; / * printing them with the new printf specifiers * /
    

21.   for ( i = ; i < count ; i + + ) {
    

22.          dfund = dfund * dinterest ;
    

23.          Dfund = Dfund * Dinterest ; / * performing maths * /
    

24.   }
    

25.  printf ( "Print final funds/n" ) ;
    

26.  printf ( "double  fund=%30.10f/n" , dfund ) ;
    

27.  printf ( "Decimal fund=%30.10Df/n" , Dfund ) ;
    

28. }

其中 atodecimal是自己寫的一個幫助函數

1. _Decimal64 atodecimal ( char * s )
   

2. {
   

3. _Decimal64 top = , bot = , result ;
   

4. int negative = , i ;
   

5.          if ( s [ ] = = ' - ' ) {
   

6.                 negative = 1 ;
   

7.                 s + + ;
   

8.          }
   

9.          if ( s [ ] = = ' + ' ) s + + ;
   

10.           for ( ; isdigit ( * s ) ; s + + ) {
    

11.                 top = top * 10 ;
    

12.                 top = top + * s - ' ' ;
    

13.          }
    

14.          if ( * s = = ' . ' ) {
    

15.                 s + + ;
    

16.                  for ( i = strlen ( s ) - 1 ; isdigit ( s [ i ] ) ; i - - ) {
    

17.                         bot = bot / 10 ;
    

18.                         bot = bot + ( _Decimal64 ) ( s [ i ] - ' ' ) / ( _Decimal64 ) 10 ;
    

19.                  }
    

20.          }
    

21.         result = top + bot ;
    

22.          if ( negative )
    

23.                 result = - result ;
    

24.         return result ;
    

25. }

這個程序用xlc 10.2編譯時，跟上參數表示使用硬件十進制浮點。不過這樣會導致編譯出來的可執行文件在power5以前的cpu上無法運行。

運行的時候輸入參數 ./dfp_hw 1 1.00000091 6000000

dfp_hw是程序的名字，1就是程序里面的 fund，1.00000091是interest，也就是利息，6000000是count，輸出結果：

value=0x22180000,0x800001b
double  fund=        1.0000000000 interest=        1.000000910000000020616539586630
Decimal fund=        1.0000000000 interest=        1.000000910000000000000000000000
Print final funds
double  fund=                235.0968403429
Decimal fund=                235.0968403137

可以看到用double存儲利息，再輸出，就不再是1.00000091了，后面有一點誤差。而用_Decimal64存儲輸入結果，再輸出，是一點誤差都沒有。

然后把interest乘6000000次，也就是1.0000091的6000000次方，輸出的結果誤差就比較明顯了。用windows自帶的計算器可以驗證，_Decimal64的結果是正確的。

現在來看看1.00000091的二進制表示。也就是0x22180000,0x800001b，注意這里這個power機器是大端的，所以前面以前是高4字節，后面是低4字節。連起來看，就是0x22180000 0800001b也就是

00100010 00011000 00000000 00000000 00001000 00000000 00000000 00011011

DPD表示方法也比較復雜，從高位開始看，第一位還是符號位0，DPD的規定如果符號位后面的兩位是00，01，或者10，那么每一位的意義如下

s 00 mmm (00)eeeeeeee (mmm)[tttttttttt][tttttttttt][tttttttttt][tttttttttt][tttttttttt]
s 01 mmm (01)eeeeeeee (mmm)[tttttttttt][tttttttttt][tttttttttt][tttttttttt][tttttttttt]
s 10 mmm (10)eeeeeeee (mmm)[tttttttttt][tttttttttt][tttttttttt][tttttttttt][tttttttttt]
其中，e是指數，e的表示方法跟前面的BID方式很像。t和m是有效數字，其中，每10位t組成一個declet，表示一個3位的十進制數。m實際的位置是在第4位到第6位，但是它邏輯上的位置是在那些t前面，所以用()表示放到e的后面。

因為2的10次方是1024，剛好能表示10的3次方。但是表示起來還是需要點技巧的，declet表示三位十進制數的規則比較復雜，這也是這個表示方法叫Densely Packed Decimal(密集十進制數)的原因。下表是編碼的方式。b9-b0代表10個二進制數，d2 d1 d0代表3個十進制數。

就我們的例子來看一下，最低的10位是0000011011，看b3b2b1，這里是101，所以就是上表第3行的情況，三位數字就是 (0000)(1001)(0001)也就是091，然后看從低位數的第3個10位二進制數，也就是00100000000，這顯然是第一種情況，也就是100，連起來就是100000091，指數部分是390，那么這個十進制的值就是 10^(390-398)*100000091 = 1.00000091.

通過這個簡單的例子，就應該對DPD方式的十進制浮點表示方式有個大概的了解了。這個方式算起來比較麻煩，所以除非有硬件支持，軟件模擬的方式都不會使用的，但是DPD轉換成十進制浮點的字符串表示就會很方便。