在较为复杂的单片机系统中,为扩大取值范围,实现复杂的计算和控制,一般都要涉及浮点数的运算。而一般单片机是没有浮点数运算指令的,必须自行编制相应软件。在进行除法计算时,通常使用的方法是比较除法[1],即利用循环移位和减法操作来得到24~32位商,效率很低。有些文献给出了一些改进方法[2],但思想不清晰,很难推广使用。这里给出一种浮点数除法运算的实用快速算法。该方法以数值计算中的预估-修正方法为指导,充分利用了16位单片机的乘除法功能,很轻易地实现了浮点数的除法。
1 浮点数格式
IEEE的浮点数标准规定了单精度(4字节)、双精度(8字节)和扩展精度(10字节)三种浮点数的格式。最常用的是单精度浮点数,格式如图1所示。但是这种格式的阶码不在同一个字节单元内,不易寻址,从而会影响运算速度。
通常在单片机上采用的是一种变形格式的浮点数,如图2所示。其中的23位尾数加上隐含的最高位1,构成一个定点原码小数,即尾数为小于1大于等于0.5的小数。有关浮点数格式的详细内容请参考有关文献[1][2]。
2 快速除法的算法原理
在16位单片机中只有16位的乘除法,而浮点数的精度(即尾数的有效位数)达24位,因此无法直接相除,但依然可以利用16位的乘除法指令来实现24位除法。不过,如果只进行一次16位的除法必定会带来很大误差,因此问题的关键在于如何消除这个误差,从而达到要求的精度。这其实就是通常数值计算中所采用的预估-修正方法。
假设两个浮点数经过预处理后,被除数和除数尾数扩展为32位(末8位为0)分别放入X和Y中。邻YL为Y的低16位,并记YH=Y-YL。显然YH≈Y,X/Y与Y/YH相差不多:
(X/Y)/(X/YH)=(YH/Y)
=YH/(YH+YL)
=1/(1+YL/YH)
≈1-YL/YH
=(YH-YL)/YL
可见只需要在X/YH的基础上再乘以一个修正因子(YH-YL)/YH,就可以得到X/Y的一次校准值。不难证明这个值已经达到了24位的精度要求。事实上,相对误差满足:
这说明这个一次校准值完全可以作为最终的结果。
3 算法的具体实现
在具体实现本算法时,主要经过下列步骤:
(1)计算预估值Q0=X/YH;
(2)计算修正因子Q1=(YH-YL)/YH;
(3)计算校准值Q=Q0×Q1,并作为最后结果。
这里的YH虽仍是32位,但其低16位已为0,计算时可以将它视为16位数,这不会影响计算精度。通过两次16位除法,就可得到精确的32位结果。例如,计算Q0时,第一次除法,X除以YH的高16位,得到的商为Q0的高16位,而16位余数末尾添0成32位,再除以YH的高16位,得到Q0的低16位(余数舍去)。由此得到了32位的Q0。
在具体运算中,X应选除以4(X左移2位),以保证Q0不会溢出(YH取高16位):
Q0'=(Y/4)/YH≤(0ffffff00H/4)/8000H=7fffh
由于X为32位(末8位为0),这一操作不影响有效数字。而Q1(YH-YL)/YH≤1,不存在溢出的问题。最后计算校准值Q时,有Q=4Q0'Q1。
在计算Q0'、Q1时,均进行了两次16位除法,使得Q0'、Q1均为精确的32位,保证了计算过程中的精度,减小了累积误差。对于YL=0即除数只有16位有效数字的特殊情况,直接有Q1=1,还能省去两次16位除法。
在计算Q时,则通过3次16位乘法实现了32位乘法,取结果的高32位,即得Q。
整个算法至多只须用4次除法、3次乘法和5次加法,就求得了浮点数商的尾数,可见计算效率是很高的,保证了运算速度。
浮点数除法流程图如图3所示。
4 程序源代码
限于篇幅,只给出源代码中的关键部分,即有效数字的计算部分。
;被除数为x,除数为y
;用yh,hl分别表示y的高16位和低16位
...
;假设x,y的有效数字部分分别在(dx,cx)和(bx,ax)中
;计算预估值Q0'=(x/4)/yh
shrl cx,#2 ;计算x/4
pu cx,bx ;计算(x/4)&pide;yh
ld fx, cx ;把商暂放入寄存器fx,即Q0'的高16位有 ;效数字
clr cx
pu cx,bx ;把余数末尾添0后面再除以yh
ld ex,cx ;把商暂放入寄存器ex,即Q0'
;的低16位有效数字
;(fx,ex)=Q0'
;计算修正因子 Q1=(yh-yl)/yh
cmp ax,0 ;判断yl是否为0
jne getQ1 ;若yl非0,计算修正因数Q1
ld ax,ex ;若yl=0,修正因数Q1=1
ld bx,fx ;(Q0'×Q1)=Q0',可以直接计算Q
sjmp getQ
getQ1:
ld hx,bx ;把yh放于寄存器hx中
neg ax
dec bx ;计算yh-yl
pu ax,hx ;计算Q1=(yh-yl)&pide;yh
ld dx,ax ;把商暂时放入寄存器dx,即Q1的高16位有;效数字
clr ax
pu ax,hx ;把余数末尾添0后再除以yh,得Q1的;低16位有效数字
ld bx dx ;(bx,ax)=Q1
;计算Q0'×Q1=(fx,ex)×(bx,ax),只取32位有效数字
ld hx,bx
mulu cx,bx,ex ;(dx,cx)=bx×ex
mulu ax,fx ;(bx,ax)=ax×fx
clr ex
add cx,ax
addc dx,bx
addc ex,0 ;(ex,dx,cx)=(dx,cx)+(bx,ax)
mulu ax,fx,hx ;(bx,ax)=fx×hx
add ax,dx ;(bx,ax)=(bx,ax)+(ex,dx)
addc bx,ex ;(bx,ax)=Q0'×Q1
;计算校准值Q=(Q0' ×Q1)×4并调整阶码
getQ:
...
代码到这里为止,浮点数商的有效数字已经全部求出。只要再执行一些调整浮点数阶码的操作,就可以得到最终结果。
在作者开发的一个80C196KC单片机系统中,涉及到了二进制-十进制数制转换、分段线性插值、数字滤波等大量浮点数的运算,都是靠加减乘除等底层函数来实现的。