关于大数阶乘的程序编写

1、大数阶乘序大数阶乘的计算是一个有趣的话题，从中学生到大学教授，许多人都投入到这个问题的探索和研究之中，并发表了他们自己的研究成果。如果你用阶乘作关键字在google上搜索，会找到许多此类文章，另外，如果你使用google学术搜索，也能找到一些计算大数阶乘的学术论文。但这些文章和论文的深度有限，并没有给出一个高速的算法和程序。 我和许多对大数阶乘感兴趣的人一样，很早就开始编制大数阶乘的程序。从2000年开始写第一个大数阶乘程序算起，到现在大约己有6-7年的时光，期间我写了多个版本的阶乘计算器，在阶乘计算器的算法探讨和程序的编写和优化上，我花费了很大的时间和精力，品尝了这一过程中的种种甘苦，我曾因

2、一个新算法的实现而带来速度的提升而兴奋，也曾因费了九牛二虎之力但速度反而不及旧的版本而懊恼，也有过因解一个bug而通宵敖夜的情形。我写的大数阶乘的一些代码片段散见于互联网络，而算法和构想则常常萦绕在我的脑海。自以为，我对大数阶乘计算器的算法探索在深度和广度上均居于先进水平。我常常想，应该写一个关于大数阶乘计算的系列文章，一为整理自己的劳动成果，更重要的是可以让同行分享我的知识和经验，也算为IT界做一点儿贡献吧。 我的第一个大数阶乘计算器始于2000年，那年夏天，我买了一台电脑，开始在家专心学习VC，同时写了我的第一个VC程序，一个仿制windows界面的计算器。该计算器的特点是高精度和高速度，

3、它可以将四则运算的结果精确到6万位以内，将三角、对数和指数函数的结果精确到300位以内，也可以计算开方和阶乘等。当时，我碰巧看到一个叫做实用计器的软件。值得称颂的是，该计算器的作者姜边竟是一个高中生，他的这个计算器功能强大，获得了各方很高的评价。该计算的功能之一是可以计算9000以内阶乘的精确值，且速度很快。在佩服之余，也激起我写一个更好更快的程序的决心，经过几次改进，终于使我的计算器在做阶乘的精确计算时 (以9000!为例)，可比实用计算器快10倍。而当精确到30多位有效数字时，可比windows自带的计算器快上7500倍。其后的2001年1月，我在csdn上看到一个贴子，题目是“有谁可以用

4、四行代码求出1000000的阶乘”，我回复了这个贴子，给出一个相对简洁的代码，这是我在网上公布的第一个大数阶乘的程序。这一时期，可以看作是我写阶乘计算器的第一个时期。 我写阶乘计算器的第二个时期始于2003年9月，在那时我写了一组专门计算阶乘的程序，按运行速度来分，分为三个级别的版本，初级版、中级版和高级版。初级版本的算法许多人都能想到，中级版则采用大数乘以大数的硬乘法，高级版本在计算大数乘法时引入分治法。期间在csdn社区发了两个贴子，“擂台赛：计算n!(阶乘)的精确值，速度最快者2000分送上”和“擂台赛：计算n!(阶乘)的精确值，速度最快者2000分送上（续）”。其高级算的版本完成于20

5、03年11月。此后，郭先强于2004年5月10日也发表了系列贴子，“擂台：超大整数高精度快速算法”、“擂台：超大整数高精度快速算法(续)”和“擂台：超大整数高精度快速算法(续2)”, 该贴重点展示了大数阶乘计算器的速度。这个贴子一经发表即引起了热列的讨论，除了我和郭先强先生外，郭雄辉也写了同样功能的程序，那时，大家都在持续改进自己的程序，看谁的程序更快。初期，郭先强的稍稍领先，中途郭子将apfloat的代码应用到阶乘计算器中，使得他的程序胜出，后期（2004年8月后）在我将程序作了进一步的改进后，其速度又稍胜于他们。在这个贴子中，arya提到一个开放源码的程序，它的大数乘法采用FNTCRT(快

6、速数论变换中国剩余定理)。郭雄辉率先采用apflot来计算大数阶乘，后来郭先强和我也参于到apfloat的学习和改进过程中。在这点上，郭先强做得非常好，他在apfloat的基础上，成功地优化和改时算法，并应用到大数阶乘计算器上，同时他也将FNT算法应用到他的超大整数高精度快速算法库中，并在2004.10.18正式推出V3.0.2.1版。此后，我在2004年9月9日也完成一个采用FNT算法的版本，但却不及郭先强的阶乘计算器快。那时，郭先强的程序是我们所知的运算速度最快的阶乘计算器，其速度超过久负盛名的数学软件Mathematica和Maple,以及通用高精度算法库GMP。我写阶乘计算器的第三个时

