摘要：本文仅讨论精度为16位有效数字以内近似计算，和上一篇文章不同，它采用一个叫做斯特林的数学公式来计算。它能够计算很大的数的阶乘，但精度较低。
在《阶乘之计算从入门到精通－近似计算之一》中，我们采用两个数来表示中间结果，使得计算的范围扩大到1千万，并可0.02秒内完成10000000！的计算。在保证接近16位有效数字的前提下，有没有更快的算法呢。很幸运，有一个叫做Stirling的公式，它可以快速计算出一个较大的数的阶乘，而且数越大，精度越高。有
http://mathworld.wolfram.com查找Stirling’s 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 PI 3.1415926535897932384626433832795
#define E 2.7182818284590452353602874713527
struct bigNum
{
double n; //科学计数法表示的尾数部分
int e; //科学计数法表示的指数部分,若一个bigNum为x,这x=n * 10^e
};
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 printfResult(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/n^2 -139/51840/n^3-571/2488320/n^4+…)
void calcFac1(struct bigNum *p,int n)
{
double logx,
s, //级数的和
item; //级数的每一项
int i;
// 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<sizeof(a1)/sizeof(double);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 -...)
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/n^3 + 1/1260/n^5)
for (item=1/(double)n,s=0.0,i=0;i<sizeof(a2)/sizeof(double);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);
}
int main(int argc, char* argv[])
{
struct bigNum r;
char buff[32];
int n;
printf("n=?");
scanf("%d",&n);
calcFac1(&r,n); //用第一种方法，计算n的阶乘
printfResult(&r,buff); //将结果转化一个字符串
printf("%d!=%s by method 1\n",n,buff);
calcFac2(&r,n); //用第二种方法，计算n的阶乘
printfResult(&r,buff); //将结果转化一个字符串
printf("%d!=%s by method 2\n",n,buff);
return 0;
}
| 网站地图 |