这次我们来学习一下如何求平方根。在计算机中很难有精确的求出数据的平方根的算法，基本都是要求一个误差可接受范围内的近似值。治理我们取误差值为1e-5。
笔者的编程语言及环境如下
- 编程语言：c++
- 编译器：Code Blocks
- 系统： windows 10 x64
求平方根的方法
- 二分法
- 递归法（牛顿法）
- 卡马克法
- 完整代码展示
二分法
二分法求平方根的思路如下：
- 设置一个初始区间[0,n]，若n<1则初始区间为[0,1]。
- 取区间中值mid的平方与数据n比较。
- 若数据n小则令区间变为[0,mid]；否则变为[mid,n(或1)]。
- 若先后两次中值的差在误差范围内，则此次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为测试数据结果如下
递归法（牛顿法）
迭代法求平方根的主要思路为：
- 令x为任意值（通常取数据n的一半）。
- 令y=（x+n/x）/2。
- 若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;xint i = *(int*)&x; // 将x值按计算机中的存储格式强转成int赋值给ii = 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 valuei = 0x5f3759df - (i >> 1); // gives initial guess//现在float最优魔值是0x5f375a86 所以上面的代码可以用 i = 0x5f375a86 - (i >> 1); 替换x = *(float*)&i; // convert bits back to floatx = x * (1.5 - xhalf * x * x); // 1st Newton stepx = x * (1.5 - xhalf * x * x); // 2nd Newton stepreturn 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 iterationreturn y;
}
