序

大数阶乘的计算是一个有趣的话题，从中学生到大学教授，许多人都投入到这个问题的探索和研究之中，并发表了他们自己的研究成果。如果你用阶乘作关键字在google上搜索，会找到许多此类文章，另外，如果你使用google学术搜索，也能找到一些计算大数阶乘的学术论文。但这些文章和论文的深度有限，并没有给出一个高速的算法和程序。

 

我和许多对大数阶乘感兴趣的人一样，很早就开始编制大数阶乘的程序。从2000年开始写第一个大数阶乘程序算起，到现在大约己有6-7年的时光，期间我写了多个版本的阶乘计算器，在阶乘计算器的算法探讨和程序的编写和优化上，我花费了很大的时间和精力，品尝了这一过程中的种种甘苦，我曾因一个新算法的实现而带来速度的提升而兴奋，也曾因费了九牛二虎之力但速度反而不及旧的版本而懊恼，也有过因解一个bug而通宵敖夜的情形。我写的大数阶乘的一些代码片段散见于互联网络，而算法和构想则常常萦绕在我的脑海。自以为，我对大数阶乘计算器的算法探索在深度和广度上均居于先进水平。我常常想，应该写一个关于大数阶乘计算的系列文章，一为整理自己的劳动成果，更重要的是可以让同行分享我的知识和经验，也算为IT界做一点儿贡献吧。

　

　 我的第一个大数阶乘计算器始于2000年，那年夏天，我买了一台电脑，开始在家专心学习VC，同时写了我的第一个VC程序，一个仿制windows界面的计算器。该计算器的特点是高精度和高速度，它可以将四则运算的结果精确到6万位以内，将三角、对数和指数函数的结果精确到300位以内，也可以计算开方和阶乘等。当时，我碰巧看到一个叫做实用计器的软件。值得称颂的是，该计算器的作者姜边竟是一个高中生，他的这个计算器功能强大，获得了各方很高的评价。该计算的功能之一是可以计算9000以内阶乘的精确值，且速度很快。在佩服之余，也激起我写一个更好更快的程序的决心，经过几次改进，终于使我的计算器在做阶乘的精确计算时 (以9000!为例)，可比实用计算器快10倍。而当精确到30多位有效数字时，可比windows自带的计算器快上7500倍。其后的2001年1月，我在csdn上看到一个贴子，题目是“有谁可以用四行代码求出1000000的阶乘”，我回复了这个贴子，给出一个相对简洁的代码，这是我在网上公布的第一个大数阶乘的程序。这一时期，可以看作是我写阶乘计算器的第一个时期。

 

　　我写阶乘计算器的第二个时期始于2003年9月，在那时我写了一组专门计算阶乘的程序，按运行速度来分，分为三个级别的版本，初级版、中级版和高级版。初级版本的算法许多人都能想到，中级版则采用大数乘以大数的硬乘法，高级版本在计算大数乘法时引入分治法。期间在csdn社区发了两个贴子，“擂台赛：计算n!(阶乘)的精确值，速度最快者2000分送上”和“擂台赛：计算n!(阶乘)的精确值，速度最快者2000分送上（续）”。其高级算的版本完成于2003年11月。此后，郭先强于2004年5月10日也发表了系列贴子，“擂台：超大整数高精度快速算法”、“擂台：超大整数高精度快速算法(续)”和“擂台：超大整数高精度快速算法(续2)”, 该贴重点展示了大数阶乘计算器的速度。这个贴子一经发表即引起了热列的讨论，除了我和郭先强先生外，郭雄辉也写了同样功能的程序，那时，大家都在持续改进自己的程序，看谁的程序更快。初期，郭先强的稍稍领先，中途郭子将apfloat的代码应用到阶乘计算器中，使得他的程序胜出，后期（2004年8月后）在我将程序作了进一步的改进后，其速度又稍胜于他们。在这个贴子中，arya提到一个开放源码的程序，它的大数乘法采用FNT＋CRT(快速数论变换＋中国剩余定理)。郭雄辉率先采用apflot来计算大数阶乘，后来郭先强和我也参于到apfloat的学习和改进过程中。在这点上，郭先强做得非常好，他在apfloat的基础上，成功地优化和改时算法，并应用到大数阶乘计算器上，同时他也将FNT算法应用到他的＜超大整数高精度快速算法库＞中，并在2004.10.18正式推出V3.0.2.1版。此后，我在2004年9月9日也完成一个采用FNT算法的版本，但却不及郭先强的阶乘计算器快。那时，郭先强的程序是我们所知的运算速度最快的阶乘计算器，其速度超过久负盛名的数学软件Mathematica和Maple,以及通用高精度算法库GMP。

　　

