代码优化-之-优化浮点数取整

代码优化一之一优化浮点数取整HouSisong@GM2007.05.19tag:浮点数转换为整数fpu,sse,sse2,读缓冲区优化，代码优化ftol，取整,f2l,ftoi,f2i,floattoint摘要：本文首先给出一个浮点数取整的需求，并使用默认的取整方式，然后通过尝试各种方法来优化它的速度；最终的浮点数取整实现速度甚至达到了初始代码的5倍（是vc6代码的18倍）！（注意：文章中的测试结果在不同的CPU和系统环境下可能有不同的结果，数据仅作参考）（2007.06.08更新：补充SSE3新增的FPU取整指令fisttp的说明）（2007.06.04更新：一些修正、补充double取整、补充FPU的RC场说明）正文：为了便于讨论，这里代码使用C++,涉及到汇编优化的时候假定为x86平台；使用的编译器为vc2005;测试使用的CPU为AMD64x24200+,测试时使用的单线程执行；为了代码的可读性，没有加入异常处理代码；A:需要优化的原始代码（使用了大量的浮点数到整数的转换）#include<stdio.h>#include<stdlib.h>#include<time.h>volatilelongtestResult;//使用一个全局域的volatile变量以避免编译器把需要测试的代码优化掉constlongtestDataCount=10000000;constlongtestCount=20;floatfSrc[testDataCount];#defineasmasmvoidftol_test_0（）{longtmp=0;for(longi=0;i<testDataCount;++i)tmp+=(long)fSrc[i];//需要优化的浮点数取整testResult=tmp;}intmain(){//intifor(longi=O;i<testDataCount;++i)fSrc[i]=(float)(rand()*(1.0/RAND_MAX)*(rand()—(RAND_MAX〉〉1))*rand()*(1.0/RAND_MAX));//testdoublestart0=(double)clock();for(longc=O;c<testCount;++c)ftol_test_0();startO=((double)clock()-startO)*(1.0/CLOCKS_PER_SEC);//outprintf("Result=%udSeconds=%8.5f",testResult,start0)；return0;}//////////////////////////////////////////////////////////////////////////////////速度测试：//===============================================================//ftol_test_01.047秒(VC6编译3.64秒)：//(使用vc2005的SSE编译选项“/arch:SSE”0.437秒)////////////////////////////////////////////////////////////////////////////////一般编译器生成的浮点数转换为整数的指令序列都比预想的速度慢很多,它的性能代价很容易被人忽略；在VC6编译器下上面的代码需要运行3.64秒，代码先修改FPU的取整模式（RC场），完成取整后在恢复RC场；VC2006生成的代码在CPU支持SSE的时候会调用使用cvttsd2si指令实现的版本，从而加快了取整的速度，达到了1.047秒，快了很多！让我们来尝试继续优化这个含有大量取整操作的函数ftol_test_0；B：最容易想到的就是用浮点协处理器（FPU）（也可以称作x87）来优化取整将设置FPU取整方式和恢复FPU的取整方式的代码放到循环体外面从而加快了速度voidftol_test_fpu（）{asm{//设置FPU的取整方式为了直接使用fistp浮点指令FNSTCWRC_Old//保存协处理器控制字，用来恢复FNSTCWRC_Edit//保存协处理器控制字，用来修改FWAITORRC_Edit,OxOFOO//改为RC=11使FPU向零取整FLDCWRC_Edit//载入协处理器控制字，RC场已经修改movecx,testDataCountxoreax,eaxtestecx,ecxjleEndLoopleaedx,[fSrc+ecx*4]negecxStartLoop:flddwordptr[edx+ecx*4]fistpisrcaddeax,isrcincecxjnzStartLoopEndLoop:movtestResult,eax;//恢复FPU的取整方式FWAITFLDCWRC_Old}}//RC场占用第11、10bit位用于控制舍入方向//RC=00向最近（或偶数）舍入RC=01向下（负无穷大）舍入//RC=10向上（正无穷大）舍入RC=11向零舍入//提示：一般的编程语言环境中，RC场都会设置为一个默认值（一般为RC=00）,//语言就可能利用这一点做快速的取整（比如Delphi中的round函数），但某些引入的//第三方库或代码可能会修改该默认值，从而造成以前运行正确的程序出现异常情况//////////////////////////////////////////////////////////////////////////////////速度测试：//===============================================================//ftol_test_fpu0.407秒////////////////////////////////////////////////////////////////////////////////SSE3增加了一条FPU取整指令fisttp,和fistp指令功能几乎相同（我的电脑上经过测试速度也相同）但默认向0取整，和RC场设置无关，所以使用fisttp的代码就可以不管RC场了，有利于简化代码和优化性能；C:利用浮点数的编码格式来''手工〃处理浮点数到整数的转换（利用了IEEE浮点编码格式）inlinelongftolieee(floatf){longa=*(long*)(&f);unsignedlongmantissa=(a&((1<<23)-1))|(1<<23);//不支持非规格化浮点数longexponent=((a&0x7fffffff)〉〉23);longr=(mantissa<<8)>>(31+127-exponent);longsign=(a〉〉31);return((r入(sign))-sign)&〜((exponent-127)〉〉31);}voidftol_test_ieee(){longtmp=0;for(longi=0;i<testDataCount;++i)tmp+=_ftol_ieee(fSrc[i]);testResult=tmp;}//////////////////////////////////////////////////////////////////////////////////速度测试：//===============================================================//ftol_test_ieee0.828秒////////////////////////////////////////////////////////////////////////////////手工实现居然超过了VC2005的SSE实现（主要是VC2005的实现函数调用开销太大）；如果能够允许存在误差的话，还有一个快速的取整算法（注意，该函数的结果和标准不完全相同）：//ftol_test_ieee_M0.438秒D:对于Double到整数的转换有一个超强的算法(利用了IEEE浮点编码格式)inlinelong_ftol_ieee_MagicNumber(doublex){staticconstdoublemagic=6755399441055744.0;//(1<<51)|(1<<52)doubletmp=x;tmp+=(x>0)?-0.499999999999:+0.499999999999;//如果需要4舍5入取整就去掉这一行tmp+=magic;return*(long*)&tmp;}voidftol_test_ieee_MagicNumber(){longtmp=0;for(longi=0;i<testDataCount;++i)tmp+=_ftol_ieee_MagicNumber(fSrc[i]);testResult=tmp;}（警告：该算法要求FPU的计算精度为高精度模式，某些程序可能为了速度而将FPU改成了低精度模式，比如在D3D中会默认调整该设置）//////////////////////////////////////////////////////////////////////////////////速度测试：//===============================================================//ftol_test_ieee_MagicNumber1.813秒////////////////////////////////////////////////////////////////////////////////如果需要4舍5入取整，速度就能快出很多，降低到0.407秒(ftol_test_ieee,ftol_test_ieee_MagicNumber的实现主要参考了：云风的《_ftol的优化》：/2005/12/_ftol_opt.html和/cgi-bin/fcarticles.cgi?show=64008这里有改动)E：借鉴vc2005的SSE实现使用cvttss2si指令voidftol_test_sse(){asm{movecx,testDataCountxoreax,eaxtestecx,ecxjleEndLoopleaedx,[fSrc+ecx*4]negecxStartLoop:cvttss2siebx,dwordptr[edx+ecx*4]addeax,ebx

