计算百万位e

shshsh_0510 · 发表于 2008-4-17 08:38:27

呵呵，liangbch已经给出了

无心人 · 发表于 2008-4-17 08:51:56

由liangbch给出的网页得到的算法
只需要长整数乘以32位小整数和长加法算法计算到足够精度

然后进行一次长除法，再移位加一了

算法复杂度显然小于他自己写的程序

具体算法俺已从他给出网页整理到19#

liangbch · 发表于 2008-4-17 10:39:46

回复 17楼 mathe.
求 x，使得 $ log10(x!)=n $有更有效的算法，只是这部分代码占用的时间很小，我在之前的代码中没有给出更有效更复杂的代码。

下面贴出代码：
函数 faclogEquation(double f) 返回方程 $log10(x!)=f $ 的解

#define PI 3.1415926535897932384626433832795
#define LOG10 2.3025850929940456840179914546844
#define LOGS2PI 0.91893853320467274178032973640562 //log(sqrt(2*pi))
typedef unsigned long DWORD;
typedef struct _logfacpair
{
unsigned long n;
double logfac;
}LOGFAC_PAIR;
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
/*
计算e,使用HugeCalc，一个简单的版本
*/
//#define DEBUG_LOG_FAC
const double logfac[]=
{
0.000000, // 0
0.000000, // 1
0.693147, // 2
1.791759, // 3
3.178054, // 4
4.787492, // 5
6.579251, // 6
8.525161, // 7
10.604603, // 8
12.801827, // 9
15.104413, // 10
17.502308, // 11
19.987214, // 12
22.552164, // 13
25.191221, // 14
27.899271, // 15
30.671860, // 16
33.505073, // 17
36.395445, // 18
39.339884, // 19
42.335616, // 20
45.380139, // 21
48.471181, // 22
51.606676, // 23
54.784729, // 24
58.003605, // 25
61.261702, // 26
64.557539, // 27
67.889743, // 28
71.257039, // 29
74.658236, // 30
78.092224, // 31
81.557959, // 32
85.054467, // 33
88.580828, // 34
92.136176, // 35
95.719695, // 36
99.330612, // 37
102.968199, // 38
106.631760, // 39
110.320640, // 40
114.034212, // 41
117.771881, // 42
121.533082, // 43
125.317271, // 44
129.123934, // 45
132.952575, // 46
136.802723, // 47
140.673924, // 48
144.565744, // 49
148.477767, // 50
152.409593, // 51
156.360836, // 52
160.331128, // 53
164.320112, // 54
168.327445, // 55
172.352797, // 56
176.395848, // 57
180.456291, // 58
184.533829, // 59
188.628173, // 60
192.739047, // 61
196.866182, // 62
201.009316, // 63
205.168199, // 64
209.342587, // 65
213.532241, // 66
217.736934, // 67
221.956442, // 68
226.190548, // 69
230.439044, // 70
234.701723, // 71
238.978390, // 72
243.268849, // 73
247.572914, // 74
251.890402, // 75
256.221136, // 76
260.564941, // 77
264.921650, // 78
269.291098, // 79
273.673124, // 80
278.067573, // 81
282.474293, // 82
286.893133, // 83
291.323950, // 84
295.766601, // 85
300.220949, // 86
304.686857, // 87
309.164194, // 88
313.652830, // 89
318.152640, // 90
322.663499, // 91
327.185288, // 92
331.717887, // 93
336.261182, // 94
340.815059, // 95
345.379407, // 96
349.954118, // 97
354.539086, // 98
359.134205, // 99
363.739376, // 100
368.354496, // 101
372.979469, // 102
377.614198, // 103
382.258589, // 104
386.912549, // 105
391.575988, // 106
396.248817, // 107
400.930948, // 108
405.622296, // 109
410.322777, // 110
415.032307, // 111
419.750806, // 112
424.478193, // 113
429.214392, // 114
433.959324, // 115
438.712914, // 116
443.475088, // 117
448.245773, // 118
453.024896, // 119
457.812388, // 120
462.608179, // 121
467.412200, // 122
472.224384, // 123
477.044665, // 124
481.872979, // 125
486.709261, // 126
491.553448, // 127
496.405478 // 128
};
const LOGFAC_PAIR logfac_tab[]=
{
{ 1, 0.000000 },
{ 2, 0.693147 },
{ 4, 3.178054 },
{ 8, 10.604603 },
{ 16, 30.671860 },
{ 32, 81.557959 },
{ 64, 205.168199 },
{ 128, 496.405478 },
{ 256, 1167.256953 },
{ 512, 2686.060309 },
{ 1024, 6078.211803 },
{ 2048, 13571.950932 },
{ 4096, 29978.648040 },
{ 8192, 65630.826536 },
{ 16384, 142613.098657 },
{ 32768, 307933.819731 },
{ 65536, 661287.962119 },
{ 131072, 1413421.993946 },
{ 262144, 3008541.898276 },
{ 524288, 6380485.734864 },
{ 1048576, 13487781.810467 },
{ 2097152, 28429191.113103 },
{ 4194304, 59765644.367806 },
{ 8388608, 125345820.522651 },
{ 16777216, 262320712.469790 },
{ 33554432, 547899575.985538 },
{ 67108864, 1142315462.606553 },
{ 134217728, 2377663555.374189 },
{ 268435456, 4941392380.307250 },
{ 536870912, 10254915309.315521 },
{ 1073741824, 21254091725.962936 },
{ 2147483648, 43996705676.866089 }
};
//使用2分查找算法在logfac中查找 x>=f 对应的下标
DWORD searchInLogFacTab1(double f)
{
DWORD low=0;
DWORD high=128;
DWORD mid;
while (true)
{
mid=(low+high)/2;
if ( mid==low)
{
while( logfac[mid]<f)
mid++;
break;
}
if ( f>=logfac[mid] )
{
low=mid;
}
else
{
high=mid;
}
}
return mid;
}
double getLogfac(DWORD n)
{
if (n<=128)
return logfac[n];
else
{
return LOGS2PI + ((double)n+0.5) * log(n)-n;
}
}
/*
求一个数x,使得满足不等式 x! > 10^n 的最小的值
解方程 x!=10^n (1)
定义
c1= 0.5log(2pi)= 0.91893853320467274178032973640562
c2= n* log(10)= n* 2.3025850929940456840179914546844
根据斯特林公式,有
x!= sqrt(2*pi*x) * (x/e)^x (pi:为圆周率，e：自然对数的底)
方程左边取对数
log(x!) = log ( sqrt(2*pi*x)* (x/e)^x)
= 0.5(log(x)+log(2pi))+ x(log(x)-1)
= (x+0.5)logx-x + 0.5log(2pi)
= (x+0.5)logx-x + c1
方程右边取对数
log(10^n)= n * log(10)2.302585
= c2
则原方程(1)化为 (x+0.5)logx-x + c1= c2 (2)
定义：f(x)=log(x!)
则有： f'(x)= 0.5+logx
若方程(2)的一个近似为x0,可导出牛顿迭代式：
x1= x0 + ( c2 -f(x0))/f'(x)
= x0 + ( c2- (x0+0.5)log(x0) +x0 - c1)/(0.5/x0+log(x0))
*/
// 解方程log(x!)=f,已知 low <= x <= high
// 返回x的值
DWORD faclogEquation(double f,int low,int high)
{
int i;
double x0,x1;
DWORD n;
if ( f <= logfac[128])
return searchInLogFacTab1(f);
x0= (low+high)/2;
#ifdef DEBUG_LOG_FAC
printf("---------\n");
printf("x0=%f\n",x0);
#endif
while(true)
{
x1= x0 + ( f- getLogfac((DWORD)x0))/(0.5/x0+log(x0));
#ifdef DEBUG_LOG_FAC
printf("x1=%f\n",x1);
#endif
if ( fabs(x1-x0)<0.75 )
break;
x0=x1;
}
#ifdef DEBUG_LOG_FAC
printf("---------\n");
#endif
if (x1<x0)
x0=x1;
n=(DWORD)x0;
while ( getLogfac(n) < f)
n++;
return n;
}
// 解方程log(x!)=f
// 返回x的值
DWORD faclogEquation(double f)
{
int i;
int low,high;
if ( f <= logfac[128])
return searchInLogFacTab1(f);
i=7;
while ( logfac_tab[i].logfac<f)
i++;
low= logfac_tab[i-1].n;
high=logfac_tab[i].n;
return faclogEquation(f,low,high);
}

