计算百万位e

shshsh_0510

呵呵，liangbch已经给出了

无心人

由liangbch给出的网页得到的算法 只需要长整数乘以32位小整数和长加法算法计算到足够精度 然后进行一次长除法，再移位加一了 算法复杂度显然小于他自己写的程序 具体算法俺已从他给出网页整理到19#

liangbch

回复 17楼 mathe. 求 x，使得 $ log10(x!)=n $有更有效的算法，只是这部分代码占用的时间很小，我在之前的代码中没有给出更有效更复杂的代码。 下面贴出代码： 函数 faclogEquation(double f) 返回 方程 $log10(x!)=f $ 的解

1. #define PI 3.1415926535897932384626433832795
2. #define LOG10 2.3025850929940456840179914546844
3. #define LOGS2PI 0.91893853320467274178032973640562 //log(sqrt(2*pi))
4. typedef unsigned long DWORD;
6. typedef struct _logfacpair
7. {
8. unsigned long n;
9. double logfac;
10. }LOGFAC_PAIR;
12. #include <stdio.h>
13. #include <stdlib.h>
14. #include <math.h>
15. #include <assert.h>
16. /*
17. 计算e,使用HugeCalc，一个简单的版本
18. */
19. //#define DEBUG_LOG_FAC
21. const double logfac[]=
22. {
23. 0.000000, // 0
24. 0.000000, // 1
25. 0.693147, // 2
26. 1.791759, // 3
27. 3.178054, // 4
28. 4.787492, // 5
29. 6.579251, // 6
30. 8.525161, // 7
31. 10.604603, // 8
32. 12.801827, // 9
33. 15.104413, // 10
34. 17.502308, // 11
35. 19.987214, // 12
36. 22.552164, // 13
37. 25.191221, // 14
38. 27.899271, // 15
39. 30.671860, // 16
40. 33.505073, // 17
41. 36.395445, // 18
42. 39.339884, // 19
43. 42.335616, // 20
44. 45.380139, // 21
45. 48.471181, // 22
46. 51.606676, // 23
47. 54.784729, // 24
48. 58.003605, // 25
49. 61.261702, // 26
50. 64.557539, // 27
51. 67.889743, // 28
52. 71.257039, // 29
53. 74.658236, // 30
54. 78.092224, // 31
55. 81.557959, // 32
56. 85.054467, // 33
57. 88.580828, // 34
58. 92.136176, // 35
59. 95.719695, // 36
60. 99.330612, // 37
61. 102.968199, // 38
62. 106.631760, // 39
63. 110.320640, // 40
64. 114.034212, // 41
65. 117.771881, // 42
66. 121.533082, // 43
67. 125.317271, // 44
68. 129.123934, // 45
69. 132.952575, // 46
70. 136.802723, // 47
71. 140.673924, // 48
72. 144.565744, // 49
73. 148.477767, // 50
74. 152.409593, // 51
75. 156.360836, // 52
76. 160.331128, // 53
77. 164.320112, // 54
78. 168.327445, // 55
79. 172.352797, // 56
80. 176.395848, // 57
81. 180.456291, // 58
82. 184.533829, // 59
83. 188.628173, // 60
84. 192.739047, // 61
85. 196.866182, // 62
86. 201.009316, // 63
87. 205.168199, // 64
88. 209.342587, // 65
89. 213.532241, // 66
90. 217.736934, // 67
91. 221.956442, // 68
92. 226.190548, // 69
93. 230.439044, // 70
94. 234.701723, // 71
95. 238.978390, // 72
96. 243.268849, // 73
97. 247.572914, // 74
98. 251.890402, // 75
99. 256.221136, // 76
100. 260.564941, // 77
101. 264.921650, // 78
102. 269.291098, // 79
103. 273.673124, // 80
104. 278.067573, // 81
105. 282.474293, // 82
106. 286.893133, // 83
107. 291.323950, // 84
108. 295.766601, // 85
109. 300.220949, // 86
110. 304.686857, // 87
111. 309.164194, // 88
112. 313.652830, // 89
113. 318.152640, // 90
114. 322.663499, // 91
115. 327.185288, // 92
116. 331.717887, // 93
117. 336.261182, // 94
118. 340.815059, // 95
119. 345.379407, // 96
120. 349.954118, // 97
121. 354.539086, // 98
122. 359.134205, // 99
123. 363.739376, // 100
124. 368.354496, // 101
125. 372.979469, // 102
126. 377.614198, // 103
127. 382.258589, // 104
128. 386.912549, // 105
129. 391.575988, // 106
130. 396.248817, // 107
131. 400.930948, // 108
132. 405.622296, // 109
133. 410.322777, // 110
134. 415.032307, // 111
135. 419.750806, // 112
136. 424.478193, // 113
137. 429.214392, // 114
138. 433.959324, // 115
139. 438.712914, // 116
140. 443.475088, // 117
141. 448.245773, // 118
142. 453.024896, // 119
143. 457.812388, // 120
144. 462.608179, // 121
145. 467.412200, // 122
146. 472.224384, // 123
147. 477.044665, // 124
148. 481.872979, // 125
149. 486.709261, // 126
150. 491.553448, // 127
151. 496.405478 // 128
152. };
154. const LOGFAC_PAIR logfac_tab[]=
155. {
156. { 1, 0.000000 },
157. { 2, 0.693147 },
158. { 4, 3.178054 },
159. { 8, 10.604603 },
160. { 16, 30.671860 },
161. { 32, 81.557959 },
162. { 64, 205.168199 },
163. { 128, 496.405478 },
164. { 256, 1167.256953 },
165. { 512, 2686.060309 },
166. { 1024, 6078.211803 },
167. { 2048, 13571.950932 },
168. { 4096, 29978.648040 },
169. { 8192, 65630.826536 },
170. { 16384, 142613.098657 },
171. { 32768, 307933.819731 },
172. { 65536, 661287.962119 },
173. { 131072, 1413421.993946 },
174. { 262144, 3008541.898276 },
175. { 524288, 6380485.734864 },
176. { 1048576, 13487781.810467 },
177. { 2097152, 28429191.113103 },
178. { 4194304, 59765644.367806 },
179. { 8388608, 125345820.522651 },
180. { 16777216, 262320712.469790 },
181. { 33554432, 547899575.985538 },
182. { 67108864, 1142315462.606553 },
183. { 134217728, 2377663555.374189 },
184. { 268435456, 4941392380.307250 },
185. { 536870912, 10254915309.315521 },
186. { 1073741824, 21254091725.962936 },
187. { 2147483648, 43996705676.866089 }
188. };
190. //使用2分查找算法在logfac中 查找 x>=f 对应的下标
191. DWORD searchInLogFacTab1(double f)
192. {
193. DWORD low=0;
194. DWORD high=128;
195. DWORD mid;
197. while (true)
198. {
199. mid=(low+high)/2;
200. if ( mid==low)
201. {
202. while( logfac[mid]<f)
203. mid++;
204. break;
205. }
206. if ( f>=logfac[mid] )
207. {
208. low=mid;
209. }
210. else
211. {
212. high=mid;
213. }
214. }
215. return mid;
216. }
218. double getLogfac(DWORD n)
219. {
220. if (n<=128)
221. return logfac[n];
222. else
223. {
224. return LOGS2PI + ((double)n+0.5) * log(n)-n;
225. }
226. }
228. /*
229. 求一个数x,使得满足 不等式 x! > 10^n 的最小的值
231. 解方程 x!=10^n (1)
232. 定义
233. c1= 0.5log(2pi)= 0.91893853320467274178032973640562
234. c2= n* log(10)= n* 2.3025850929940456840179914546844
237. 根据斯特林公式,有
238. x!= sqrt(2*pi*x) * (x/e)^x (pi:为圆周率，e：自然对数的底)
240. 方程左边取对数
241. log(x!) = log ( sqrt(2*pi*x)* (x/e)^x)
242. = 0.5(log(x)+log(2pi))+ x(log(x)-1)
243. = (x+0.5)logx-x + 0.5log(2pi)
244. = (x+0.5)logx-x + c1
245. 方程右边取对数
246. log(10^n)= n * log(10)2.302585
247. = c2
249. 则原方程(1)化为 (x+0.5)logx-x + c1= c2 (2)
251. 定义：f(x)=log(x!)
252. 则有： f'(x)= 0.5+logx
254. 若 方程(2)的一个近似为x0,可导出牛顿迭代式：
255. x1= x0 + ( c2 -f(x0))/f'(x)
256. = x0 + ( c2- (x0+0.5)log(x0) +x0 - c1)/(0.5/x0+log(x0))
257. */
258. // 解方程log(x!)=f,已知 low <= x <= high
259. // 返回x的值
260. DWORD faclogEquation(double f,int low,int high)
261. {
262. int i;
264. double x0,x1;
265. DWORD n;
267. if ( f <= logfac[128])
268. return searchInLogFacTab1(f);
269. x0= (low+high)/2;
270. #ifdef DEBUG_LOG_FAC
271. printf("---------\n");
272. printf("x0=%f\n",x0);
273. #endif
274. while(true)
275. {
276. x1= x0 + ( f- getLogfac((DWORD)x0))/(0.5/x0+log(x0));
277. #ifdef DEBUG_LOG_FAC
278. printf("x1=%f\n",x1);
279. #endif
280. if ( fabs(x1-x0)<0.75 )
281. break;
282. x0=x1;
283. }
284. #ifdef DEBUG_LOG_FAC
285. printf("---------\n");
286. #endif
287. if (x1<x0)
288. x0=x1;
289. n=(DWORD)x0;
291. while ( getLogfac(n) < f)
292. n++;
293. return n;
294. }
296. // 解方程log(x!)=f
297. // 返回x的值
298. DWORD faclogEquation(double f)
299. {
300. int i;
301. int low,high;
303. if ( f <= logfac[128])
304. return searchInLogFacTab1(f);
306. i=7;
307. while ( logfac_tab[i].logfac<f)
308. i++;
309. low= logfac_tab[i-1].n;
310. high=logfac_tab[i].n;
311. return faclogEquation(f,low,high);
312. }

