无符号数开方源程序代码

0

一．如果对一个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;
}

下面的这段程序是一种运算速度更快的书写形式：
/*****************************************/
/*Function: 开根号处理                                */
/*入口参数：被开方数，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;
}//在VC6下通过测试

/*****************************************/
/*Function: 开根号处理                                 */
/*入口参数：被开方数，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;
}//在VC6下通过测试
二．
float insqrt(float x)
{
float y=x;
float xhalf=0.5f*x;
long  i=*(long*)&x;
i=0x5f375a86-(i>>1);
x=*(float*)&i;
x=x*(1.5f-xhalf*x*x);
x=x*(1.5f-xhalf*x*x);
return y*x;
}

0