复制代码

liangbch · 发表于 2008-4-17 10:50:32

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

liangbch · 发表于 2008-4-17 10:52:36

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

无心人 · 发表于 2008-4-17 11:03:28

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

无心人 · 发表于 2008-4-17 11:09:52

另外
任何FFT均慢于最优化的NTT

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

liangbch · 发表于 2008-4-17 11:10:09

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

文件 fft.h 的内容

#ifndef FFT_h
#define FFT_h
typedef double Real;
extern double FFTSquareWorstError;
extern long AllocatedMemory;
void InitializeFFT(long MaxFFTLength);
void MulWithFFT(Real * ACoef, long ASize,
Real * BCoef, long BSize,
Real * CCoef);
#endif

复制代码

文件 BigInt.h的内容

#ifndef BigInt_h
#define BigInt_h
#include "FFT.h"
typedef struct BigIntTag {
Real * Coef;
long Size, SizeMax;
} BigInt;
/*
* Decrease the base if the worst error in FFT become too large
*/
#define BASE 10000.
#define invBASE 0.0001
#define NBDEC_BASE 4
void InitializeBigInt(BigInt * A, long MaxSize);
void FreeBigInt(BigInt * A);
void PrintBigInt(BigInt * A);
void UpdateBigInt(BigInt * A);
void AddBigInt(BigInt * A, BigInt * B, BigInt * C);
void MulBigInt(BigInt * A, BigInt * B, BigInt * C);
void Inverse(BigInt * A, BigInt * B, BigInt * tmpBigInt);
#endif

