1. 程式人生 > >牛頓迭代法(Newton's Method)

牛頓迭代法(Newton's Method)

簡介

牛頓迭代法(簡稱牛頓法)由英國著名的數學家牛頓爵士最早提出。但是,這一方法在牛頓生前並未公開發表。


牛頓法的作用是使用迭代的方法來求解函式方程的根。簡單地說,牛頓法就是不斷求取切線的過程。

對於形如f(x)=0的方程,首先任意估算一個解x0,再把該估計值代入原方程中。由於一般不會正好選擇到正確的解,所以有f(x)=a。這時計算函式在x0處的斜率,和這條斜率與x軸的交點x1。

f(x)=0中精確解的意義是,當取得解的時候,函式值為零(即f(x)的精確解是函式的零點)。因此,x1比x0更加接近精確的解。只要不斷以此方法更新x,就可以取得無限接近的精確的解。

但是,有可能會遇到牛頓迭代法無法收斂的情況。比如函式有多個零點,或者函式不連續的時候。

牛頓法舉例

下面介紹使用牛頓迭代法求方根的例子。牛頓迭代法是已知的實現求方根最快的方法之一,只需要迭代幾次後就能得到相當精確的結果。

首先設x的m次方根為a。

下面程式使用牛頓法求解平方根。

const float EPS = 0.00001; 
int sqrt(double x) { 
	if(x == 0) return 0; 
	double result = x; /*Use double to avoid possible overflow*/ 
	double lastValue; 
	do{ 
		lastValue = result; 
		result = result / 2.0f + x / 2.0f / result; 
	}while(abs(result - lastValue) > EPS);
 return (double)result;
 }

更快的方法

文獻2提到了比上述程式更快的求解平方根的非典型牛頓迭代法。介紹如下。

1999年12月,美國id Software公司釋出了名為“雷神之錘III”的電子遊戲。它是第一個支援軟體加速的遊戲,取得了極大成功。(由於影響力過大,文化部於2004年將它列入了非法遊戲名單)


雷神之錘III並不是id Software公司的第一次成功。早在1993年開始,這家公司就以“毀滅戰士”系列遊戲名聞天下。1995年,“毀滅戰士”的安裝數超過了當年微軟的windows 95。據傳比爾蓋茨才曾經考慮買下id software。(id software公司後來被推出過“上古卷軸”系列的Bethesda公司買下)

id Software所取得的成功很大程度上要歸功於它的創始人約翰·卡馬克。馬克爾也是一個著名的程式設計師,他是id Software遊戲引擎的主要負責人。 回到剛才提到的雷神之錘,馬克爾是開源軟體的積極推動者,他於2005年公佈了雷神之錘III的原始碼。至此人們得以通過研究這款遊戲引擎的原始檔來檢視它成功的祕密。

在其中一個名字為q_math.c的檔案中發現瞭如下程式碼段。
float Q_rsqrt( float number ) { 
	long i; float x2, y; const float threehalfs = 1.5F;
	x2 = number * 0.5F; 
	y = number; 
	i = * ( long * ) &y; // evil floating point bit level hacking 
	i = 0x5f3759df - ( i >> 1 ); // what the fuck? 
	y = * ( float * ) &i; 
	y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration 
	// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
	#ifndef Q3_VM #
	ifdef __linux__ assert( !isnan(y) ); // bk010122 - FPE?
	#endif
	#endif return y; 
}

這段程式碼的作用就是求number的平方根,並且返回它的倒數。

經過測試,它的效率比上述牛頓法程式要快幾十倍。也比c++標準庫的sqrt()函式要快好幾倍。此段程式碼有一個奇怪的句子:

i = 0x5f3759df - ( i >> 1 ); // what the fuck? 

這句話的註釋是“what the fuck?”,翻譯過來就是“我靠?”

任何受過程式訓練的人看到這句大概都會在想,這句話到底在搞什麼鳥?

之所以會出現這種奇怪的註釋,要麼是此段程式的作者(可能是馬克爾)根本不知道該如何解釋清楚,或者是維護這段程式的程式設計師完全看不懂這句話,所以有點兒抓毛。而實際上,它的作用(再加上y = y * ( threehalfs - ( x2 * y * y ) )這句牛頓迭代)就是求平方根。

至於是為什麼,本博主也不知道。

以雷神之錘III程式為藍本可以寫出比sqrt()更強大的求平方根函式:
int sqrt(float x) { 
	if(x == 0) return 0; 
	float result = x; 
	float xhalf = 0.5f*result; 
	int i = *(int*)&result; 
	i = 0x5f375a86- (i>>1); // what the fuck? 
	result = *(float*)&i; 
	result = result*(1.5f-xhalf*result*result); // Newton step, repeating increases accuracy 
	result = result*(1.5f-xhalf*result*result); 
	return 1.0f/result; 
}


參考文獻:
1.wikipedia.org
2.http://www.2cto.com/kf/201206/137256.html