很牛的牛頓迭代法

在MIT公開課《計算機科學與編程導論》的第五講中，講到編寫求解平方根的函數sqrt時，提到了牛頓迭代法。今天仔細一查，發現這是一個用途很廣、很牛的計算方法。

首先，考慮如何編寫一個開平方根的函數sqrt(float num, float e)。參數num是要求開平方根的實數，參數e是計算結果可以達到多大誤差。這是一個無法得到精確解，只能求出近似解的問題。該如何編寫呢？

1. 傳統的二分法

我們可以先猜測一個值作為解，看這個值是否在誤差范圍內。如果沒有達到誤差要求，就構造一個更好的猜測值，繼續迭代。猜測值可以簡單地初始化為num/2，但怎樣在下一輪迭代前構造一個更好的猜測值呢？我們不妨參照二分查找算法的思路，解的初始范圍是[0, num]，用二分法逐漸縮小范圍。

     private static float sqrt(float num, float e) {

              

               float  low = 0F;

               float  high = num;

               float  guess, e0;

               int  count = 0;

              

               do  {

                     guess = (low + high) / 2;

                      if  (guess*guess > num) {

                           high = guess;

                           e0 = guess*guess - num;

                     }  else  {

                           low = guess;

                           e0 = num - guess*guess;

                     }

                     

                     count++;

                     System. out .printf( "Try %f, e: %f\n" , guess, e0);

              }  while  (e0 > e);

              System. out .printf( "Try %d times, result: %f\n" , count, guess);

              

               return  guess;

       }

在區間[low, high]內取中點(low+high)/2作為猜測值。如果guess*guess大於num，說明猜測值偏大，則在區間[low, guess]進行下一輪迭代，否則在區間[guess, high]中繼續。當誤差值e0小於能夠接受的誤差值e時停止迭代，返回結果。

取num=2, e=0.01進行測試，運行結果如下：

Try 1.000000, e: 1.000000

Try 1.500000, e: 0.250000

Try 1.250000, e: 0.437500

Try 1.375000, e: 0.109375

Try 1.437500, e: 0.066406

Try 1.406250, e: 0.022461

Try 1.421875, e: 0.021729

Try 1.414063, e: 0.000427

Try 8 times, result: 1.414063

可見嘗試了八次才達到0.01的誤差。

2. 神奇的牛頓法

仔細思考一下就能發現，我們需要解決的問題可以簡單化理解。

從函數意義上理解：我們是要求函數f(x) = x²，使f(x) = num的近似解，即x² - num = 0的 近似解。

從幾何意義上理解：我們是要求拋物線g(x) = x² - num與x軸交點（g(x) = 0）最接近的點。

我們假設g(x0)=0，即x0是正解，那么我們要做的就是讓近似解x不斷逼近x0，這是函數導數的定義：

可以由此得到

從幾何圖形上看，因為導數是切線，通過不斷迭代，導數與x軸的交點會不斷逼近x0。

3. 牛頓法的實現與測試

     public  static  void  main(String[] args) {

               float  num = 2;

               float  e = 0.01F;

              sqrt(num, e);

              sqrtNewton(num, e);

              

              num = 2;

              e = 0.0001F;

              sqrt(num, e);

              sqrtNewton(num, e);

              

              num = 2;

              e = 0.00001F;

              sqrt(num, e);

              sqrtNewton(num, e);

       }

     private  static  float  sqrtNewton( float  num,  float  e) {

              

               float  guess = num / 2;

               float  e0;

               int  count = 0;

              

               do  {

                     guess = (guess + num / guess) / 2;

                     e0 = guess*guess - num;

                     

                     count++;

                     System. out .printf( "Try %f, e: %f\n" , guess, e0);

              }  while  (e0 > e);

              System. out .printf( "Try %d times, result: %f\n" , count, guess);

              

               return  guess;

       }

與二分法的對比測試結果：

Try 1.000000, e: 1.000000

Try 1.500000, e: 0.250000

Try 1.250000, e: 0.437500

Try 1.375000, e: 0.109375

Try 1.437500, e: 0.066406

Try 1.406250, e: 0.022461

Try 1.421875, e: 0.021729

Try 1.414063, e: 0.000427

Try 8 times, result: 1.414063

Try 1.500000, e: 0.250000

Try 1.416667, e: 0.006945

Try 2 times, result: 1.416667

Try 1.000000, e: 1.000000

Try 1.500000, e: 0.250000

Try 1.250000, e: 0.437500

Try 1.375000, e: 0.109375

Try 1.437500, e: 0.066406

Try 1.406250, e: 0.022461

Try 1.421875, e: 0.021729

Try 1.414063, e: 0.000427

Try 1.417969, e: 0.010635

Try 1.416016, e: 0.005100

Try 1.415039, e: 0.002336

Try 1.414551, e: 0.000954

Try 1.414307, e: 0.000263

Try 1.414185, e: 0.000082

Try 14 times, result: 1.414185

Try 1.500000, e: 0.250000

Try 1.416667, e: 0.006945

Try 1.414216, e: 0.000006

Try 3 times, result: 1.414216

Try 1.000000, e: 1.000000

Try 1.500000, e: 0.250000

Try 1.250000, e: 0.437500

Try 1.375000, e: 0.109375

Try 1.437500, e: 0.066406

Try 1.406250, e: 0.022461

Try 1.421875, e: 0.021729

Try 1.414063, e: 0.000427

Try 1.417969, e: 0.010635

Try 1.416016, e: 0.005100

Try 1.415039, e: 0.002336

Try 1.414551, e: 0.000954

Try 1.414307, e: 0.000263

Try 1.414185, e: 0.000082

Try 1.414246, e: 0.000091

Try 1.414215, e: 0.000004

Try 16 times, result: 1.414215

Try 1.500000, e: 0.250000

Try 1.416667, e: 0.006945

Try 1.414216, e: 0.000006

Try 3 times, result: 1.414216

可以看到隨着對誤差要求的更加精確，二分法的效率很低下，而牛頓法的確非常高效。

可在兩三次內得到結果。

如果搞不清牛頓法的具體原理，可能就要像我一樣復習下數學知識了。極限、導數、泰勒展開式、單變量微分等。

4. 更快的方法

在Quake源碼中有段求sqrt的方法，大概思路是只進行一次牛頓迭代，得到能夠接受誤差范圍內的結果。

因此肯定是更快的。

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

參考文章

Quake3源碼中的sqrt  http://www.matrix67.com/blog/archives/362

牛頓迭代方程的近似解  http://blueve.me/archives/369