复制代码

liangbch

回复15楼， 我刚才用刘楚雄的计算器计算了一下，计算100万位e值需要3.50秒，可能是他用的 ooura FFT 快于楼主的FNT.

liangbch

回复22楼，如果频繁的使用 长整数乘以32位小整数，绝不会比9#算法更快。因为 大整数 乘以 单精度 整数没有快速算法，而大数乘以大数有快速算法。

无心人

你竟然有这么奇怪的论调啊？

无心人

另外 任何FFT均慢于最优化的NTT FFT快， 说明HugeCalc的NTT算法不行 GxQ的NTT就是NTT+CRT 没利用好其他好的NTT 因为NTT算法有很多种 在不同情况下，各有优势

liangbch

http://numbers.computation.free.fr/Constants 网站上给出一个计算任意位 e 的程序，其大整数乘法用FFT实现，FFT的代码非常简单（当然效率差点儿）。整个程序的代码量并不大，下面我贴出这个程序完整的代码。 文件 fft.h 的内容

1. #ifndef FFT_h
2. #define FFT_h
4. typedef double Real;
6. extern double FFTSquareWorstError;
7. extern long AllocatedMemory;
9. void InitializeFFT(long MaxFFTLength);
11. void MulWithFFT(Real * ACoef, long ASize,
12. Real * BCoef, long BSize,
13. Real * CCoef);
15. #endif