　　我写阶乘计算器的第三个时间约开始于2006年，在2005年8月收到北大刘楚雄老师的一封e-mail，他提到了他写的一个程序在计算阶乘时比我们的更快。这使我非常吃惊，在询问后得知，其核心部分使用的是ooura写的FFT函数。ooura的FFT代码完全公开，是世界上运行的最快的FFT程序之一，从这点上，再次看到了我们和世界先进水平的差距。佩服之余，我决定深入学习FFT算法，看看能否写出和ooura速度相当或者更快的程序，同时一个更大计划开始形成，即写一组包括更多算法的阶乘计算器，包括使用FFT算法的终极版和使用无穷级数的stirling公式来计算部分精度的极速版，除此之外，我将重写和优化以前的版本，力争使速度更快，代码更优。这一计划的进展并不快，曾一度停止。

　　

　　目前，csdn上blog数量正在迅速地增加，我也萌生了写blog的计划，借此机会，对大数阶乘之计算作一个整理，用文字和代码详述我的各个版本的算法和实现，同时也可能分析一些我在网上看到的别人写的程序，当然在这一过程中，我会继续编写未完成的版本或改写以前己经实现的版本，争取使我公开的第一份代码都是精品，这一过程可能是漫长的，但是我会尽力做下去。

菜鸟篇

程序1，一个最直接的计算阶乘的程序

 

#include "stdio.h"

#include "stdlib.h"

 

int main(int argc, char* argv[])

{

long i,n,p;

printf("n=?");

scanf("%d",&n);

p=1;

for (i=1;i<=n;i++)

p*=i;

printf("%d!=%d/n",n,p);

return 0;

}

 

程序2，稍微复杂了一些，使用了递归，一个c++初学者写的程序

 

#include <iostream.h>
long int fac(int n);
void main()
{
int n;
cout<<"input a positive integer:";
cin>>n;
long fa=fac(n);
cout<<n<<"! ="<<fa<<endl;
}
long int fac(int n)
{
long int p;
if(n==0) p=1;
else
p=n*fac(n-1);
return p;
}

程序点评，这两个程序在计算12以内的数是正确，但当n>12,程序的计算结果就完全错误了，单从算法上讲，程序并没有错，可是这个程序到底错在什么地方呢？看来程序作者并没有意识到，一个long型整数能够表示的范围是很有限的。当n>=13时，计算结果溢出，在C语言，整数相乘时发生溢出时不会产生任何异常，也不会给出任何警告。既然整数的范围有限，那么能否用范围更大的数据类型来做运算呢？这个主意是不错，那么到底选择那种数据类型呢？有人想到了double类型，将程序1中long型换成double类型，结果如下:

 

#include "stdio.h"

#include "stdlib.h"

 

int main(int argc, char* argv[])

{

double i,n,p;

printf("n=?");

scanf("%lf",&n);

p=1.0;

for (i=1;i<=n;i++)

p*=i;

printf("%lf!=%.16g/n",n,p);

return 0;

}

 

运行这个程序，将运算结果并和windows计算器对比后发现，当于在170以内时，结果在误差范围内是正确。但当N>=171,结果就不能正确显示了。这是为什么呢？和程序1类似，数据发生了溢出，即运算结果超出的数据类型能够表示的范围。看来C语言提供的数据类型不能满足计算大数阶乘的需要，为此只有两个办法。1.找一个能表示和处理大数的运算的类库。2.自己实现大数的存储和运算问题。方法1不在本文的讨论的范围内。本系列的后续文章将围绕方法2来展开。

 

