求开平方

定义

  所谓整数平方根即。

算法

  算法1.猜试法

利用等差级数公式:        

                
  这样的话, 从1开始一直算到数列的前项和第一次大于x的时候,即是所求。下面给出source code(C): 

unsigned linear_search(unsigned long x)
{
    unsigned long sum_n = 1;
    unsigned n = 1;

    if(x <= 1)
    {
        return x;
    }

    while(sum_n <= x)
    {
        n++;
        sum_n += (n<<1) - 1;
    }

    return (n-1);

}

 

  这种方法无异于穷举法，其唯一的优点是：每次的迭代用到了前面迭代的结果，所以会有一些效率的增益。对于该算法的改进就是不穷举，改用我们熟悉的二分查找法来做。

 

unsigned bi_search(unsigned long x)
{
    unsigned long sum_n = 0;
    unsigned n = (x >> 1);
    unsigned top = x;
    unsigned bottom = 0;

    if (x <= 1)
    {
        return x;
    }

    for (;;)
    {
        sum_n = n * n;
        if (sum_n < x)
        {
          bottom = n;
          n += ((top - bottom) >> 1);
          if (n == bottom)
              return n;
        }
        else 
        if (sum_n > x)
        {
          top = n;
          n -= ((top - bottom) >>1);
          if (n == top)
              return n-1;
        }
        else
        {
            return n;
        }
    }

}

 

算法2 Newton 法

把这个问题转换为方程求根问题，即:，求x。

 而方程求根的问题可以用Newton 法来解决。现在的问题有一点不同，即所求的根必须是整数。通过证明，我们可以发现，Newton迭代公式是适用于整数情况的，于是有：

                  
至于是怎么证明的，可以参考hacker’s delight。

另外，初值的选择也是很重要的一环，这里我们选择大于等于的最小的2的幂次数。

OK，下面给出程序：

 

unsigned newton_method(unsigned long x)
{
    unsigned long x1 = x - 1;
    unsigned s = 1;
    unsigned g0,g1;

    /* 初值设定 */
    if (x1 > 65535) {s += 8; x1 >>= 16;}
    if (x1 > 255)   {s += 4; x1 >>= 8;}
    if (x1 > 15)    {s += 2; x1 >>= 4;}
    if (x1 > 3)     {s += 1; x1 >>= 2;}

    /*迭代*/
    g0 = 1 << s;
    g1 = (g0 + (x >> s)) >> 1;
    while(g1 < g0)
    {
       g0 = g1;
       g1 = (g0 + x/g0) >> 1;
    }
    return g0;
}

浮点数：

首先是最普通的CRT里自带的sqrt，只需要引用math.h就可以使用了：

[cpp] view plaincopyprint?

1. #include <math.h>  
2. result = sqrt(number);  

 

 

接下来是传统的牛顿迭代法，我们计算开方的时候就是手工不断尝试每一位最合适的数字，然后一步步收敛求得更精确的答案。牛顿迭代法就是通过折半来快速收敛，每迭代一次就得到更精确的结果，下面这段程序采用的是迭代10次的牛顿迭代法：

[cpp] view plaincopyprint?

1. float newtonSqrt(float a){  
2.   int i;  
3.   float x;  
4.   x=a;  
5.   for(i=1;i<=10;i++)  
6.     x=(x+a/x)/2;  
7.   return x;  
8. }