原创 转贴：开方运算

如果对一个32位无符号数开方，那么结果一定是一个16位无符号数。

现设被开方的数为a，开方结果：

b = b[15] * 2^15 + b[14] * 2^14 + ... + b[0] * 2^0

上式中的b[15]表示16位无符号数b的第15位（最高位），2^15表示2的15次方；以此类推。

开方的基本原理如下：
首先假设b[15] = 1，此时如果(b[15] * 2^15)^2 <= a成立，则b[15]=1，否则b[15]=0；然后再假设b[14]=1，此时如果(b[15] * 2^15 + b[14] * 2^14)^2 <= a成立，则b[14]=1，否则b[14]=0。以此类推，经过16次比较之后可以得出最终的开方结果。

程序如下：

unsigned int fixed_sqrt(unsigned long a)
{
    unsigned char i;
    unsigned int b="0",c;
    for(i=0; i<16; i++)
    {
        c = b | (1 << (15 - i);
        if(c * c <= a) b = c;
    }
    return b;
}

此种方法虽然简单，但是要进行16次16位乘16位的乘法运算，对于没有乘法器的单片机来讲速度太慢.

我们继续对上面的求解方法分析，目标是把不等式中的乘法运算消除。

第一次运算时，已经假设b[15]=1，所以不等式(b[15] * 2^15)^2 <= a可以转化为以下形式：
2^30 <= a

因为2^30可以通过移位运算来求解，所以成功的消除了不等式中的乘法。

第二次运算时，b[15]已知，已假设b[14]=1，所以不等式(b[15] * 2^15 + b[14] * 2^14)^2 <= a可以转化为以下形式：

(b[15] * 2^15)^2 + 2 * b[15] * 2^15 * b[14] * 2^14 + (2^14)^2 <= a

2 * b[15] * 2^15 * b[14] * 2^14 + (2^14)^2 <= a - (b[15] * 2^15)^2

2 * b[15] * 2^15 * 2^14 + 2^28 <= a - (b[15] * 2^15)^2

(b[15] * 2^15) * 2^15 + 2^28 <= a - (b[15] * 2^15)^2

注意观察不等式右端可知，当b[15]=1时，不等式右端的值可通过上一次运算时不等式两端的值相减求得。

第三次运算时，b[15]、b[14]已知，已假设b[13]=1 ，所以不等式(b[15] * 2^15 + b[14] * 2^14 + b[13] * 2^13)^2 <= a可以转化为以下形式：

((b[15] * 2^15 + b[14] * 2^14) + b[13] * 2^13)^2 <= a

(b[15] * 2^15 + b[14] * 2^14)^2 + 2 * (b[15] * 2^15 + b[14] * 2^14) * b[13] * 2^13 + (b[13] * 2^13)^2 <= a

2 * (b[15] * 2^15 + b[14] * 2^14) * b[13] * 2^13 + (b[13] * 2^13)^2 <= a - (b[15] * 2^15 + b[14] * 2^14)^2

2 * (b[15] * 2^15 + b[14] * 2^14) * 2^13 + (2^13)^2 <= a - (b[15] * 2^15 + b[14] * 2^14)^2

(b[15] * 2^15 + b[14] * 2^14）* 2^14 + (2^13)^2 <= a - (b[15] * 2^15 + b[14] * 2^14)^2

注意观察不等式右端可知，当b[14]=1时，不等式右端的值可通过上一次运算时不等式两端的值相减求得。

以此类推，便可得到以下函数：

unsigned int fixed_sqrt(unsigned long a)
{
    unsigned char i;
    unsigned int b = 0;
    unsigned long c;
    for(i=0; i<16; i++)
    {
        c = ((unsigned long)1 << ((15-i) << 1)) + ((unsigned long)b << (16 - i));
        if(c <= a )
        {
            a -= c;

            b += (1 << (15-i));
        }
    }
    return b;
}

下面的这段程序是一种运算速度更快的书写形式：
/*****************************************/
/*: 开根号处理                                */
/*入口参数：被开方数，32位整型           */
/*出口参数：开方结果，16位整型           */
/****************************************/
unsigned int sqrt_fixed(unsigned long a)
{
    unsigned long i,c;
    unsigned long b = 0;
    for(i = 0x40000000; i != 0; i >>= 2)
    {
        c = i + b;
        b >>= 1;
        if(c <= a)
        {
            a -= c;
            b += i;
        }
    }
    return (unsigned int)b;
}

/*****************************************/
/*: 开根号处理                                 */
/*入口参数：被开方数，16位整型            */
/*出口参数：开方结果，8位整型             */
/****************************************/
unsigned short sqrt_fixed(unsigned short a)
{
  unsigned short i,c;
  unsigned short b = 0;
  for(i = 0x4000; i!= 0; i >>= 2)
  {
    c = i + b;
    b >>= 1;
    if(c <= a)
 {
      a -= c;
      b += i;
 }
  }
  return (unsigned short)b;
}