近似计算之一
在<阶乘之计算从入门到精通-菜鸟篇>中提到,使用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,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",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次方
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。