7、间约开始于2006年，在2005年8月收到北大刘楚雄老师的一封e-mail，他提到了他写的一个程序在计算阶乘时比我们的更快。这使我非常吃惊，在询问后得知，其核心部分使用的是ooura写的FFT函数。ooura的FFT代码完全公开，是世界上运行的最快的FFT程序之一，从这点上，再次看到了我们和世界先进水平的差距。佩服之余，我决定深入学习FFT算法，看看能否写出和ooura速度相当或者更快的程序，同时一个更大计划开始形成，即写一组包括更多算法的阶乘计算器，包括使用FFT算法的终极版和使用无穷级数的stirling公式来计算部分精度的极速版，除此之外，我将重写和优化以前的版本，力争使速度更快，代码更

8、优。这一计划的进展并不快，曾一度停止。目前，csdn上blog数量正在迅速地增加，我也萌生了写blog的计划，借此机会，对大数阶乘之计算作一个整理，用文字和代码详述我的各个版本的算法和实现，同时也可能分析一些我在网上看到的别人写的程序，当然在这一过程中，我会继续编写未完成的版本或改写以前己经实现的版本，争取使我公开的第一份代码都是精品，这一过程可能是漫长的，但是我会尽力做下去。 菜鸟篇程序1，一个最直接的计算阶乘的程序#include "stdio.h"#include "stdlib.h" int main(int argc, char* argv)

9、long i,n,p; printf("n=?"); scanf("%d",&n); p=1; for (i=1;i<=n;i+) p*=i; printf("%d!=%dn",n,p); return 0; 程序2，稍微复杂了一些，使用了递归，一个c+初学者写的程序#include <iostream.h> long int fac(int n); void main() int n; cout<<"input a positive integer:" cin>>

10、n; long fa=fac(n); cout<<n<<"! ="<<fa<<endl; long int fac(int n) long int p; if(n=0) p=1; else p=n*fac(n-1); return p; 程序点评，这两个程序在计算12以内的数是正确，但当n>12,程序的计算结果就完全错误了，单从算法上讲，程序并没有错，可是这个程序到底错在什么地方呢？看来程序作者并没有意识到，一个long型整数能够表示的范围是很有限的。当n>=13时，计算结果溢出，在C语言，整数相乘时发生溢出时不会

11、产生任何异常，也不会给出任何警告。既然整数的范围有限，那么能否用范围更大的数据类型来做运算呢？这个主意是不错，那么到底选择那种数据类型呢？有人想到了double类型，将程序1中long型换成double类型，结果如下:#include "stdio.h"#include "stdlib.h" int main(int argc, char* argv) double i,n,p; printf("n=?"); scanf("%lf",&n); p=1.0; for (i=1;i<=n;i+) p*=i

12、; printf("%lf!=%.16gn",n,p); return 0; 运行这个程序，将运算结果并和windows计算器对比后发现，当于在170以内时，结果在误差范围内是正确。但当N>=171,结果就不能正确显示了。这是为什么呢？和程序1类似，数据发生了溢出，即运算结果超出的数据类型能够表示的范围。看来C语言提供的数据类型不能满足计算大数阶乘的需要，为此只有两个办法。1.找一个能表示和处理大数的运算的类库。2.自己实现大数的存储和运算问题。方法1不在本文的讨论的范围内。本系列的后续文章将围绕方法2来展开。大数的表示 1.大数，这里提到的大数指有效数字非常多的数，

13、它可能包含少则几十、几百位十进制数，多则几百万或者更多位十进制数。有效数字这么多的数只具有数学意义，在现实生活中，并不需要这么高的精度，比如银河系的直径有10万光年，如果用原子核的直径来度量，31位十进制数就可使得误差不超过一个原子核。 2.大数的表示： 2.1定点数和浮点数我们知道，在计算机中，数是存贮在内存(RAM)中的。在内存中存储一个数有两类格式，定点数和浮点数。定点数可以精确地表示一个整数，但数的范围相对较小，如一个32比特的无符号整数可表示04294967295之间的数，可精确到910位数字(这里的数字指10进制数字，如无特别指出，数字一律指10进制数字),而一个8字节的无符号整数

