牛頓迭代公式

2019-11-08 18:51:40

字體：大中小

供稿：網(wǎng)友

牛頓迭代法（Newton's method）又稱為牛頓-拉夫遜方法（Newton-Raphson method），它是牛頓在17世紀(jì)提出的一種在實(shí)數(shù)域和復(fù)數(shù)域上近似求解方程的方法。多數(shù)方程不存在求根公式，因此求精確根非常困難，甚至不可能，從而尋找方程的近似根就顯得特別重要。方法使用函數(shù)f(x)的泰勒級(jí)數(shù)的前面幾項(xiàng)來尋找方程f(x) = 0的根。牛頓迭代法是求方程根的重要方法之一，其最大優(yōu)點(diǎn)是在方程f(x) = 0的單根附近具有平方收斂，而且該法還可以用來求方程的重根、復(fù)根。另外該方法廣泛用于計(jì)算機(jī)編程中。

設(shè)r是f(x) = 0的根，選取x0作為r初始近似值，過點(diǎn)（x0,f(x0)）做曲線y = f(x)的切線L，L的方程為y = f(x0)+f'(x0)(x-x0)，求出L與x軸交點(diǎn)的橫坐標(biāo) x1 = x0-f(x0)/f'(x0)，稱x1為r的一次近似值。

過點(diǎn)（x1,f(x1)）做曲線y = f(x)的切線，并求該切線與x軸交點(diǎn)的橫坐標(biāo) x2 = x1-f(x1)/f'(x1)，稱x2為r的二次近似值。重復(fù)以上過程，得r的近似值序列，其中x(n+1)=x(n)－f(x(n))/f'(x(n))，稱為r的n+1次近似值，上式稱為牛頓迭代公式。

根據(jù)牛頓迭代的原理，可以得到以下的迭代公式：X(n+1)=[X(n)+p/Xn]/2

來看看WIKI：

其操作過程：

