负数的阶乘

c语⾔⼤整数阶乘计算器,⼤数阶乘_yuanmercu_oxxdl_新浪博

客

菜鸟篇

程序1，⼀个最直接的计算阶乘的程序

#include"stdio.h"

#include"stdlib.h"

intmain(intargc,char*argv[])

{

longi,n,p;

printf("n=?");

scanf("%d",&n);

p=1;

for(i=1;i<=n;i++)

p*=i;

printf("%d!=%dn",n,p);

return0;

}

程序2，稍微复杂了⼀些，使⽤了递归，⼀个c++初学者写的程序

#includelongintfac(intn);voidmain(){intn;cout<>n;longfa=fac(n);cout<

程序点评，这两个程序在计算12以内的数是正确，但当n>12,程序的计算结果就完全错误了，单从算法上讲，程序并没有错，可是这个程序

到底错在什么地⽅呢？看来程序作者并没有意识到，⼀个long型整数能够表⽰的范围是很有限的。当n>=13时，计算结果溢出，在C语⾔，

整数相乘时发⽣溢出时不会产⽣任何异常，也不会给出任何警告。既然整数的范围有限，那么能否⽤范围更⼤的数据类型来做运算呢？这个

主意是不错，那么到底选择那种数据类型呢？有⼈想到了double类型，将程序1中long型换成double类型，结果如下:

#include"stdio.h"

#include"stdlib.h"

intmain(intargc,char*argv[])

{

doublei,n,p;

printf("n=?");

scanf("%lf",&n);

p=1.0;

for

(i=1;i<=n;i++)

p*=i;

printf("%lf!=%.16gn",n,p);

return0;

}

运⾏这个程序，将运算结果并和windows计算器对⽐后发现，当于在170以内时，结果在误差范围内是正确。但当N>=171,结果就不能正

确显⽰了。这是为什么呢？和程序1类似，数据发⽣了溢出，即运算结果超出的数据类型能够表⽰的范围。看来C语⾔提供的数据类型不能

满⾜计算⼤数阶乘的需要，为此只有两个办法。1.找⼀个能表⽰和处理⼤数的运算的类库。2.⾃⼰实现⼤数的存储和运算问题。⽅法1不在

本⽂的讨论的范围内。本系列的后续⽂章将围绕⽅法2来展开。

⼤数的表⽰

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处理器⼤多数使⽤⼩尾(LittleEndian)字节序。

Address[0]:0x78

Address[1]:0x56

Address[2]:0x34

Address[3]:0x12

⼤尾(Big

Endian)就是⾼位字节排放在内存的低端，低位字节排放在内存的⾼端。例如对于⼀个4字节的整数0x12345678,将在内存中按照如下顺

序排放,

Motorola处理器⼤多数使⽤⼤尾(BigEndian)字节序。

Address[0]:0x12

Address[1]:0x34

Address[2]:0x56

Address[3]:0x78

类似的，⼀个⼤数的各个元素的排列⽅式既可以采⽤低位在前的⽅式，也可以采⽤⾼位在前的⽅式，说不上那个更好，各有利弊吧。我习惯

使⽤⾼位在前的⽅式。

2.5不完全精度的⼤数表⽰：

尽管以上的表⽰法可准确的表⽰⼀个整数，但有时可能只要求计算结果只精确到有限的⼏位。如⽤

windows⾃带的计算器计算1000的阶乘时，只能得到⼤约32位的数字，换名话说，windows计算器的精度为32位。1000的阶乘是⼀个

整数，但我们只要它的前⼏位有效数字，象windows计算器这样，只能表⽰部分有效数字的表⽰法叫不完全精度，不完全精度不但占⽤空间

省，更重要的是，在只要求计算结果为有限精度的情况下，可⼤⼤减少计算量。⼤数的不完全精度的表⽰法除了需要⽤数组存储有数数字

外，还需要⼀个数来表⽰第⼀个有效数字的权，10的阶乘约等于4.937e+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;

charn;

};

当编译器采⽤默认设置，在通常的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,毫⽆疑问，前者的执⾏速度更快。

近似计算之⼀

在＜阶乘之计算从⼊门到精通－菜鸟篇＞中提到，使⽤double型数来计算阶乘，当n>170,计算结果就超过double数的最⼤范围⽽发⽣了溢

出，故当n＞170时，就不能⽤这个⽅法来计算阶乘了，果真如此吗？No,只要肯动脑筋,办法总是有的。

通过windows计算器，我们知道，171！＝1.24241031e+309，虽然这个数不能直接⽤double型的

数来表⽰，但我们可以⽤别的⽅法来表⽰。通过观察这个数，我们发现，这个数的表⽰法为科学计算法，它⽤两部分组成，⼀是尾数部分

1.24241031，另⼀个指数部分309。不妨我们⽤两个数来表⽰这个超⼤的数，⽤double型的数来表

⽰尾数部分，⽤⼀个long型的数来表⽰指数部分。这会涉及两个问题：其⼀是输出，这好说，在输出时将这两个部分合起来就可以了。另⼀

