【轉】代碼優化之-優化除法（牛頓迭代法）

tag:代碼優化,除法,牛頓迭代,減法代替除法,除法優化

   說明：文章中的不少數據可能在不一樣的CPU或不一樣的系統環境下有不一樣的結果，數據僅供參考
   x86系列的CPU對於位運算、加、減等基本指令都能在1個CPU週期內完成(如今的CPU還能亂序執行，從而使指令的平均CPU週期更小)；如今的CPU,作乘法也是很快的(須要幾個CPU週期，每一個週期可能啓動一個新的乘指令(x87))，但做爲基本指令的除法卻超出不少人的預料，它是一條很慢的操做，整數和浮點的除法都慢；我測試的英特爾P5賽揚CPU浮點數的除法差很少是37個CPU週期，整數的除法是80個CPU週期，AMD2200+浮點數的除法差很少是21個CPU週期，整數的除法是40個CPU週期。(改變FPU運算精度對於除法無效) (SSE指令集的低路單精度數除法指令DIVPS 18個CPU週期,四路單精度數除法指令DIVSS 36個CPU週期)   (x86求餘運算和除法運算是用同一條CPU指令實現的； 聽說，不少CPU的整數除法都是用數學協處理器的浮點除法器完成的；有一個推論就是，浮點除法和整數除法不能並行執行. (ps:intel的p4 imul指令可能有14週期(或15-18)的延遲才能獲得結果)

本文將給出一些除法的優化方法或替代算法 (警告:某些替代算法並不能保證徹底等價!)

1.儘可能少用除法
   好比： if (x/y>z) ...
   改爲： if ( ((y>0)&&(x>y*z))||((y<0)&&(x<y*z)) ) ...

   (少用求餘)
   好比： ++index; if (index>=count) index=index % count; //assert(index<count*2);
   改爲： ++index; if (index>=count) index=index - count;

2.用減法代替除法
   若是知道被除數是除數的很小的倍數，那麼能夠用減法來代替除法
   好比： uint32 x=200;
       uint32 y=70;
           uint32 z=x/y;
   改爲： uint z=0;
           while (x>=y)
           {
              x-=y; ++z;
           }

   一個用減法和移位完成的除法 (若是你沒有除法指令可用
     uint32 div(uint64 u,uint32 z) //return u/z   
     {
         uint32 x=(uint32)(u>>32);
         uint32 y=(uint32)u;
         //y保存商 x保存餘數

         for (int i=0;i<32;++i)
         {
             uint32 t=((int32)x) >> 31; 
             x=(x<<1)|(y>>31);
             y=y<<1;
             if ((x|t)>=z)
             {
                 x-=z;
                 ++y;
             }
         }
         return y;
     }
     (該函數通過了測試；z==0須要本身處理；對於有符號除法，能夠用取絕對值的方法(固然也不是輕鬆就能
寫出徹底等價的有符號除法的； 若是不需s的64bit長度，僅須要32bit,那麼能夠化簡這個函數，但改進很少)
   
3.用移位代替除法 (不少編譯器能自動作好這個優化)
   要求除數是2的次方的常量； (同理：對於某些應用，能夠優先選取這樣的數來作除數)
   好比： uint32 x=213432575;
       uint32 y=x/8;
   改爲： y=x>>3;

   對於有符號的整數；
   好比： int32 x=213432575;
       int32 y=x/8;
   改爲： if (x>=0) y=x>>3; 
           else y=(x+(1<<3-1))>>3;

4.合併除法   (替代方法不等價，不少編譯器都不會幫你作這種優化)
   適用於不能用其它方法避免除法的時候；
   好比： double x=a/b/c;
   改爲： double x=a/(b*c);

   好比： double x=a/b+c/b;
   改爲： double x=(a+c)/b;

   好比： double x=a/b;
    double y=c/d;
    double z=e/f;
   改爲： double tmp=1.0/(b*d*f);
           double x=a*tmp*d*f;
           double y=c*tmp*b*f;
           double z=e*tmp*b*d;    
   
5.把除法佔用的時間充分利用起來 
   CPU在作除法的時候，能夠不用等待該結果(也就是後面插入的指令不使用該除法結果)，而插入多條簡單整數
指令(不包含整數除法，並且結果不能是一個全局或外部變量等)，把除法佔用的時間節約出來； 
   (當除法不可避免的時候，這個方法頗有用)

6.用查表的方法代替除法
    （適用於除數和被除數的可能的取值範圍較小的狀況，不然空間消耗太大）
    好比   uint8 x; uint8 y; 
       uint8 z=x/y;
    改爲   uint8 z=table[x][y];    //其中table是預先計算好的表，table[j]=i/j; 
           //對於除零的狀況須要根據你的應用來處理
    或者：uint8 z=table[x<<8+y]; //其中table=(i>>8)/(i&(1<<8-1));

    好比   uint8 x; 
       uint8 z=x/17;
    改爲   uint8 z=table[x];    //其中table=i/17;

7.用乘法代替除法
     (替代方法不等價，不少編譯器都不會幫你作這種優化)
     好比： double x=y/21.3432575;
     改爲： double x=y*(1/21.3432575); //若是編譯器不能優化掉(1/21.3432575),請預先計算出該結果

   對於整數，可使用與定點數類似的方法來處理倒數
   (該替代方法不等價)
   好比： uint32 x,y; x=y/213;
   改爲： x=y*((1<<16)/213) >> 16;
   某些應用中y*((1<<16)/213)可能會超出值域，這時候能夠考慮使用int64來擴大值域
           uint32 x=((uint64)y)*((1<<31)/213) >> 31;
      也可使用浮點數來擴大值域
           uint32 x=(uint32)(y*(1.0/213)); (警告: 浮點數強制類型轉換到整數在不少高級語言裏都是
一條很慢的操做，在下幾篇文章中將給出其優化的方法)

   對於這種方法，某些除法是有與之徹底等價的優化方法的：
     無符號數除以3： uint32 x,y;   x=y/3; 
     推理:                       y               y           y               (1<<33)+1   y      
             x=y/3 => x=[-] => x=[- + ---------] => x=[--------- * -------] // []表示取整 
                                   3               3 3*(1<<33)            3       (1<<33) 
                                        y      1              
              (可證實: 0 <= --------- < - )
                              3*(1<<33)   3                                           
             即： x=(uint64(y)*M)>>33; 其中魔法數M=((1<<33)+1)/3=2863311531=0xAAAAAAAB;

     無符號數除以5： uint32 x,y;   x=y/5; (y<(1<<31))
     推理:                       y               y      3*y            (1<<33)+3     y      
             x=y/5 => x=[-] => x=[- + ---------] => x=[--------- * -------] 
                                   5               5   5*(1<<33)             5       (1<<33) 
                                      3*y      1              
              (可證實: 0 <= --------- < - )
                               5*(1<<33)   5                                                           
             即： x=(uint64(y)*M)>>33; 其中魔法數M=((1<<33)+3)/5=1717986919=0x66666667;

     無符號數除以7： uint32 x,y;   x=y/7; 
     推理 ：略                                                         
          即：x=((uint64(y)*M)>>33+y)>>3; 其中魔法數M=((1<<35)+3)/7-(1<<32)=613566757=0x24924925;

     對於這種徹底等價的替代,還有其餘替代公式再也不討論，對於有符號除法的替代算法沒有討論，不少數都有徹底等價的替代算法，那些魔法數也是有規律可循的；有「進取心」的編譯器應該幫助用戶處理掉這個優化(vc6中就已經見到! 偷懶的辦法是直接看vc6生成的彙編碼；

8.對於已知被除數是除數的整數倍數的除法，可以獲得替代算法；或改進的算法；
    這裏只討論除數是常數的狀況；
    好比對於(32位除法)：x=y/a; //已知y是a的倍數,並假設a是奇數
    (若是a是偶數,先轉化爲a=a0*(1<<k); y/a==(y/a0)>>k;a0爲奇數)
    求得a',使 (a'*a) % (1<<32) = 1；
    那麼: x=y/a => x=(y/a)*((a*a') %(1<<32)) => x=(y*a') % (1<<32); //這裏並不須要實際作一個求餘指令   
    (該算法能夠同時支持有符號和無符號除法)
    好比   uint32 y; 
       uint32 x=y/7;   //已知y是7的倍數
    改爲   uint32 x=(uint32)(y*M);   //其中M=(5*(1<<32)+1)/7 

9.近似計算除法 (該替代方法不等價)
     優化除數255的運算(257也能夠,或者1023,1025等等)(1026也能夠,推導出的公式略有不一樣)
     好比顏色處理中 ： uint8 color=colorsum/255;
     改爲：             uint8 color=colorsum/256; 也就是 color=colorsum>>8;
     這個偏差在顏色處理中不少時候都是能夠接受的
     若是要減少偏差能夠改爲:     uint8 color=(colorsum+(colorsum>>8))>>8;　
       推導: x/255=(x+x/255)/(255+1)=(x+A)>>8; A=x/255;
       把A改爲A=x>>8 (引入小的偏差);帶入公式就獲得了:　x/255約等於(x+(x>>8))>>8的公式
       同理能夠有x/255約等於(x+(x>>8)+(x>>16))>>8等其它更精確的公式(請推導出偏差項已肯定是否精度足夠)
     
10. 牛頓迭代法實現除法
     (不少CPU的內部除法指令就是用該方法或相似的方法實現的) 
     假設有一個函數y=f(x); 求0=f(x)時，x的值;(這裏不討論有多個解的狀況或逃逸或陷入死循環或陷入混沌狀態的問題)
     
                                 (參考圖片)
     求函數的導函數爲 y=f'(x);   //能夠這樣來理解這個函數的意思：x點處函數y=f(x)的斜率；
     a.取一個合適的x初始值x0; 並獲得y0=f(x0);
     b.過(x0,y0)做一條斜率爲f'(x0)的直線，與x軸交於x1;
     c.而後用獲得的x1做爲初始值，進行迭代；
     當進行足夠屢次的迭代之後，認爲x1將會很是接近於方程0=f(x)的解，這就是牛頓迭代；
     把上面的過程化簡,獲得牛頓迭代公式: x(n+1)=x(n)-y(x(n))/y'(x(n))
     
     這裏給出利用牛頓迭代公式求倒數的方法； (用倒數獲得除法: y = x/a = x* (1/a) )
        求1/a,   
        令a=1/x; 有方程 y=a-1/x;
        求導得y'=1/x^2;
        代入牛頓迭代方程 x(n+1)=x(n)-y(x(n))/y'(x(n));
        有迭代式 x_next=(2-a*x)*x; //可證實：該公式爲2階收斂公式； 也就是說計算出的解的有效精度成倍增加；
       
     證實收斂性:令x=(1/a)+dx; //dx爲一個很小的量
     則有x_next-(1/a)=(2-a*(1/a+dx))*(1/a+dx)-1/a
                     =(-a)*dx^2 //^表示指數運算符
      證畢.
    
    程序能夠用該方法來實現除法，並按照本身的精度要求來決定迭代次數；
    (對於迭代初始值，可使用查表法來獲得,或者利用IEEE浮點格式獲得一個近似計算的表達式;在SSE指令集中有一條RCPPS(RCPSS)指令也能夠用來求倒數的近似值；有了初始值就能夠利用上面的迭代公式獲得更精確的結果)

    附錄: 用牛頓迭代法來實現開方運算   
       //開方運算能夠表示爲 y=x^0.5=1/(1/x^0.5)； 先求1/x^0.5
       求1/a^0.5,   
       令a=1/x^2; 有方程y=a-1/x^2;
       求導得y'=2/x^3;
       代入牛頓方程 x(n+1)=x(n)-y(x(n))/y'(x(n));
       有迭代式 x_next=(3-a*x*x)*x*0.5; //可證實：該公式爲2階收斂公式 //方法同上，證實過程略！

 

 

 參考來源：當程序愛上游戲的博客->代碼優化之-優化除法，http://micheal19840929.iteye.com/blog/592085