大数的表示

 

 

1.大数，这里提到的大数指有效数字非常多的数，它可能包含少则几十、几百位十进制数，多则几百万或者更多位十进制数。有效数字这么多的数只具有数学意义，在现实生活中，并不需要这么高的精度，比如银河系的直径有10万光年，如果用原子核的直径来度量，31位十进制数就可使得误差不超过一个原子核。

 

2.大数的表示：

2.1定点数和浮点数

　我们知道，在计算机中，数是存贮在内存(RAM)中的。在内存中存储一个数有两类格式，定点数和浮点数。定点数可以精确地表示一个整数，但数的范围相对较小，如一个32比特的无符号整数可表示0－4294967295之间的数，可精确到9－10位数字(这里的数字指10进制数字，如无特别指出，数字一律指10进制数字),而一个8字节的无符号整数则能精确到19位数字。浮点数能表示更大的范围，但精度较低。当表示的整数很大的，则可能存在误差。一个8字节的双精度浮点数可表示2.22*10^-308到 1.79*10^308之间的数,可精确到15－16位数字.

　

　2.2日常生活中的数的表示：

　对于这里提到的大数，上文提到的两种表示法都不能满足需求。为此，必需设计一种表示法来存储大数。我们以日常生活中的十进制数为例，看看是如何表示的。如一个数N被写成"12345",则这个数可以用一个数组a来表示，a[0]=1,　a[1]=2,　a[2]=3,　a[3]=4,　a[4]=5,这时数N=　a[4]*10^0　+a[3]*10^1　+a[2]*10^2　+a[1]*10^3　+a[0]*10^4,　(10^4表示10的4次方，下同),10^i可以叫做权，在日常生活中，a[0]被称作万位，也说是说它的权是10000，类似的，a[1]被称作千位，它的权是1000。

 

2.3　大数在计算机语言表示：

　　在日常生活中，我们使用的阿拉伯数字只有0－9共10个，按照书写习惯，一个字符表示1位数字。计算机中，我们常用的最小数据存储单位是字节，C语言称之为char,多个字节可表示一个更大的存储单位。习惯上，两个相邻字节组合起来称作一个短整数，在32位的C语言编译器中称之为short,汇编语语言一般记作word,4个相邻的字节组合起来称为一个长整数，在32位的C语言编译器中称之为long,汇编语言一般记作DWORD。在计算机中，按照权的不同，数的表示可分为两种，2进制和10进制，严格说来，应该是2^k进制和10^K进制，前者具占用空间少，运算速度快的优点。后者则具有容易显示的优点。我们试举例说明：

例1：若一个大数用一个长为len的short型数组A来表示，并采用权从大到小的顺序依次存放，数N表示为A[0] * 65536^(len-1)+A[1] * 65536^(len-2)+...A[len-1] * 65536^0,这时65536称为基，其进制2的16次方。

例2：若一个大数用一个长为len的short型数组A来表示并采用权从大到小的顺序依次存放，数N=A[0] * 10000^(len-1)+A[1] * 10000^(len-2)+...A[len-1] * 10000^0,这里10000称为基，其进制为10000，即：10^4，数组的每个元素可表示4位数字。一般地，这时数组的每一个元素为小于10000的数。类似的，可以用long型数组，基为2＾32＝4294967296来表示一个大数; 当然可以用long型组，基为1000000000来表示，这种表示法，数组的每个元素可表示9位数字。当然，也可以用char型数组，基为10。最后一种表示法，在新手写的计算大数阶乘程序最为常见，但计算速度却是最慢的。使用更大的基，可以充分发挥CPU的计算能力，计算量将更少，计算速度更快，占用的存储空间也更少。

 

