阶乘之计算从入门到精通－近似计算之二

   摘要：本文仅讨论精度为16位有效数字以内近似计算，和上一篇文章不同，它采用一个叫做斯特林的数学公式来计算。它能够计算很大的数的阶乘，但精度较低。

在《阶乘之计算从入门到精通－近似计算之一》中，我们采用两个数来表示中间结果，使得计算的范围扩大到1千万，并可0.02秒内完成10000000！的计算。在保证接近16位有效数字的前提下，有没有更快的算法呢。很幸运，有一个叫做Stirling的公式，它可以快速计算出一个较大的数的阶乘，而且数越大，精度越高。有 http://mathworld.wolfram.com查找Stirling’s Series可找到2个利用斯特林公式求阶乘的公式，一个是普通形式，一个是对数形式。前一个公式包含一个无穷级数和的形式，包含级数前5项的公式为n !=sqrt(2*PI*n)*(n/e)ⁿ*(1+ 1/12/n+ 1/288/n² –139/51840/n³-571/2488320/n⁴+…),这里的PI为圆周率，e为自然对数的底。后一个公式也为一个无穷级数和的形式，一般称这个级数为斯特林(Stirling)级数,包括级数前3项的n!的对数公式为： ln(n!)=0.5*ln(2*PI)+(n+0.5)*ln(n)-n+(1/12/n -1/360/n³ + 1/1260/n⁵)-…,下面我们分别用这两个公式计算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;
}

 

liangbch@263.net ，版权所有，转载请注明出处。