阶乘之计算从入门到精通―入门篇之三汇编的威力

 

摘要： 本文提供了2个计算阶乘的程序。第1个程序采用在C中嵌入汇编代码的方法，改进上篇中了程序2的瓶颈部分，使速度提高到原先的3倍多。第2个程序进一步改进了算法，在计算1万的阶乘时，比上一篇中的程序2快5-6倍，计算10000的阶乘， 在迅驰1.7G的仅需0.25秒。

     

      在上一篇文章《大数阶乘之计算从入门到精通 -入门篇之二》中，我们给出两个计算大数阶乘的程序，其中第2个程序由于用到64位整数的除法，速度反而更慢。在本文中，我们采用在C语言中嵌入汇编代码的方法，改进瓶颈部分，使计算速度提高到原先3倍多。

我们首先看一下计算大数阶乘的核心代码(见下)，可以看到，在每次循环中，需要计算1次64位的乘法，1次64位的加法，2次64位的除法( 求余视为除法)。我们知道，除法指令是慢于乘法指令的，对于64位的除法更是如此。在 vc中，两个64位数的除法是调 aulldiv函数来实现的，两个64位数 的求余运算是调用 aullrem来实现的。通过分析这两个函数的源代码(汇编代码)得知，在做除法运算时，首先检查除数，如果它小于2的32次方,则需两次除法以得到一个64位商，如果除数大于等于2 ³²，运算将更为复杂。同样在 做求余运算时，如果除数小于2 ³²，为了得到一个32位的余数，也需2次除法。在我们这个例子中，除数为1000000000，小于2 ³²，因此，这段代码需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,*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+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",n,*pHead);
	for (p=pHead+1;p<=pTail;p++)
		printf("%09d",*p);
	printf("/n");
	
	free(buff);//释放分配的内存 
}

liangbch@263.net ，版权所有，转载请注明出处。