单片机快速开平方的算法

C 语言中开平方的算法中要开平方的话,可以在头文件中加 #include <math.h> 。然后调 sqrt(n) 函数即可。但在低性能单片机中要开平方,可以用到下面算法。

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

我们先来看看十进制下是如何手工计算开方的。一般的数字是十进制表示法,形如 12 的字符串 ab 表示数字 10a+b, 开平方的过程就是找到一个数字(用某字符串表示)使它的平方逼近已知数字。先将数字每两位分一节,之所以两位一节,是因为 100 (两位)是 10 (一位)的平方。例如 1296 ,我们将其拆成两节, 12 看成一节, 96 看成一节。假设两位数字符 ab 的平方是 1296 ,我们现在要确定出字符 ab 。

笔算格式和除法相似,因为 3 的平方最接近 12 (就是 1296 的第一节),从而确定 a=3 ,然后和除法一样下移数字,只是这里下移两位,移下来组成新数 396 ,由于根的第一位已经确定,也就是 10a+b 中的 a 已经确定,现在要确定 b 。因为

\[(10a+b)^2=100a^2+20ab+b^2=100a^2+(20a+b)b\]

而 100a2 已经被从最高节中减掉了。

因此对于余下的新数 396 ,只需用 (20a+b)b 来逼近 396 ,使得 (20a+b)b<396 且 (20a+b+1)(b+1)>396 从而确定根的第二位 b 。具体过程如下操作,可以确定 1296 的平方根是 36 。

                   3  6       平方根
                 -------------------
                  12,96       平方数
                   9
                 -------
20a+b=20x3+6=66    3 96
                   3 96
                   -----
                      0

我们再举个更普遍的例子,求 15239.91 的平方根。先将 15239.91 以小数点为基点每两位分一节 1,52,39.91 ;分成四节对应平方根至少是 4 位有效数字的小数(用字符 abc.d 表示),最高节 1 ,对应根的最高位 1 (a=1) 。现在确定根的第二位、第三位和第四位(十位、个位和十分位),按照上例的方法确定根的第二位为 2 (b=2) 。确定了 b ,接下来再次相减,即把刚确定的 (20a+b)b 减掉,从而再次获得新余数 8 ,然后下移两位,移下来组成新数 839 。如此循环,得到平方根为 123.4 。注意到最后的余数不为零,我们可以继续每次循环补两个零以求得更精确平方根。如补两个零后求得平方根为 123.45 。

\[\begin{split}&(100a+10b+c+0.1d)^2 \\ &=(100a+10b+c)^2+\frac{(20(100a+10b+c) + d)d}{100} \\ &=100(10a+b)^2+(20(10a+b)+c)c+\frac{(20(100a+10b+c) + d)d}{100} \\\end{split}\]
                              1  2  3. 4     平方根
                         --------------------------
                              1,52,39.91     平方数
                              1
                            -------------
20a+b=20x1+2=22                 52
                                44
                              -----------
20(10a+b)+c=20x12+3=243          8 39
                                 7 29
                               ----------
20(100a+10b+c)+d=20x123+4=2464   1 10 91
                                   98 56
                                 --------
                                   12 35

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

\[(2a+b)^2=4a^2+4ab+b^2=4a^2+(4a+b)b\]

我们来看一个例子,计算 100 ( 二进制 1100100) 的平方根。这里每一步不再是把 a 乘以 20 ,而是把 a 乘以 4 ,也就是把 a 右移两位,而由于 b 的值只能为 0 或者 1 ,所以我们只需要判断余数和 4a+1 的大小关系,如果余数大于等于 4a+1 那么该节的商是 1 ,否则该节的商是 0 。

              1  0  1  0    平方根
             ---------------------
              1,10,01,00    平方数
              1
             ------------
100x1+0=100     10
                 0
             ------------
100x10+1=1001   10 01
                10 01
             ------------
                      00

下面给出完成的 C 语言程序,其中 root 表示 p , rem 表示每步计算之后的余数, divisor 表示 4p+1 ,通过 value>>30value 的最高 2 位,通过 value<<=2 将计算后的最高 2 位剔除。其中 root 的两次 <<1 相当于 4p 。程序完全是按照手工计算改写的,应该不难理解。做完 16 次循环,就可以求出 32 位无符号整数的平方根的整数部分。值得注意的是,在 16 次循环之后还多做了一步求根操作,这是为了求出平方根小数点后第一个数,然后四舍五入到整数。网络上一般的求平方根算法是直接舍弃平方根小数部分取整,并不会四舍五入取整。

uint16_t sqrt(uint32_t value)
{
    uint32_t rem = 0;
    uint32_t divisor = 0;
    uint32_t root = 0;
    for(uint32_t i=0; i<16; i++)
    {
        root <<= 1;
        rem = ((rem << 2) + (value >> 30));
        value <<= 2;
        divisor = (root << 1) + 1;
        if(divisor <= rem)
        {
            rem -= divisor;
            root++;
        }
    }

    rem <<= 2;
    divisor = (root << 2) + 1;
    if(divisor <= rem)
    {
        root++;
    }
    return (uint16_t)(root);
}