开放运算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;//93ms

	start=clock();
	sqrtNewton(A,B,N);
	end=clock();
	cout<<"牛顿迭代法："<<end-start<<"ms"<<endl;//1474ms

	start=clock();
	sqrtImproveNewton(A,B,N);
	end=clock();
	cout<<"改进的牛顿迭代法："<<end-start<<"ms"<<endl;//539ms

	system("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 y0
	x = *(float*)&i; // convert bits BACK to float
	x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
	x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
	x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy

	return 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++的缘故？这个不清楚，反正是证实了库函数的强大！！！