• 回复
  • 收藏
  • 点赞
  • 分享
  • 发新帖

最快的开平方算法

最快的开平方算法(中值定理法)



作者:李义 2006





关键词:最快 开平方根 算法 中值定理 开方




整数平方数中值定理:
设a、b、c为顺序排列间距为P的3个整数,A、B、C是它们的平方
则有:b2=(a 2+c2)/2-R,即:B=(A+C)/2-R
其中:修正值R=P2
特别地,如果间隔P=1、2、 4、 8、 16、…2 n (或Pn=2Pn-1)时
则:    修正值R=1、4、16、64、256、…22n (或Rn=4Rn-1)

证明:
已知:a=b+P
c=b-P
有:a 2=(b+P)2=b2+2Pb+P2
c2=(b-P)2=b2-2Pb+P2
则:a2+c2=2b2+2P2
即:b2=(a2+c2)/2-P2
特别地
当:间隔   P=2 n=2*2 n -1=2 Pn-1时(n为自然数)
则:修正值 R=P2=22n=(2 Pn-1)2=4(P n-1)2=4Rn-1
(证明完)


根据以上定理,可以实现整数快速开平方根计算:

先构建一个长度为N的数组1:

数组长 N=Ni+1 1 2 3 4 5 …
间隔 P=2Pi 2 4 8 16 32 …
修正值 R=4Ri 4 16 64 256 1024 …

以及一个对应2PN(这里N=4、2PN=32)的典型数和它的平方数组2:
按N=4间隔
排列的数 d=di+2PN 32 64 96 128 160 192 224 256 …
该数的平方 D=d2 1024 4096 9216 16384 25600 36864 50176 65536 …
显然,N值越大则数组2越小、程序代码效率越高、用时(插值次数)越多.

以2字节整数开方为例的计算流程如下:
其中,被开方数D(范围0~65536),其平方根d(范围0~256)


注:1、查表可以从任何位置开始,对计算速度影响不大.其中D=0、D=1、D=Di、 D>65280判断可以省去.

2、此算法完全没有乘法试算,其1/2、1/4除法运算可由二进制移位简单实现,且为完全补偿后的精确插值,所以递归速度非常快(这里4次).

3、最后运算已经包括了小数部分的精确4舍5入算法.

4、此算法略加改动,即可实现更长字节整数或定长浮点数平方根精确解,其逆运算也可以实现乘方运算.

一个C语言实例:

//  sqrt_2  中值定理法开平方程序(直接查表-插值)
//  输入D (两字节无符号整数)
//  输出d (一字节无符号整数)

