关于sqrt函数的算法

我们平时经常会有一些数据运算的操作,需要调用sqrt,exp,abs等函数,那么时候你有没有想过:这个些函数系统是如何实现的?就拿最常用的sqrt函数来说吧,系统怎么来实现这个经常调用的函数呢?

虽然有可能你平时没有想过这个问题,不过正所谓是“临阵磨枪,不快也光”,你“眉头一皱,计上心来”,这个不是太简单了嘛,用二分的方法,在一个区间中,每次拿中间数的平方来试验,如果大了,就再试左区间的中间数;如果小了,就再拿右区间的中间数来试。比如求sqrt(16)的结果,你先试(0+16)/2=8,8*8=64,64比16大,然后就向左移,试(0+8)/2=4,4*4=16刚好,你得到了正确的结果sqrt(16)=4。然后你三下五除二就把程序写出来了:

float SqrtByBisection(float n) //用二分法
{
	if(n<0) //小于0的按照你需要的处理
		return n;
	float mid,last;
	float low,up;
	low=0,up=n;
	mid=(low+up)/2;
	do
	{
		if(mid*mid>n)
			up=mid;
		else
			low=mid;
		last=mid;
		mid=(up+low)/2;
	}while(abs(mid-last) > eps);//精度控制
	return mid;
} 

然后看看和系统函数性能和精度的差别(其中时间单位不是秒也不是毫秒,而是CPU Tick,不管单位是什么,统一了就有可比性) 

从图中可以看出,二分法和系统的方法结果上完全相同,但是性能上整整差了几百倍。为什么会有这么大的区别呢?难道系统有什么更好的办法?难道。。。。哦,对了,回忆下我们曾经的高数课,曾经老师教过我们“牛顿迭代法快速寻找平方根”,或者这种方法可以帮助我们,具体步骤如下:

求出根号a的近似值:首先随便猜一个近似值x,然后不断令x等于x和a/x的平均数,迭代个六七次后x的值就已经相当精确了。  例如,我想求根号2等于多少。假如我猜测的结果为4,虽然错的离谱,但你可以看到使用牛顿迭代法后这个值很快就趋近于根号2了:  (       4  + 2/4        ) / 2 = 2.25  (     2.25 + 2/2.25     ) / 2 = 1.56944..  ( 1.56944..+ 2/1.56944..) / 2 = 1.42189..  ( 1.42189..+ 2/1.42189..) / 2 = 1.41423..  .... 这种算法的原理很简单,我们仅仅是不断用(x,f(x))的切线来逼近方程x^2-a=0的根。根号a实际上就是x^2-a=0的一个正实根,这个函数的导数是2x。也就是说,函数上任一点(x,f(x))处的切线斜率是2x。那么,x-f(x)/(2x)就是一个比x更接近的近似值。代入 f(x)=x^2-a得到x-(x^2-a)/(2x),也就是(x+a/x)/2。

相关的代码如下:

float SqrtByNewton(float x)
{
	float val = x;//最终
	float last;//保存上一个计算的值
	do
	{
		last = val;
		val =(val + x/val) / 2;
	}while(abs(val-last) > eps);
	return val;
}

然后我们再来看下性能测试:

哇塞,性能提高了很多,可是和系统函数相比,还是有这么大差距,这是为什么呀?想啊想啊,想了很久仍然百思不得其解。突然有一天,我在网上看到一个神奇的方法,于是就有了今天的这篇文章,废话不多说,看代码先:

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;
}

然后我们最后一次来看下性能测试:

这次真的是质变了,结果竟然比系统的还要好。。。

时间: 2024-12-12 15:09:40

关于sqrt函数的算法的相关文章

15、【常见算法】实现sqrt函数

问题:手写代码实现sqrt函数,即求一个整数的平方根 分析:二分查找思想 1 #include <iostream> 2 3 using namespace std; 4 5 int sqrt(int x) 6 { 7 long left = 0; 8 9 if(x == 1) 10 return 1; 11 long right = x; 12 13 long mid = left + (right - left)/2; 14 while(left + 1 < right) 15 {

