中值定理法(开平方算法)

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

作者:李义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来讲,查表环节也可省去,程序代码大幅减少,计算流程如下:

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;