WV.24-大数阶乘算法4-近似计算之二

近似计算之二

 

在《阶乘之计算从入门到精通－近似计算之一》中，我们采用两个数来表示中间结果，使得计算的范围扩大到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;
}