复制代码

FFT.c 的内容

/*
* Xavier Gourdon : Sept. 99 (xavier.gourdon.free.fr)
*
* FFT.c : Very basic FFT file.
* The FFT encoded in this file is short and very basic.
* It is just here as a teaching file to have a first
* understanding of the FFT technique.
*
* A lot of optimizations could be made to save a factor 2 or 4 for time
* and space.
* - Use a 4-Step (or more) FFT to avoid data cache misses.
* - Use an hermitian FFT to take into account the hermitian property of
* the FFT of a real array.
* - Use a quad FFT (recursion N/4->N instead of N/2->N) to save 10 or
* 15% of the time.
*
* Informations can be found on
* http://xavier.gourdon.free.fr/Constants/constants.html
*/
#include "FFT.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define PI 3.1415926535897932384626
typedef struct ComplexTag {
Real R,I;
} Complex;
long FFTLengthMax;
Complex * OmegaFFT;
Complex * ArrayFFT0, * ArrayFFT1;
Complex * ComplexCoef;
double FFTSquareWorstError;
long AllocatedMemory;
void InitializeFFT(long MaxLength)
{
long i;
Real Step;
FFTLengthMax = MaxLength;
OmegaFFT = (Complex *) malloc(MaxLength/2*sizeof(Complex));
ArrayFFT0 = (Complex *) malloc(MaxLength*sizeof(Complex));
ArrayFFT1 = (Complex *) malloc(MaxLength*sizeof(Complex));
ComplexCoef = (Complex *) malloc(MaxLength*sizeof(Complex));
Step = 2.*PI/(double) MaxLength;
for (i=0; 2*i<MaxLength; i++) {
OmegaFFT[i].R = cos(Step*(double)i);
OmegaFFT[i].I = sin(Step*(double)i);
}
FFTSquareWorstError=0.;
AllocatedMemory = 7*MaxLength*sizeof(Complex)/2;
}
void RecursiveFFT(Complex * Coef, Complex * FFT, long Length, long Step,
long Sign)
{
long i, OmegaStep;
Complex * FFT0, * FFT1, * Omega;
Real tmpR, tmpI;
if (Length==2) {
FFT[0].R = Coef[0].R + Coef[Step].R;
FFT[0].I = Coef[0].I + Coef[Step].I;
FFT[1].R = Coef[0].R - Coef[Step].R;
FFT[1].I = Coef[0].I - Coef[Step].I;
return;
}
FFT0 = FFT;
RecursiveFFT(Coef ,FFT0,Length/2,Step*2,Sign);
FFT1 = FFT+Length/2;
RecursiveFFT(Coef+Step,FFT1,Length/2,Step*2,Sign);
Omega = OmegaFFT;
OmegaStep = FFTLengthMax/Length;
for (i=0; 2*i<Length; i++, Omega += OmegaStep) {
/* Recursion formula for FFT :
FFT[i] <- FFT0[i] + Omega*FFT1[i]
FFT[i+Length/2] <- FFT0[i] - Omega*FFT1[i],
Omega = exp(2*I*PI*i/Length) */
tmpR = Omega[0].R*FFT1[i].R-Sign*Omega[0].I*FFT1[i].I;
tmpI = Omega[0].R*FFT1[i].I+Sign*Omega[0].I*FFT1[i].R;
FFT1[i].R = FFT0[i].R - tmpR;
FFT1[i].I = FFT0[i].I - tmpI;
FFT0[i].R = FFT0[i].R + tmpR;
FFT0[i].I = FFT0[i].I + tmpI;
}
}
/* Compute the complex Fourier Transform of Coef into FFT */
void FFT(Real * Coef, long Length, Complex * FFT, long NFFT)
{
long i;
/* Transform array of real coefficient into array of complex */
for (i=0; i<Length; i++) {
ComplexCoef[i].R = Coef[i];
ComplexCoef[i].I = 0.;
}
for (; i<NFFT; i++)
ComplexCoef[i].R = ComplexCoef[i].I = 0.;
RecursiveFFT(ComplexCoef,FFT,NFFT,1,1);
}
/* Compute the inverse Fourier Transform of FFT into Coef */
void InverseFFT(Complex * FFT, long NFFT, Real * Coef, long Length)
{
long i;
Real invNFFT = 1./(Real) NFFT, tmp;
RecursiveFFT(FFT, ComplexCoef, NFFT, 1, -1);
for (i=0; i<Length; i++) {
/* Closest integer to ComplexCoef[i].R/NFFT */
tmp = invNFFT*ComplexCoef[i].R;
Coef[i] = floor(0.5+tmp);
if ((tmp-Coef[i])*(tmp-Coef[i])>FFTSquareWorstError)
FFTSquareWorstError = (tmp-Coef[i])*(tmp-Coef[i]);
}
}
void Convolution(Complex * A, Complex * B, long NFFT, Complex * C)
{
long i;
Real tmpR, tmpI;
for (i=0; i<NFFT; i++) {
tmpR = A[i].R*B[i].R-A[i].I*B[i].I;
tmpI = A[i].R*B[i].I+A[i].I*B[i].R;
C[i].R = tmpR;
C[i].I = tmpI;
}
}
void MulWithFFT(Real * ACoef, long ASize,
Real * BCoef, long BSize,
Real * CCoef)
{
long NFFT = 2;
while (NFFT<ASize+BSize)
NFFT *= 2;
if (NFFT>FFTLengthMax) {
printf("Error, FFT Size is too big in MulWithFFT\n");
return;
}
FFT(ACoef, ASize, ArrayFFT0, NFFT);
FFT(BCoef, BSize, ArrayFFT1, NFFT);
Convolution(ArrayFFT0,ArrayFFT1,NFFT,ArrayFFT0);
InverseFFT(ArrayFFT0,NFFT,CCoef, ASize+BSize-1);
}

