计算百万位e
讨论深入具体e和pi一样，是数学分析中使用非常广泛的无理数之一。千百年来，有无数人投入到计算e和pi的计算工作中。相对于计算pi，计算e的算法较为简单，而且复杂度也低于计算pi。以Pifast 4.3为例，在我的电脑上（PIV2.6）计算100万位e值，仅需1.23秒，而计算100万pi则需要3.8秒。本文首先借助HugeCalc，探讨计算e的各种算法。当然大家也可以自己写个计算e的程序，比拼一下谁的程序很快，谁的算法更优。
楼下使用的HugeCacl配置文件定义为：
NumOfCores= 1
SSE2Support = 1
CacheL1size = 8
CacheL2size = 512
CacheL3size = 0
Cacheburst= 64
[ 本帖最后由 liangbch 于 2008-4-24 17:24 编辑 ] 我讨论一下最简单，最容易理解的算法。
$e=1+ 1/1! + 1/2! + 1/3! + ... 1/n! $
欲达到K位有效数字，需使得 $n!> 10^K$
算法描述
$ e=1+ 1/{1!} + 1/{2!} + 1/{3!} + ... 1/{n!} $
以n=5为例，我们将各项通分
$e= 1/1$
$+1/(1*2) $
$+ 1/(1*2*3) $
$+ 1/(1*2*3*4) $
$+ 1/(1*2*3*4*5) $
容易得到算法
n=5,i=1, item=1, s=1;
for (i=1;i<=5;i++)
{
item /=i ;
s+=item;
} 简单算法也不很慢吧？
不过是个O(n^2)算法 本帖给出使用HugeCalc库，使用楼上的算法，计算e的一个实现。该算法是一个O(n^2)的算法，在我的电脑上，计算10万位e需要7.61秒。#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "../../HugeCalc/HugeCalc_API/CppAPI/Include/HugeCalc.h" // 公共接口
#include "../../HugeCalc/HugeCalc_API/CppAPI/Include/HugeInt.h" // 10进制系统
#pragma message( "automatic link to ../../HugeCalc/HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
#pragma comment( lib, "../../HugeCalc/HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
#define integer CHugeInt
//返回计算d位有效数字所需的最大的n
static UINT32 maxNForCalcE(int d)
{
double f=d*log(10);
double r=0;
UINT32 i=0;
do
{
i++; r += log(i);
}while (r<f);
return i;
}
void calc_e_1(integer *p,int len)
{
integer item;
integer r;
UINT32 n=maxNForCalcE(len);
UINT32 i;
item=1L;
item.DecLShift(len); //增加有效数字
r=1L;
for (i=1;i<=n;i++)
{
item /= i;
r += item;
}
*p=r;
}
void calc_e1_main()
{
int d,len;
integer e;
char *pBuff=NULL;
const char *pTmp;
char fileName;
FILE *fp=NULL;
printf("calculate e (simple version)\n");
printf("How many digital is needed?");
scanf("%d",&d);
len=(d+17)/9*9; //取9的倍数，可能稍微快一点
HugeCalc::EnableTimer( TRUE );
HugeCalc::ResetTimer();
calc_e_1(&e,len);
// 记录计算耗时
UINT32 u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_us );
printf("计算用时 %.6f 秒\n", (double)u32Timer_Calc/1000000.0);
//----------------------------------------
HugeCalc::EnableTimer( TRUE );
HugeCalc::ResetTimer();
pBuff=new char;
if (pBuff==NULL)
goto thisExit;
pTmp= e.GetStr(FS_NORMAL);
strcpy(pBuff+1,pTmp);
pBuff=pBuff;
pBuff='.'; //插入小数点
pBuff=0; //去掉多余的数字
sprintf(fileName,"e1_%d.txt",d);
fp=fopen(fileName,"wb");
fwrite(pBuff,d+2,1,fp);
fclose(fp); fp=NULL;
u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_us );
printf("输出到文件用时 %.6f 秒\n", (double)u32Timer_Calc/1000000.0);
thisExit:
if (fp!=NULL)
{
fclose(fp); fp=NULL;
}
if (pBuff!=NULL)
{
delete[] pBuff; pBuff=NULL;
}
} 本帖将给出计算e的另一个算法，这个算法复杂度同4#，只是4楼主要使用大数除以单精度整数的算法，而该算法主要使用大数乘以单精度整数的算法。该算法比4#更快，在我的电脑上计算10万e值，需要5.88秒。
算法描述
$e=1+ 1/{1!} + 1/{2!} + 1/{3!} + ... 1/{n!}$
以n=5为例，我们将各项通分
$ e= {1xx2xx3xx4xx5}/{5!}+{2xx3xx4xx5}/{5!}+{3xx4xx5}/{5!}+{4xx5}/{5!}+5/{5!} + 1/{5!}$
调整各项顺序，并提取公分母得
$e=(1/{5!})xx ( 1+5+ 5xx4 + 5xx4xx3 + 5xx4xx3xx2 + 5xx4xx3xx2xx1)$
若 e 表示为分数p/q,则
$p = 1 +5 +5xx4 + 5xx4xx3 + 5xx4xx3xx2 + 5xx4xx3xx2xx1$
$q=5xx4xx3xx2 xx1$
容易得到算法
n=5
step1:q=1,p=1
step2:q *= n,p +=q ,n--
step3:如果n>0,转step2,否则转step4
step4:将分数 p/q化为小数，得到e 本帖给出使用HugeCalc，使用楼上算法计算e的一个实现。#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "../../HugeCalc/HugeCalc_API/CppAPI/Include/HugeCalc.h" // 公共接口
#include "../../HugeCalc/HugeCalc_API/CppAPI/Include/HugeInt.h" // 10进制系统
#pragma message( "automatic link to ../../HugeCalc/HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
#pragma comment( lib, "../../HugeCalc/HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
#define integer CHugeInt
//计算e的稍复杂的版本，用到大数加法，以及大数乘法和大数除法
//该程序快于第一个版本，特别是要求精度较低时，优势更为明显。
//欲达到n位有效数字，需使得 n!> 10^n
//返回计算d位有效数字所需的最大的n
static UINT32 maxNForCalcE(int d)
{
double f=d*log(10);
double r=0;
UINT32 i=0;
do
{
i++; r += log(i);
}while (r<f);
return i;
}
void calc_e_2(integer *r,int len)
{
integer p,q;
UINT32 n=maxNForCalcE(len);
p=1; q=1;
do
{
q*=n;
p+=q;
n--;
}while (n>0);
p.DecLShift( q.GetDigits()-1); //增加有效数字
*r = p / q;
}
void calc_e2_main()
{
int d,len;
integer e;
char *pBuff=NULL;
const char *pTmp;
char fileName;
FILE *fp=NULL;
printf("calculate e (simple version2)\n");
printf("How many digital is needed?");
scanf("%d",&d);
len=(d+17)/9*9; //取9的倍数，可能稍微快一点
HugeCalc::EnableTimer( TRUE );
HugeCalc::ResetTimer();
calc_e_2(&e,len);
// 记录计算耗时
UINT32 u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_us );
printf("计算用时 %.6f 秒\n", (double)u32Timer_Calc/1000000.0);
//----------------------------------------
HugeCalc::EnableTimer( TRUE );
HugeCalc::ResetTimer();
pBuff=new char;
if (pBuff==NULL)
goto thisExit;
pTmp= e.GetStr(FS_NORMAL);
strcpy(pBuff+1,pTmp);
pBuff=pBuff;
pBuff='.'; //插入小数点
pBuff=0; //去掉多余的数字
sprintf(fileName,"e2_%d.txt",d);
fp=fopen(fileName,"wb");
fwrite(pBuff,d+2,1,fp);
fclose(fp); fp=NULL;
u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_us );
printf("输出到文件用时 %.6f 秒\n", (double)u32Timer_Calc/1000000.0);
thisExit:
if (fp!=NULL)
{
fclose(fp); fp=NULL;
}
if (pBuff!=NULL)
{
delete[] pBuff; pBuff=NULL;
}
} 本帖将给出一个快速的算法，这个算法能够充分利用大数乘法，计算速度比前两个算法快的多（在PIV2.6上，计算10万位e需0.66秒，计算100万位e需要8.35秒）
算法详述：
$e=1+1/{1!} + 1/{2!} + 1/{3!}+ 1/{4!} + ... 1/{N!}$ (一)
现在我们将讨论如何化解
$ 1/{1!} + 1/{2!} + 1/{3!}+ 1/{4!} + ... 1/{N!}$ (二)
为$p/q$
以计算 $ 1/{1!} + 1/{2!} + 1/{3!} + 1/{4!} +1/{5!} + 1/{6!} + 1/{7!} + 1/{8!}$ 为例
将这个多项式分成2部分得
第一部分为：$( 1/{0!}) xx( 1/1 + 1/{1xx2} + 1/{1xx2xx3}+1/{1xx2xx3xx4})$
第二部分为：$1/{4!}(1/5+1/{5xx6}+1/{5xx6xx7}+1/{5xx6xx7xx8})$
暂不考虑系数$1/{0!}$和$1/{4!}$，则问题转化为求
$1/n+ 1/{nxx(n+1)}+1/{n(n+1)(n+2)} ...+1/{n(n+1)(n+2)...m}$ (三)
的形式。易看出，若这个多项的的计算结果为$ p/q $,则 $q=nxx(n+1)xx(n+2)xx...(m) $
当m-n较小时，我们可以按照calc_Frac_Sum中的方法直接计算,
当m-n较大时，我们可以将这个多项式再分成2段，
一段为从n 到 (n+m)/2,另一端为一段为从(n+m)/2+1 到 m
令r=(n+m)/2、s=(n+m)/2+1
为了计算(三)，我们先分别计算
$1/n + 1/{(n(n+1))} + ...1/{(n(n+1)...r)}$ (四)
$1/s + 1/{(s(s+1))} + ...1/{(s(s+1)...m)} $(五)
令(三)的计算结果为$p/q$，
(四)的计算结果为${p1}/{q1}$，
(五)的计算结果为${p2}/{q2}$,
则：$1/n+ 1/{nxx(n+1)}+1/{n(n+1)(n+2)}+...+1/{n(n+1)(n+2)...r}$
$+1/{n(n+1)(n+2)...s} + ... + 1/{n(n+1)(n+2)...m}$
$= p/q$
$= {p1}/{q1} + (1/{nxx(n+1)xx(n+2)xx..r}) xx {p2}/{q2}$
$= {p1}/{q1} + 1/{q1} xx {p2}/{q2}$
$= {p1xxq2}/{q1xxq2} + {p2}/{q1xxq2}$
$= {p1xxq2+p2}/{q1xxq2}$
所以：$p=p1xxq2+p2, q=q1xxq2$
下面，我们以$(1/5+1/{5xx6} + 1/{5xx6xx7} + 1/{5xx6xx7xx8})$ 说明计算过程：
$ p/q= 1/5 + 1/{5xx6} + 1/{5xx6xx7} + 1/{5xx6xx7xx8}$
$= (1/5 + 1/{5xx6}) + (1/{5xx6xx7} + 1/{5xx6xx7xx8})$
$= ((6+1)/(5xx6)) + 1/(5xx6) xx ( 1/7 + 1/{7xx8} )$
$= ((6+1)/(5xx6)) + 1/(5xx6) xx ( {8+1}/{7xx8} ) $
$= (7/30) + 1/(30) xx ( 9/56)$
$= (7 xx 56 + 9 )/(30xx56)$
按此法计算，使用递归方式，所有的多项式终将合并，这时n=1,m=N,因此(二)将变为(三)
则有e
$=1+1/{1!}+1/{2!}+1/{3!} + ...1/{N!}$
$=1+p/q$
$ =(p+q)/q$
最后我们作一次除法，将分数(p+q)/q化为小数即可 看不懂哦 本帖给出使用HugeCalc，使用7#楼算法计算e的一个实现,使用该程序，在PIV2.6上，计算10万位e需0.412秒，计算100万位e需要5.06秒。#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "../../../HugeCalc_API/CppAPI/Include/HugeCalc.h" // 公共接口
#include "../../../HugeCalc_API/CppAPI/Include/HugeInt.h" // 10进制系统
#pragma message( "automatic link to ../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
#pragma comment( lib, "../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
#define integer CHugeInt
#define USE_HUGE_CALC_DIV
#if 0
# define MODIFIED_BY_GxQ
#endif
#define MIN_SPLIT 32
static DWORD g_divTime=0;
// 计算e的一个快速的版本，采用2分法，尽量使用以及大数乘法，
// 收尾阶段采用大数除法
//欲达到n位有效数字，需使得 n!> 10^n
//返回计算d位有效数字所需的最大的n
static UINT32 maxNForCalcE(int d)
{
double f=d*log(10);
double r=0;
UINT32 i=0;
//loop untill log(n!)>log(10^d)
do
{
i++; r += log(i);
}while (r<f);
return i;
}
// calc 1/n1 + 1/(n1*(n1+1)) + 1/(n1*(n1+1)+(n1+2))+ 1/(n1*(n1+1)+(n1+2)+...m)
// 计算结果为 p/q
// for example:
// 1/6 + 1/6*7 + 1/6*7*8 + 1/6*7*8*9
// = 7*8*9/6*7*8*9 + 8*9/6*7*8*9 + 9/6*7*8*9 + 1/ 6*7*8*9
// calculate
// q=6*7*8*9, p=1 + 9+ 8*9 + 7*8*9
// steps:
// step1: p=1, q=1;
// step2: q*=9, p+=q ,q=9, q=1+9
// step3: q*=8, p+=q ,q=9*8, q=1+8*9
// step4: q*=7, p+=q ,q=9*8*7, q=1+ 8*9 + 7*8*9
// step5: q*=6
static void calc_Frac_Sum(UINT32 n1,UINT32 n2, integer &p,integer &q)
{
p=(UINT32)1;
q=(UINT32)1;
for ( UINT32 i=n2;i>n1;i--)
{
q*=i;
p+=q;
}
q*=i;
}
#ifndef MODIFIED_BY_GxQ
void splitCalc(UINT32 n1,UINT32 n2, integer &p,integer &q)
{
if (n2-n1<=MIN_SPLIT)
{
calc_Frac_Sum(n1,n2,p,q);
}
else
{
integer p1,q1,p2,q2;
splitCalc( n1, (n1+n2)/2, p1,q1);
splitCalc( (n1+n2+2)/2 ,n2, p2,q2);
p=p1*q2 + p2;
q=q1*q2;
}
}
#else
integer tmp;
void splitCalc(UINT32 n1,UINT32 n2, integer &p,integer &q)
{
if (n2-n1<=MIN_SPLIT)
{
calc_Frac_Sum(n1,n2,p,q);
}
else
{
integer q1,p2;
splitCalc( n1, (n1+n2)/2, p, q1);
splitCalc( (n1+n2+2)/2 ,n2, p2,q );
// ( p *= q ) += p2;
tmp.Swap( p );
p.Mul( tmp, q ) += p2;
// q *= q1;
tmp.Swap( q );
q.Mul( tmp, q1 );
}
}
#endif
void calc_e_3(integer *r, int len )
{
UINT32 n;
int shift;
integer p,q;
n= maxNForCalcE(len)+2;
splitCalc(1,n,p,q);// p/q= 1/1!+1/2!+1/3!+1/4!+...
p+=q;// p/q +=1 // e= p/q= 1 + 1/1!+1/2!+1/3!+1/4!+...
shift=(len+8)/9*9;//get upbounds get times of 9
p.DecLShift(shift); //增加有效数字
g_divTime=GetTickCount();
#ifdef USE_HUGE_CALC_DIV
#ifndef MODIFIED_BY_GxQ
p=p /q;
#else
p /= q;
#endif
#else
{
void inverse( const integer *s, integer *pResult,int len);
integer t;
inverse(&q,&t,len);
p *=t;
}
#endif
g_divTime=GetTickCount()-g_divTime;
*r=p;
}
void calc_e3_main()
{
int len;
integer e;
char *pBuff=NULL;
const char *pTmp;
char fileName;
FILE *fp=NULL;
printf("calculate e (fast version1)\n");
#ifdef MODIFIED_BY_GxQ
printf("defined 'MODIFIED_BY_GxQ'\n");
#endif
while (true)
{
printf("How many digital is needed(0 will exit)?");
scanf("%d",&len);
if (len==0)
break;
HugeCalc::EnableTimer( TRUE );
HugeCalc::ResetTimer();
calc_e_3(&e,len);
// 记录计算耗时
UINT32 u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_us );
printf("计算用时 %.6f 秒\n", (double)u32Timer_Calc/1000000.0);
if (g_divTime>=2)
{
printf("其中除法用时 %f秒\n", (double)g_divTime/1000.00);
}
//----------------------------------------
HugeCalc::EnableTimer( TRUE );
HugeCalc::ResetTimer();
pBuff=new char;
if (pBuff==NULL)
goto thisExit;
pTmp= e.GetStr(FS_NORMAL);
strcpy(pBuff+1,pTmp);
pBuff=pBuff;
pBuff='.'; //插入小数点
pBuff=0; //去掉多余的数字
sprintf(fileName,"e3_%d.txt",len);
FILE *fp=fopen(fileName,"wb");
fwrite(pBuff,len+2,1,fp);
fclose(fp); fp=NULL;
delete[] pBuff; pBuff=NULL;
u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_us );
printf("输出到文件用时 %.6f 秒\n", (double)u32Timer_Calc/1000000.0);
}
thisExit:
if (fp!=NULL)
{
fclose(fp); fp=NULL;
}
if (pBuff!=NULL)
{
delete[] pBuff; pBuff=NULL;
}
}
[ 本帖最后由 liangbch 于 2008-4-17 14:34 编辑 ] :)
我想
我想
是不是有更精彩更复杂的算法
比如类似阶乘的？