代碼最佳化－－最佳化除法

轉自：http://blog.csdn.net/housisong/article/details/1116423

 

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[i][j]=i/j; 
           //對於除零的情況需要根據你的應用來處理
    或者：uint8  z=table[x<<8+y];  //其中table[i]=(i>>8)/(i&(1<<8-1));

    比如   uint8  x;
         uint8  z=x/17;
    改成   uint8  z=table[x];    //其中table[i]=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等其它更精確的公式(請推匯出誤差項已確定是否精度足夠)

       比如: x/255 寫成 (x+(x>>8)+1)>>8 ,當0<=x<65535的時候都完全等價    

  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階收斂公式  //方法同上 證明過程略