14、则能精确到19位数字。浮点数能表示更大的范围，但精度较低。当表示的整数很大的，则可能存在误差。一个8字节的双精度浮点数可表示2.22*10-308到 1.79*10308之间的数,可精确到1516位数字.2.2日常生活中的数的表示：对于这里提到的大数，上文提到的两种表示法都不能满足需求。为此，必需设计一种表示法来存储大数。我们以日常生活中的十进制数为例，看看是如何表示的。如一个数N被写成"12345",则这个数可以用一个数组a来表示，a0=1,a1=2,a2=3,a3=4,a4=5,这时数N=a4*100+a3*101+a2*102+a1*103+a0*104,(104表示

15、10的4次方，下同),10i可以叫做权，在日常生活中，a0被称作万位，也说是说它的权是10000，类似的，a1被称作千位，它的权是1000。 2.3大数在计算机语言表示：在日常生活中，我们使用的阿拉伯数字只有09共10个，按照书写习惯，一个字符表示1位数字。计算机中，我们常用的最小数据存储单位是字节，C语言称之为char,多个字节可表示一个更大的存储单位。习惯上，两个相邻字节组合起来称作一个短整数，在32位的C语言编译器中称之为short,汇编语语言一般记作word,4个相邻的字节组合起来称为一个长整数，在32位的C语言编译器中称之为long,汇编语言一般记作DWORD。在计算机中，按照权的不

16、同，数的表示可分为两种，2进制和10进制，严格说来，应该是2k进制和10K进制，前者具占用空间少，运算速度快的优点。后者则具有容易显示的优点。我们试举例说明： 例1：若一个大数用一个长为len的short型数组A来表示，并采用权从大到小的顺序依次存放，数N表示为A0 * 65536(len-1)+A1 * 65536(len-2)+.Alen-1 * 655360,这时65536称为基，其进制2的16次方。 例2：若一个大数用一个长为len的short型数组A来表示并采用权从大到小的顺序依次存放，数N=A0 * 10000(len-1)+A1 * 10000(len-2)+.Alen-1 *

17、100000,这里10000称为基，其进制为10000，即：104，数组的每个元素可表示4位数字。一般地，这时数组的每一个元素为小于10000的数。类似的，可以用long型数组，基为2324294967296来表示一个大数; 当然可以用long型组，基为1000000000来表示，这种表示法，数组的每个元素可表示9位数字。当然，也可以用char型数组，基为10。最后一种表示法，在新手写的计算大数阶乘程序最为常见，但计算速度却是最慢的。使用更大的基，可以充分发挥CPU的计算能力，计算量将更少，计算速度更快，占用的存储空间也更少。 2.4 大尾序和小尾序，我们在书写一个数时，总是先写权较大的数字，

18、后写权较小的数字，但计算机中的数并不总是按这个的顺序存放。小尾（Little Endian）就是低位字节排放在内存的低端，高位字节排放在内存的高端。例如对于一个4字节的整数0x12345678,将在内存中按照如下顺序排放, Intel处理器大多数使用小尾(Little Endian)字节序。 Address0: 0x78 Address1: 0x56Address2: 0x34 Address3:0x12大尾（Big Endian）就是高位字节排放在内存的低端，低位字节排放在内存的高端。例如对于一个4字节的整数0x12345678,将在内存中按照如下顺序排放, Motorola处理器大多数使用

19、大尾(Big Endian)字节序。Address0: 0x12 Address1: 0x34Address2: 0x56 Address3:0x78 类似的，一个大数的各个元素的排列方式既可以采用低位在前的方式，也可以采用高位在前的方式，说不上那个更好，各有利弊吧。我习惯使用高位在前的方式。2.5 不完全精度的大数表示：尽管以上的表示法可准确的表示一个整数，但有时可能只要求计算结果只精确到有限的几位。如用 windows自带的计算器计算1000的阶乘时，只能得到大约32位的数字，换名话说，windows计算器的精度为32位。1000的阶乘是一个整数，但我们只要它的前几位有效数字，象windo