incecxjnzStartLoopincecxjnzStartLoopEndLoop:movtestResult,eax;//////////////////////////////////////////////////////////////////////////////////速度测试：//===============================================================//ftol_test_sse0.422秒////////////////////////////////////////////////////////////////////////////////F:cvttss2si是一个单指令单数据流的指令，我们可以使用它的单指令多数据流的版本:cvttps2dq指令；它能同时将4个float取整！StartLoop://—次循环处理16个floatcvttps2dqxmm2,xmmwordptr[edx+ecx*4]cvttps2dqxmm3,xmmwordptr[edx+ecx*4+16]cvttps2dqxmm4,xmmwordptr[edx+ecx*4+16*2]cvttps2dqxmm5,xmmwordptr[edx+ecx*4+16*3]padddxmm2,xmm3padddxmm4,xmm5addecx,16padddxmm0,xmm2padddxmm1,xmm4jnzStartLoopEndLoop:padddxmm0,xmm1movapsxmm1,xmm0movhlpsxmm1,xmm0padddxmm0,xmm1movapsxmm2,xmm0shufpsxmm2,xmm0,1padddxmm0,xmm2movdeax,xmm0movresult,eax}returnresult;}voidftol_test_sse_expand16(){longtmp=0;for(longi=0;i<testDataCount;i+=2000){tmp+=ftol_sse_expand16(&fSrc[i],2000);//2000=16*125}//todo:因为testDataCount是2000的倍数，所以这里不用处理边界了testResult=tmp;}//////////////////////////////////////////////////////////////////////////////////速度测试：//===============================================================//ftol_test_sse_expand160.281秒////////////////////////////////////////////////////////////////////////////////G:由于函数需要读取大量的数据来处理，所以可以考虑优化读缓冲区（也可以考虑使用显式预读指令）