看一段牛人的程序：

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;   正常人寫的程序：
double sqr(double n) {     double k=1.0;     while(abs(k*k-n)>1e-9) {         k=(k+n/k)/2;     }     return k; }差別夠大的，但它們效果是一樣的，看來數(shù)學(xué)不好，還真不行啊。
在3D圖形編程中，經(jīng)常要求平方根或平方根的倒數(shù)，例如：求向量的長度或?qū)⑾蛄繗w一化。C數(shù)學(xué)函數(shù)庫中的sqrt具有理想的精度，但對(duì)于3D游戲程式來說速度太慢。我們希望能夠在保證足夠的精度的同時(shí)，進(jìn)一步提高速度。
Carmack在QUAKE3中使用了下面的算法，它第一次在公眾場(chǎng)合出現(xiàn)的時(shí)候，幾乎震住了所有的人。據(jù)說該算法其實(shí)并不是Carmack發(fā)明的，它真正的作者是Nvidia的Gary Tarolli（未經(jīng)證實(shí)）。－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－//// 計(jì)算參數(shù)x的平方根的倒數(shù)//float InvSqrt (float x){float xhalf = 0.5f*x;int i = *(int*)&x;i = 0x5f3759df - (i >> 1); // 計(jì)算第一個(gè)近似根x = *(float*)&i;x = x*(1.5f - xhalf*x*x); // 牛頓迭代法return x;}－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－該算法的本質(zhì)其實(shí)就是牛頓迭代法（Newton-Raphson Method，簡(jiǎn)稱NR），而NR的基礎(chǔ)則是泰勒級(jí)數(shù)（Taylor Series）。NR是一種求方程的近似根的方法。首先要估計(jì)一個(gè)與方程的根比較靠近的數(shù)值，然后根據(jù)公式推算下一個(gè)更加近似的數(shù)值，不斷重復(fù)直到可以獲得滿意的精度。其公式如下：－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－函數(shù)：y=f(x)其一階導(dǎo)數(shù)為：y'=f'(x)則方程：f(x)=0 的第n+1個(gè)近似根為x[n+1] = x[n] - f(x[n]) / f'(x[n])－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ NR最關(guān)鍵的地方在于估計(jì)第一個(gè)近似根。如果該近似根與真根足夠靠近的話，那么只需要少數(shù)幾次迭代，就可以得到滿意的解。
現(xiàn)在回過頭來看看如何利用牛頓法來解決我們的問題。求平方根的倒數(shù)，實(shí)際就是求方程1/(x^2)-a=0的解。將該方程按牛頓迭代法的公式展開為：x[n+1]=1/2*x[n]*(3-a*x[n]*x[n]) 將1/2放到括號(hào)里面，就得到了上面那個(gè)函數(shù)的倒數(shù)第二行。
接著，我們要設(shè)法估計(jì)第一個(gè)近似根。這也是上面的函數(shù)最神奇的地方。它通過某種方法算出了一個(gè)與真根非常接近的近似根，因此它只需要使用一次迭代過程就獲得了較滿意的解。它是怎樣做到的呢？所有的奧妙就在于這一行：i = 0x5f3759df - (i >> 1); // 計(jì)算第一個(gè)近似根
超級(jí)莫名其妙的語句，不是嗎？但仔細(xì)想一下的話，還是可以理解的。我們知道，IEEE標(biāo)準(zhǔn)下，float類型的數(shù)據(jù)在32位系統(tǒng)上是這樣表示的（大體來說就是這樣，但省略了很多細(xì)節(jié)，有興趣可以GOOGLE）：－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－bits：31 30 ... 031：符號(hào)位30-23：共8位，保存指數(shù)（E）22-0：共23位，保存尾數(shù)（M）－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－所以，32位的浮點(diǎn)數(shù)用十進(jìn)制實(shí)數(shù)表示就是：M*2^E。開根然后倒數(shù)就是：M^(-1/2)*2^(-E/2)?，F(xiàn)在就十分清晰了。語句i> >1其工作就是將指數(shù)除以2，實(shí)現(xiàn)2^(E/2)的部分。而前面用一個(gè)常數(shù)減去它，目的就是得到M^(1/2)同時(shí)反轉(zhuǎn)所有指數(shù)的符號(hào)。
至于那個(gè)0x5f3759df，呃，我只能說，的確是一個(gè)超級(jí)的Magic Number。
那個(gè)Magic Number是可以推導(dǎo)出來的，但我并不打算在這里討論，因?yàn)閷?shí)在太繁瑣了。簡(jiǎn)單來說，其原理如下：因?yàn)镮EEE的浮點(diǎn)數(shù)中，尾數(shù)M省略了最前面的1，所以實(shí)際的尾數(shù)是1+M。如果你在大學(xué)上數(shù)學(xué)課沒有打瞌睡的話，那么當(dāng)你看到(1+M)^(-1/2)這樣的形式時(shí)，應(yīng)該會(huì)馬上聯(lián)想的到它的泰勒級(jí)數(shù)展開，而該展開式的第一項(xiàng)就是常數(shù)。下面給出簡(jiǎn)單的推導(dǎo)過程：－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－對(duì)于實(shí)數(shù)R>0，假設(shè)其在IEEE的浮點(diǎn)表示中，指數(shù)為E，尾數(shù)為M，則：
R^(-1/2)= (1+M)^(-1/2) * 2^(-E/2)
將(1+M)^(-1/2)按泰勒級(jí)數(shù)展開，取第一項(xiàng)，得：
原式= (1-M/2) * 2^(-E/2)= 2^(-E/2) - (M/2) * 2^(-E/2)
如果不考慮指數(shù)的符號(hào)的話，(M/2)*2^(E/2)正是(R>>1)，而在IEEE表示中，指數(shù)的符號(hào)只需簡(jiǎn)單地加上一個(gè)偏移即可，而式子的前半部分剛好是個(gè)常數(shù)，所以原式可以轉(zhuǎn)化為：
原式 = C - (M/2)*2^(E/2) = C - (R>>1)，其中C為常數(shù)
所以只需要解方程：R^(-1/2)= (1+M)^(-1/2) * 2^(-E/2)= C - (R>>1)求出令到相對(duì)誤差最小的C值就可以了－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－上面的推導(dǎo)過程只是我個(gè)人的理解，并未得到證實(shí)。而Chris Lomont則在他的論文中詳細(xì)討論了最后那個(gè)方程的解法，并嘗試在實(shí)際的機(jī)器上尋找最佳的常數(shù)C。有興趣的朋友可以在文末找到他的論文的鏈接。
所以，所謂的Magic Number，并不是從N元宇宙的某個(gè)星系由于時(shí)空扭曲而掉到地球上的，而是幾百年前就有的數(shù)學(xué)理論。只要熟悉NR和泰勒級(jí)數(shù)，你我同樣有能力作出類似的優(yōu)化。
在http://GameDev.net上有人做過測(cè)試，該函數(shù)的相對(duì)誤差約為0.177585%，速度比C標(biāo)準(zhǔn)庫的sqrt提高超過20%。如果增加一次迭代過程，相對(duì)誤差可以降低到e-004 的級(jí)數(shù)，但速度也會(huì)降到和sqrt差不多。據(jù)說在DOOM3中，Carmack通過查找表進(jìn)一步優(yōu)化了該算法，精度近乎完美，而且速度也比原版提高了一截（正在努力弄源碼，誰有發(fā)我一份）。
值得注意的是，在Chris Lomont的演算中，理論上最優(yōu)秀的常數(shù)（精度最高）是0x5f37642f，并且在實(shí)際測(cè)試中，如果只使用一次迭代的話，其效果也是最好的。但奇怪的是，經(jīng)過兩次NR后，在該常數(shù)下解的精度將降低得非常厲害（天知道是怎么回事?。?。經(jīng)過實(shí)際的測(cè)試，Chris Lomont認(rèn)為，最優(yōu)秀的常數(shù)是0x5f375a86。如果換成64位的double版本的話，算法還是一樣的，而最優(yōu)常數(shù)則為0x5fe6ec85e7de30da（又一個(gè)令人冒汗的Magic Number - -b）。
這個(gè)算法依賴于浮點(diǎn)數(shù)的內(nèi)部表示和字節(jié)順序，所以是不具移植性的。如果放到Mac上跑就會(huì)掛掉。如果想具備可移植性，還是乖乖用sqrt好了。但算法思想是通用的。大家可以嘗試推算一下相應(yīng)的平方根算法。
下面給出Carmack在QUAKE3中使用的平方根算法。Carmack已經(jīng)將QUAKE3的所有源代碼捐給開源了，所以大家可以放心使用，不用擔(dān)心會(huì)受到律師信。－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－//// Carmack在QUAKE3中使用的計(jì)算平方根的函數(shù)//float CarmSqrt(float x){union{int intPart;float floatPart;} convertor;union{int intPart;float floatPart;} convertor2;convertor.floatPart = x;convertor2.floatPart = x;convertor.intPart = 0x1FBCF800 + (convertor.intPart >> 1);convertor2.intPart = 0x5f3759df - (convertor2.intPart >> 1);return 0.5f*(convertor.floatPart + (x * convertor2.floatPart));}－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
故事還沒完，普渡大學(xué)的數(shù)學(xué)家Chris Lomont看了以后也覺得很有趣，決定要研究一下卡馬克弄出來的這個(gè)猜測(cè)值有什么奧秘。Lomont也是個(gè)牛人，在精心研究之后從理論上也推導(dǎo)出一個(gè)最佳猜測(cè)值，和卡馬克的數(shù)字非常接近, 0x5f37642f。Lomont計(jì)算出結(jié)果以后非常滿意，于是拿自己計(jì)算出的起始值和卡馬克的神秘?cái)?shù)字做比賽，看看誰的數(shù)字能夠更快更精確的求得平方根。結(jié)果居然還是卡馬克贏了... 
最后Lomont怒了，采用暴力方法一個(gè)數(shù)字一個(gè)數(shù)字試過來，終于找到一個(gè)比卡馬克數(shù)字要好上那么一丁點(diǎn)的數(shù)字，雖然實(shí)際上這兩個(gè)數(shù)字所產(chǎn)生的結(jié)果非常近似，這個(gè)暴力得出的數(shù)字是0x5f375a86。
數(shù)學(xué)的重要性不言而喻了！
1. C++標(biāo)準(zhǔn)模板庫從入門到精通 
http://edu.csdn.net/course/detail/3324
2.跟老菜鳥學(xué)C++
http://edu.csdn.net/course/detail/29013. 跟老菜鳥學(xué)python
http://edu.csdn.net/course/detail/25924. 在VC2015里學(xué)會(huì)使用tinyxml庫
http://edu.csdn.net/course/detail/25905. 在Windows下SVN的版本管理與實(shí)戰(zhàn) 
 http://edu.csdn.net/course/detail/2579
6.Visual Studio 2015開發(fā)C++程序的基本使用 
http://edu.csdn.net/course/detail/2570
7.在VC2015里使用PRotobuf協(xié)議
http://edu.csdn.net/course/detail/25828.在VC2015里學(xué)會(huì)使用MySQL數(shù)據(jù)庫
http://edu.csdn.net/course/detail/2672