入门篇之一

 

在《大数阶乘之计算从入门到精通-大数的表示》中,我们学习了如何表示和存储一个大数。在这篇文章中,我们将讨论如何对大数做乘法运算,并给出一个可以求出一个大整数阶乘的所有有效数字的程序。

大整数的存储和表示已经在上一篇文章做了详细的介绍。其中最简单的表示法是:大数用一个字符型数组来表示,数组的每一个元素表示一位十进制数字,高位在前,低位在后。那么,用这种表示法,如何做乘法运算呢?其实这个问题并不困难,可以使用模拟手算的方法来实现。回忆一下我们在小学时,是如何做一位数乘以多位数的乘法运算的。例如: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+1)/2的n次方,下面是推导过程。

Caes 1: n是奇数,则中间的那个数mid= (n+1)/2, 除了这个数外,我们可以将1到n之间的数分成n/2组,每组的两个数为 mid-i和mid+i (i=1到mid-1),如1,2,3,4,5,6,7 可以分为数4,和3对数,它们是(3,5),(2,6)和(1,7),容易知道,每对数的积都于小mid*mid,故n!小于(n+1)/2 的n的次方。

 

Case 2: n 是个偶数,则中间的两个数(n-1)/2和(n+1)/2, 我们将(n+1)/2记做mid,则其它的几对数是(mid-2,mid+1),(mid-3)(mid+2)等等,容易看出,n!小于mid 的n次方。

由以上两种情况可知,对于任意大于1的正整数n, n!<(n+1)/2的n次方。

如果想求出n!更准确的上限,可以使用司特林公式,参见该系列文章《阶乘之计算从入门到精通-近似计算之二》。

 

到此,我们已经解决大数阶乘之计算的主要难题,到了该写出一个完整程序的时候了,下面给出一个完整的代码。

 

#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((n+1)/2)); //计算n!有数数字的个数

blkLen+=4; //保险起见,多加4位

if (n<=1)

{ printf("%d!=0/n",n); return;}

 

char *arr=(char *)malloc(blkLen);

if (arr==NULL)

{ printf("alloc memory fail/n"); 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("/n");

 

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;

}