浮点算术

转载自 http://www.aogosoft.com/downpage.asp?mode=viewtext&id=186 作者 ham。

作者：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

        ;把两指数相加再减去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)

因此：

其中，

那

因此乘法的结合律的相对误差大约在以内。

程序员在做浮点运算时应该避免测试两个被计算的值是否相等，虽然这极不可能出现。

好，关于浮点算术，这期的内容就写到这里，下期还有更精彩的，待续！