20、ws计算器这样，只能表示部分有效数字的表示法叫不完全精度，不完全精度不但占用空间省，更重要的是，在只要求计算结果为有限精度的情况下，可大大减少计算量。大数的不完全精度的表示法除了需要用数组存储有数数字外，还需要一个数来表示第一个有效数字的权，10的阶乘约等于4.023872600770937e+2567，则第一个有效数字的权是102567，这时我们把2567叫做阶码。在这个例子中，我们可以用一个长为16的char型数组和一个数来表示，前者表示各位有效数字，数组的各个元素依次为:4,0,2,3,8,7,2,6,0,0,7,7,0,9,3,7，后者表示阶码，值为2567。 2.6 大数的链式存储法

21、 如果我们搜索大数阶乘的源代码，就会发现，有许多程序采用链表存储大数。尽管这种存储方式能够表示大数，也不需要事先知道一个特定的数有多少位有效数字，可以在运算过程中自动扩展链表长度。但是，如果基于运算速度和内存的考虑，强烈不建议采用这种存储方式，因为：1. 这种存储方式的内存利用率很低。基于大数乘法的计算和显示，一般需要定义双链表，假如我们用1个char表示1位十进制数，则可以这样定义链表的节点： struct _node struct _node* pre;struct _node* next;char n; ;当编译器采用默认设置，在通常的32位编译器，这个结构体将占用12字节。但这并不等于

22、说，分配具有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命中率。数组的各

23、个元素的地址是连续的，而链表的各个节点在内存中的地址是不连续的，而且具有更大的数据量。因此前者的cache的命中率高于后者，从而导致运行速度高于后者。2.3对数组的顺序访问也比链表快，如p1表示数组当前元素的地址，则计算数组的下一个地址时一般用p1+,而对链表来说则可能是p2=p2->next,毫无疑问，前者的执行速度更快。 近似计算之一 在阶乘之计算从入门到精通菜鸟篇中提到，使用double型数来计算阶乘，当n>170,计算结果就超过double数的最大范围而发生了溢出，故当n170时，就不能用这个方法来计算阶乘了，果真如此吗？No,只要肯动脑筋,办法总是有的。通过windows

24、计算器，我们知道，171！1.2410180702176678234248405241031e+309，虽然这个数不能直接用double型的数来表示，但我们可以用别的方法来表示。通过观察这个数，我们发现，这个数的表示法为科学计算法，它用两部分组成，一是尾数部分1.2410180702176678234248405241031，另一个指数部分309。不妨我们用两个数来表示这个超大的数，用double型的数来表示尾数部分，用一个long型的数来表示指数部分。这会涉及两个问题：其一是输出，这好说，在输出时将这两个部分合起来就可以了。另一个就是计算部分了，这是难点所在（其实也不难）。下面我们分析一下，

25、用什么方法可以保证不会溢出呢？我们考虑170！，这个数约等于7.257415e+306，可以用double型来表示，但当这个数乘以171就溢出了。我们看看这个等式：7.257415e+3067.257415e+306 * 100 (注1)（如用两个数来表示，则尾数部分7.257415e+306，指数部分0)(7.257415e+306 / 10300 )* (100*10300)=（7.257415e6）*（10 300）（如用两个数来表示，则尾数部分7.257415e+6，指数部分300） 依照类似的方法，在计算过程中，当尾数很大时，我们可以重新调整尾数和指数，缩小尾数，同时相应地增大指数，

26、使其表示的数的大小不变。这样由于尾数很小，再乘以一个数就不会溢出了，下面给出完整的代码。 程序3.#include "stdafx.h"#include "math.h"#define MAX_N 10000000.00 /能够计算的最大的n值，如果你想计算更大的数对数，可将其改为更大的值#define MAX_MANTISSA (1e308/MAX_N) /最大尾数 struct bigNum double n1; /表示尾数部分 int n2; /表示指数部分 ; void calcFac(struct bigNum *p,int n) int i;