复制代码

BigInt.c的内容

/*
* Xavier Gourdon : Sept. 99 (xavier.gourdon.free.fr)
*
* BigInt.c : Basic file for manipulation of Large Integers.
* It rely on FFT.c to multiply number with the FFT
* technique.
*
* Several optimizations can be made :
* - Implementation of a classic multiplication for small BigInts.
* - Avoid using UpdateBigInt for the sum of BigInts.
* - Use a larger base for the internal data storage of BigInts to
* save space.
* ...
*
* Informations can be found on
* http://xavier.gourdon.free.fr/Constants/constants.html
*/
#include "BigInt.h"
#include "FFT.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
void InitializeBigInt(BigInt * A, long MaxSize)
{
A->Coef = (Real *) malloc(MaxSize*sizeof(Real));
AllocatedMemory += MaxSize*sizeof(Real);
if (A->Coef==0)
printf("Not enough memory\n");
A->Size = 0;
A->SizeMax = MaxSize;
}
void FreeBigInt(BigInt * A)
{
free(A->Coef);
}
void PrintBigInt(BigInt * A)
{
long i,j,Digit=0,Dec;
double Pow, Coef;
printf("%.0lf\n",A->Coef[A->Size-1]);
for (i=A->Size-2; i>=0; i--) {
Pow = BASE*0.1;
Coef = A->Coef[i];
for (j=0; j<NBDEC_BASE; j++) {
Dec = (long) (Coef/Pow);
Coef -= Dec*Pow;
Pow *= 0.1;
printf("%ld",Dec);
Digit++;
if (Digit%10==0) printf(" ");
if (Digit%50==0) printf(": %ld\n",Digit);
}
}
}
/*
* Update A to the normal form (0<= A.Coef[i] < BASE)
*/
void UpdateBigInt(BigInt * A)
{
long i;
Real carry=0., x;
for (i=0; i<A->Size; i++) {
x = A->Coef[i] + carry;
carry = floor(x*invBASE);
A->Coef[i] = x-carry*BASE;
}
if (carry!=0) {
while (carry!=0.) {
x = carry;
carry = floor(x*invBASE);
A->Coef[i] = x-carry*BASE;
i++;
A->Size=i;
}
if (A->Size > A->SizeMax) {
printf("Error in UpdateBigInt, Size>SizeMax\n");
}
}
else {
while (i>0 && A->Coef[i-1]==0.) i--;
A->Size=i;
}
}
/*
* Compute C = A + B
*/
void AddBigInt(BigInt * A, BigInt * B, BigInt * C)
{
long i;
if (A->Size<B->Size) {
AddBigInt(B,A,C);
return;
}
/* We now have B->Size<A->Size */
for (i=0; i<B->Size; i++)
C->Coef[i] = A->Coef[i] + B->Coef[i];
for (; i<A->Size; i++)
C->Coef[i] = A->Coef[i];
C->Size = A->Size;
UpdateBigInt(C);
}
/*
* Compute C = A*B
*/
void MulBigInt(BigInt * A, BigInt * B, BigInt * C)
{
MulWithFFT(A->Coef,A->Size,B->Coef,B->Size,C->Coef);
C->Size = A->Size+B->Size-1;
UpdateBigInt(C);
}
/*
* Compute the inverse of A in B
*/
void Inverse(BigInt * A, BigInt * B, BigInt * tmpBigInt)
{
double x;
long i, N,NN,Delta;
int Twice=1, Sign;
BigInt AA;
/* Initialization */
x = A->Coef[A->Size-1]+invBASE*(A->Coef[A->Size-2]
+ invBASE*A->Coef[A->Size-3]);
x = BASE/x;
B->Coef[1] = floor(x);
B->Coef[0] = floor((x-B->Coef[1])*BASE);
B->Size = 2;
/* Iteration used : B <- B+(1-AB)*B */
N=2;
while (N<A->Size) {
/* Use only the first 2*N limbs of A */
NN = 2*N;
if (NN>A->Size)
NN = A->Size;
AA.Coef = A->Coef+A->Size-NN;
AA.Size = NN;
MulBigInt(&AA,B,tmpBigInt);
Delta = NN+B->Size-1;
/* Compute BASE^Delta-tmpBigInt in tmpBigInt */
if (tmpBigInt->Size==Delta) {
Sign=1;
for (i=0; i<Delta; i++)
tmpBigInt->Coef[i] = BASE-1 - tmpBigInt->Coef[i];
UpdateBigInt(tmpBigInt);
}
else {
Sign = -1;
tmpBigInt->Coef[Delta] = 0.;
}
MulBigInt(tmpBigInt,B,tmpBigInt);
for (i=0; i<tmpBigInt->Size-2*N+1; i++)
tmpBigInt->Coef[i] = tmpBigInt->Coef[i+2*N-1];
tmpBigInt->Size -= 2*N-1;
for (i=B->Size-1; i>=0; i--)
B->Coef[i+NN-N] = B->Coef[i];
for (i=NN-N-1; i>=0; i--)
B->Coef[i]=0.;
B->Size += NN-N;
if (Sign==-1) {
for (i=0; i<tmpBigInt->Size; i++)
tmpBigInt->Coef[i] = -tmpBigInt->Coef[i];
}
AddBigInt(B,tmpBigInt,B);
/* Once during the process, redo one iteration with same precision */
if (8*N>A->Size && Twice) {
Twice=0;
B->Size = N;
for (i=0; i<N; i++)
B->Coef[i] = B->Coef[i+N];
}
else
N *= 2;
}
}

