/***********************************************************************/ /* */ /* */ /* Program: */ /* -------- */ /* */ /* sseii.c, Calculates e = exp(1), Long long integer (64 bit) method. */ /* */ /* */ /* Version: */ /* -------- */ /* */ /* First version: 1998-02-19 */ /* 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 the 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 sseii sseii.c -lm */ /* */ /* The time wraps around after 4295 seconds on some computers. */ /* */ /* */ /* Output: */ /* ------- */ /* */ /* Output is to the file "sseii.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. */ /* */ /* The program use long long integer (64 bit) division and gives */ /* probably correct result up to 2*10^9 digits of e (slightly less */ /* than maxint = 2^31 digits). But it will take a VERY long time ! */ /* The result has however only been verified up to 10^6 digits. */ /* */ /* The base for calculations is 10^8 (integer 32 bit arrays). */ /* */ /* The program code is optimized for speed by the author. */ /* */ /* */ /* For 64 bit CPUs: */ /* ---------------- */ /* */ /* This program might be very fast on a computer with a 64 bit */ /* integer calculation unit in the CPU. */ /* */ /* long long int, Means a 64 bit integer. */ /* int, Means a 32 bit integer. */ /* */ /* If necessary, change this to appropriate declarations in */ /* the code below. See comments in the code ! */ /* */ /* */ /* 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 18.33 gcc -O3 */ /* na Celeron 1200 MHz 17.69 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 sseii sseii.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) */ /*---------------------------------------------------------------------*/ /* sseii 59.87 s 5095.05 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. */ /* */ /* */ /* 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 +O3). */ /* */ /* The CPU-times for calculations of e follows in the table below. */ /* Observe that the times can vary a few percent depending load etc. */ /* */ /* Observe that these are 32 bit computers which means that */ /* long long int (64 bit) integer calculations are slow ! */ /* */ /* ------------------------------------------------------------------ */ /* P 180 MMX OD HP 735 */ /* Digits (s) (h m s) (s) (h m s) */ /* ------------------------------------------------------------------- */ /* 10^5 193 00 03 13 816 00 13 36 */ /* 10^6 16954 04 42 34 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]="sseii.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; /* 32 (or 64) bit integer */ int imin,il,ll,ll2; /* 32 (or 64) bit integer */ int remi,qi,basi; /* 32 (or 64) bit integer */ long long int qii,sii,basii,jii,remii; /* Must be 64 bit integer ! */ double lt,l2f,sum,jd,lld; clock_t time1,timediff1; int *v1,*v2; /* Arrays of 32 bit integer */ 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)); /* 32 bit integer array */ v2=(int *)calloc(l2+1,sizeof(int)); /* 32 bit integer array */ 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("Long long integer (64 bit) method, numerically stable.\n"); printf("\n"); printf("Factorial=%8d; Length=%8d; Decimals=%8d\n",fac,il,n); printf("\n"); printf("Calculating ...\n"); basi=100000000; basii=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); lld=lld+lt*log(jd); ll=(int)(0.125*lld)-3; if (ll > 0) { ll2=l2-ll; } else { ll2=l2; } /* if */ jii=(long long int)(j); /* jii is a 64 bit integer */ remii=0; /* remii is a 64 bit integer */ for (i=ll2; i>=1; i--) { qii=(long long int)(v1[i])+basii*remii; /* 64 bit integer calculation */ sii=qii / jii; /* 64 bit integer calculation */ v1[i]=(int)(sii); /* 64 bit integer converted to 32 bit */ remii=qii-jii*sii; /* 64 bit integer calculation */ } /* 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 */