pxorxmm0,xmm0pxorxmm1,xmm1movecx,count16negecxStartLoop:[edx+ecx*4][edx+ecx*4+16][edx+ecx*4+16*2][edx+ecx*4+16*3][edx+ecx*4][edx+ecx*4+16][edx+ecx*4+16*2][edx+ecx*4+16*3]cvttps2dqxmm3,xmmwordptrcvttps2dqxmm4,xmmwordptrcvttps2dqxmm5,xmmwordptrpadddxmm2,xmm3padddxmm4,xmm5addecx,16padddxmm0,xmm2padddxmm1,xmm4jnzStartLoopEndLoop:padddxmm0,xmm1movapsxmm1,xmm0movhlpsxmm1,xmm0padddxmm0,xmm1movapsxmm2,xmm0shufpsxmm2,xmm0,1padddxmm0,xmm2movdeax,xmm0movresult,eax}returnresult;}voidftol_test_sse_expand16_prefetch(){longtmp=0;//////////////////////////////////////////////////////////////////////////////////速度测试：//===============================================================//ftol_test_sse_expand16_prefetch0.219秒////////////////////////////////////////////////////////////////////////////////H:补充Double的取整，完整测试源代码#include<stdio.h>#include<stdlib.h>#include<time.h>volatilelongtestResult;//使用一个全局域的volatile变量以避免编译器把需要测试的代码优化掉constlongtestDataCount=10000000;constlongtestCount=20;doublefSrc[testDataCount];#defineasmasmvoiddftol_test_0(){longtmp=0;for(longi=0;i<testDataCount;++i){tmp+=(long)fSrc[i];//需要优化的浮点取整}testResult=tmp;}voiddftol_test_fpu(){longisrc;asm//设置FPU的取整方式为了直接使用fistp浮点指令{FNSTCWRC_Old//保存协处理器控制字，用来恢复FNSTCWRC_Edit//保存协处理器控制字，用来修改FWAITORRC_Edit,OxOFOO//改为RC=11使FPU向零取整FLDCWRCEdit//载入协处理器控制字，RC场已经修改//}//asm//{movecx,testDataCountxoreax,eaxtestecx,ecxjleEndLoopleaedx,[fSrc+ecx*8]negecxStartLoop:fldqwordptr[edx+ecx*8]fistpisrcaddeax,isrcincecxjnzStartLoop

EndLoop:movtestResult,eax;//}//asm//恢复FPU的取整方式//{FWAITFLDCWRC_Old}}inlinelongdftol_ieee_MagicNumber(doublex)(1<<5//如果需要(1<<5//如果需要4staticconstdoublemagic=6755399441055744.0;//1)|(1<<52)doubletmp=x;tmp+=(x>0)?-0.499999999999:+0.499999999999;舍5入取整就去掉这一行tmp+=magic;return*(long*)&tmp;}voiddftol_test_ieee_MagicNumber(){longtmp=0;for(longi=0;i<testDataCount;++i)tmp+=dftol_ieee_MagicNumber(fSrc[i]);testResult=tmp;}voiddftol_test_sse2(){asm{movecx,testDataCountxoreax,eaxtestecx,ecxjleEndLoopleaedx,[fSrc+ecx*8]negecxStartLoop:cvttsd2siebx,qwordptr[edx+ecx*8]addeax,ebxincecxjnzStartLoopEndLoop:movtestResult,eax;}}longdftol_sse2_expand8(double*psrc,longcount8){longresult;asm{movecx,count8testecx,ecxjleEndLooppxorxmm0,xmm0pxorxmm1,xmm1movedx,psrcleaedx,[edx+ecx*8]negecxStartLoop://—次循环处理8个doublecvttpd2dqxmm2,xmmwordptr[edx+ecx*8]cvttpd2dqxmm3,xmmwordptr[edx+ecx*8+16]cvttpd2dqxmm4,xmmwordptr[edx+ecx*8+16*2]cvttpd2dqxmm5,xmmwordptr[edx+ecx*8+16*3]padddxmm2,xmm3padddxmm4,xmm5addecx,8padddxmm0,xmm2padddxmm1,xmm4jnzStartLoopEndLoop:padddxmm0,xmm1movapsxmm1,xmm0shufpsxmm1,xmm0,1padddxmm0,xmm1movdeax,xmm0movresult,eax}returnresult;}voiddftol_test_sse2_expand8(){longtmp=0;for(longi=0;i<testDataCount;i+=2000){tmp+=dftol_sse2_expand8(&fSrc[i],2000);//2000=8*256}//todo:因为testDataCount是2000的倍数，所以这里不用处理边界了testResult=tmp;}longdftol_sse2_expand8_prefetch(double*psrc,longcount8){longresult;asm{movecx,count8testecx,ecxjleEndLoop//预读movedx,psrcleaedx,[edx+ecx*8]negecxReadStartLoop:moveax,dwordptr[edx+ecx*8]addecx,8jnzReadStartLooppxorxmm0,xmm0pxorxmm1,xmm1movecx,count8negecxStartLoop:cvttpd2dqxmm2,xmmwordptr[edx+ecx*8]cvttpd2dqxmm3,xmmwordptr[edx+ecx*8+16]cvttpd2dqxmm4,xmmwordptr[edx+ecx*8+16*2]cvttpd2dqxmm5,xmmwordptr[edx+ecx*8+16*3]padddxmm2,xmm3padddxmm4,xmm5addecx,8padddxmm0,xmm2padddxmm1,xmm4jnzStartLoopEndLoop:padddxmm0,xmm1movapsxmm2,xmm0shufpsxmm2,xmm0,1padddxmm0,xmm2movdeax,xmm0movresult,eax}returnresult;}voiddftol_test_sse2_expand8_prefetch(){longtmp=0;for(longi=0;i<testDataCount;i+=2000){tmp+=dftol_sse2_expand8_prefetch(&fSrc[i],2000);}testResult=tmp;}intmain(){//intifor(longi=O;i<testDataCount;++i)fSrc[i]=(float)(rand()*(1.0/RAND_MAX)*(rand()—(RAND_MAX〉〉1))*rand()*(1.0/RAND_MAX));//testdoublestart0=(double)clock();for(longc=O;c