计算百万位e

liangbch

e和$pi$一样，是数学分析中使用非常广泛的无理数之一。千百年来，有无数人投入到计算e和$pi$的计算工作中。相对于计算$pi$，计算e的算法较为简单，而且复杂度也低于计算$pi$。以Pifast 4.3为例，在我的电脑上（PIV2.6）计算100万位e值，仅需1.23秒，而计算100万$pi$则需要3.8秒。
   本文首先借助HugeCalc，探讨计算e的各种算法。当然大家也可以自己写个计算e的程序，比拼一下谁的程序很快，谁的算法更优。

  楼下使用的HugeCacl配置文件定义为：
[UserSetCPUID]
NumOfCores  =        1
SSE2Support =        1
CacheL1size =        8
CacheL2size =        512
CacheL3size =        0
Cacheburst  =        64

[ 本帖最后由 liangbch 于 2008-4-24 17:24 编辑 ]

liangbch

我讨论一下最简单，最容易理解的算法。

  $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)算法

liangbch

本帖给出使用HugeCalc库，使用楼上的算法，计算e的一个实现。该算法是一个O(n^2)的算法，在我的电脑上，计算10万位e需要7.61秒。

1. #include <stdlib.h>
2. #include <stdio.h>
3. #include <math.h>
5. #include "../../HugeCalc/HugeCalc_API/CppAPI/Include/HugeCalc.h" // 公共接口
6. #include "../../HugeCalc/HugeCalc_API/CppAPI/Include/HugeInt.h" // 10进制系统
8. #pragma message( "automatic link to ../../HugeCalc/HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
9. #pragma comment( lib, "../../HugeCalc/HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
11. #define integer CHugeInt
13. //返回计算d位有效数字所需的最大的n
14. static UINT32 maxNForCalcE(int d)
15. {
16. double f=d*log(10);
17. double r=0;
18. UINT32 i=0;
19. do
20. {
21. i++; r += log(i);
22. }while (r<f);
23. return i;
24. }
26. void calc_e_1(integer *p,int len)
27. {
28. integer item;
29. integer r;
30. UINT32 n=maxNForCalcE(len);
31. UINT32 i;
33. item=1L;
34. item.DecLShift(len); //增加有效数字
36. r=1L;
37. for (i=1;i<=n;i++)
38. {
39. item /= i;
40. r += item;
41. }
42. *p=r;
43. }
45. void calc_e1_main()
46. {
47. int d,len;
48. integer e;
49. char *pBuff=NULL;
50. const char *pTmp;
51. char fileName[32];
52. FILE *fp=NULL;
54. printf("calculate e (simple version)\n");
55. printf("How many digital is needed?");
56. scanf("%d",&d);
57. len=(d+17)/9*9; //取9的倍数，可能稍微快一点
59. HugeCalc::EnableTimer( TRUE );
60. HugeCalc::ResetTimer();
62. calc_e_1(&e,len);
64. // 记录计算耗时
65. UINT32 u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_us );
66. printf("计算用时 %.6f 秒\n", (double)u32Timer_Calc/1000000.0);
68. //----------------------------------------
69. HugeCalc::EnableTimer( TRUE );
70. HugeCalc::ResetTimer();
72. pBuff=new char[e.GetDigits()+3];
73. if (pBuff==NULL)
74. goto thisExit;
76. pTmp= e.GetStr(FS_NORMAL);
77. strcpy(pBuff+1,pTmp);
79. pBuff[0]=pBuff[1];
80. pBuff[1]='.'; //插入小数点
81. pBuff[d+2]=0; //去掉多余的数字
83. sprintf(fileName,"e1_%d.txt",d);
84. fp=fopen(fileName,"wb");
85. fwrite(pBuff,d+2,1,fp);
86. fclose(fp); fp=NULL;
88. u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_us );
89. printf("输出到文件用时 %.6f 秒\n", (double)u32Timer_Calc/1000000.0);
91. thisExit:
92. if (fp!=NULL)
93. {
94. fclose(fp); fp=NULL;
95. }
96. if (pBuff!=NULL)
97. {
98. delete[] pBuff; pBuff=NULL;
99. }
100. }

复制代码

liangbch

本帖将给出计算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

liangbch

本帖给出使用HugeCalc，使用楼上算法计算e的一个实现。

