[从零开始学算法]求平方根

这次我们来学习一下如何求平方根。在计算机中很难有精确的求出数据的平方根的算法，基本都是要求一个误差可接受范围内的近似值。治理我们取误差值为1e-5。
笔者的编程语言及环境如下

1. 编程语言：c++
2. 编译器：Code Blocks
3. 系统： windows 10 x64
   

二分法

二分法求平方根的思路如下：

1. 设置一个初始区间[0,n]，若n<1则初始区间为[0,1]。
2. 取区间中值mid的平方与数据n比较。
3. 若数据n小则令区间变为[0,mid]；否则变为[mid,n(或1)]。
4. 若先后两次中值的差在误差范围内，则此次mid为平方根的近似值。

代码如下：

double sqrt(double number)//二分法
{
    double low=0.0,high=number>1?number:1;
    double mid=0,last=0;
    do{
            last=mid;
        mid=(low+high)/2;
    if(mid*mid>number)high=mid;
    else low=mid;
    }while(fabs(mid-last)>1e-5);
    return mid;
}

以123456789012为测试数据结果如下

递归法（牛顿法）

迭代法求平方根的主要思路为：

1. 令x为任意值（通常取数据n的一半）。
2. 令y=（x+n/x）/2。
3. 若x与y的差在误差范围内，则y为平方根近似值；否则令x=y执行第2步

代码如下：

double sqrt1(double number,double x)
{
    double y=(x+number/x)/2;
    return fabs(y-x)<0.001?y:sqrt1(number,y);
}

以123456789012测试结果如下：

卡马克法

可以在上文看到牛顿迭代算法的初值是任意的，一个好的初值可以大大降低迭代次数。卡马克算法所求的是平方根的倒数（在实际应用中平方根的倒数比平方根用的多）。卡马克算法可以理解为优化的迭代算法，其选择了一个极好的初值，将算法优化到一次或两次迭代就可以达到误差范围内。因为其涉及到浮点数在计算机内的存储方式及Magic Number，这里将不再赘述，有兴趣的可以自行百度。
因为其涉及到浮点数的存储，而系统中float与double的存储不一样，所以有float与double两种版本。
float版代码如下：

float sqrt(float x)
{
    float xhalf = 0.5 * x;x
    int i = *(int*)&x; // 将x值按计算机中的存储格式强转成int赋值给i
    i = 0x5f3759df - (i >> 1); 
    //现在float最优魔值是0x5f375a86 所以上面的代码可以用  i = 0x5f375a86 - (i >> 1); 替换
    x = *(float*)&i; 
    x = x * (1.5 - xhalf * x * x); // 一次牛顿迭代
    //x = x * (1.5 - xhalf * x * x); // 二次牛顿迭代
    return x;
}

说实话，笔者对于 int i = (int)&x;这种写法也就见过这么一次，可以说卡马克算法无论是从思想上还是代码上都极具魅力，你能想象到求平方根算法连一个循环一个迭代都没有么！世上没有取巧的代码，只有强大的理论。
同样以123456789012测试，其效果如下：

我们可以看到数据同一二种方法有了数据差异，但是数据差异仅有0.1%，可以说是完全可以接受。如果迭代两次的话效果如下图：

误差来到了0.0002%，可太可怕了。
下面附上处理double型数据的代码，就不进行效果展示了：

double sqrt(double number)
{
    __int64 i;
    double x2, y;
    const double threehalfs = 1.5F;

    x2 = number * 0.5F;
    y  = number;
    i  = * ( __int64 * ) &y;
    i  = 0x5fe6ec85e7de30da - ( i >> 1 );
    y  = * ( double * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) ); // 一次牛顿迭代
    //y  = y * ( threehalfs - ( x2 * y * y ) ); // 二次迭代
    return y;
}

完整代码展示

现在对以上三种方法进行整体代码演示：

#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
double sqrt0(double number);//二分法
double sqrt1(double number,double x);//牛顿迭代法
float sqrt2(float x);//卡马克快速平方根
double sqrt3(double number);//double版
int main()
{
    freopen("in.txt","r",stdin);
    cout.setf(ios::fixed,ios::floatfield);
    cout.precision(7);
    double number,once,twice,others;
    cin>>number;
    cout<<"测试数据："<<number<<endl;
    cout<<"二分法求平方根："<<sqrt0(number)<<endl;
    cout<<"迭代法求平方根: "<<sqrt1(number,number/2)<<endl;
    others=sqrt0(number);
    once=1/sqrt3(number);//卡马克法得出的是平方根倒数
    twice=1/sqrt2(number);//在一定数据范围内float和double相互转化可以不损失精度
    cout<<"一次迭代卡马克法求平方根: "<<once<<"  相差比"<<fabs(once-others)/others<<endl;
    cout<<"二次迭代卡马克法求平方根: "<<twice<<"  相差比"<<fabs(twice-others)/others<<endl;
    return 0;
}
double sqrt0(double number)
{
    double low=0.0,high=number>1?number:1;
    double mid=0,last=0;
    do
    {
        last=mid;
        mid=(low+high)/2;
        if(mid*mid>number)high=mid;
        else low=mid;
    }
    while(fabs(mid-last)>1e-5);
    return mid;
}
double sqrt1(double number,double x)
{
    double y=(x+number/x)/2;
    return fabs(y-x)<0.001?y:sqrt1(number,y);
}
float sqrt2(float x)
{
    float xhalf = 0.5 * x;
    int i = *(int*)&x; // get bits for floating value
    i = 0x5f3759df - (i >> 1); // gives initial guess
    //现在float最优魔值是0x5f375a86 所以上面的代码可以用  i = 0x5f375a86 - (i >> 1); 替换
    x = *(float*)&i; // convert bits back to float
    x = x * (1.5 - xhalf * x * x); // 1st Newton step
    x = x * (1.5 - xhalf * x * x); // 2nd Newton step
    return x;
}
double sqrt3(double number)
{
    __int64 i;
    double x2, y;
    const double threehalfs = 1.5F;

    x2 = number * 0.5F;
    y  = number;
    i  = * ( __int64 * ) &y;
    i  = 0x5fe6ec85e7de30da - ( i >> 1 );
    y  = * ( double * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
    //y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration
    return y;
}