浮点算术.txt

会变语言实现的一些程序

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