阶乘之计算从入门到精通－任意阶乘计算

摘要：本文讨论如何使用一个简单的算法计算一个大整数n的阶乘，大数采用char数组存储，一个元素表示1位10进制数。本中给出一个完整的计算大数阶乘的程序,该程序在迅驰1.7G笔记本上计算10000的阶乘大约2.7秒。
 
    在《大数阶乘之计算从入门到精通－大数的表示》中，我们学习了如何表示和存储一个大数。在这篇文章中，我们将讨论如何对大数做乘法运算，并给出一个可以求出一个整数n的阶乘的所有有效数字的程序。大整数的存储和表示已经在上一篇文章做了详细的介绍。其中最简单的表示法是:大数用一个字符型数组来表示，数组的每一个元素表示一位十进制数字，高位在前，低位在后。那么，用这种表示法，如何做乘法运算呢？其实这个问题并不困难，可以使用模拟手算的方法来实现。回忆一下我们在小学时，是如何做一位数乘以多位数的乘法运算的。例如：2234*8。
  
 
我们将被乘数表示为一个数组A[], a[1],a[2],a[3],a[4]分别为2,2,3,4,a[0]置为0。
Step1: 从被乘数的个位a[4]起，取出一个数字4.
Step2: 与乘数8相乘，其积是两位数32，其个位数2作为结果的个位,存入a[4], 十位数3存入进位c。
Step3: 取被乘数的上一位数字a[3]与乘数相乘，并加上上一步计算过程的进位C,得到27，将这个数的个位7作为结果的倒数第二位，存入a[3]，十位数2存入进位c。
Step4：重复Step3,取a[i]（i依次为4,3,2,1）与乘数相乘并加上c,其个位仍存入a[i],　十位数字存入c，直到i等于1为止。
Step5:将最后一步的进位c作为积的最高位a[0]。
 
这一过程其实和小学学过的多位数乘以一位数的珠算乘法一模一样，学过珠算乘法的朋友还有印象吗？
 
在计算大数阶乘的过程中，乘数一般不是一位数字，那么如何计算呢？我们可以稍作变通，将上次的进位加上本次的积得到数P,　将P除以10的余数做为结果的本位，将P除以10的商作为进位。当被乘数的所有数字都和乘数相乘完毕后，将进位C放在积的最前面即可。下面给出C语言代码。
一个m位数乘以n位数，其结果为m+n-1,或者m+n位，所以需首先定义一个至少m+n个元素的数组，并置前n位为0。
 
计算一个m位的被乘数乘以一个n位的整数k，积仍存储于数组a

void mul(unsigned char a[],unsigned long k,int m,int n)
{
	int i;
	unsigned long p;
	unsigned long c=0;
	
	for ( i=m+n-1; i>=n;i--)
	{
		p= a[i] * k +c;
		a[i]=(unsigned char)( p % 10);
		c= p / 10;
	}
	
	while (c>0)
	{
		a[i]=(unsigned char)( c % 10);
		i--;
		c /=10;
	}
}
 
int main(int argc, char* argv[])
{
	int i;
	unsigned char a[]={0,0,0,2,3,4,5};
	mul(a,678,4,3);
	i=0;
	while ( a[i]==0)
		i++;
	for (;i<4+3;i++)
		printf("%c",a[i]+’0’); //由于数a[i]（0<=a[i] <=9）对应的可打印字任符为’0’到’9’,所以显示为i＋’0’ 
	return 0;
}

　　

从上面的例子可知，在做乘法之前，必须为数组保留足够的空间。具体到计算n!的阶乘时，必须准备一个能容纳的n!的所有位数的数组或者内存块。即数组采有静态分配或者动态分配。前者代码简洁，但只适应于n小于一个固定的值，后者灵活性强，只要有足够的内存，可计算任意n的阶乘，我们这里讨论后一种情况，如何分配一块大小合适的内存。

n!有多少位数呢？我们给出一个近似的上限值：n!<(n+12)nn!<(n+12)n，下面是推导过程。
Caes 1: n是奇数，则中间的那个数m=(n+1)/2, 除了这个数外，我们可以将1到n之间的数分成n/2组，每组的两个数为 m-i 和m+i (i=1到m-1),如1,2,3,4,5,6,7 可以分为数4，和3对数，它们是(3,5),(2,6)和(1,7),容易知道，每对数的积都于小 m2m2，故 n!<mn=(n+12)nn!<mn=(n+12)n的次方。
 
 
     Case 2: n 是个偶数，则中间的两个数(n/2)和(n/2+1),其它的几对儿数是(n/2-1,n/2+2),(n/2-2,n/2+3)等等。容易证明，这几对儿数中，中间的那对数n/2和n/2+1的乘积最大，因为共有n/2对数，故
    n!<(n2(n2+1))n/2=(n2+2n4)n2<((n+12)2)n2=(n+12)nn!<(n2(n2+1))n/2=(n2+2n4)n2<((n+12)2)n2=(n+12)n

     由以上两种情况可知，对于任意大于1的正整数n, n!<(n+12)nn!<(n+12)n。如果想求出n!更准确的上限，可以使用司特林公式，参见该系列文章《阶乘之计算从入门到精通－近似计算之二》。
 
到此，我们已经解决大数阶乘之计算的主要难题，到了该写出一个完整程序的时候了，下面给出一个完整的代码。

程序运行

新建一个VS工程

选择空项目

添加一个cpp文件

把代码贴进去

遇到 scanf（）函数报错问题

右键项目-属性

图示选否

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

void calcFac(unsigned long n)
{
	unsigned long i, j, head, tail;
	int blkLen = (int)(n*log10(double(n + 1) / 2)); //计算n!有数数字的个数 
	blkLen += 4;  //保险起见，多加4位 

	if (n <= 1)
	{
		printf("%d!=0
", n);     return;
	}

	char *arr = (char *)malloc(blkLen);
	if (arr == NULL)
	{
		printf("alloc memory fail
"); return;
	}

	memset(arr, 0, sizeof(char)*blkLen);
	head = tail = blkLen - 1;
	arr[tail] = 1;

	for (i = 2; i <= n; i++)
	{
		unsigned long c = 0;
		for (j = tail; j >= head; j--)
		{
			unsigned long prod = arr[j] * i + c;
			arr[j] = (char)(prod % 10);
			c = prod / 10;
		}
		while (c>0)
		{
			head--;
			arr[head] = (char)(c % 10);
			c /= 10;
		}
	}
	printf("%d!=
", n);
	for (i = head; i <= tail; i++)
		printf("%c", arr[i] + '0');
	printf("
");

	free(arr);
}

void testCalcFac()
{
	int n;
	while (1)
	{
		printf("n=?
");
		scanf("%ld", &n);
		if (n == 0)
			break;
		calcFac(n);
	}
}

int main(int argc, char* argv[])
{
	testCalcFac();
	return 0;
}