如何在单片机中实现快速开平方的算法

C语言中开平方的算法中要开平方的话，可以在头文件中加#include 。然后调sqrt（n）;函数即可。但在单片机中要开平方。可以用到下面算法：

算法1：

本算法只采用移位、加减法、判断和循环实现，因为它不需要浮点运算，也不需要乘除运算，因此可以很方便地运用到各种芯片上去。

我们先来看看10进制下是如何手工计算开方的。

先看下面两个算式，

x = 10*p + q （1）

公式（1）左右平方之后得：

x^2 = 100*p^2 + 20pq + q^2 （2）

现在假设我们知道x^2和p，希望求出q来，求出了q也就求出了x^2的开方x了。

我们把公式（2）改写为如下格式：

q = （x^2 - 100*p^2）/（20*p+q） （3）

这个算式左右都有q，因此无法直接计算出q来，因此手工的开方算法和手工除法算法一样有一步需要猜值。

我们来一个手工计算的例子：计算1234567890的开方

首先我们把这个数两位两位一组分开，计算出最高位为3。也就是（3）中的p，最下面一行的334为余数，也就是公式（3）中的（x^2 - 100*p^2）近似值

3 --------------- | 12 34 56 78 90 9 --------------- | 3 34

下面我们要找到一个0-9的数q使它最接近满足公式（3）。我们先把p乘以20写在334左边：

3 q --------------- | 12 34 56 78 90 9 --------------- 6q| 3 34

我们看到q为5时（60+q*q）的值最接近334，而且不超过334。于是我们得到：

3 5 --------------- | 12 34 56 78 90 9 --------------- 65| 3 34 | 3 25 --------------- 9 56

接下来就是重复上面的步骤了，这里就不再啰嗦了。

这个手工算法其实和10进制关系不大，因此我们可以很容易的把它改为二进制，改为二进制之后，公式（3）就变成了：

q = （x^2 - 4*p^2）/（4*p+q） （4）

我们来看一个例子，计算100（二进制1100100）的开方：

1 0 1 0 --------------- | 1 10 01 00 1 --------------- 100| 0 10 | 0 00 --------------- | 10 011001| 10 01 --------------- 0 00

这里每一步不再是把p乘以20了，而是把p乘以4，也就是把p右移两位，而由于q的值只能为0或者1，所以我们只需要判断余数（x^2 - 4*p^2）和（4*p+1）的大小关系，如果余数大于等于（4*p+q）那么该上一个1，否则该上一个0。

下面给出完成的C语言程序，其中root表示p，rem表示每步计算之后的余数，divisor表示（4*p+1），通过a》》30取a的最高 2位，通过a《《=2将计算后的最高2位剔除。其中root的两次《《1相当于4*p。程序完全是按照手工计算改写的，应该不难理解。

unsigned short sqrt（unsigned long a）{

unsigned long rem = 0;

unsigned long root = 0;

unsigned long divisor = 0;

for（int i=0; i《16; i++）{

root 《《= 1;

rem = （（rem 《《 2） + （a 》》 30））;

a 《《= 2;

divisor = （root《《1） + 1;

if（divisor 《= rem）{

rem -= divisor;

root++;

}

}

return （unsigned short）（root）;

}

算法2 ：单片机开平方的快速算法

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

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

绍给大家，希望会有些帮助。

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 《 pow（2， 1/2+k） 《 pow（2， k+1），

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. 实现代码

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

-------------------------------------------------------------------------------

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

/*Function： 开根号处理 */

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

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

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

unsigned int sqrt_16（unsigned long M）

{

unsigned int N， i;

unsigned long tmp， ttp; // 结果、循环计数

if （M == 0） // 被开方数，开方结果也为0

return 0;

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 ++;

}

}

return N;

}

以上都是网络查找的资料，有些晦涩难懂，不过在实际运用中可以使用这些算法。