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

来源：互联网发布：java中compare2 编辑：程序博客网时间：2024/04/28 15:40

在＜阶乘之计算从入门到精通－菜鸟篇＞中提到，使用double型数来计算阶乘，当n>170,计算结果就超过double数的最大范围而发生了溢出，故当n＞170时，就不能用这个方法来计算阶乘了，果真如此吗？No,只要肯动脑筋,办法总是有的。

通过windows计算器，我们知道，171！＝1.2410180702176678234248405241031e+309，虽然这个数不能直接用double型的数来表示，但我们可以用别的方法来表示。通过观察这个数，我们发现，这个数的表示法为科学计算法，它用两部分组成，一是尾数部分1.2410180702176678234248405241031，另一个指数部分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" #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;     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(struct 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 buff[32];     int n;          printf("n=?");      scanf("%d",&n);     calcFac(&r,n);           //计算n的阶乘     printfResult(&r,buff);   //将结果转化一个字符串     printf("%d!=%s/n",n,buff);     return 0;}

以上代码中的数的表示方式为：数的值等于尾数乘以 10 ^ 指数部分，当然我们也可以表示为：尾数乘以 2 ^ 指数部分，这们将会带来这样的好处，在调整尾数部分和指数部分时，不用除法，可以依据浮点数的格式直读取阶码和修改阶码(上文提到的指数部分的标准称呼)，同时也可在一定程序上减少误差。为了更好的理解下面的代码，有必要了解一下浮点数的格式。浮点数主要分为32bit单精度和64bit双精度两种。本方只讨论64bit双精度(double型)浮点数的格式，一个double型浮点数包括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这个16bit数，最高位清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; struct 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浮点数格式，取得尾数，仅仅适用于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);      }      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[32];int n;printf("n=?"); scanf("%d",&n);calcFac(&r,n);           //计算n的阶乘printfResult(&r,buff);   //将结果转化一个字符串printf("%d!=%s/n",n,buff);return 0;}

程序优化的威力：

程序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的部分代码：

inline short GetExpBase2(double a) // 获得 a 的阶码{    // 按照IEEE 754浮点数格式，取得阶码，仅仅适用于Intel 系列 cpu     WORD *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(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.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 += GetExpBase2(p->n1);        p->n1 = GetMantissa(p->n1);        p->n1 *=(double)(i);   }}

注1：10＾0，表示10的0次方