根号怎么开

因为工作的需要，要在单片机上实现开根号的操作。目前开平方的方法大部分是用牛顿迭代

法。我在查了一些资料以后找到了一个比牛顿迭代法更加快速的方法。不敢独享，介绍给大

家，希望会有些帮助。

1.原理

因为排版的原因，用pow(X,Y)表示X的Y次幂，用B[0]，B[1]，...，B[m-1]表示一个序列，

其中[x]为下标。

假设：

B[x],b[x]都是二进制序列,取值0或1。

M=B[m-1]*pow(2,m-1)+B[m-2]*pow(2,m-2)+...+B[1]*pow(2,1)+B[0]*pow

(2,0)

N=b[n-1]*pow(2,n-1)+b[n-2]*pow(2,n-2)+...+b[1]*pow(2,1)+n[0]*pow

(2,0)

pow(N,2)=M

(1)N的最高位b[n-1]可以根据M的最高位B[m-1]直接求得。

设m已知,因为pow(2,m-1)<=M<=pow(2,m)，所以pow(2,(m-1)/2)<=N<=

pow(2,m/2)

如果m是奇数，设m=2*k+1,

那么pow(2,k)<=N

n-1=k,n=k+1=(m+1)/2

如果m是偶数，设m=2k,

那么pow(2,k)>N>=pow(2,k-1/2)>pow(2,k-1),

n-1=k-1,n=k=m/2

所以b[n-1]完全由B[m-1]决定。

余数M[1]=M-b[n-1]*pow(2,2*n-2)

(2)N的次高位b[n-2]可以采用试探法来确定。

因为b[n-1]=1，假设b[n-2]=1，则pow(b[n-1]*pow(2,n-1)+b[n-1]*pow(2,n-2),

2)=b[n-1]*pow(2,2*n-2)+(b[n-1]*pow(2,2*n-2)+b[n-2]*pow(2,2*n-4)),

然后比较余数M[1]是否大于等于(pow(2,2)*b[n-1]+b[n-2])*pow(2,2*n-4)。这种

比较只须根据B[m-1]、B[m-2]、...、B[2*n-4]便可做出判断，其余低位不做比较。

若M[1]>=(pow(2,2)*b[n-1]+b[n-2])*pow(2,2*n-4),则假设有效，b[n-2]=

1；

余数M[2]=M[1]-pow(pow(2,n-1)*b[n-1]+pow(2,n-2)*b[n-2],2)=M[1]-

(pow(2,2)+1)*pow(2,2*n-4)；

若M[1]<(pow(2,2)*b[n-1]+b[n-2])*pow(2,2*n-4),则假设无效，b[n-2]=

0；余数M[2]=M[1]。

(3)同理，可以从高位到低位逐位求出M的平方根N的各位。

使用这种算法计算32位数的平方根时最多只须比较16次，而且每次比较时不必把M的各

位逐

一比较，尤其是开始时比较的位数很少，所以消耗的时间远低于牛顿迭代法。

2.流程图

(制作中，稍候再上)

3.实现代码

这里给出实现32位无符号整数开方得到16位无符号整数的C语言代码。

//--------------------------------------------------------------------------------

/****************************************/

/*Function:开根号处理*/

/*入口参数：被开方数，长整型*//*出口参数：开方结果，整型*/

/****************************************/

unsignedintsqrt_16(unsignedlongM)

{

unsignedintN,i;

unsignedlongtmp,ttp;//结果、循环计数

if(M==0)//被开方数，开方结果也为0

return0;

N=0;

tmp=(M>>30);//获取最高位：B[m-1]

M<<=2;

if(tmp>1)//最高位为1

{

N++;//结果当前位为1，否则为默认的0

tmp-=N;

}

for(i=15;i>0;i--)//求剩余的15位

{

N<<=1;//左移一位

tmp<<=2;

tmp+=(M>>30);//假设

ttp=N;

ttp=(ttp<<1)+1;

M<<=2;

if(tmp>=ttp)//假设成立

{

tmp-=ttp;N++;

}

}

returnN;

}

原来上面说的叫做“卡马克算法”：

求平方根的函数

更新：有人问这个算法的原理。其实原理很简单。就是牛顿迭代求根。卡马克算法牛X的

地方就是他选了一个常数作为起始值。而这个起始值让他只用一次迭代就够了。

从这里看来的。QuakeIII自然就是传奇高手卡马克的杰作之一了。在有的CPU上，这个函

数比普通的(float)(1.0/sqrt(x)快4倍！快的原因之一是用了一个神秘常数，0x5f3759df。普渡

大学的ChrisLomont在这篇论文里讨论了这个常数的意义，尝试用严格的方法推导出这个

常数（他还提到有人认为这个函数是在NVidia工作过的GaryTarolli写的）。Chris推出的

常数是0x5f37642f)，和Q_rsqrt里的稍有不同，而且实际表现也稍有不如。卡马克到底怎么

推出这个常数的就是谜了。这种高手不写书，实在可惜。

floatQ_rsqrt(floatnumber)

{

longi;

floatx2,y;

constfloatthreehalfs=1.5F;

x2=number*0.5F;

y=number;

i=*(long*)&y;//evilfloatingpointbitlevelhacking

i=0x5f3759df-(i>>1);//whatthefunc?

y=*(float*)&i;

y=y*(threehalfs-(x2*y*y));//1stiteration

//y=y*(threehalfs-(x2*y*y));//2nditeration,thiscanbe

removed

#ifndefQ3_VM

#ifdef__linux__

asrt(!isnan(y));//bk010122-FPE?

#endif[/#]

#endif[/#]

returny;

}