【转】代码优化之-优化除法（牛顿迭代法）

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