Page 21 - 《软件学报》2020年第12期

赵世忠 等:循环迭代程序的一种可信计算算法 3687

数字.这时,自变量的误差可能传播到函数值 f(x 0 ),从而使得后者也含有错误有效数字.不妨设 N x 与 N f 分别为 x 0
与 f(x 0 )各自近似值中含有的正确有效数字个数.我们称 N x −N f 为 f(x 0 )(或 f(x)关于 x 0 )的错数.一般情形下,若错数
大于 0,则它代表计算结果中错误有效数字的个数.
比如,设:
4
f(x)=(x−1) ,x 0 =0.9993 (4)
在 53 位的双精度下,x 0 变成了二进制下的 0.11111111110100100001111111110010111001001000111010001,对应
的 十进制数为 0.99929999999999996607158436745521612465381622314453125. 这样 ,N x =16( 理论

上,3= 2.9 = 2.999... .因此,与 3 相比,3.000023 与 2.999923 均具有 5 位正确有效数字.即连续 9 或 0 均为正确数字).
即 x 0 的表示值具有 16 位正确有效数字.对于 f(x 0 )来说,若在 Ubuntu 15.04 的 Gcc 4.9.2 或 Windows 7 的 Visual
Studio 2010 下执行语句 “printf(“%1.28f”,pow((double)0.9993−1,4));”, 则可得相同结果 , 为
4
2.401000000000466e−013.这样,与正确结果 2.401e−013 相比,它有 N f =13 位正确有效数字.因此,(x−1) 在 0.9993
处的错数为
N x −N f =16−13=3 (5)
由于 53 位的双精度运算相当于十进制下的 16 位有效位数下的运算,因此实质上,上述错数 3 表示 16 位(或双精
度下)的结果中包含的错误有效数字的个数.
m
已知实数 x≠0.若使用科学记数法,不妨设 x=±0.t 1 t 2 …×10 (t 1 ≠0,t i 为 0~9 之一),则整数 m 为其指数.比如,
0.9993 的指数为 0,2.401e−013 的指数为−12.
若 f(x)可导,不妨设 x 0 ,f(x 0 )以及 f′(x 0 )的指数分别为 m 1 ,m 2 与 m 0 ,则 f(x 0 )的错数为 m 1 −m 2 +m 0 .其中,误差为 1.
对于式(4)中的函数与自变量来说,不难算得 m 0 =−8.因此,对应错数为 0−(−12)+(−8)=4,它表示函数值的错误
数字比自变量的错误数字多 4 位或 3 位.或更进一步,它表示双精度下,函数值的 16 位计算结果中含有 4 位或 3
位错误数字(实际是 3 位).这也正是式(5)中错数为 3 的根本原因.
对于一个迭代来说,若其具有收敛值(即极限值),则迭代可看作为一个函数,并且存在一个自然数 N,当 n>N
2
时,m 1 =m 2 .这时,错数等于 m 1 −m 2 +m 0 =m 0 .比如,对于迭代式(1)来说,它可被近似看作为函数 111−1130/x+ 3000/x ,
并且当自变量取 6 时,函数值也为 6.因此 m 1 =m 2 =1.不难验证,上述函数在 6 的导数值约为 3.6.因此,错数为 m 0 =1.
这样,16 次迭代就可能将 16 位正确有效位数消耗殆尽,这也正是图 1 中程序的 u 17 不再含有正确有效数字的原
因.另外,容易验证:不论是在 Gcc 4.9.2 下还是 Visual Studio 2010 下,均是从 u 17 开始,不再含有正确有效的数字.
同样,类似的,容易求得式(2)对应函数的错数也为 1.因此,几乎任何系统任何精度下,迭代计算均出错.
由上可知,构造这种计算机计算错误的循环迭代并不是件难事.首先,给出一个函数 f(x),然后令 x=f(x),再进
一步获得满足上式的一个解 x 0 ,最后选择合适的初值(以便使得极限值是 x 0 ).这样,只要 f′(x 0 )的指数大于 0,那么
迭代必定会出错.
比如,下面这两个迭代:
⎧ x = 12.3 ⎧ y = 0.5
⎨ 1 与 ⎨ 1 (6)
−
⎩ x = n 212.3 2460/ x n− 1 ⎩ y = n sin(121*arcsin(y n− 1 ))
前者沿用 Muller 与 Kahan 的思路,是一个更为简单的迭代.容易验证,其每一次迭代的值均等于初值 x 1 = 12.3.然
而容易算得,其对应函数的错数为 2.因此,双精度下,前 7 次迭代结果(从 x 2 到 x 8 )一般含有正确有效数字,而第 16
次迭代的结果(即 x 17 )一定没有正确有效数字.事实上,Gcc 4.9.2 与 Visual Studio 2010 的结果中,分别从 x 14 与 x 15
开始,不再含有正确有效的数字.
对于后者来说,同样容易知道,其每一个 y n 均为 0.5,对应错数为 3.这样,每次迭代的结果中,即使仅有 2 位错
误数字,那么含有正确有效数字的也仅仅是 y 1 到 y 8 .因此,y 9 中一定没有正确有效数字.事实上,无论 Gcc 4.9.2 还
是 Visual Studio 2010,均是从 y 9 开始,结果中不再含有正确有效数字.