开放运算sqrt(《编程珠玑(续)》)

下面程序是对欧式距离的计算，其中第一个程序调用的是库函数，而第二个程序是用牛顿迭代法进行计算sqrt，第三个程序是对第二个程序的优化，减少了循环和收敛测试的开销。

#include<iostream>#include<algorithm>#include<ctime>#include<iomanip>using namespace std;void mySqrt(double A[],double B[],int N){int K=2;//每K组数算一次欧式距离double sum=0;double temp=0;for(int j=0;j<N-K;j++){for(int i=j;i<K;i++){temp=abs(A[i]-B[i]);sum+=temp*temp;}sqrt(sum);}}void sqrtNewton(double A[],double B[],int N){int K=2;//每K组数算一次欧式距离double sum=0;double temp=0;double maxValue=0;for(int j=0;j<N-K;j++){for(int i=j;i<K;i++){temp=abs(A[i]-B[i]);maxValue=max(maxValue,temp);sum+=temp*temp;}double eps=1.0e-7;double x=maxValue;//初始值为10个值中最大的double preResult=maxValue;while(1){x=0.5*(x+sum/x);if(abs(x-preResult)<eps)break;preResult=x;}}}void sqrtImproveNewton(double A[],double B[],int N){int K=2;//每K组数算一次欧式距离double sum=0;double temp=0;double maxValue=0;for(int j=0;j<N-K;j++){for(int i=j;i<K;i++){temp=abs(A[i]-B[i]);maxValue=max(maxValue,temp);sum+=temp*temp;}maxValue*=2;maxValue=0.5*(maxValue+sum/maxValue);maxValue=0.5*(maxValue+sum/maxValue);maxValue=0.5*(maxValue+sum/maxValue);maxValue=0.5*(maxValue+sum/maxValue);}}int main(){const int N=10000000;double *A=new double[N],*B=new double[N];generate(A,A+N,rand);generate(B,B+N,rand);clock_t start,end;start=clock();mySqrt(A,B,N);end=clock();cout<<"调用库函数sqrt："<<end-start<<"ms"<<endl;//93msstart=clock();sqrtNewton(A,B,N);end=clock();cout<<"牛顿迭代法："<<end-start<<"ms"<<endl;//1474msstart=clock();sqrtImproveNewton(A,B,N);end=clock();cout<<"改进的牛顿迭代法："<<end-start<<"ms"<<endl;//539mssystem("pause");return 0;}

运行结果：

可见，还是调用库函数的效率高，这一点和《编程珠玑》中的不一样，我也没有找出为什么书上用牛顿迭代法会比库函数快，库函数中的初始值选择也是很有讲究的，在《C标准库》中sprt的实现如下：

没有找到这里选择初始值的依据，但是结果是比我写的有效。

为此，我专门拿出sqrt函数来和库函数进行对比：其中神奇的方法来自：http://blog.csdn.net/stormbjm/article/details/8191737

#include<iostream>#include<algorithm>#include<ctime>#include<iomanip>using namespace std;double sqrtNewton(double a){double eps=1.0e-7;double x=a;//初始值double preResult=x;while(1){x=0.5*(x+a/x);if(abs(x-preResult)<eps)break;preResult=x;//cout<<setprecision(20)<<x<<endl;}return x;}float InvSqrt(float x){float xhalf = 0.5f*x;int i = *(int*)&x; // get bits for floating VALUE i = 0x5f375a86- (i>>1); // gives initial guess y0x = *(float*)&i; // convert bits BACK to floatx = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracyx = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracyx = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracyreturn 1/x;}int main(){<span style="white-space:pre"></span>const int N=10000000;<span style="white-space:pre"></span>clock_t start,end;<span style="white-space:pre"></span>start=clock();<span style="white-space:pre"></span>for(int i=0;i<N;i++)<span style="white-space:pre"></span>sqrtNewton(65536.0);<span style="white-space:pre"></span>end=clock();<span style="white-space:pre"></span>cout<<"简单牛顿迭代："<<end-start<<"ms"<<endl;//6763ms<span style="white-space:pre"></span>start=clock();<span style="white-space:pre"></span>for(int i=0;i<N;i++)<span style="white-space:pre"></span>sqrt(65536.0);<span style="white-space:pre"></span>end=clock();<span style="white-space:pre"></span>cout<<"库函数："<<end-start<<"ms"<<endl;//79ms<span style="white-space:pre"></span>start=clock();<span style="white-space:pre"></span>for(int i=0;i<N;i++)<span style="white-space:pre"></span>InvSqrt(65536.0);<span style="white-space:pre"></span>end=clock();<span style="white-space:pre"></span>cout<<"所谓的神奇方法："<<end-start<<"ms"<<endl;//657ms<span style="white-space:pre"></span>system("pause");<span style="white-space:pre"></span>return 0;}

结果：

经过对多个不同的数的开方运算，简单牛顿迭代法是受不同数字影响最大的，因为它的初始值就是值本身，库函数的运行时间几乎和要开方的数的值没有关系，可见其所选取的初始值的效率是很高的。

而且发现库函数比”神奇的方法“更有效，也许是我用的C++的缘故？这个不清楚，反正是证实了库函数的强大！！！