浮点算术.txt

会变语言实现的一些程序

浮点算术  作者：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