个就是计算部分了，这是难点所在(其实也不难)。下⾯我们分析⼀下，⽤什么⽅法可以保证不会溢出呢？

我们考虑170！，这个数约等于7.257415e+306，可以⽤double型来表⽰，但当这个数乘以171就溢出了。我们看看这个等式：

7.257415e+306

＝7.257415e+306*10^0

(注1)(如⽤两个数来表⽰，则尾数部分7.257415e+306，指数部分0)

＝(7.257415e+306/10^300)*(10^0*10^300)

=(7.257415e6)*(10^300) (如⽤两个数来表⽰，则尾数部分7.257415e+6，指数部分300)

依照类似的⽅法，在计算过程中，当尾数很⼤时，我们可以重新调整尾数和指数，缩⼩尾数，同时相应地增⼤指数，使其表⽰的数的⼤⼩不

变。这样由于尾数很⼩，再乘以⼀个数就不会溢出了，下⾯给出完整的代码。

程序3.

#include"stdafx.h"

#include"math.h"

#defineMAX_N

10000000.00//能够计算的最⼤的n值，如果你想计算更⼤的数对数，可将其改为更⼤的值

#defineMAX_MANTISSA(1e308/MAX_N)//最⼤尾数

structbigNum

{

double

n1;//表⽰尾数部分

intn2;//表⽰指数部分

};

voidcalcFac(structbigNum*p,intn)

{

inti;

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;

}

}

voidprintfResult(structbigNum*p,charbuff[])

{

while(p->n1>=10.00)

{p->n1/=10.00;p->n2++;}

sprintf(buff,"%.14fe%d",p->n1,p->n2);

}

intmain(intargc,char*argv[])

{

structbigNumr;

charbuff[32];

intn;

printf("n=?");

scanf("%d",&n);

calcFac(&r,n);//计算n的阶乘

printfResult(&r,buff);//将结果转化⼀个字符串

printf("%d!=%sn",n,buff);

return0;

}

以上代码中的数的表⽰⽅式为：数的值等于尾数乘以10^指数部分，当然我们也可以表⽰为：尾数乘以2^

指数部分，这们将会带来这样的好处，在调整尾数部分和指数部分时，不⽤除法，可以依据浮点数的格式直读取阶码和修改阶码(上⽂提到

的指数部分的标准称呼)，同时也可在⼀定程序上减少误差。为了更好的理解下⾯的代码，有必要了解⼀下浮点数的格式。浮点数主要分为

32bit单精度和64bit双精度两种。本⽅只讨论64bit双精度(double型)浮点数的格式，⼀个double型浮点数包括8个字节(64bit)，我们把

最低位记作bit0,最⾼位记作bit63,则⼀个浮点数各个部分定义为:第⼀部分尾数：bit0⾄bit51，共计52bit，第⼆部分阶码:bit52-bit62,共

计11bit，第三部分符号位:bit63，0:表⽰正数，1表⽰负数。如⼀个数为

*2^

exp,则exp表⽰指数部分，范围为－1023到1024，实际存储时采⽤移码的表⽰法，即将exp的值加上0x3ff，使其变为⼀个0到2047范围

内的⼀个值。函数GetExpBa2

中各语句含义如下：1.“(*pWord&

0x7fff)”，得到⼀个bit48-bit63这个16bit数，最⾼位清0。2.“>>4”,将其右移4位以清除最低位的4bit尾数，变成⼀个11bit的数(最⾼

位5位为零)。3(rank-0x3ff)”,

减去0x3ff还原成真实的指数部分。以下为完整的代码。

程序4：

#include"stdafx.h"

#include"math.h"

#defineMAX_N

10000000.00//能够计算的最⼤的n值，如果你想计算更⼤的数对数，可将其改为更⼤的值

#defineMAX_MANTISSA(1e308/MAX_N)//最⼤尾数

typedefunsignedshortWORD;

structbigNum

{

double

n1;//表⽰尾数部分

intn2;//表⽰阶码部分

};

shortGetExpBa2(doublea)//获得a的阶码

{

//按照IEEE754浮点数格式，取得阶码，仅仅适⽤于Intel系列cpu

WORD*pWord=(WORD*)(&a)+3;

WORDrank

=((*pWord&0x7fff)

>>4);

return

(short)(rank-0x3ff);

}

doubleGetMantissa(doublea)//获得a的尾数

{

//按照IEEE754浮点数格式，取得尾数，仅仅适⽤于Intel系列cpu

WORD*pWord=(WORD*)(&a)+3;

*pWord&=0x800f;//清除阶码

*pWord|=0x3ff0;//重置阶码

returna;

}

voidcalcFac(structbigNum*p,intn)

{

inti;

p->n1=1;

p->n2=0;

for(i=1;i<=n;i++)

{

if

(p->n1>=MAX_MANTISSA) //继续相乘可能溢出，调整之

{

p->n2+=GetExpBa2(p->n1);

p->n1=GetMantissa(p->n1);

}

p->n1*=(double)i;

}

}

voidprintfResult(structbigNum*p,charbuff[])

