摘  要: 介绍一种在8096/98系列单片机上实现的单精度浮点数快速除法。该算法采用了预估-修正的数值计算方法,并充分利用了16位CPU中的乘除法指令,计算速度快、精度高,有很强的实用性。

　　关键词: 浮点数  除法 尾数  预估-修正  误差  精度

　　在较为复杂的单片机系统中,为扩大取值范围,实现复杂的计算和控制,一般都要涉及浮点数的运算。而一般单片机是没有浮点数运算指令的,必须自行编制相应软件。在进行除法计算时,通常使用的方法是比较除法^[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位，并记Y_H=Y-Y_L。显然Y_H≈Y，X/Y与X/Y_H相差不多:

　　可见只需要在X/YH的基础上再乘以一个修正因子(YH-YL)/YH,就可以得到X/Y的一次校准值。不难证明这个值已经达到了24位的精度要求。事实上,相对误差满足:

　　这说明这个一次校准值完全可以作为最终的结果。

3 算法的具体实现

　　这里的YH虽仍是32位,但其低16位已为0,计算时可以将它视为16位数,这不会影响计算精度。通过两次16位除法,就可得到精确的32位结果。例如,计算Q₀时,第一次除法,X除以YH的高16位,得到的商为Q₀的高16位,而16位余数末尾添0成32位,再除以YH的高16位,得到Q₀的低16位(余数舍去)。由此得到了32位的Q₀。

　　在具体运算中,X应先除以4(X右移2位),以保证Q₀不会溢出(Y_H取高16位)：

　　在计算Q₀′、Q₁时,均进行了两次16位除法,使得Q₀′、Q₁均为精确的32位,保证了计算过程中的精度,减小了累积误差。对于Y_L=0即除数只有16位有效数字的特殊情况,直接有Q₁=1,还能省去两次16位除法。

　　在计算Q时,则通过3次16位乘法实现了32位乘法,取结果的高32位,即得Q。

　　整个算法至多只须用4次除法、3次乘法和5次加法,就求得了浮点数商的尾数,可见计算效率是很高的,保证了运算速度。

　　浮点数除法流程图如图3所示。

4 程序源代码

　　限于篇幅,只给出源代码中的关键部分,即有效数字的计算部分。

；被除数为x，除数为y

；用yh，yl分别表示y的高16位和低16位

…

；假设x，y的有效数字部分分别在(d_x，c_x)和(b_x，a_x)中

；计算预估值Q₀′=(x/4)/y_h

shrl  cx， #2         ；计算x/4

divu cx， bx          ；计算(x/4)÷yh

ld fx， cx            ；把商暂放入寄存器f_x，即Q₀′的高16位有

                      ；效数字

clr   cx

divu cx， bx          ；把余数末尾添0后再除以yh

ld  ex， cx     　　　；把商暂放入寄存器e_x， 即Q₀′

                      ；的低16位有效数字

                      ；(f_x，e_x) = Q₀′

；计算修正因子 Q₁=(yh-y_l)/yh

cmp ax， 0            ；判断yl是否为0

jne   getQ1           ；若yl非0，计算修正因数Q1

ld   ax， ex          ；若yl=0， 修正因数Q1=1

ld   bx， fx          ；(Q₀′×Q₁)=Q₀′，可以直接计算Q

sjmp  getQ

getQ1：

ld  hx， bx           ；把yh放于寄存器hx中

neg  ax

dec  bx               ；计算yh-yl

divu ax， hx          ；计算Q1=(yh-yl)÷yh

ld  dx， ax           ；把商暂时放入寄存器d_x，即Q₁的高16位有

                      ；效数字

clr   ax

divu ax， hx          ；把余数末尾添0后再除以yh,得Q1的

                      ；低16位有效数字

ld  bx dx　           ；(b_x，a_x) = Q₁

；计算Q₀′×Q₁=(f_x，e_x)×(b_x，a_x)，只取32位有效数字

ld     hx， bx

mulu   cx， bx， ex   ；(d_x，c_x) = b_x×e_x

mulu   ax， fx        ；(b_x，a_x) = a_x×f_x

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 = (Q₀′×Q₁)×4并调整阶码

getQ：

　　…　

　　代码到这里为止,浮点数商的有效数字已经全部求出。只要再执行一些调整浮点数阶码的操作,就可以得到最终结果。

　　在作者开发的一个80C196KC单片机系统中,涉及到了二进制-十进制数制转换、分段线性插值、数字滤波等大量浮点数的运算,都是靠加减乘除等底层函数来实现的。

　　此外,本算法思路清晰,因此很容易加以推广。例如,为了得到更高的精度，可取修正因子:

参考文献

1 复旦大学计算机系微机开发研究室.十六位单片机8096的原理和设计方法.重庆：科学技术文献出版社重

  庆分社，1988

2 涂时亮，姚志石.单片微机MCS-96/98实用子程序.上海：复旦大学出版社， 1991

3 李庆扬.数值分析. 武汉：华中工学院出版社， 1986