Algorithm 版 (精华区)

发信人: sino (茶水博士), 信区: Theory
标  题: [转载] 新型的快速高准确度开方算法及程序设计
发信站: 哈工大紫丁香 (Sun Aug 27 14:53:19 2000), 转信

发信人: jpeg (雨季), 信区: algorithm
发信站: NJU Lily BBS (Tue Jul 13 19:37:31 1999), 转信

【 以下文字转载自 EEdepartment 讨论区 】
【 原文由 jpeg 所发表 】
 电子技术应用
APPLICATION OF ELECTRONIC TECHNIQUE
1999年 第3期 No.3 1999
  

--------------------------------------------------------------------------------
 
新型的快速高准确度开方算法及程序设计

李芙英　王恒福

　　摘　要：介绍一种新型的快速高准确度开方算法，特别适用于需要用计算机进行型式
开方运算场合。算法巧妙地将开方变量由两个减少为一个，将变量变化区间由整个实数轴
缩小为［0,1］区间，进而采用查表与插值相结合的方法，实现了高准确度、快速开方运
算。在单片机80c196kb
上，利用PL/M96语言编程进行了运算，效果良好。
　　关键词：开方算法　整型数运算　程序设计

　　在测量系统中，经常会遇到计算有效值问题，而计算有效值不可避免地要在计算机上
进行开方运算。如在电子系统中，对畸变的电压、电流波形用付立叶算法进行谐波分析时
，计算各次谐波分量幅值公式为Ak=，其中ak、bk为k次谐波的实部和虚部。对这种计算既
要求快速实时，又要求
准确度高。本文介绍了一种适用于计算机上进行的快速、高准确度的形式的开方算法，具
有较高的使用价值。

1　算法描述
1.1　问题的转化
　　在计算机上进行a2+b2形式的开方运算，通常采用的是迭代法，但迭代法运算速度慢
，不能满足工程上实时快速运算的要求；为了提高运算速度，本算法首先将a2+b2形式加
以变换，使变量区间由整个实轴缩小为一个有限的区间，然后采用查表和插值的方法实现
快速高准确的开方运算。
　　为了叙述方便起见，设被开方数y为a、b的平方和开方，且｜a｜>｜b｜，则
　　

　　其中是很容易求得的，因此只要得到对应于x的值，则值也就知道了。由此可见，计
算机的运算就转化成计算在区间［0,1］上某一点的值了，而且将变量由原来的两个a、b
减少成为一个x，把变量变化区间由整个实轴转化为区间［0,1］。这样的转变为快速高准
确度的开方算法(查表法)
提供了一种可能。
1.2　用插值法实现开方运算
　　在计算机里提前制作一张根据x就可以查得值的表，这样已知x就能得到1+x2的值，但
问题在于x变化区间虽然缩小到［0,1］，为了确保计算的准确度，这张表将是非常庞大的
。
　　例如：假设将区间［0,1］等分N段，步长为，则查表的相对误差为：
　　

(2)

　　当x=1时，，取得最大值；
　　若要保证万分之一的准确度，则，N≥5000。即这张表至少要存入5000个数，才能保
证在区间［0,1］查表得到的的值达万分之一准确度，而在单片机测量系统中，这张表要
占用这么大的存储区是不合适的［1］。由于开方函数是一条以x为自变量的曲线，因此我
们采用折线段逼近曲线方
法来实现。这样的开方运算在同样准确要求下，在计算机中存入的数字量将会大大地减少。
　　设：逼近所需折线段数为m，即将区间［0,1］等分为m段，将区间等分点对应曲线点
依次连成直线段逼近曲线。
　　此种方法带来相对误差：
　　

(3)

　　当x=0时，取得最大值
　　若要达到万分之一的准确度，则，N≥35.5，则N可取36，因此只需在计算机中存入7
2个数(每点存二个数：值及斜率值)即可。
1.3　算法在整型数运算的实现
　　以上讨论了计算的数学方法，那么在计算机上实施时又怎样能得到高准确度的运算结
果。在这里数型选取直接会影响到计算准确度。由于考虑到速度的要求，不宜采用浮点数
运算，根据硬件和软件环境，在这里我们选用整型数运算，为此还需研究由此带来的一些
特殊问题。
　　运算首先要解决用整型数表示问题。
　　例：｜b｜=｜a｜，则x=1，其他情况x=0，这种表示法显然是不行的，必须采用先将
｜b｜乘以2m倍(即在双字型变量中左移若干位，再除以｜a｜)，则得到n=2mx，即x=，在
程序中n也是用一个整型数表示，所以取m=16。
　　在程序中是用128个折线段来逼近曲线的，根据n的二进制数值的前7位来选择基点处
｜x=xi的值以及斜率值，然后再由实际的x和基点处xi的差来计算出x处1+x2的值。

