牛頓迭代法(Newton's Method)
阿新 • • 發佈:2018-12-26
簡介
牛頓迭代法(簡稱牛頓法)由英國著名的數學家牛頓爵士最早提出。但是,這一方法在牛頓生前並未公開發表。
對於形如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年將它列入了非法遊戲名單)
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