1. #include <stdlib.h>
2. #include <stdio.h>
3. #include <math.h>
5. #include "../../HugeCalc/HugeCalc_API/CppAPI/Include/HugeCalc.h" // 公共接口
6. #include "../../HugeCalc/HugeCalc_API/CppAPI/Include/HugeInt.h" // 10进制系统
8. #pragma message( "automatic link to ../../HugeCalc/HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
9. #pragma comment( lib, "../../HugeCalc/HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
12. #define integer CHugeInt
14. //计算e的稍复杂的版本，用到大数加法，以及大数乘法和大数除法
15. //该程序快于第一个版本，特别是要求精度较低时，优势更为明显。
17. //欲达到n位有效数字，需使得 n!> 10^n
18. //返回计算d位有效数字所需的最大的n
19. static UINT32 maxNForCalcE(int d)
20. {
21. double f=d*log(10);
22. double r=0;
23. UINT32 i=0;
24. do
25. {
26. i++; r += log(i);
27. }while (r<f);
28. return i;
29. }
31. void calc_e_2(integer *r,int len)
32. {
33. integer p,q;
34. UINT32 n=maxNForCalcE(len);
36. p=1; q=1;
37. do
38. {
39. q*=n;
40. p+=q;
41. n--;
42. }while (n>0);
44. p.DecLShift( q.GetDigits()-1); //增加有效数字
45. *r = p / q;
46. }
48. void calc_e2_main()
49. {
50. int d,len;
51. integer e;
52. char *pBuff=NULL;
53. const char *pTmp;
54. char fileName[32];
55. FILE *fp=NULL;
57. printf("calculate e (simple version2)\n");
58. printf("How many digital is needed?");
59. scanf("%d",&d);
60. len=(d+17)/9*9; //取9的倍数，可能稍微快一点
62. HugeCalc::EnableTimer( TRUE );
63. HugeCalc::ResetTimer();
65. calc_e_2(&e,len);
67. // 记录计算耗时
68. UINT32 u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_us );
69. printf("计算用时 %.6f 秒\n", (double)u32Timer_Calc/1000000.0);
71. //----------------------------------------
72. HugeCalc::EnableTimer( TRUE );
73. HugeCalc::ResetTimer();
75. pBuff=new char[e.GetDigits()+3];
76. if (pBuff==NULL)
77. goto thisExit;
79. pTmp= e.GetStr(FS_NORMAL);
80. strcpy(pBuff+1,pTmp);
82. pBuff[0]=pBuff[1];
83. pBuff[1]='.'; //插入小数点
84. pBuff[d+2]=0; //去掉多余的数字
86. sprintf(fileName,"e2_%d.txt",d);
88. fp=fopen(fileName,"wb");
89. fwrite(pBuff,d+2,1,fp);
90. fclose(fp); fp=NULL;
92. u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_us );
93. printf("输出到文件用时 %.6f 秒\n", (double)u32Timer_Calc/1000000.0);
95. thisExit:
96. if (fp!=NULL)
97. {
98. fclose(fp); fp=NULL;
99. }
100. if (pBuff!=NULL)
101. {
102. delete[] pBuff; pBuff=NULL;
103. }
104. }

复制代码

liangbch

本帖将给出一个快速的算法，这个算法能够充分利用大数乘法，计算速度比前两个算法快的多（在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化为小数即可

无心人

看不懂哦

liangbch

本帖给出使用HugeCalc，使用7#楼算法计算e的一个实现,使用该程序，在PIV2.6上，计算10万位e需0.412秒，计算100万位e需要5.06秒。