2.4 大尾序和小尾序，我们在书写一个数时，总是先写权较大的数字，后写权较小的数字，但计算机中的数并不总是按这个的顺序存放。小尾（Little Endian）就是低位字节排放在内存的低端，高位字节排放在内存的高端。例如对于一个4字节的整数0x12345678,将在内存中按照如下顺序排放, Intel处理器大多数使用小尾(Little Endian)字节序。

Address[0]: 0x78

Address[1]: 0x56

Address[2]: 0x34

Address[3]:0x12

大尾（Big Endian）就是高位字节排放在内存的低端，低位字节排放在内存的高端。例如对于一个4字节的整数0x12345678,将在内存中按照如下顺序排放, Motorola处理器大多数使用大尾(Big Endian)字节序。

Address[0]: 0x12

Address[1]: 0x34

Address[2]: 0x56

Address[3]:0x78

　类似的，一个大数的各个元素的排列方式既可以采用低位在前的方式，也可以采用高位在前的方式，说不上那个更好，各有利弊吧。我习惯使用高位在前的方式。　　　

 

　2.5 不完全精度的大数表示：

　尽管以上的表示法可准确的表示一个整数，但有时可能只要求计算结果只精确到有限的几位。如用 windows自带的计算器计算1000的阶乘时，只能得到大约32位的数字，换名话说，windows计算器的精度为32位。1000的阶乘是一个整数，但我们只要它的前几位有效数字，象windows计算器这样，只能表示部分有效数字的表示法叫不完全精度，不完全精度不但占用空间省，更重要的是，在只要求计算结果为有限精度的情况下，可大大减少计算量。大数的不完全精度的表示法除了需要用数组存储有数数字外，还需要一个数来表示第一个有效数字的权，10的阶乘约等于4.023872600770937e+2567，则第一个有效数字的权是10^2567，这时我们把2567叫做阶码。在这个例子中，我们可以用一个长为16的char型数组和一个数来表示，前者表示各位有效数字，数组的各个元素依次为:4,0,2,3,8,7,2,6,0,0,7,7,0,9,3,7，后者表示阶码，值为2567。

 

 

2.6 大数的链式存储法

如果我们搜索大数阶乘的源代码，就会发现，有许多程序采用链表存储大数。尽管这种存储方式能够表示大数，也不需要事先知道一个特定的数有多少位有效数字，可以在运算过程中自动扩展链表长度。但是，如果基于运算速度和内存的考虑，强烈不建议采用这种存储方式，因为：

 

 

1. 这种存储方式的内存利用率很低。基于大数乘法的计算和显示，一般需要定义双链表，假如我们用1个char表示1位十进制数，则可以这样定义链表的节点：

struct _node

{

struct _node* pre;

struct _node* next;

char n;

};

当编译器采用默认设置，在通常的32位编译器，这个结构体将占用12字节。但这并不等于说，分配具有1000个节点的链表需要1000*12字节。不要忘记，操作系统或者库函数在从内存池中分配和释放内存时，也需要维护一个链表。实验表明，在VC编译的程序，一个节点总的内存占用量为 sizeof(struct _node) 向上取16的倍数再加8字节。也就是说，采用这种方式表示n位十进制数需要 n*24字节，而采用1个char型数组仅需要n字节。

 

 

2采用链表方式表示大数的运行速度很慢.

2.1如果一个大数需要n个节点，需要调用n次malloc(C)或new(C++)函数，采用动态数组则不要用调用这么多次malloc.

 

2.2 存取数组表示的大数比链表表示的大数具有更高的cache命中率。数组的各个元素的地址是连续的，而链表的各个节点在内存中的地址是不连续的，而且具有更大的数据量。因此前者的cache的命中率高于后者，从而导致运行速度高于后者。

 

2.3对数组的顺序访问也比链表快，如p1表示数组当前元素的地址，则计算数组的下一个地址时一般用p1++,而对链表来说则可能是p2=p2->next,毫无疑问，前者的执行速度更快。

 

 

近似计算之一

 

在＜阶乘之计算从入门到精通－菜鸟篇＞中提到，使用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