在
＜阶乘之计算从入门到精通－菜鸟篇＞中提到，使用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次方
liangbch@263.net ，版权所有，转载请注明出处。
8万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
