单片机快速开平方的算法 ====================== C 语言中开平方的算法中要开平方的话，可以在头文件中加 ``#include `` 。然后调 ``sqrt(n)`` 函数即可。但在低性能单片机中要开平方，可以用到下面算法。 本算法只采用移位、加减法、判断和循环实现，因为它不需要浮点运算，也不需要乘除运算，因此可以很方便地运用到各种芯片上去。 我们先来看看十进制下是如何手工计算开方的。一般的数字是十进制表示法，形如 12 的字符串 ab 表示数字 10a+b，开平方的过程就是找到一个数字（用某字符串表示）使它的平方逼近已知数字。先将数字每两位分一节，之所以两位一节，是因为 100（两位）是 10（一位）的平方。例如 1296，我们将其拆成两节，12 看成一节，96 看成一节。假设两位数字符 ab 的平方是 1296，我们现在要确定出字符 ab。 笔算格式和除法相似，因为 3 的平方最接近 12（就是 1296 的第一节），从而确定 a=3，然后和除法一样下移数字，只是这里下移两位，移下来组成新数 396，由于根的第一位已经确定，也就是 10a+b 中的 a 已经确定，现在要确定 b。因为 .. code:: none (10+b)² = 100a² + (20a+b)b 而 100a\ :sup:`2`\ 已经被从最高节中减掉了。 因此对于余下的新数 396，只需用 (20a+b)b 来逼近 396，使得 (20a+b)b < 396 且 (20a+b+1)(b+1) > 396 从而确定根的第二位 b。具体过程如下操作，可以确定 1296 的平方根是 36。 .. code:: none 3 6 平方根 ------------------- 12,96 平方数 9 ------- 20a+b=20×3+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。 .. code:: none (100a+10b+c+0.1d)² = (100a+10b+c)²+0.01(20(100a+10b+c)+d)d = (100a+10b)²+(20(10a+1b)+c)c+0.01(20(100a+10b+c)+d)d = 10000a²+100(20a+b)b+(20(10a+1b)+c)c+0.01(20(100a+10b+c)+d)d .. code:: none 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 这个手工算法其实和十进制关系不大，因此我们可以很容易的把它改为二进制，改为二进制之后，公式就变成了： .. code:: none (2a+b)² = 4a²+(4a+b)b 我们来看一个例子，计算 100（二进制 1100100）的平方根。这里每一步不再是把 a 乘以 20，而是把 a 乘以 4，也就是把 a 右移两位，而由于 b 的值只能为 0 或者 1，所以我们只需要判断余数和 4a+1 的大小关系，如果余数大于等于 4a+1 那么该节的商是 1，否则该节的商是 0。 .. code:: none 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>>30`` 取 ``value`` 的最高 2 位，通过 ``value<<=2`` 将计算后的最高 2 位剔除。其中 ``root`` 的两次 ``<<1`` 相当于 4p。程序完全是按照手工计算改写的，应该不难理解。做完 16 次循环，就可以求出 32 位无符号整数的平方根的整数部分。值得注意的是，在 16 次循环之后还多做了一步求根操作，这是为了求出平方根小数点后第一个数，然后四舍五入到整数。 .. code:: none 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); }