阶乘之计算从入门到精通－大数的表示
http://blog.csdn.net/liangbch/article/details/1562014
1.大数，这里提到的大数指有效数字非常多的数，它可能包含少则几十、几百位十进制数，多则几百万或者更多位十进制数。有效数字这么多的数只具有数学意义，在现实生活中，并不需要这么高的精度，比如银河系的直径有10万光年，如果用原子核的直径来度量，31位十进制数就可使得误差不超过一个原子核。
2.大数的表示：
2.1定点数和浮点数
我们知道，在计算机中，数是存贮在内存(RAM)中的。在内存中存储一个数有两类格式，定点数和浮点数。定点数可以精确地表示一个整数，但数的范围相对较小，如一个32比特的无符号整数可表示0－4294967295之间的数，可精确到9－10位数字(这里的数字指10进制数字，如无特别指出，数字一律指10进制数字),而一个8字节的无符号整数则能精确到19位数字。浮点数能表示更大的范围，但精度较低。当表示的整数很大的，则可能存在误差。一个8字节的双精度浮点数可表示2.22*10^-308到 1.79*10^308之间的数,可精确到15－16位数字.
2.2日常生活中的数的表示：
对于这里提到的大数，上文提到的两种表示法都不能满足需求。为此，必需设计一种表示法来存储大数。我们以日常生活中的十进制数为例，看看是如何表示的。如一个数N被写成"12345",则这个数可以用一个数组a来表示，a[0]=1, a[1]=2, a[2]=3, a[3]=4, a[4]=5,这时数N= a[4]*10^0 +a[3]*10^1 +a[2]*10^2 +a[1]*10^3 +a[0]*10^4, (10^4表示10的4次方，下同),10^i可以叫做权，在日常生活中，a[0]被称作万位，也说是说它的权是10000，类似的，a[1]被称作千位，它的权是1000。
2.3 大数在计算机语言表示：
在日常生活中，我们使用的阿拉伯数字只有0－9共10个，按照书写习惯，一个字符表示1位数字。计算机中，我们常用的最小数据存储单位是字节，C语言称之为char,多个字节可表示一个更大的存储单位。习惯上，两个相邻字节组合起来称作一个短整数，在32位的C语言编译器中称之为short,汇编语语言一般记作word,4个相邻的字节组合起来称为一个长整数，在32位的C语言编译器中称之为long,汇编语言一般记作DWORD。在计算机中，按照权的不同，数的表示可分为两种，2进制和10进制，严格说来，应该是2^k进制和10^K进制，前者具占用空间少，运算速度快的优点。后者则具有容易显示的优点。我们试举例说明：
例1：若一个大数用一个长为len的short型数组A来表示，并采用权从大到小的顺序依次存放，数N表示为A[0] * 65536^(len-1)+A[1] * 65536^(len-2)+...A[len-1] * 65536^0,这时65536称为基，其进制2的16次方。
例2：若一个大数用一个长为len的short型数组A来表示并采用权从大到小的顺序依次存放，数N=A[0] * 10000^(len-1)+A[1] * 10000^(len-2)+...A[len-1] * 10000^0,这里10000称为基，其进制为10000，即：10^4，数组的每个元素可表示4位数字。一般地，这时数组的每一个元素为小于10000的数。类似的，可以用long型数组，基为2＾32＝4294967296来表示一个大数; 当然可以用long型组，基为1000000000来表示，这种表示法，数组的每个元素可表示9位数字。当然，也可以用char型数组，基为10。最后一种表示法，在新手写的计算大数阶乘程序最为常见，但计算速度却是最慢的。使用更大的基，可以充分发挥CPU的计算能力，计算量将更少，计算速度更快，占用的存储空间也更少。
2.4 大尾序和小尾序，我们在书写一个数时，总是先写权较大的数字，后写权较小的数字，但计算机中的数并不总是按这个的顺序存放。小尾（Little Endian）就是低位字节排放在内存的低端，高位字节排放在内存的高端。例如对于一个4字节的整数0x12345678,将在内存中按照如下顺序排放, Intel处理器大多数使用小尾(Little Endian)字节序。
Address[0]: 0x78
Address[1]: 0x56
Address[2]: 0x34
Address[3]:0x12
大尾（Big Endian）就是高位字节排放在内存的低端，低位字节排放在内存的高端。例如对于一个4字节的整数0x12345678,将在内存中按照如下顺序排放, Motorola处理器大多数使用大尾(Big Endian)字节序。
Address[0]: 0x12
Address[1]: 0x34
Address[2]: 0x56
Address[3]:0x78
类似的，一个大数的各个元素的排列方式既可以采用低位在前的方式，也可以采用高位在前的方式，说不上那个更好，各有利弊吧。我习惯使用高位在前的方式。
2.5 不完全精度的大数表示：
尽管以上的表示法可准确的表示一个整数，但有时可能只要求计算结果只精确到有限的几位。如用 windows自带的计算器计算1000的阶乘时，只能得到大约32位的数字，换名话说，windows计算器的精度为32位。1000的阶乘是一个整数，但我们只要它的前几位有效数字，象windows计算器这样，只能表示部分有效数字的表示法叫不完全精度，不完全精度不但占用空间省，更重要的是，在只要求计算结果为有限精度的情况下，可大大减少计算量。大数的不完全精度的表示法除了需要用数组存储有数数字外，还需要一个数来表示第一个有效数字的权，1000的阶乘约等于4.023872600770937e+2567，则第一个有效数字的权是10^2567，这时我们把2567叫做阶码。在这个例子中，我们可以用一个长为16的char型数组和一个数来表示，前者表示各位有效数字，数组的各个元素依次为:4,0,2,3,8,7,2,6,0,0,7,7,0,9,3,7，后者表示阶码，值为2567。
2.6 大数的链式存储法
如果我们搜索大数阶乘的源代码，就会发现，有许多程序采用链表存储大数。尽管这种存储方式能够表示大数，也不需要事先知道一个特定的数有多少位有效数字，可以在运算过程中自动扩展链表长度。但是，如果基于运算速度和内存的考虑，强烈
不建议采用这种存储方式，因为：
1. 这种存储方式的内存利用率很低。基于大数乘法的计算和显示，一般需要定义双链表，假如我们用1个char表示1位十进制数，则可以这样定义链表的节点：
struct _node
{
struct _node* pre;
struct _node* next;
char n;
};
当编译器采用默认设置，在通常的32位编译器，这个结构体将占用12字节。但这并不等于说，分配具有1000个节点的链表需要1000*12字节。不要忘记，操作系统或者库函数在从内存池中分配和释放内存时，也需要维护一个链表。实验表明，在VC编译的程序，一个节点总的内存占用量为 sizeof(struct _node) 向上取16的倍数再加8字节。也就是说，采用这种方式表示n位十进制数需要 n*24字节，而采用1个char型数组仅需要n字节。
2采用链表方式表示大数的运行速度很慢.
2.1如果一个大数需要n个节点，需要调用n次malloc(C)或new(C++)函数，采用动态数组则不要用调用这么多次malloc.
2.2 存取数组表示的大数比链表表示的大数具有更高的cache命中率。数组的各个元素的地址是连续的，而链表的各个节点在内存中的地址是不连续的，而且具有更大的数据量。因此前者的cache的命中率高于后者，从而导致运行速度高于后者。
2.3对数组的顺序访问也比链表快，如p1表示数组当前元素的地址，则计算数组的下一个地址时一般用p1++,而对链表来说则可能是p2=p2->next,毫无疑问，前者的执行速度更快。
摘要：本文讨论如何使用一个简单的算法计算一个大整数的阶乘，大数采用char数组存储，一个元素表示1位10进制数。本中给出一个完整的计算大数阶乘的程序,该程序在迅驰1.7G笔记本上计算10000的阶乘大约2.7秒。
在《大数阶乘之计算从入门到精通－大数的表示》中，我们学习了如何表示和存储一个大数。在这篇文章中，我们将讨论如何对大数做乘法运算，并给出一个可以求出一个大整数阶乘的所有有效数字的程序。
大整数的存储和表示已经在上一篇文章做了详细的介绍。其中最简单的表示法是:大数用一个字符型数组来表示，数组的每一个元素表示一位十进制数字，高位在前，低位在后。那么，用这种表示法，如何做乘法运算呢？其实这个问题并不困难，可以使用模拟手算的方法来实现。回忆一下我们在小学时，是如何做一位数乘以多位数的乘法运算的。例如：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
计算一个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;
}
//由于数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次方。
摘要：本文采用和《大数阶乘之计算从入门到精通―入门篇之一》几乎相同的算法思想计算阶乘，和上篇不同，本文给出的程序采用一个数组元素表示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中，long long可表示一个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 ，版权所有，转载请注明出处。
619
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