27、 double MAX_POW10_LOG=(floor(log10(1e308/MAX_N); /最大尾数的常用对数的整数部分, double MAX_POW10= (pow(10.00, MAX_POW10_LOG); / 10 MAX_POW10_LOG p->n1=1; p->n2=0; for (i=1;i<=n;i+) if (p->n1>=MAX_MANTISSA) p->n1 /= MAX_POW10; p->n2 += MAX_POW10_LOG; p->n1 *=(double)i; void printfResult(str

28、uct bigNum *p,char buff) while (p->n1 >=10.00 ) p->n1/=10.00; p->n2+; sprintf(buff,"%.14fe%d",p->n1,p->n2); int main(int argc, char* argv) struct bigNum r; char buff32; int n; printf("n=?"); scanf("%d",&n); calcFac(&r,n); /计算n的阶乘 printfResult(&

29、amp;r,buff); /将结果转化一个字符串 printf("%d!=%sn",n,buff); return 0; 以上代码中的数的表示方式为：数的值等于尾数乘以 10 指数部分，当然我们也可以表示为：尾数 乘以 2 指数部分，这们将会带来这样的好处，在调整尾数部分和指数部分时，不用除法，可以依据浮点数的格式直读取阶码和修改阶码(上文提到的指数部分的标准称呼)，同时也可在一定程序上减少误差。为了更好的理解下面的代码，有必要了解一下浮点数的格式。浮点数主要分为32bit单精度和64bit双精度两种。本方只讨论64bit双精度(double型)浮点数的格式，一个doubl

30、e型浮点数包括8个字节(64bit)，我们把最低位记作bit0,最高位记作bit63,则一个浮点数各个部分定义为:第一部分尾数：bit0至bit51，共计52bit，第二部分阶码:bit52-bit62,共计11bit，第三部分符号位:bit63，0:表示正数，1表示负数。如一个数为0.xxxx * 2 exp,则exp表示指数部分，范围为1023到1024，实际存储时采用移码的表示法，即将exp的值加上0x3ff，使其变为一个0到2047范围内的一个值。函数GetExpBase2 中各语句含义如下：1.“(*pWord & 0x7fff)”，得到一个bit48-bit63这个16bi

31、t数，最高位清0。2.“>>4”,将其右移4位以清除最低位的4bit尾数，变成一个11bit的数（最高位5位为零）。3(rank-0x3ff)”, 减去0x3ff还原成真实的指数部分。以下为完整的代码。 程序4：#include "stdafx.h"#include "math.h"#define MAX_N 10000000.00 /能够计算的最大的n值，如果你想计算更大的数对数，可将其改为更大的值#define MAX_MANTISSA (1e308/MAX_N) /最大尾数typedef unsigned short WORD; str

32、uct bigNum double n1; /表示尾数部分int n2; /表示阶码部分;short GetExpBase2(double a) / 获得 a 的阶码 / 按照IEEE 754浮点数格式，取得阶码，仅仅适用于Intel 系列 cpu WORD *pWord=(WORD *)(&a)+3; WORD rank = ( (*pWord & 0x7fff) >>4 ); return (short)(rank-0x3ff); double GetMantissa(double a) / 获得 a 的 尾数/ 按照IEEE 754浮点数格式，取得尾数，仅仅适

33、用于Intel 系列 cpu WORD *pWord=(WORD *)(&a)+3;*pWord &= 0x800f; /清除阶码*pWord |= 0x3ff0; /重置阶码return a; void calcFac(struct bigNum *p,int n)int i;p->n1=1;p->n2=0;for (i=1;i<=n;i+)if (p->n1>=MAX_MANTISSA)/继续相乘可能溢出，调整之 p->n2 += GetExpBase2(p->n1); p->n1 = GetMantissa(p->n1

34、); p->n1 *=(double)i; void printfResult(struct bigNum *p,char buff)double logx=log10(p->n1)+ p->n2 * log10(2);/求计算结果的常用对数int logxN=(int)(floor(logx); /logx的整数部分sprintf(buff,"%.14fe%d",pow(10,logx-logxN),logxN);/转化为科学计算法形式的字符串 int main(int argc, char* argv)struct bigNum r;char buff

35、32;int n;printf("n=?"); scanf("%d",&n);calcFac(&r,n); /计算n的阶乘printfResult(&r,buff); /将结果转化一个字符串printf("%d!=%sn",n,buff);return 0; 程序优化的威力：程序4是一个很好的程序，它可以计算到1千万（当n更大时，p->n2可能溢出）的阶乘，但是从运行速度上讲，它仍有潜力可挖，在采用了两种优化技术后，（见程序5），速度竟提高5倍，甚至超出笔者的预计。第一种优化技术，将频繁调用的函数定义成i

36、nline函数，inline是C+关键字，如果使用纯C的编译器，可采用MACRO来代替。第二种优化技术，将循环体内的代码展开，由一个循环步只做一次乘法，改为一个循环步做32次乘法。实际速度：计算10000000！，程序4需要0.11667秒，程序5只需要0.02145秒。测试环境为迅驰1.7G,256M RAM。关于程序优化方面的内容，将在后续文章专门讲述。下面是被优化后的部分代码，其余代码不变。 程序5的部分代码： inline short GetExpBase2(double a) / 获得 a 的阶码/ 按照IEEE 754浮点数格式，取得阶码，仅仅适用于Intel 系列 cpu WOR

37、D *pWord=(WORD *)(&a)+3; WORD rank = ( (*pWord & 0x7fff) >>4 ); return (short)(rank-0x3ff);inline double GetMantissa(double a) / 获得 a 的 尾数 / 按照IEEE 754浮点数格式，取得尾数，仅仅适用于Intel 系列 cpu WORD *pWord=(WORD *)(&a)+3; *pWord &= 0x800f; /清除阶码 *pWord |= 0x3ff0; /重置阶码 return a; void calcFac

38、(struct bigNum *p,int n) int i,t; double f,max_mantissa; p->n1=1;p->n2=0;t=n-32; for (i=1;i<=t;i+=32) p->n2 += GetExpBase2(p->n1); p->n1 = GetMantissa(p->n1); f=(double)i; p->n1 *=(double)(f+0.0); p->n1 *=(double)(f+1.0); p->n1 *=(double)(f+2.0); p->n1 *=(double)(f+3

39、.0); p->n1 *=(double)(f+4.0); p->n1 *=(double)(f+5.0); p->n1 *=(double)(f+6.0); p->n1 *=(double)(f+7.0); p->n1 *=(double)(f+8.0); p->n1 *=(double)(f+9.0); p->n1 *=(double)(f+10.0); p->n1 *=(double)(f+11.0); p->n1 *=(double)(f+12.0); p->n1 *=(double)(f+13.0); p->n1 *=

40、(double)(f+14.0); p->n1 *=(double)(f+15.0); p->n1 *=(double)(f+16.0); p->n1 *=(double)(f+17.0); p->n1 *=(double)(f+18.0); p->n1 *=(double)(f+19.0); p->n1 *=(double)(f+20.0); p->n1 *=(double)(f+21.0); p->n1 *=(double)(f+22.0); p->n1 *=(double)(f+23.0); p->n1 *=(double)(f

41、+24.0); p->n1 *=(double)(f+25.0); p->n1 *=(double)(f+26.0); p->n1 *=(double)(f+27.0); p->n1 *=(double)(f+28.0); p->n1 *=(double)(f+29.0); p->n1 *=(double)(f+30.0); p->n1 *=(double)(f+31.0); for (;i<=n;i+) p->n2 += GetExpBase2(p->n1); p->n1 = GetMantissa(p->n1); p-

42、>n1 *=(double)(i); 注1：100，表示10的0次方近似计算之二 在阶乘之计算从入门到精通近似计算之一中，我们采用两个数来表示中间结果，使得计算的范围扩大到1千万，并可0.02秒内完成10000000！的计算。在保证接近16位有效数字的前提下，有没有更快的算法呢。很幸运，有一个叫做Stirling的公式，它可以快速计算出一个较大的数的阶乘，而且数越大，精度越高。有 :/mathworld.wolfram 查找Stirlings Series可找到2个利用斯特林公式求阶乘的公式，一个是普通形式，一个是对数形式。前一个公式包含一个无穷级数和的形式，包含级数前5项的公式为n!=

43、sqrt(2*PI*n)*(n/e)n*(1+ 1/12/n+ 1/288/n2 139/51840/n3-571/2488320/n4+),这里的PI为圆周率，e为自然对数的底。后一个公式也为一个无穷级数和的形式，一般称这个级数为斯特林(Stirling)级数,包括级数前3项的n!的对数公式为：ln(n!)=0.5*ln(2*PI)+(n+0.5)*ln(n)-n+(1/12/n -1/360/n3 + 1/1260/n5)-,下面我们分别用这两个公式计算n!.完整的代码如下：#include "stdafx.h"#include "math.h"#d

44、efine PI 3.1415926535897932384626433832795#define E 2.7182818284590452353602874713527struct bigNum double n; /科学计数法表示的尾数部分 int e; /科学计数法表示的指数部分,若一个bigNum为x,这x=n * 10e;const double a1= 1.00, 1.0/12.0, 1.0/288, -139/51840,-571.0/2488320.0 ;const double a2= 1.0/12.0, -1.0/360, 1.0/1260.0 ;void printfRe

45、sult(struct bigNum *p,char buff) while (p->n >=10.00 ) p->n/=10.00; p->e+; sprintf(buff,"%.14fe%d",p->n,p->e);/ n!=sqrt(2*PI*n)*(n/e)n*(1+ 1/12/n+ 1/288/n2 -139/51840/n3-571/2488320/n4+)void calcFac1(struct bigNum *p,int n) double logx, s, /级数的和 item; /级数的每一项 int i; / xn=

46、 pow(10,n * log10(x); logx= n* log10(double)n/E); p->e= (int)(logx); p->n= pow(10.0, logx-p->e); p->n *= sqrt( 2* PI* (double)n); /求(1+ 1/12/n+ 1/288/n2 -139/51840/n3-571/2488320/n4+) for (item=1.0,s=0.0,i=0;i<sizeof(a1)/sizeof(double);i+) s+= item * a1i; item /= (double)n; /item= 1/(

47、ni) p->n *=s;/ln(n!)=0.5*ln(2*PI)+(n+0.5)*ln(n)-n+(1/12/n -1/360/n3 + 1/1260/n5 -.)void calcFac2(struct bigNum *p,int n) double logR; double s, /级数的和 item; /级数的每一项 int i; logR=0.5*log(2.0*PI)+(double)n+0.5)*log(n)-(double)n; /s= (1/12/n -1/360/n3 + 1/1260/n5) for (item=1/(double)n,s=0.0,i=0;i<

48、sizeof(a2)/sizeof(double);i+) s+= item * a2i; item /= (double)(n)* (double)n; /item= 1/(n(2i+1) logR+=s; /根据换底公式，log10(n)=ln(n)/ln(10) p->e = (int)( logR / log(10); p->n = pow(10.00, logR/log(10) - p->e);int main(int argc, char* argv) struct bigNum r; char buff32; int n; printf("n=?&qu

49、ot;); scanf("%d",&n); calcFac1(&r,n); /用第一种方法，计算n的阶乘 printfResult(&r,buff); /将结果转化一个字符串 printf("%d!=%s by method 1n",n,buff); calcFac2(&r,n); /用第二种方法，计算n的阶乘 printfResult(&r,buff); /将结果转化一个字符串 printf("%d!=%s by method 2n",n,buff); return 0;程序运行时间的测量 算

50、法的好坏有好多评价指标，其中一个重要的指标是时间复杂度。如果两个程序完成一个同样的任务，即功能相同，处理的数据相同，那么运行时间较短者为优。操作系统和库函数一般都提供了对时间测量的函数，这么函数一般都会返回一个代表当前时间的数值，通过在运行某个程序或某段代码之前调用一次时间函数来得到一个数值，在程序或者代码段运行之后再调用一次时间函数来得到另一个数值，将后者减去前者即为程序的运行时间。在windwos平台(指windwow95及以后的版本，下同),常用的用于测量时间的函数或方法有三种：1.API函数GetTickCount或C函数clock,2.API函数QueryPerformanceCounter，3:汇编指令RDSTC1API函数GetTickCount：函数原形：DWORD GetTickCount(VOID)；该函数取回从电脑开机至现在的毫秒数，即每个时间单位为1毫秒。他的分辨率比较低，常用在测量用时较长程序，如果你的程序用时为100毫秒以上，可以使用这个函数.另一个和GetTickCount类似的函数是clock，该函数的回的时间的单位的是CLOCKS_PER_SEC，在windows95/2000操作系统，该值是1000,也就是说，在windows平台，这两个函数的功能几乎完全相同。2.API函数QueryPerformanceCounter：函数原形：BOOL Q