复制代码

文件 BigInt.h的内容

1. #ifndef BigInt_h
2. #define BigInt_h
4. #include "FFT.h"
6. typedef struct BigIntTag {
7. Real * Coef;
8. long Size, SizeMax;
9. } BigInt;
11. /*
12. * Decrease the base if the worst error in FFT become too large
13. */
14. #define BASE 10000.
15. #define invBASE 0.0001
16. #define NBDEC_BASE 4
18. void InitializeBigInt(BigInt * A, long MaxSize);
19. void FreeBigInt(BigInt * A);
21. void PrintBigInt(BigInt * A);
23. void UpdateBigInt(BigInt * A);
24. void AddBigInt(BigInt * A, BigInt * B, BigInt * C);
25. void MulBigInt(BigInt * A, BigInt * B, BigInt * C);
26. void Inverse(BigInt * A, BigInt * B, BigInt * tmpBigInt);
29. #endif

复制代码

FFT.c 的内容

1. /*
2. * Xavier Gourdon : Sept. 99 (xavier.gourdon.free.fr)
3. *
4. * FFT.c : Very basic FFT file.
5. * The FFT encoded in this file is short and very basic.
6. * It is just here as a teaching file to have a first
7. * understanding of the FFT technique.
8. *
9. * A lot of optimizations could be made to save a factor 2 or 4 for time
10. * and space.
11. * - Use a 4-Step (or more) FFT to avoid data cache misses.
12. * - Use an hermitian FFT to take into account the hermitian property of
13. * the FFT of a real array.
14. * - Use a quad FFT (recursion N/4->N instead of N/2->N) to save 10 or
15. * 15% of the time.
16. *
17. * Informations can be found on
18. * http://xavier.gourdon.free.fr/Constants/constants.html
19. */
20. #include "FFT.h"
22. #include <stdlib.h>
23. #include <stdio.h>
24. #include <math.h>
26. #define PI 3.1415926535897932384626
28. typedef struct ComplexTag {
29. Real R,I;
30. } Complex;
32. long FFTLengthMax;
33. Complex * OmegaFFT;
34. Complex * ArrayFFT0, * ArrayFFT1;
35. Complex * ComplexCoef;
36. double FFTSquareWorstError;
37. long AllocatedMemory;
39. void InitializeFFT(long MaxLength)
40. {
41. long i;
42. Real Step;
44. FFTLengthMax = MaxLength;
45. OmegaFFT = (Complex *) malloc(MaxLength/2*sizeof(Complex));
46. ArrayFFT0 = (Complex *) malloc(MaxLength*sizeof(Complex));
47. ArrayFFT1 = (Complex *) malloc(MaxLength*sizeof(Complex));
48. ComplexCoef = (Complex *) malloc(MaxLength*sizeof(Complex));
49. Step = 2.*PI/(double) MaxLength;
50. for (i=0; 2*i<MaxLength; i++) {
51. OmegaFFT[i].R = cos(Step*(double)i);
52. OmegaFFT[i].I = sin(Step*(double)i);
53. }
54. FFTSquareWorstError=0.;
55. AllocatedMemory = 7*MaxLength*sizeof(Complex)/2;
56. }
58. void RecursiveFFT(Complex * Coef, Complex * FFT, long Length, long Step,
59. long Sign)
60. {
61. long i, OmegaStep;
62. Complex * FFT0, * FFT1, * Omega;
63. Real tmpR, tmpI;
65. if (Length==2) {
66. FFT[0].R = Coef[0].R + Coef[Step].R;
67. FFT[0].I = Coef[0].I + Coef[Step].I;
68. FFT[1].R = Coef[0].R - Coef[Step].R;
69. FFT[1].I = Coef[0].I - Coef[Step].I;
70. return;
71. }
73. FFT0 = FFT;
74. RecursiveFFT(Coef ,FFT0,Length/2,Step*2,Sign);
75. FFT1 = FFT+Length/2;
76. RecursiveFFT(Coef+Step,FFT1,Length/2,Step*2,Sign);
78. Omega = OmegaFFT;
79. OmegaStep = FFTLengthMax/Length;
80. for (i=0; 2*i<Length; i++, Omega += OmegaStep) {
81. /* Recursion formula for FFT :
82. FFT[i] <- FFT0[i] + Omega*FFT1[i]
83. FFT[i+Length/2] <- FFT0[i] - Omega*FFT1[i],
84. Omega = exp(2*I*PI*i/Length) */
85. tmpR = Omega[0].R*FFT1[i].R-Sign*Omega[0].I*FFT1[i].I;
86. tmpI = Omega[0].R*FFT1[i].I+Sign*Omega[0].I*FFT1[i].R;
87. FFT1[i].R = FFT0[i].R - tmpR;
88. FFT1[i].I = FFT0[i].I - tmpI;
89. FFT0[i].R = FFT0[i].R + tmpR;
90. FFT0[i].I = FFT0[i].I + tmpI;
91. }
92. }
94. /* Compute the complex Fourier Transform of Coef into FFT */
95. void FFT(Real * Coef, long Length, Complex * FFT, long NFFT)
96. {
97. long i;
98. /* Transform array of real coefficient into array of complex */
99. for (i=0; i<Length; i++) {
100. ComplexCoef[i].R = Coef[i];
101. ComplexCoef[i].I = 0.;
102. }
103. for (; i<NFFT; i++)
104. ComplexCoef[i].R = ComplexCoef[i].I = 0.;
106. RecursiveFFT(ComplexCoef,FFT,NFFT,1,1);
107. }
109. /* Compute the inverse Fourier Transform of FFT into Coef */
110. void InverseFFT(Complex * FFT, long NFFT, Real * Coef, long Length)
111. {
112. long i;
113. Real invNFFT = 1./(Real) NFFT, tmp;
115. RecursiveFFT(FFT, ComplexCoef, NFFT, 1, -1);
116. for (i=0; i<Length; i++) {
117. /* Closest integer to ComplexCoef[i].R/NFFT */
118. tmp = invNFFT*ComplexCoef[i].R;
119. Coef[i] = floor(0.5+tmp);
120. if ((tmp-Coef[i])*(tmp-Coef[i])>FFTSquareWorstError)
121. FFTSquareWorstError = (tmp-Coef[i])*(tmp-Coef[i]);
122. }
123. }
125. void Convolution(Complex * A, Complex * B, long NFFT, Complex * C)
126. {
127. long i;
128. Real tmpR, tmpI;
130. for (i=0; i<NFFT; i++) {
131. tmpR = A[i].R*B[i].R-A[i].I*B[i].I;
132. tmpI = A[i].R*B[i].I+A[i].I*B[i].R;
133. C[i].R = tmpR;
134. C[i].I = tmpI;
135. }
136. }
138. void MulWithFFT(Real * ACoef, long ASize,
139. Real * BCoef, long BSize,
140. Real * CCoef)
141. {
142. long NFFT = 2;
144. while (NFFT<ASize+BSize)
145. NFFT *= 2;
147. if (NFFT>FFTLengthMax) {
148. printf("Error, FFT Size is too big in MulWithFFT\n");
149. return;
150. }
151. FFT(ACoef, ASize, ArrayFFT0, NFFT);
152. FFT(BCoef, BSize, ArrayFFT1, NFFT);
153. Convolution(ArrayFFT0,ArrayFFT1,NFFT,ArrayFFT0);
154. InverseFFT(ArrayFFT0,NFFT,CCoef, ASize+BSize-1);
155. }

