入门篇之三汇编的威力
在上一篇文章《大数阶乘之计算从入门到精通-入门篇之二》中,我们给出两个计算大数阶乘的程序,其中第2个程序由于用到64位整数的除法,速度反而更慢。在本文中,我们采用在C语言中嵌入汇编代码的方法,改进瓶颈部分,使计算速度提高到原先3倍多。
我们首先看一下计算大数阶乘的核心代码(见下),可以看到,在每次循环中,需要计算1次64位的乘法,1次64位的加法,2次64位的除法(求余视为除法)。我们知道,除法指令是慢于乘法指令的,对于64位的除法更是如此。在vc中,两个64位数的除法是调aulldiv函数来实现的,两个64位数的求余运算是调用aullrem来实现的。通过分析这两个函数的源代码(汇编代码)得知,在做除法运算时,首先检查除数,如果它小于2的32次方,则需两次除法以得到一个64位商,如果除数大于等于232,运算将更为复杂。同样在做求余运算时,如果除数小于232,为了得到一个32位的余数,也需2次除法。在我们这个例子中,除数为1000000000,小于232,因此,这段代码需4次除法。
for (carry=0,p=pTail;p>=pHead;p--) { prod=(UINT64)(*p) * (UINT64)i +carry; *p=(DWORD)(prod % TEN9); carry=prod / TEN9; }
我们注意到,在这段代码中,在进行除法运算时,其商总是小于2的32次方,因此,这个64位数的除法和求余运算可以用一条除法指令来实现。下面我们用C函数中嵌入汇编代码的方法,执行一次除法指令同时得到商和余数。函数原形为:DWORD Div_TEN9_2 (UINT64 x,DWORD *pRemainder );这个函数返加x 除以1000000000的商,除数存入pRemainder。下面是函数Div_TEN9_2和计算阶乘的代码,用C中嵌入式汇编实现。关于C代码中嵌入汇编的用法,详见MSDN。
inline DWORD Div_TEN9_2(UINT64 x,DWORD *pRemainder ) { _asm { mov eax,dword ptr [x] // low DWORD mov edx,dword ptr [x+4] // high DWORD mov ebx,TEN9 div ebx mov ebx,pRemainder mov [ebx],edx // remainder-> *remainder // eax,return value } } void calcFac2(DWORD n) { DWORD *buff,*pHead,*pTail,*p; DWORD t,i,len,carry; UINT64 prod; if (n==0) { printf("%d!=1",n); return; } //---计算并分配所需的存储空间 t=GetTickCount(); len=calcResultLen(n,TEN9); buff=(DWORD*)malloc( sizeof(DWORD)*len); if (buff==NULL) return ; //以下代码计算n! pHead=pTail=buff+len-1; for (*pTail=1,i=2;i<=n;i++) { for (carry=0,p=pTail;p>=pHead;p--) { prod=(UINT64)(*p) * (UINT64)i +(UINT64)carry; carry=Div_TEN9_2(prod,p ); } while (carry>0) { pHead--; *pHead=(DWORD)(carry % TEN9); carry /= TEN9; } } t=GetTickCount()-t; //显示计算结果 printf("It take %d ms/n",t); printf("%d!=%d",n,*pHead); for (p=pHead+1;p<=pTail;p++) printf("%09d",*p); printf("/n"); free(buff); //释放分配的内存 }
注意,本文题名虽为汇编的威力,用汇编语言并不总是那么有效,一般情况下,将关键代码用汇编改写后,性能提升的幅度并不像上面的例子那么明显,可能至多提高30%,本程序是特例。
最后的改进:如果我们分析一下这几个计算阶乘的函数,就会发现,计算阶乘其实际上是一个二重循环,内循环部分计算出(i-1)! * i, 外循环则依次计算2!,3!,4!,直到n!,假如们己计算出来r=(i-1)!,可否先算出 prod= i*(i+1)* …m, 使得i*(i+1)* …刚好小于2^32,而i*(i+1)* …m*(m+1)则 >=2^32,再计算r * prod,如此一来,可减少外循环的次数,从而提高速度。理论和测试结果都表明,当计算30000以下的阶乘时,速度可提高1倍以上。下面给出代码。
void calcFac3(DWORD n) { DWORD *buff,*p; DWORD t,carry; UINT64 prod; if (n==0) { printf("%d!=1",n); return; } //---计算并分配所需的存储空间 t=GetTickCount(); len=calcResultLen(n,TEN9); buff=(DWORD*)malloc( sizeof(DWORD)*len); if (buff==NULL) return ; //以下代码计算n! pHead=pTail=buff+len-1; for (*pTail=1,i=2;i+15<n;) { UINT64 t=i++; while (t<4294967296I64) { t *= (UINT64)i; i++; } i--; t/=i; for (carry=0,p=pTail;p>=pHead;p--) { prod=(UINT64)(*p) * t +(UINT64)carry; carry=Div_TEN9_2(prod,p ); } while (carry>0) { pHead--; *pHead=(DWORD)(carry % TEN9); carry /= TEN9; } } for (;i<=n;i++) { for (carry=0,p=pTail;p>=pHead;p--) { prod=(UINT64)(*p) * (UINT64)i +(UINT64)carry; carry=Div_TEN9_2(prod,p ); } while (carry>0) { pHead--; *pHead=(DWORD)(carry % TEN9); carry /= TEN9; } } printf("It take %d ms/n",GetTickCount()-t); //显示计算结果 printf("%d!=%d",*pHead); for (p=pHead+1;p<=pTail;p++) printf("%09d",*p); printf("/n"); free(buff);//释放分配的内存 }