计算机浮点数运算的舍入误差分析

计算机浮点数运算的舍入误差分析

计算机中的字数应表示为浮点形式。由于计算机浮点数中尾数的位数是有限的,在进行科学计算的时候,对一般实数总是要按舍入原则表示为浮点数。这样,计算机浮点数就是该实数的近似值。因此,我们有必要对计算机浮点数和浮点数的算术运算的舍入误差进行系统的理论分析。事实上,这正是在科学计算中对各种数值计算方法作误差分析的基础。由算术运算的误差界限,可以建立更复杂运算的误差界限,因而可对各种算法的舍入误差进行分析。1计算机浮点系统的构成设F表示系统的计算机浮点数全体所构成的集合,即F={0}∪{f∶f=±βJt∑k=1dkβ-k,0≤dk≤β,d1≠0,L≤j≤U}(1)其中每个dk都是整数。β是计算机浮点数的基数,当β=10时,就是通常使用的十进制数;当β=2,8,16时,分别表示二进制数、八进制数和十六进制数。本文就一般基数情形进行统一讨论。又±t∑k=1dkβ-k=±0.d1d2⋯dt是尾数,自然数t是字长,整数j是阶,阶有固定的下限L和上限U,L和U均为正整数或零;因此,F是由四元数组(β,t,L,U)来描述的,即F=F(β,t,L,U),称为计算机浮点数系统,显然它是包含2(βt-βt-1)(U-L+1)+1个数的有限集。F中的数对称地分布在区间[-M,m]和[m,M]中,其中{Μ=βU(β-1)(t∑k=1β-k)=βU(1-β-t)m=β-1βL=βL-1(2)以下将计算机浮点数简称为浮点数,并把实数x的浮点表示为fl(x),若x=0时,fl(0)=0;当m≤|x|≤M时,对舍入法有fl(x)={0.d1d2⋯dt×βJ,当0≤dt+1<β/2时(0.d1d2⋯dt+β-t)×βJ,当dt+1≥β/2时(3)对截断法,有fl(x)=0.d1d2…dt×βJ。例如对二进制系统的舍入法有fl(x)={0.d1d2⋯dt×2J,当dt+1=0(0.d1d2⋯dt+2-t)×2J,当dt+1=1(4)截断法有fl(x)=0.d1d2…dt×2J。两个浮点数相加或相减,首先要对阶,即将两数的小数点对齐,使它的阶相等。具体做法是把阶小的浮点数的尾数右移,每移动一位其阶就加1,直到两浮点数的阶相等。然后相加或相减。下面,我们将用“。”来表示加、减、乘、除四则运算中任意一种运算。在运算中,若|a。b|>M是上溢;若0<|a。b|<m是下溢。以下我们将在不发生溢出的情况下进行讨论。2浮点数的间距实数表示成浮点数之后,我们可得到它的相对误差限,这是进行误差分析的基础。误差分析的基本定理是:设m<|x|<M,其中m和M由(2)式给出,则有fl(x)=x(1+δ),|δ|≤u(5)其中u为浮点数的相对误差限(机器精度),即u={12β1-t,舍入法β1-t,截入法(6)证明若fl(x)=0,(5)式显然成立。现不妨假定x<0,同理可记x>0的情形。设x∈(-βJ,-βJ-1),故x的阶为J,于是在该区间中的所有t位浮点数以间距βJ-t分部。对舍入法,有|fl(x)-x|≤12βJ-t,12βJ-t是fl(x)的绝对误差限,又|fl(x)-x|≤12βJ-t=12βJ-1β1-t≤12β1-t从而|fl(x)-xx|≤12β1-t(7)这样我们就得到对浮点数使用舍入法的相对误差限。令δ=fl(x)-xx,由(7)知|δ|≤u=12β1-t,从而有fl(x)=x(1+δ),|δ|≤u=12β1-t(8)对截断法,有|fl(x)-x|≤βJ-t,βJ-t为使用截断法的绝对误差限,因为|fl(x)-x|≤βJ-t=βJ-1β1-t≤|x|β1-t从而有|fl(x)-xx|≤u=β1-t(9)令δ=fl(x)-xx,得|fl(x)-x|=x(1+δ),|δ|≤u=β1-t(10)联合(8)和(10)式,即得(5)式和(6)式,证毕。特别对10进制系统,(6)式为u={5×10-t,舍入法10×10-t,截入法对二进制系统为u={2-t,舍入法2×2-t,截入法在某些情形,为使用便利,我们也将fl(x)的相对误差表为fl(x)=x1+δ,|δ|≤u(11)3约化法证明由浮点数表示的特征可以导出算术运算的舍入误差。设x,y∈F,fl(x。y)的意义是先让x,y进行运算,运算的结果是精确的实数x。y,再按舍入规则表示成F中的浮点数。x。y与fl(x。y)的误差就是舍入误差。关于浮点数算术运算的舍入误差,据(5)和(6)式得如下结论:fl(x˚y)=(x˚y)(1+δ),|δ|≤u(12)其中u由(6)式确定。利用这些算术运算的相对误差界,可以得到更复杂运算的相对误差界。例如s1=fl(x1×y1)=x1×y1(1+γ1)=x1×y1(1+γ1)(1+δ1),|γ1|≤u(13)式中约定δ1=0。s2=fl(x1×y1+x2×y2)=fl(fl(x1×y1)+fl(x2×y2))=fl(s1+fl(x2×y2))=fl(x1×y2)(1+γ1)+(x2×y2)(1+γ2)=[(x1×y1)(1+γ1)+(x2×y2)(1+γ2)]×(1+δ2)=(x1×y1)(1+γ1)(1+δ1)×(1+δ2)+(x2×y2)(1+γ2)(1+δ2)=2∑i=1xi×yi(1+γi)2∏j=i(1+δj),|γi|,|δi|≤u,δ1=0(14)一般地,有sn=fl(n∑i=1xi×yi)=n∑i=1xi×yi(1+γi)n∏j=i(1+δj)|γi|,|δj|≤u,δ1=0(15)现用数学归纳法证明之。当n=1时,据(13)式,上式显然成立。设sk-1=flk-1∑i=1xi×yi=k-1∑i=1xi×yi(1+γi)k-1∏j=i(1+δj)成立,则sk=fl(sk-1+fl(xk×yk))=[sk-1+xk×yk(1+γk)](1+δk)=[k-1∑i=1xi×yi(1+γi)k-1∏j=i(1+δj)+xk×yk(1+γk)]×(1+δk)=k-1∑i=1xi×yi(1+γi)k∏j=i(1+δj)+xk×yk(1+γk)(1+δk)=k∑i=1xi×yi(1+γi)k∏j=i(1+δj)|γk|,|δk|≤u;δ1=0所以对k成立,根据数学归纳法,对任何自然数n,(15)式均成立。后面还将继续对这个问题进行进一步分析。4数学上的加减运算常用的误差分析方法主要有两种,即向前误差分析方法和向后误差分析方法。人们注意到上面误差分析的特点是,将浮点数含舍入误差的浮点运算,表示为近似数的精确数学运算。例如fl(x±y)=(x±y)(1±δ)=x(1+δ)±y(1±δ)浮点数的x,y的加(减)运算,可以看做是近似数x(1+δ)和y(1+δ)精确的数学上的加(减)运算。这种把计算过程产生的误差返回到原始资料的误差的误差分析方法称为向后误差分析方法。直接估计计算解与精确解之间的误差,这种误差分析方法称为向前误差分析方法或直接误差分析法。例如对加法运算就是要估计|fl(x+y)-(x+y)|的界限。向后误差分析将浮点数的实际计算转化为实数的精确运算,所以在误差分析过程中就可毫无困难地使用实数的代数运算法则。向前误差分析就没有这个优点;向前误差分析另一缺点是所得到的误差界常常十分保守,且误差分析方法非常麻烦。此外,在实际问题中,原始数据本身常带有误差,所以使用向后误差分析能将计算过程的误差与原始数据的误差比较而确定计算过程是否稳定。总之,向后误差分析方法是一种十分有效的、功能齐全的误差分析方法,它广泛地应用在线性方程组求解的误差分析中。事实上,比如用Gauss列主元消去法求解线性方程组时,由于引进了舍入误差,其计算解可归结为系数矩阵作扰动而得到的扰动方程组的精确解。5e0.0.11.1x19式在误差估计式中常常会遇到∏(1+δk)形式的乘积,现在我们对它进行估计。结论是:若|δi|≤u,且nu≤0.01,则1-nu≤n∏i=1(1+δi)≤1+1.01nu|δi|≤u(16)证明:因|δi|≤u,所以-u≤δi≤u,从而有1-u≤1+δi≤1+u,得(1-u)n≤n∏i=1(1+δi)≤(1+u)n(17)由于0<u<1,成立不等式1-nu≤(1-u)n,立即得到左端(1-u)n的下界1-1.01nu≤1-nu≤(1-u)n(18)另方面,由ex的幂级数展开式,当x>0时,1+x≤ex≤1+x+x2xex,若0≤x≤0.01,则1+x≤ex≤1+x+0.012xe0.01≤1+1.01x(19)上式利用了e0.01=1.01005<2的结果。令x=u,由(19)式左端知(1+u)n≤enu(20)再令x=nu,据假设nu≤0.01,由不等式(19)右端得enu≤1+1.01nu(21)联合(20)和(21)式,得(17)式右端上界(1+u)n≤1+1.01nu(22)这就证明了,当nu≤0.01时,乘积n∏i=1(1+δi)的估计式(16)。由(18)式,1-1.01nu≤1-nu,(16)式可以写为1-1.01nu≤n∏i=1(1+δi)≤1+0.1nu(23)从而得-1.01nu≤n∏i=1(1+δi)-1≤1.01nu(24)令n∏i=1(1+δi)-1=δ,由上式得|δ|≤1.01nu,故(24)式,从而(16)式又可写成n∏i=1(1+δi)=1+δ,|δ|≤1.01nu(25)或n∏i=1(1+δi)=1+1.01nθu,|θ|≤1有了重要的不等式(16)等,再继续对(15)式进行分析。令1+εi=(1+γi)n∏j=i(1+δi),(δ1=0)由(23)式有1-1.01nu≤1+εi≤1+1.01nu,从而-1.01nu≤εi≤1.01nu,即|εi|≤1.01nu(26)因sn=fl(n∑i=1xi×yi)=n∑i=1xi×yi(1+γi)×n∏j=i(1+δj)=n∑i=1(1+εi)xi×yi由(26)式得|fl(n∑i=1xi×yi)-n∑i=1xi×yi|≤n∑i=1|εi|⋅|xi×yi|≤1.01nun∑i=1|xi×yi|(27)于是fln∑i=1xi×yi的相对误差限为|fl(∑i=1nxi×yi)-∑i=1nxi×yi∑i=1nxi×yi|≤1.01nu∑i=1n|xi×yi|∑i=1nxi×yi(28)根据此式,若|∑i=1nxi×yi|≪∑i=1n|xi×yi|,则fl∑i=1n|xi×yi|的相对误差限可能会很大,因此为提高精度,人们先用双精度计算,然后把计算结果舍入为单精度。6浮点数n的性质边形定义利用以上得到的各结论,我们不难获得矩阵运算的舍入误差。定义|E|=[|eij|],(E∈Rn×n)设|G|=|gij|,且|eij|≤|gij|,(i,j=1,2,⋯,n),则规定|E|≤|G|,(i,j=1,2,⋯,n)如A,B是由F中的浮点数构成的n×n矩阵,又浮点数α∈F,则由(12)式易得fl(αA)=αA+E,|E|≤u|αA|fl(A+B)=(A+B)+E,|E|≤u|A+B|其中E为误差矩阵,又据(27)式,得fl(AB)=AB+E,|E|≤1.01nu|A|⋅|B|由(28)式可知,若|AB|≪|A|·|B|,fl(AB)的相对误差可能较