复制代码

主程序(e.c)的内容

/*
* Xavier Gourdon : Sept. 99 (xavier.gourdon.free.fr)
*
* e.c : Basic file for computation of e=2.7181... with the binary
* splitting method.
* It relies on the files FFT.c and BigInt.c to handle
* Large Integers.
*
* Several optimizations can be made :
* - Use a specialized routine to handle binary splitting when the
* indexes are close.
* ...
*
* Informations can be found on
* http://xavier.gourdon.free.fr/Constants/constants.html
*/
#include "FFT.h"
#include "BigInt.h"
#include <stdio.h>
#include <math.h>
#include <time.h>
BigInt tmpBigInt; /* Temporary BigInt used in the binary splitting */
/*
* Compute the numerator P and denominator Q of
* P/Q = 1/(N0+1) + 1/(N0+1)/(N0+2) + ... + 1/(N0+1)/.../N1
*/
void BinarySplittingE(long N0, long N1, BigInt * P, BigInt * Q)
{
BigInt PP, QQ;
long NMid;
if (N1-N0==1) {
P->Size = Q->Size = 1;
P->Coef[0] = 1.;
Q->Coef[0] = (double) N1;
UpdateBigInt(P);
UpdateBigInt(Q);
return;
}
NMid = (N0+N1)/2;
BinarySplittingE(N0,NMid,P,Q);
/* To save memory, take the non used coefficients of P and Q
for coefficient of the BigInt used in the second splitting part */
PP.Coef = P->Coef + P->Size;
PP.SizeMax = P->SizeMax-P->Size;
QQ.Coef = Q->Coef + Q->Size;
QQ.SizeMax = Q->SizeMax-Q->Size;
BinarySplittingE(NMid,N1,&PP,&QQ);
MulBigInt(P,&QQ,&tmpBigInt);
AddBigInt(&tmpBigInt,&PP,P);
MulBigInt(Q,&QQ,Q);
}
/*
* Returns as a BigInt the constant E with NbDec decimal digits
*/
BigInt ECompute(long NbDec)
{
BigInt P, Q, tmp;
long i, MaxSize, MaxFFTSize, SeriesSize;
double logFactorial, logMax;
/* MaxSize should be more than NbDec/NBDEC_BASE (see BinarySplitting) */
MaxSize = NbDec/NBDEC_BASE + 10 + (long) (2.*log((double)NbDec)/log(2.));
InitializeBigInt(&P,MaxSize);
InitializeBigInt(&Q,MaxSize);
MaxFFTSize = 2;
while (MaxFFTSize < MaxSize)
MaxFFTSize *= 2;
MaxFFTSize *= 2;
InitializeFFT(MaxFFTSize);
/* Temporary BigInts are needed */
InitializeBigInt(&tmpBigInt,MaxFFTSize);
InitializeBigInt(&tmp,MaxFFTSize);
printf("Total Allocated memory = %ld K\n",AllocatedMemory/1024);
/* Compute the number of term SeriesSize of the series required.
SeriesSize is such that log(SeriesSize !) > NbDec*log(10.) */
SeriesSize=1;
logFactorial=0.;
logMax = (double)NbDec * log(10.);
while (logFactorial<logMax) {
logFactorial += log((double) SeriesSize);
SeriesSize++;
}
printf("Starting series computation\n");
/* Compute the numerator P and the denominator Q of
sum_{i=0}^{SeriesSize} 1/i! */
BinarySplittingE(0,SeriesSize,&P,&Q);
AddBigInt(&P,&Q,&P);
printf("Starting final division\n");
/* Compute the inverse of Q in tmpBigInt */
Inverse(&Q,&tmpBigInt,&tmp);
/* Comute P/Q in tmpBigInt */
MulBigInt(&P,&tmpBigInt,&tmpBigInt);
/* Put the number of required digits in P */
P.Size = 1+NbDec/NBDEC_BASE;
for (i=1; i<=P.Size; i++)
P.Coef[P.Size-i] = tmpBigInt.Coef[tmpBigInt.Size-i];
FreeBigInt(&Q);
FreeBigInt(&tmpBigInt);
return P;
}
void main()
{
double StartTime;
BigInt E;
long NbDec;
printf("*** E computation *** \n\n");
printf("Enter the number of decimal digits : ");
scanf("%ld",&NbDec);
printf("\n");
StartTime = (double) clock();
E = ECompute(NbDec);
printf("\n");
printf("Total time : %.2lf seconds\n",((double) clock() - StartTime)/CLOCKS_PER_SEC);
printf("Worst error in FFT (should be less than 0.25): %.10lf\n",sqrt(FFTSquareWorstError));
printf("\n");
printf("E = ");
//PrintBigInt(&E);
printf("\n");
}

复制代码

shshsh_0510 · 发表于 2008-4-17 11:12:59

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

无心人 · 发表于 2008-4-17 11:14:20

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

账号		自动登录	找回密码
密码			欢迎注册

[擂台] 计算百万位e