/***********************************************************************/ /* */ /* */ /* Program: */ /* -------- */ /* */ /* sseff.c, Calculates e = exp(1), Floating point method. */ /* */ /* */ /* Version: */ /* -------- */ /* */ /* First version: 1998-02-16 */ /* Latest update: 2003-02-19 */ /* */ /* */ /* Author of the program: */ /* ---------------------- */ /* */ /* Stefan Spaennare, Lund, Sweden */ /* */ /* E-mail: stefan@spaennare.se */ /* */ /* If you have problems with the program please let me know! */ /* */ /* */ /* Copy-right: */ /* ----------- */ /* */ /* This program is free to copy if this information follows the */ /* program. */ /* */ /* */ /* IMPORTANT ! */ /* ----------- */ /* */ /* Please don't use the result from this program for important */ /* purposes without doing an independent check of the result! */ /* See also notice below ! */ /* */ /* */ /* Function: */ /* --------- */ /* */ /* This C program calculates e = exp(1) to many decimal places. */ /* */ /* The parameter ni is the number of decimals. */ /* */ /* */ /* Compilation: */ /* ------------ */ /* */ /* The program works in UNIX and Linux and I recommend that you */ /* optimize the program (i.e. use cc -O2 option or similar option) */ /* when you compile (gcc if you use a GNU compiler). */ /* */ /*> cc -O2 -o sseff sseff.c -lm */ /* */ /* The time wraps around after 4295 seconds on some computers. */ /* */ /* */ /* Output: */ /* ------- */ /* */ /* Output is to the file "sseff.dat". This file name can easily be */ /* changed. */ /* */ /* */ /* Algorithm and accuracy: */ /* ----------------------- */ /* */ /* The algorithm is from a student exercise in Pascal programming */ /* at Lund Institute of Technology in February 1998. */ /* */ /* This version probably gives correct result if less than */ /* 45*10^6 ! (factorial) fits into the n digits. This is because */ /* (2^52) / (10^8) = 45*10^6 approximatley. A double precision */ /* floating point number is here supposed to have 52 bits in */ /* the fraction part following the IEEE standard. */ /* */ /* This means that e probably can be calculated with up to 320*10^6 */ /* correct digits with this program. But it will take a very long */ /* time ! The program has however only been verified to 10^6 */ /* (one million) digits. */ /* */ /* The base for calculations is 10^8 (integer arrays). */ /* */ /* The program code is optimized for speed by the author. */ /* */ /* */ /* Benchmark */ /* ========= */ /* */ /* The test was run on a Intel Celeron (Tualatin core) computer */ /* with 256 kbyte cache memory and 512 Mbyte primary memory. */ /* The operating system was Red Hat 7.2 Linux. Compiler tested */ /* was GNU GCC 2.96 (gcc) for Linux. Under Windows XP the compiler */ /* GNU DJGPP 3.0.4 (gcc, djgpp) for DOS was used. */ /* */ /* Numbers tested: */ /* --------------- */ /* */ /* na = 100000 = 10^5 = 100*10^3 */ /* */ /* Results: */ /* -------- */ /* */ /* n: Computer: CPU-time (s): Compiler options: */ /* ================================================================ */ /* */ /* na Celeron 1200 MHz 28.93 gcc -O3 */ /* na Celeron 1200 MHz 118.41 gcc -O2 (djgpp) */ /* */ /* */ /* Timing, benchmark and memory: */ /* ----------------------------- */ /* */ /* The CPU-times are given for a 350 MHz Pentium II computer */ /* (100 MHz memory buss) and 256 Mbyte memory. The cache memory is */ /* 512 kbyte at half CPU speed. The operating system was Red Hat */ /* Linux 6.0 with the egcs 2.91.66 (egcs-1.1.2 release) C-compiler. */ /* The CPU load was on average 1.00 (no other programs running). */ /* */ /* The programs were compiled using the following line: */ /* */ /*>egcs -O2 -fomit-frame-pointer -funroll-loops -o sseff sseff.c -lm */ /* */ /* How to compile see also the file "compall". TC below means */ /* Time Complexity. The CPU-times are given for n=10^5=100000 and */ /* and n=10^6=1000000 decimal digits. */ /* */ /* */ /* Program 10^5 10^6 TC Memory/(10^6) */ /*---------------------------------------------------------------------*/ /* sseff 93.43 s 7809.49 s O(n^2) 1.4 Mbyte */ /*---------------------------------------------------------------------*/ /* */ /* */ /* Old timing and benchmarks: */ /* -------------------------- */ /* */ /* The time depends on the factorial that fits into the n digits. */ /* Empirically this means that the computing time approximately */ /* increases with a factor of four if the number of digits doubles. */ /* */ /* Unfortunately the fmod function is very slow on some (HP) */ /* computers. */ /* */ /* The program is well adapted for benchmark testing and */ /* has been tested on two different computers: */ /* */ /* 1: A DELL DIMENSION XPS P90 computer with a 180 MHz Pentium MMX */ /* overdrive processor, 96 Mbyte memory and 256 kbyte cache. */ /* The operating system was Linux Red Hat 4.1 and the code was */ /* fully optimized (gcc -O). */ /* */ /* 2: A HP 735 computer (99 MHz) with 144 Mbyte memory and 256 kbyte */ /* cache. The operating system was HP UNIX and the code was fully */ /* optimized (cc +O4). */ /* */ /* The CPU-times for calculations of e follows in the table below. */ /* Observe that the times can vary a few percent depending load etc. */ /* */ /* ------------------------------------------------------------------ */ /* P 180 MMX OD HP 735 */ /* Digits (s) (h m s) (s) (h m s) */ /* ------------------------------------------------------------------- */ /* 10^5 176 00 02 56 1225 00 20 25 */ /* 10^6 15747 04 22 27 Not calc. */ /* ------------------------------------------------------------------- */ /* */ /***********************************************************************/ /********************************************************************************/ /* */ /* Notice */ /* ====== */ /* */ /* I make no warranties that this program is (1) free of errors, (2) consistent */ /* with any standard merchantability, or (3) meeting the requirements of a */ /* particular application. This software shall not, partly or as a whole, */ /* participate in a process, whose outcome can result in injury to a person or */ /* loss of property. It is solely designed for analytical work. Permission to */ /* use, copy and distribute is hereby granted without fee, providing that the */ /* header above including this notice appears in all copies. */ /* */ /* Stefan Spaennare */ /* */ /********************************************************************************/ #include #include #include #include #include double timetic=(double)(CLOCKS_PER_SEC); char fname[80]="sseff.dat"; void fppp(t,outfile) int t; FILE *outfile; { if (t > 9999999) { fprintf(outfile,"%8d",t); } else if (t > 999999) { fprintf(outfile,"0%7d",t); } else if (t > 99999) { fprintf(outfile,"00%6d",t); } else if (t > 9999) { fprintf(outfile,"000%5d",t); } else if (t > 999) { fprintf(outfile,"0000%4d",t); } else if (t > 99) { fprintf(outfile,"00000%3d",t); } else if (t > 9) { fprintf(outfile,"000000%2d",t); } else { fprintf(outfile,"0000000%1d",t); } /* if */ } /* fppp */ int main(argc,argv) int argc; char *argv[]; { FILE *outfile; int i,j,n,l1,l2,fac; int imin,il,ll,ll2; int remi,qi,basi; double lt,l2f,sum,jd,jdinv,lld; double rem,q,bas; clock_t time1,timediff1; int *v1,*v2; if (argc != 2) { printf("Usage: %s ni\n",argv[0]); exit(0); } /* if */ n=atoi(argv[1]); l1=n/8; l2=(int)(1.005*(double)(l1)+0.5)+5; n=8*l1; v1=(int *)calloc(l2+1,sizeof(int)); v2=(int *)calloc(l2+1,sizeof(int)); lt=1.0/log(10.0); l2f=8.0*(double)(l2); i=1; sum=0.0; while (sum < l2f) { sum=sum+lt*log((double)(i)); i++; } /* while */ fac=i; il=(int)(sum); printf("\n"); printf("e = exp(1)\n"); printf("==========\n"); printf("\n"); printf("Floating point method, probably numerically stable.\n"); printf("\n"); printf("Factorial=%8d; Length=%8d; Decimals=%8d\n",fac,il,n); printf("\n"); printf("Calculating ...\n"); bas=100000000.0; basi=100000000; for (i=1; i<=l2; i++) { v1[i]=0; v2[i]=0; } /* for i */ v1[l2]=1; v2[l2]=1; lld=0.0; time1=clock(); for (j=1; j<=fac; j++) { jd=(double)(j); jdinv=1.0/jd; lld=lld+lt*log(jd); ll=(int)(0.125*lld)-3; if (ll > 0) { ll2=l2-ll; } else { ll2=l2; } /* if */ rem=0.0; for (i=ll2; i>=1; i--) { q=(double)(v1[i])+bas*rem; rem=(double)((int)(fmod(q,jd)+0.5)); v1[i]=(int)(jdinv*(q-rem)+0.5); } /* for i */ remi=0; for (i=1; i<=ll2; i++) { qi=v1[i]+v2[i]+remi; if (qi >= basi) { v2[i]=qi-basi; remi=1; } else { v2[i]=qi; remi=0; } /* if */ } /* for i */ if (remi > 0) { for (i=ll2+1; i<=l2; i++) { qi=v1[i]+v2[i]+remi; if (qi >= basi) { v2[i]=qi-basi; remi=1; } else { v2[i]=qi; remi=0; } /* if */ } /* for i */ } /* if */ } /* for j */ timediff1=clock()-time1; printf("\n"); printf("CPU time=%10.2f s\n",(double)(timediff1)/timetic); printf("\n"); printf("Calculation finished, output to %s.\n\n",fname); outfile=fopen(fname,"w"); fprintf(outfile,"%d.\n",v2[l2]); imin=l2-l1; j=1; for (i=l2-1; i>=imin; i--) { fppp(v2[i],outfile); if ((j % 9) == 0) { fprintf(outfile,"\n"); } /* if */ j++; } /* for i */ fclose(outfile); } /* End */