1. #include <stdlib.h>
2. #include <stdio.h>
3. #include <math.h>
5. #include "../../../HugeCalc_API/CppAPI/Include/HugeCalc.h" // 公共接口
6. #include "../../../HugeCalc_API/CppAPI/Include/HugeInt.h" // 10进制系统
8. #pragma message( "automatic link to ../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
9. #pragma comment( lib, "../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
12. #define integer CHugeInt
13. #define USE_HUGE_CALC_DIV
15. #if 0
16. # define MODIFIED_BY_GxQ
17. #endif
19. #define MIN_SPLIT 32
20. static DWORD g_divTime=0;
22. // 计算e的一个快速的版本，采用2分法，尽量使用以及大数乘法，
23. // 收尾阶段采用大数除法
26. //欲达到n位有效数字，需使得 n!> 10^n
27. //返回计算d位有效数字所需的最大的n
29. static UINT32 maxNForCalcE(int d)
30. {
31. double f=d*log(10);
32. double r=0;
33. UINT32 i=0;
34. //loop untill log(n!)>log(10^d)
35. do
36. {
37. i++; r += log(i);
38. }while (r<f);
39. return i;
40. }
42. // calc 1/n1 + 1/(n1*(n1+1)) + 1/(n1*(n1+1)+(n1+2))+ 1/(n1*(n1+1)+(n1+2)+...m)
43. // 计算结果为 p/q
44. // for example:
46. // 1/6 + 1/6*7 + 1/6*7*8 + 1/6*7*8*9
47. // = 7*8*9/6*7*8*9 + 8*9/6*7*8*9 + 9/6*7*8*9 + 1/ 6*7*8*9
48. // calculate
49. // q=6*7*8*9, p=1 + 9+ 8*9 + 7*8*9
50. // steps:
51. // step1: p=1, q=1;
52. // step2: q*=9, p+=q , q=9, q=1+9
53. // step3: q*=8, p+=q , q=9*8, q=1+8*9
54. // step4: q*=7, p+=q , q=9*8*7, q=1+ 8*9 + 7*8*9
55. // step5: q*=6
56. static void calc_Frac_Sum(UINT32 n1,UINT32 n2, integer &p,integer &q)
57. {
58. p=(UINT32)1;
59. q=(UINT32)1;
60. for ( UINT32 i=n2;i>n1;i--)
61. {
62. q*=i;
63. p+=q;
64. }
65. q*=i;
66. }
68. #ifndef MODIFIED_BY_GxQ
69. void splitCalc(UINT32 n1,UINT32 n2, integer &p,integer &q)
70. {
71. if (n2-n1<=MIN_SPLIT)
72. {
73. calc_Frac_Sum(n1,n2,p,q);
74. }
75. else
76. {
77. integer p1,q1,p2,q2;
78. splitCalc( n1, (n1+n2)/2, p1,q1);
79. splitCalc( (n1+n2+2)/2 ,n2, p2,q2);
80. p=p1*q2 + p2;
81. q=q1*q2;
82. }
83. }
84. #else
85. integer tmp;
86. void splitCalc(UINT32 n1,UINT32 n2, integer &p,integer &q)
87. {
88. if (n2-n1<=MIN_SPLIT)
89. {
90. calc_Frac_Sum(n1,n2,p,q);
91. }
92. else
93. {
94. integer q1,p2;
95. splitCalc( n1, (n1+n2)/2, p, q1);
96. splitCalc( (n1+n2+2)/2 ,n2, p2,q );
98. // ( p *= q ) += p2;
99. tmp.Swap( p );
100. p.Mul( tmp, q ) += p2;
102. // q *= q1;
103. tmp.Swap( q );
104. q.Mul( tmp, q1 );
105. }
106. }
107. #endif
109. void calc_e_3(integer *r, int len )
110. {
111. UINT32 n;
112. int shift;
113. integer p,q;
115. n= maxNForCalcE(len)+2;
116. splitCalc(1,n,p,q); // p/q= 1/1!+1/2!+1/3!+1/4!+...
118. p+=q; // p/q +=1 // e= p/q= 1 + 1/1!+1/2!+1/3!+1/4!+...
119. shift=(len+8)/9*9; //get upbounds get times of 9
120. p.DecLShift(shift); //增加有效数字
122. g_divTime=GetTickCount();
123. #ifdef USE_HUGE_CALC_DIV
124. #ifndef MODIFIED_BY_GxQ
125. p=p /q;
126. #else
127. p /= q;
128. #endif
129. #else
130. {
131. void inverse( const integer *s, integer *pResult,int len);
132. integer t;
133. inverse(&q,&t,len);
134. p *=t;
135. }
136. #endif
137. g_divTime=GetTickCount()-g_divTime;
138. *r=p;
139. }
141. void calc_e3_main()
142. {
143. int len;
144. integer e;
145. char *pBuff=NULL;
146. const char *pTmp;
147. char fileName[32];
148. FILE *fp=NULL;
150. printf("calculate e (fast version1)\n");
151. #ifdef MODIFIED_BY_GxQ
152. printf("defined 'MODIFIED_BY_GxQ'\n");
153. #endif
154. while (true)
155. {
156. printf("How many digital is needed(0 will exit)?");
157. scanf("%d",&len);
158. if (len==0)
159. break;
161. HugeCalc::EnableTimer( TRUE );
162. HugeCalc::ResetTimer();
164. calc_e_3(&e,len);
166. // 记录计算耗时
167. UINT32 u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_us );
168. printf("计算用时 %.6f 秒\n", (double)u32Timer_Calc/1000000.0);
169. if (g_divTime>=2)
170. {
171. printf("其中除法用时 %f秒\n", (double)g_divTime/1000.00);
172. }
173. //----------------------------------------
174. HugeCalc::EnableTimer( TRUE );
175. HugeCalc::ResetTimer();
177. pBuff=new char[e.GetDigits()+3];
178. if (pBuff==NULL)
179. goto thisExit;
181. pTmp= e.GetStr(FS_NORMAL);
182. strcpy(pBuff+1,pTmp);
184. pBuff[0]=pBuff[1];
185. pBuff[1]='.'; //插入小数点
186. pBuff[len+2]=0; //去掉多余的数字
188. sprintf(fileName,"e3_%d.txt",len);
189. FILE *fp=fopen(fileName,"wb");
190. fwrite(pBuff,len+2,1,fp);
191. fclose(fp); fp=NULL;
192. delete[] pBuff; pBuff=NULL;
194. u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_us );
195. printf("输出到文件用时 %.6f 秒\n", (double)u32Timer_Calc/1000000.0);
196. }
198. thisExit:
199. if (fp!=NULL)
200. {
201. fclose(fp); fp=NULL;
202. }
203. if (pBuff!=NULL)
204. {
205. delete[] pBuff; pBuff=NULL;
206. }
207. }

复制代码

[ 本帖最后由 liangbch 于 2008-4-17 14:34 编辑 ]

无心人

我想 我想 是不是有更精彩更复杂的算法 比如类似阶乘的？