分享
 
 
 

代码优化-之-优化除法

王朝other·作者佚名  2006-08-31
窄屏简体版  字體: |||超大  

<<代码优化-之-优化除法>> 作者:HouSisong@263.net

(说明:文章中的很多数据可能在新的CPU或不同的CPU或不同的系统环境下有不同的结果,可能不能面面俱到)

x86系列的CPU对于位运算、加、减等基本指令都能在1个CPU周期内完成(现在的CPU还能乱序执行,从而使指令的平均CPU周期更小);现在的CPU,做乘法也是很快的(1个CPU周期左右,或者是需要两\三个周期,但每个周期能启动一个新的乘指令),但作为基本指令的除法却超出很多人的预料,它是一条很慢的操作,整数和浮点的除法都慢;我测试的英特尔P5赛扬CPU浮点数的除法差不多是37个CPU周期,整数的除法是80个CPU周期,AMD2200+浮点数的除法差不多是21个CPU周期,整数的除法是40个CPU周期。(改变FPU运算精度对于除法无效) (SSE指令集的低路单精度数除法指令DIVPS 18个CPU周期,四路单精度数除法指令DIVSS 36个CPU周期) (x86求余运算和除法运算是用同一条CPU指令实现的; 据说,很多CPU的整数除法都是用数学协处理器的浮点除法器完成的;有一个推论就是,浮点除法和整数除法不能并行执行)

本文将给出一些除法的优化方法或替代算法 (警告:某些替代算法并不能保证完全等价!)

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 s,uint32 z) //return u/z

{

uint32 x=(uint32)(s>>32);

uint32 y=(uint32)s;

//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;//利用a为奇数,y是a的倍数可以推导出a'的存在性

那么: 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 //^表示指数运算符

证毕.

程序可以用该方法来实现除法,并按照自己的精度要求来决定迭代次数;

(初始值的选取不再讨论)

附录: 用牛顿迭代法来实现开方运算

//开方运算可以表示为 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阶收敛公式 //方法同上 证明过程略

11. 用快速乘法器实现快速除法

这是存在的一种快速除法算法,当前最好用硬件实现,请查阅相关资料

 
 
 
免责声明:本文为网络用户发布,其观点仅代表作者个人观点,与本站无关,本站仅提供信息存储服务。文中陈述内容未经本站证实,其真实性、完整性、及时性本站不作任何保证或承诺,请读者仅作参考,并请自行核实相关内容。
2023年上半年GDP全球前十五强
 百态   2023-10-24
美众议院议长启动对拜登的弹劾调查
 百态   2023-09-13
上海、济南、武汉等多地出现不明坠落物
 探索   2023-09-06
印度或要将国名改为“巴拉特”
 百态   2023-09-06
男子为女友送行,买票不登机被捕
 百态   2023-08-20
手机地震预警功能怎么开?
 干货   2023-08-06
女子4年卖2套房花700多万做美容:不但没变美脸,面部还出现变形
 百态   2023-08-04
住户一楼被水淹 还冲来8头猪
 百态   2023-07-31
女子体内爬出大量瓜子状活虫
 百态   2023-07-25
地球连续35年收到神秘规律性信号,网友:不要回答!
 探索   2023-07-21
全球镓价格本周大涨27%
 探索   2023-07-09
钱都流向了那些不缺钱的人,苦都留给了能吃苦的人
 探索   2023-07-02
倩女手游刀客魅者强控制(强混乱强眩晕强睡眠)和对应控制抗性的关系
 百态   2020-08-20
美国5月9日最新疫情:美国确诊人数突破131万
 百态   2020-05-09
荷兰政府宣布将集体辞职
 干货   2020-04-30
倩女幽魂手游师徒任务情义春秋猜成语答案逍遥观:鹏程万里
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案神机营:射石饮羽
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案昆仑山:拔刀相助
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案天工阁:鬼斧神工
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案丝路古道:单枪匹马
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案镇郊荒野:与虎谋皮
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案镇郊荒野:李代桃僵
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案镇郊荒野:指鹿为马
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案金陵:小鸟依人
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案金陵:千金买邻
 干货   2019-11-12
 
推荐阅读
 
 
 
>>返回首頁<<
 
靜靜地坐在廢墟上,四周的荒凉一望無際,忽然覺得,淒涼也很美
© 2005- 王朝網路 版權所有