char a,b,c,p;
int A,B,C,D,R,K;
  void main()
  {D=11111;    // 被开方数
  if(D>50176){A=0;   a=0;  C=50176;c=224;break;}; // 查表
  if(D>36864){A=50176;a=224;C=36864;c=192;break;};      
  if(D>25600){A=36864;a=192;C=25600;c=160;break;};
  if(D>16384){A=25600;a=160;C=16384;c=128;break;};
  if(D>9216) {A=16384; a=128;C= 9216; c= 96; break;};
  if(D>4096) {A= 9216; a= 96; C= 4096; c= 64; break;};
  if(D>1024) {A= 4096; a= 64; C= 1024; c= 32; break;};
            A= 1024; a= 32; C= 0;    c= 0;  break;
p=16;R=256;                // 初始化数据
  do{ b=c+p;B=C;B>>=1;       // 插值计算循环
     if(A!=0){K=A;K>>=1;}
        else K=0x8000;        // 65536>>=1的数
     B+=K;B-=R;
     if(D>B){C=B;c=b;}
        else{A=B;a=b;}
     p>>=1;R>>=2;
     }while(p!=1);             // 循环4次结束
   K=A-C;K>>=2;A-=K; C+=K;  // 小数部分四舍五入
   if(D      else{
      if(D         else b=a;}
}                        //开方结束

进一步研究表明,由于循环内所有运算都是加、减、位移、比较等简单运算,所花费的时间很少,可以适当加大循环次数.

特别地,如果把间隔P加大到128,对应修正值R=13684,则循环次数N=7,对应数组2就简化到:

按N=7间隔排列的数 d=di+Pn 0 256 512 …
该数的平方         Di=d*d 0 65536 262144 …

这时,对于两字节数被开方数D来讲,查表环节也可省去,程序代码大幅减少,计算流程如下:

500) {this.resized=true; this.width=500; this.alt='这是一张缩略图,点击可放大。\n按住CTRL,滚动鼠标滚轮可自由缩放';this.style.cursor='hand'}" onclick="if(!this.resized) {return true;} else {window.open('http://u.dianyuan.com/bbs/u/50/1172121266.jpg?x-oss-process=image/watermark,g_center,image_YXJ0aWNsZS9wdWJsaWMvd2F0ZXJtYXJrLnBuZz94LW9zcy1wcm9jZXNzPWltYWdlL3Jlc2l6ZSxQXzQwCg,t_20');}" onmousewheel="return imgzoom(this);">



C语言程序的一个例子:
unsigned char a,b,p=0x80;
unsigned int K,A,B,C,R=0x4000,D=60000;


sqt1(){
do{
   b=a-p;B=C;B>>=1;
   if(A){K=A;K>>=1;}
     else K=0x8000;
   B+=K;B-=R;
   if(D>B)C=B;
     else{A=B;a=b;}
   p>>=1;R>>=2;
  }while(p!=1);     //循环7次结束
   p=(A-C)>>2;A-=p; C+=p;
   b=a;
   if(D < A)b--;
   if(D < C)b--;
   }         //输出方根b

程序里只用了一个特别的数128(及其它的平方数16384),就能够把两字节数0~65535范围内的任意整数的整数平方根精确(小数部分严格4舍5入)求解.

程序思想还可以继续延伸到更长字节无符号整数的开方,只需要修改对应的初始值p、R就行了.



结论 :

本文首先提出并证明了整数平方数中值定理,进而提出一种基于此定理的的快速开方算法,并给出了具体的计算流程和C语言程序实例.由于全部运算不使用乘法运算或幂运算,只使用加、减、移位、逻辑等简单运算,只引入极少的初始变量,在经过有限次循环后即可迅速逼近任意有限大整数的整数平方根的精确解(小数部分严格4舍5入),从而把整数开平方运算的迄今最快速度提高了一个数量级.

此算法对于由没有集成硬件乘法器的芯片组成的微处理(控制)系统、同时要担负大量开方运算而时间特别紧张的应用——如大型游戏程序、图形处理系统中两坐标点的距离计算(方根)、交流电压有效值等控制计算(均方根)——具有一针见血的效果,对于任何高级程序语言、游戏机乃至计算器、民用或工业控制系统中的开平方函数代码的优化具有显著积极的意义.
全部回复(8)
正序查看
倒序查看
xlylyl
LV.1
2
2007-02-22 15:25
我是刚学单片机需要各位给予帮助:
    我手上有一个仿真板,P1口上接了八个LED,P3.2,P3.3对地各接了一个轻触开关,
    1.我想按一下P3.2上开关,P1口上的八个灯就向上流动,按一下就亮一个灯.
    2.按一下P3.3上开关,P1口上的八个灯就向下流动,按一下就亮一个灯.
0
回复
sdjufeng
LV.6
3
2007-02-23 23:15
C51自带的运算程序远不是最佳的,如果想快速的开平方,我建议采用外挂asm文件,用汇编语言来实现,不仅如此,在我的程序中,就是双字节乘除法我都是利用自己编写的汇编子程序,而不用系统自带的,速度提高数倍.
0
回复
ldfa
LV.4
4
2007-02-24 22:11
查表是最快的了
0
回复
da2007
LV.2
5
2007-02-28 20:45
李哥:会欣赏的人还是少数的!!
呵呵~~~!
0
回复
nc965
LV.6
6
2007-03-02 00:30
@da2007
李哥:会欣赏的人还是少数的!!呵呵~~~!
你把它发到自动化坛子或者单片机坛子里去看看
0
回复
luck5263
LV.1
7
2011-12-01 16:33

能给出32位无符号开平方算法在一般的51单片机上运行时间的数据吗?我在72MHZ的ARM7(16位模式)上运行竟然要20多微秒

0
回复
zq2007
LV.11
8
2011-12-01 16:58
@luck5263
能给出32位无符号开平方算法在一般的51单片机上运行时间的数据吗?我在72MHZ的ARM7(16位模式)上运行竟然要20多微秒
不会C语言怎么办?有没有更贴合实际的方法?
0
回复
cheng111
LV.11
9
2011-12-01 21:49
@zq2007
不会C语言怎么办?有没有更贴合实际的方法?
有,以前书上是有估算公式的,大学的高等数学书上。
0
回复