复制代码

BigInt.c的内容

1. /*
2. * Xavier Gourdon : Sept. 99 (xavier.gourdon.free.fr)
3. *
4. * BigInt.c : Basic file for manipulation of Large Integers.
5. * It rely on FFT.c to multiply number with the FFT
6. * technique.
7. *
8. * Several optimizations can be made :
9. * - Implementation of a classic multiplication for small BigInts.
10. * - Avoid using UpdateBigInt for the sum of BigInts.
11. * - Use a larger base for the internal data storage of BigInts to
12. * save space.
13. * ...
14. *
15. * Informations can be found on
16. * http://xavier.gourdon.free.fr/Constants/constants.html
17. */
19. #include "BigInt.h"
20. #include "FFT.h"
22. #include <math.h>
23. #include <stdlib.h>
24. #include <stdio.h>
26. void InitializeBigInt(BigInt * A, long MaxSize)
27. {
28. A->Coef = (Real *) malloc(MaxSize*sizeof(Real));
29. AllocatedMemory += MaxSize*sizeof(Real);
31. if (A->Coef==0)
32. printf("Not enough memory\n");
33. A->Size = 0;
34. A->SizeMax = MaxSize;
35. }
37. void FreeBigInt(BigInt * A)
38. {
39. free(A->Coef);
40. }
42. void PrintBigInt(BigInt * A)
43. {
44. long i,j,Digit=0,Dec;
45. double Pow, Coef;
47. printf("%.0lf\n",A->Coef[A->Size-1]);
49. for (i=A->Size-2; i>=0; i--) {
50. Pow = BASE*0.1;
51. Coef = A->Coef[i];
52. for (j=0; j<NBDEC_BASE; j++) {
53. Dec = (long) (Coef/Pow);
54. Coef -= Dec*Pow;
55. Pow *= 0.1;
56. printf("%ld",Dec);
57. Digit++;
58. if (Digit%10==0) printf(" ");
59. if (Digit%50==0) printf(": %ld\n",Digit);
60. }
61. }
62. }
64. /*
65. * Update A to the normal form (0<= A.Coef[i] < BASE)
66. */
67. void UpdateBigInt(BigInt * A)
68. {
69. long i;
70. Real carry=0., x;
72. for (i=0; i<A->Size; i++) {
73. x = A->Coef[i] + carry;
74. carry = floor(x*invBASE);
75. A->Coef[i] = x-carry*BASE;
76. }
77. if (carry!=0) {
78. while (carry!=0.) {
79. x = carry;
80. carry = floor(x*invBASE);
81. A->Coef[i] = x-carry*BASE;
82. i++;
83. A->Size=i;
84. }
85. if (A->Size > A->SizeMax) {
86. printf("Error in UpdateBigInt, Size>SizeMax\n");
87. }
88. }
89. else {
90. while (i>0 && A->Coef[i-1]==0.) i--;
91. A->Size=i;
92. }
93. }
95. /*
96. * Compute C = A + B
97. */
98. void AddBigInt(BigInt * A, BigInt * B, BigInt * C)
99. {
100. long i;
102. if (A->Size<B->Size) {
103. AddBigInt(B,A,C);
104. return;
105. }
106. /* We now have B->Size<A->Size */
107. for (i=0; i<B->Size; i++)
108. C->Coef[i] = A->Coef[i] + B->Coef[i];
109. for (; i<A->Size; i++)
110. C->Coef[i] = A->Coef[i];
112. C->Size = A->Size;
113. UpdateBigInt(C);
114. }
116. /*
117. * Compute C = A*B
118. */
119. void MulBigInt(BigInt * A, BigInt * B, BigInt * C)
120. {
121. MulWithFFT(A->Coef,A->Size,B->Coef,B->Size,C->Coef);
122. C->Size = A->Size+B->Size-1;
123. UpdateBigInt(C);
124. }
126. /*
127. * Compute the inverse of A in B
128. */
129. void Inverse(BigInt * A, BigInt * B, BigInt * tmpBigInt)
130. {
131. double x;
132. long i, N,NN,Delta;
133. int Twice=1, Sign;
134. BigInt AA;
136. /* Initialization */
137. x = A->Coef[A->Size-1]+invBASE*(A->Coef[A->Size-2]
138. + invBASE*A->Coef[A->Size-3]);
139. x = BASE/x;
140. B->Coef[1] = floor(x);
141. B->Coef[0] = floor((x-B->Coef[1])*BASE);
142. B->Size = 2;
144. /* Iteration used : B <- B+(1-AB)*B */
145. N=2;
146. while (N<A->Size) {
147. /* Use only the first 2*N limbs of A */
148. NN = 2*N;
149. if (NN>A->Size)
150. NN = A->Size;
151. AA.Coef = A->Coef+A->Size-NN;
152. AA.Size = NN;
153. MulBigInt(&AA,B,tmpBigInt);
154. Delta = NN+B->Size-1;
155. /* Compute BASE^Delta-tmpBigInt in tmpBigInt */
156. if (tmpBigInt->Size==Delta) {
157. Sign=1;
158. for (i=0; i<Delta; i++)
159. tmpBigInt->Coef[i] = BASE-1 - tmpBigInt->Coef[i];
160. UpdateBigInt(tmpBigInt);
161. }
162. else {
163. Sign = -1;
164. tmpBigInt->Coef[Delta] = 0.;
165. }
167. MulBigInt(tmpBigInt,B,tmpBigInt);
168. for (i=0; i<tmpBigInt->Size-2*N+1; i++)
169. tmpBigInt->Coef[i] = tmpBigInt->Coef[i+2*N-1];
170. tmpBigInt->Size -= 2*N-1;
171. for (i=B->Size-1; i>=0; i--)
172. B->Coef[i+NN-N] = B->Coef[i];
173. for (i=NN-N-1; i>=0; i--)
174. B->Coef[i]=0.;
175. B->Size += NN-N;
176. if (Sign==-1) {
177. for (i=0; i<tmpBigInt->Size; i++)
178. tmpBigInt->Coef[i] = -tmpBigInt->Coef[i];
179. }
180. AddBigInt(B,tmpBigInt,B);
181. /* Once during the process, redo one iteration with same precision */
182. if (8*N>A->Size && Twice) {
183. Twice=0;
184. B->Size = N;
185. for (i=0; i<N; i++)
186. B->Coef[i] = B->Coef[i+N];
187. }
188. else
189. N *= 2;
190. }
191. }

