阶乘之计算从入门到精通―入门篇之二

 

摘要：本文采用和《大数阶乘之计算从入门到精通―入门篇之一》几乎相同的算法思想计算阶乘，和上篇不同，本文给出的程序采用一个数组元素表示 4位或者9位10进制数的方法，使得计算速度更快，占用内存更省。本文给出两个计算阶乘的函数，程序代码简洁，速度也不慢。其中第一个程序在计算1万的阶乘时需约18.5K的内存，在迅驰1.7G笔记本用时为 0.86秒 。

 

       在《大数阶乘之计算从入门到精通―入门篇之一》中，我们给出一个计算阶乘的程序，它采用 char型数组存贮大数，1个元素表示1位十进制数字，在计算时，一次乘法可计算一位数字和一个整数的乘积。该算法具有简单直观的优点，但缺点也是明显的，速度不快，占用内存空间也较多，本文将给出一个改后的程序，有效的克服了这些缺点。

学过 80x86汇编的人都知道，8086/8088的CPU可对两个16比特的数相乘，其结果为32比特,80386及其后的32位CPU可对两个32比特的数相乘，结果为64比特（以下写作bit）。8086 CPU等16位CPU已完全淘汰，这是不去讨论它。在32位c语言编译器中，unsigned long(DWORD)型变量可以表示一个32bit的整数，unsigned short(WORD)型变量可表示一个16bit的整数。两个65535以内的数相乘，其结果完全可以存贮在一个unsigned long变量中。另外，在好多32位编译器中，也提供了64bit的整数类型（如在VC中，unsigned __int64可表示一个64bit的整数,在GCC中，longlong可表示一个64位的整数）。同理两个40亿以内的数相乘，其结果可以用一个unsigned __int64 型的变量来存储。让一个具有一次可计算两个32bit数乘法能力的CPU一次只计算1个1位10进制数和一个整数的乘法，实在是一种浪费。下面我们提出两种大数的表示法和运算方法。

第一种方法：用 WORD型数组表示大数，用DWORD型变量表示两个WORD型变量的乘积。数组的每个元素表示4位十进制数。在运算时，从这个数的最后一个元素开始，依次乘以当前乘数并加上上次的进位，其和的低4位数依然存在原位置，高几位向前进位。当乘数i小于42万时，其乘积加上进位可以用一个DWORD型变量表示，故这个程序可以计算上到42万的阶乘，当计算42万的阶乘时，占用内存空间小于1.1兆字节。至于速度，在计算1000/10000的阶乘时，在迅驰1.7G电脑约需0.015/0.86秒。

 

#include "stdlib.h" 
#include "stdio.h" 
#include "math.h" 

#define PI 3.1415926535897932384626433832795 
#define RAD 10000 
typedef unsigned long DWORD;
typedef unsigned short WORD;

//用stirling公式估算结果长度,稍微算得大一点 
DWORD calcResultLen(DWORD n,DWORD rad)
{
	double r=0.5*log(2*PI) + (n+0.5)*log(n)-n;
	return (DWORD)(r/log(rad)+2);
}

void calcFac1(DWORD n)
{
	DWORD i,carry,prod,len;
	WORD *buff,*pHead,*pTail,*p;
	
	if (n==0)
	{ printf("%d!=1",n); return; }
	
	//---计算并分配所需的存储空间 
	len=calcResultLen(n,RAD);
	buff=(WORD*)malloc( sizeof(WORD)*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=(DWORD)(*p) * i +carry;
			*p=(WORD)(prod % RAD);
			carry=prod / RAD;
		}
		while (carry>0)
		{
			pHead--;
			*pHead=(WORD)(carry % RAD);
			carry /= RAD;
		}
	}
	
	//显示计算结果 
	printf("%d!=%d",n,*pHead);
	for (p=pHead+1;p<=pTail;p++)
		printf("%04d",*p);
	printf("/n");
	free(buff);//释放分配的内存 
}

int main(int argc, char* argv[])
{
	DWORD n;
	printf("please input n:");
	scanf("%d",&n);
	calcFac1(n);
	return 0;
}

第二种方法：用 DWORD型数组表示大数，用unsigned __int64 表示两个DWORD型变量的乘积。数组的每个元素表示9位十进制数。在运算时，从这个数的最后一个元素开始，依次乘以当前乘数i(i=1..n)并加上上次的进位,其和的低9位数依然存在原位置，高几位向前进位。从算法上讲，这个程序能够计算到40亿的阶乘，在实际计算过程中，仅受限于内存的大小。至于速度，比前一个程序要慢一些，原因在于unsigned __int64的除法较慢。我们将在下一篇文章给出解决方法，下面是采用该方法计算阶乘的代码。

#define TEN9 1000000000 
void calcFac2(DWORD n)
{
	DWORD *buff,*pHead,*pTail,*p;
	DWORD t,i,len;
	UINT64 carry,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 +carry;
			*p=(DWORD)(prod % TEN9);
			carry=prod / TEN9;
		}
		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);//释放分配的内存 
}

 

 

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