{

doublelogx=log10(p->n1)+p->n2*

log10(2);//求计算结果的常⽤对数

intlogxN=(int)(floor(logx));//logx的整数部分

sprintf(buff,"%.14fe%d",pow(10,logx-logxN),logxN);//转化为科学计算法形式的字符串

}

intmain(intargc,char*argv[])

{

structbigNumr;

charbuff[32];

intn;

printf("n=?");

scanf("%d",&n);

calcFac(&r,n);//计算n的阶乘

printfResult(&r,buff);//将结果转化⼀个字符串

printf("%d!=%sn",n,buff);

return0;

}

程序优化的威⼒：

程序4是⼀个很好的程序，它可以计算到1千万(当n更⼤时，p->n2可能溢出)的阶乘，但是从运⾏速度上讲，它仍有潜⼒可挖，在采⽤了两

种优化技术后，(见程序5)，速度竟提⾼5倍，甚⾄超出笔者的预计。

第⼀种优化技术，将频繁调⽤的函数定义成inline函数，inline是C++关键字，如果使⽤纯C的编译器，可采⽤MACRO来代替。

第⼆种优化技术，将循环体内的代码展开，由⼀个循环步只做⼀次乘法，改为⼀个循环步做32次乘法。

实际速度：计算10000000！，程序4需要0.11667秒，程序5只需要0.02145秒。测试环境为迅驰1.7G,256M

RAM。关于程序优化⽅⾯的内容，将在后续⽂章专门讲述。下⾯是被优化后的部分代码，其余代码不变。

程序5的部分代码：

inlineshortGetExpBa2(doublea)//获得a

的阶码

{

//按照IEEE754浮点数格式，取得阶码，仅仅适⽤于Intel系列cpu

WORD

*pWord=(WORD*)(&a)+3;

WORDrank

=((*pWord&0x7fff)

>>4);

return

(short)(rank-0x3ff);

}

inlinedoubleGetMantissa(doublea)//获得a的尾数

{

//按照IEEE754浮点数格式，取得尾数，仅仅适⽤于Intel系列cpu

WORD

*pWord=(WORD*)(&a)+3;

*pWord&=

0x800f;//清除阶码

*pWord|=

0x3ff0;//重置阶码

return

a;

}

voidcalcFac(structbigNum*p,intn)

{

inti,t;

doublef,max_mantissa;

p->n1=1;p->n2=0;t=n-32;

for

(i=1;i<=t;i+=32)

{

p->n2+=GetExpBa2(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.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

*=(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+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+=GetExpBa2(p->n1);

p->n1=GetMantissa(p->n1);

p->n1*=(double)(i);

}

}

注1：10＾0，表⽰10的0次⽅

近似计算之⼆

Series可找到2个利⽤斯特林公式求阶乘的公式，⼀个是普通形式，⼀个是对数形式。前⼀个公式包含⼀个⽆穷级数和的形式，包含级数前

5项的公式为n!=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"

#define

PI3.97932384626433832795

#defineE2.74713527

structbigNum

{

doublen;//科学计数法表⽰的尾数部分

inte;

//科学计数法表⽰的指数部分,若⼀个bigNum为x,这x=n*10^e

};

constdoublea1[]=

{1.00,1.0/12.0,1.0/288,-139/51840,-571.0/2488320.0};

constdoublea2[]=

{1.0/12.0,-1.0/360,1.0/1260.0};

voidprintfResult(structbigNum*p,charbuff[])

{

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/n^2

-139/51840/n^3-571/2488320/n^4+…)

voidcalcFac1(structbigNum*p,intn)

{

doublelogx,

s,//级数的和

item;//级数的每⼀项

inti;

//x^n=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/n^2-139/51840/n^3-571/2488320/n^4+…)

for

(item=1.0,s=0.0,i=0;i

{

s+=item*a1[i];

item/=(double)n;//item=1/(n^i)

}

p->n*=s;

}

//ln(n!)=0.5*ln(2*PI)+(n+0.5)*ln(n)-n+(1/12/n-1/360/n^3+

1/1260/n^5-...)

voidcalcFac2(structbigNum*p,intn)

{

doublelogR;

doubles,//级数的和

item;//级数的每⼀项

int

i;

logR=0.5*log(2.0*PI)+((double)n+0.5)*log(n)-(double)n;

// s=(1/12/n-1/360/n^3+1/1260/n^5)

for

(item=1/(double)n,s=0.0,i=0;i

{

s+=item*a2[i];

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);

}

intmain(intargc,char*argv[])

{

structbigNumr;

charbuff[32];

intn;

printf("n=?");

scanf("%d",&n);

calcFac1(&r,n);//⽤第⼀种⽅法，计算n的阶乘

printfResult(&r,buff);//将结果转化⼀个字符串

printf("%d!=%sbymethod1n",n,buff);

calcFac2(&r,n);//⽤第⼆种⽅法，计算n的阶乘

printfResult(&r,buff);//将结果转化⼀个字符串

printf("%d!=%sbymethod2n",n,buff);

return0;

}