一个Sqrt函数引发的血案

我们平时经常会有一些数据运算的操作,需要调用sqrt,exp,abs等函数,那么时候你有没有想过:这个些函数系统是如何实现的?就拿最常用的sqrt函数来说吧,系统怎么来实现这个经常调用的函数呢? 虽然有可能你平时没有想过这个问题,不过正所谓是"临阵磨枪,不快也光",你"眉头一皱,计上心来",这个不是太简单了嘛,用二分的方法,在一个区间中,每次拿中间数的平方来试验,如果大了,就再试左区间的中间数:如果小了,就再拿右区间的中间数来试.比如求sqrt(16)的结果,你先试

【转载】一个Sqrt函数引发的血案

转自:http://www.cnblogs.com/pkuoliver/archive/2010/10/06/sotry-about-sqrt.html 源码下载地址:http://diducoder.com/sotry-about-sqrt.html 好吧,我承认我标题党了,不过既然你来了,就认真看下去吧,保证你有收获. 我们平时经常会有一些数据运算的操作,需要调用sqrt,exp,abs等函数,那么时候你有没有想过:这个些函数系统是如何实现的?就拿最常用的sqrt函数来说吧,系统怎么来实现这

sqrt函数实现之卡马克方法

sqrt函数的实现主要有三种方式: 二分法 牛顿法 卡马克方法 这里主要介绍高效的卡马克方法.卡马克方法起源于<雷神之锤III竞技场>中使用的平方根倒数速算法,下列代码是平方根倒数速算法在<雷神之锤III竞技场>源代码中的应用实例.示例剥离了C语言预处理器的指令,但附上了原有的注释: float Q_rsqrt( float number ) { long i; float x2, y; const float threehalfs = 1.5F; x2 = number * 0.

二叉树链表中一个结论的证明以及 Successor() 函数中算法的解释

首先, 先明确 “结点 x 的 successor” 的概念. 这句话的概念是, successor 是整个二叉树中, key 仅比 x 的 key 大的结点. 求证: 若 x 有两个孩子, 那么其 successor 没有左孩子. 证明: 若 x 有两个孩子, 则其右子树存在.在二叉树链表中, 对于任意 x 的右子树都是 key 比 x 的 key 大的结点的集合, 即{y|y->key > x->key}.根据 successor  的定义, 若 x->rightChild ≠

四舍五入的一些简单写法(利用内置函数,算法2种写法)

?       //内置函数的写法        //网上零售价和折扣价在计算结束需要进行进位,规则如下:         //个位为1,2,3,4进位到5,例如计算后的价格为1201,则价格为1205:         //个位为6,7,8,9进位到0,例如计算后的价格为1209.则价格为1210:         public static string ChangePrice(double price)         {             int changed = 0;     

自己实现的一个sqrt函数

看了MIT的6.00的计算机科学导论课的第一节,其中,老师讲的一个sqrt的实现算法记忆犹新. 代码如下: #include <stdio.h> #define E 1e-10 double Abs(double a,double b){ if(a-b<0){ return b - a; } return a - b; } double Sqrt(double value){ double i,last=0,cur=0; for(i=0;i*i<value;i++); i--; d

C语言计算程序中某一个函数或算法的执行时间

计算程序中某一个函数或算法的执行时间 #include <stdio.h> #include <time.h> #include <stdlib.h> int main() { long i = 10000000L; clock_t start, finish; double duration; printf( "Time to do %ld empty loops is ", i) ; start = clock(); while( i-- );

STL进阶--成员函数 vs 算法

容器的成员函数 vs 算法 容器中同名的函数 List: void remove(const T); template<class Comp> void remove_if(Comp); void unique(); template<class Comp> void unique(Comp); void sort(); template<class Comp> void sort(Comp); void merge(list&); template<cl