2　开方程序的编制
　　开方程序在80c196kb上，利用PL/M96语言编程进行了运算。具体程序如下［2］：
SQRT：DO；
DECLARE SQT(128) WORD DATA (… …)；
/* SQT(i)=(1+((i)/(128))2-1)*216的取整值*/
DECLARE LOP(128) WORD DATA (… …)；
/* LOP(i)=SQT(i+1)-SQT(i)*/
DECLARE (a,b) WORD；
DECLARE (BIGGER $ ONE，SMALLER $ ONE)WORD；
DECLARE NUMBER WORD；
DECLARE BASE $ POINT WORD；
DELCARE BASE $ NUMBER WORD；
DELCARE MUL $ NUMBER WORD；
DELCARE SQRTION WORD；
IF　a＞=32768 THEN　　　/* 求｜a｜和 ｜b｜*/
a=NOT a+1；
IFb>=32768 THEN
b=NOT b+1；
IF a=b THEN
SORTION=a+HIGH(DOUBLE(a)*27146)；
ELSE DO；
　IF a>b THEN
　　DO；
　　　BIGGER $ ONE=a；
　　　SMALLER $ ONE=b；
　　END；
　ELSE
　　DO；
　　　BIGGER $ ONE=b；
　　　SMALLER $ ONE=a；
　　END；
/*该部分作用是将max(｜a｜,｜b｜)存入BIGGER $ ONE中 ，min(｜a｜,｜b｜)存入SMAL
LER $ ONE中；*/
NUMBER=LOW (SHL (DOUBLE (SMALLER $ ONE)，16) /BIGGER$ONE)；
/*NUMBER=216×SMALLER $ ONE/BIGGER $ ONE*/
BASE $ POINT=SHR (NUMBER,9);
/*取NUMBER高7位值，便得到了基点*/
BASE $ NUMBER=SHR (NUMBER AND 01FFH,2);
/*该值表示x和xi间的距离，之所以右移2，是为了防止后面运算溢出*/
MUL $ NUMER = SQT (BASE $ POINT) + SHR (BASE $ NUMBER*LOP (BASE*POINT),7);
/*因为LOP(i)=SQT(i+1)-SQT(i)，而xi+1和xi之间分成了27份，故要现除以27，即右移7
位，当然运算中要先乘后除*/
SQRTION=BIGGER $ ONE + HIGH (DOUBLE (BIGER $ ONE) *MUL $ NUMBER);
/*之所以把BIGGER $ ONE扩展为双字，是因为在PL/M-96中，字乘字结果仍为字，高位信
息被丢失，故先扩展为双字，乘完后再截取高位即可*/
RETURN SQRTION；
END；/*ELSE块结束*/
END SQRT；(过程结束)

3　测试结果
　　 测试结果表明该程序误差极小，例如a=10000，b=2000时，a2+b2=10198.04，用整数
表示则为10198，而用该程序得到的结果也是10198。而它的运算速度却非常快，主要运算
为两次乘一次除，完成一次开方运算只需要150个时钟周期，因此是一个非常适合于数字
测量系统的开方程序。
　　该算法应用于光电式VIP测量系统中，取得了良好的效果。

作者单位：北京清华大学(100084)

参考文献

　[1]孙涵芳.Intel 16位单片机.北京：北京航空航天大学出版社，1998.3
　[2]袁涛，孙腾谌.PL/M-96程序设计语言及其应用.北京：清华大学出版社，1990.8

收稿日期：1998-11-05
 

--
※ 来源:．NJU Lily BBS dii.nju.edu.cn．[FROM: e248jxy.nju.edu]

--
※ 修改:．fib 於 Aug 27 13:27:32 修改本文．[FROM: bbs.hit.edu.cn]
--
※ 转寄:．南京大学小百合 bbs.nju.edu.cn．[FROM: bbs.hit.edu.cn]

--
☆ 来源:．哈工大紫丁香 bbs.hit.edu.cn．[FROM: fib.bbs@bbs.nju.edu.]