开平方数的快速算法
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
开平方数的快速算法.txt我爸说过的最让我感动的一句话:“孩子,好好学习吧,爸以前玩麻将都玩儿10块的,现在为了供你念书,改玩儿1块的了。”整数平方数中值定理:
设a、b、c为顺序排列间距为P的3个整数,A、B、C是它们的平方
则有:b^2=(a^2+c^2)/2-R,即:B=(A+C)/2-R
其中:修正值R=P^2
特别地,如果间隔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=b^2+2Pb+P^2
c2=(b-P)^2=b^2-2Pb+P^2
则:a^2+c^2=2b^2+2P^2
即:b^2=(a^2+c^2)/2-P^2
特别地
当:间隔 P=2 n=2*2 n -1=2 Pn-1时(n为自然数)
则:修正值 R=P^2=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
if(Delse 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;
if(Dif(D
程序里只用了一个特别的数128(及其它的平方数16384),就能够把两字节数0~65535范围内的任意整数的整数平方根精确(小数部分严格4舍5入)求解。
程序思想还可以继续延伸到更长字节无符号整数的开方,只需要修改对应的初始值p、R就行了。
结论 :
本文首先提出并证明了整数平方数中值定理,进而提出一种基于此定理的的快速开方算法,并给出了具体的计算流程和C语言程序实例。由于全部运算不使用乘法运算或幂运算,只使用加、减、移位、逻辑等简单运算,只引入极少的初始变量,在经过有限次循环后即可迅速逼近任意有限大整数的整数平方根的精确解(小数部分严格4舍5入。
以下是大明轮王点评:
说自己最快是要证明的:) 否定他只要举出一个更快的,更简单的。中国古代就有地道的开方术,原理是一位一位的尝试,写出来就是程序sqrt_1,很容易理解,但是用了乘法。通过二项式定理进行一系列优化,可以得到程序sqrt_5,只用了加减和移位,对结果四舍五入,供你参考比较! (程序是我摘录的,非原创) 如果需要软件实现的浮点库, 开源的Softfloat很好。
#define N_BITS 32
#define MAX_BIT ((N_BITS + 1) / 2 - 1)
unsigned long sqrt_1(unsigned long x)
{
register int i; register unsigned long m, r, root;
root = 0;
m = 1 << MAX_BIT;
for (i = MAX_BIT; i >= 0; i--)
{ r = root + m;
if (r * r <= x) root = r;
m >>= 1;
} return root;
}
unsigned long sqrt_5(unsigned long x)
{
register unsigned long xroot, m2, x2;
xroot = 0;
m2 = 1 << MAX_BIT * 2;
do
{
x2 = xroot + m2;
xroot
>>= 1;
if (x2 <= x) { x -= x2; xroot += m2; }
} while (m2 >>= 2);
if (xroot < x) return xroot + 1;
return xroot;
}