复制代码

主程序(e.c)的内容

1. /*
2. * Xavier Gourdon : Sept. 99 (xavier.gourdon.free.fr)
3. *
4. * e.c : Basic file for computation of e=2.7181... with the binary
5. * splitting method.
6. * It relies on the files FFT.c and BigInt.c to handle
7. * Large Integers.
8. *
9. * Several optimizations can be made :
10. * - Use a specialized routine to handle binary splitting when the
11. * indexes are close.
12. * ...
13. *
14. * Informations can be found on
15. * http://xavier.gourdon.free.fr/Constants/constants.html
16. */
18. #include "FFT.h"
19. #include "BigInt.h"
21. #include <stdio.h>
22. #include <math.h>
23. #include <time.h>
25. BigInt tmpBigInt; /* Temporary BigInt used in the binary splitting */
27. /*
28. * Compute the numerator P and denominator Q of
29. * P/Q = 1/(N0+1) + 1/(N0+1)/(N0+2) + ... + 1/(N0+1)/.../N1
30. */
31. void BinarySplittingE(long N0, long N1, BigInt * P, BigInt * Q)
32. {
33. BigInt PP, QQ;
34. long NMid;
36. if (N1-N0==1) {
37. P->Size = Q->Size = 1;
38. P->Coef[0] = 1.;
39. Q->Coef[0] = (double) N1;
40. UpdateBigInt(P);
41. UpdateBigInt(Q);
42. return;
43. }
44. NMid = (N0+N1)/2;
45. BinarySplittingE(N0,NMid,P,Q);
46. /* To save memory, take the non used coefficients of P and Q
47. for coefficient of the BigInt used in the second splitting part */
48. PP.Coef = P->Coef + P->Size;
49. PP.SizeMax = P->SizeMax-P->Size;
50. QQ.Coef = Q->Coef + Q->Size;
51. QQ.SizeMax = Q->SizeMax-Q->Size;
53. BinarySplittingE(NMid,N1,&PP,&QQ);
54. MulBigInt(P,&QQ,&tmpBigInt);
55. AddBigInt(&tmpBigInt,&PP,P);
57. MulBigInt(Q,&QQ,Q);
58. }
60. /*
61. * Returns as a BigInt the constant E with NbDec decimal digits
62. */
63. BigInt ECompute(long NbDec)
64. {
65. BigInt P, Q, tmp;
66. long i, MaxSize, MaxFFTSize, SeriesSize;
67. double logFactorial, logMax;
69. /* MaxSize should be more than NbDec/NBDEC_BASE (see BinarySplitting) */
70. MaxSize = NbDec/NBDEC_BASE + 10 + (long) (2.*log((double)NbDec)/log(2.));
71. InitializeBigInt(&P,MaxSize);
72. InitializeBigInt(&Q,MaxSize);
74. MaxFFTSize = 2;
75. while (MaxFFTSize < MaxSize)
76. MaxFFTSize *= 2;
77. MaxFFTSize *= 2;
78. InitializeFFT(MaxFFTSize);
79. /* Temporary BigInts are needed */
80. InitializeBigInt(&tmpBigInt,MaxFFTSize);
81. InitializeBigInt(&tmp,MaxFFTSize);
83. printf("Total Allocated memory = %ld K\n",AllocatedMemory/1024);
85. /* Compute the number of term SeriesSize of the series required.
86. SeriesSize is such that log(SeriesSize !) > NbDec*log(10.) */
87. SeriesSize=1;
88. logFactorial=0.;
89. logMax = (double)NbDec * log(10.);
90. while (logFactorial<logMax) {
91. logFactorial += log((double) SeriesSize);
92. SeriesSize++;
93. }
95. printf("Starting series computation\n");
96. /* Compute the numerator P and the denominator Q of
97. sum_{i=0}^{SeriesSize} 1/i! */
98. BinarySplittingE(0,SeriesSize,&P,&Q);
99. AddBigInt(&P,&Q,&P);
101. printf("Starting final division\n");
102. /* Compute the inverse of Q in tmpBigInt */
103. Inverse(&Q,&tmpBigInt,&tmp);
104. /* Comute P/Q in tmpBigInt */
105. MulBigInt(&P,&tmpBigInt,&tmpBigInt);
107. /* Put the number of required digits in P */
108. P.Size = 1+NbDec/NBDEC_BASE;
109. for (i=1; i<=P.Size; i++)
110. P.Coef[P.Size-i] = tmpBigInt.Coef[tmpBigInt.Size-i];
112. FreeBigInt(&Q);
113. FreeBigInt(&tmpBigInt);
115. return P;
116. }
119. void main()
120. {
121. double StartTime;
122. BigInt E;
123. long NbDec;
125. printf("*** E computation *** \n\n");
126. printf("Enter the number of decimal digits : ");
127. scanf("%ld",&NbDec);
128. printf("\n");
129. StartTime = (double) clock();
130. E = ECompute(NbDec);
131. printf("\n");
132. printf("Total time : %.2lf seconds\n",((double) clock() - StartTime)/CLOCKS_PER_SEC);
133. printf("Worst error in FFT (should be less than 0.25): %.10lf\n",sqrt(FFTSquareWorstError));
134. printf("\n");
135. printf("E = ");
136. //PrintBigInt(&E);
137. printf("\n");
138. }

复制代码

shshsh_0510

呵呵，借用李白观黄鹤楼的感慨：眼前有景道不得，liangbch程序在上头 只剩学习了

无心人

我有一种感觉 计算e的算法效率远不如计算pi