/***********************************************************************/ /* */ /* */ /* Program: */ /* -------- */ /* */ /* hex2dec.c */ /* */ /* This program is included in a package, consisting of 6 fast */ /* programs, to calculate e=exp(1) hexadecimal. Programs to */ /* convert from hexadecimal to decimal are also included. */ /* */ /* */ /* Version: */ /* -------- */ /* */ /* First version: November 1998 */ /* Latest update: 2003-06-15 */ /* */ /* */ /* Author of the program: */ /* ---------------------- */ /* */ /* Stefan Spaennare, Lund, Sweden */ /* */ /* E-mail: stefan@spaennare.se */ /* */ /* Home page: http://www.spaennare.se/index.html */ /* Program page: http://www.spaennare.se/ssprog.html */ /* */ /* If you have problems with the program please let me know! */ /* */ /* */ /* References: */ /* ----------- */ /* */ /* 1. The program now use a very fast FFT (fftsg_h.c) by Takuya Ooura: */ /* */ /* http://momonga.t.u-tokyo.ac.jp/~ooura/ */ /* */ /* */ /* 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: */ /* --------- */ /* */ /* How to calculate e=exp(1) hexadecimal and decimal: */ /* */ /* First calculate the hexadecimal numbers a and b using one of the */ /* ssehex.c, ssehex2.c or ssehex3.c. Here n=100000 is the number */ /* of desired DECIMAL digits. Type the following: */ /* */ /*>ssehex 100000 */ /* */ /* Then calculate e=exp(1) hexadecimal by e=a/b: */ /* */ /*>hexediv ehex_a.dat ehex_b.dat e5.hex */ /* */ /* Then convert e=exp(1) hexadecimal to decimal by: */ /* */ /*>hex2dec e5.hex e5.dec */ /* */ /* or */ /* */ /*>fhex2dec e5.hex e5.dec */ /* */ /* */ /* Compilation: */ /* ------------ */ /* */ /* On Linux computers just run "compall" to compile all six programs: */ /* */ /*>compall */ /* */ /* The program works in UNIX, Linux and Windows (DOS using DJGPP) */ /* It is recommended to optimize the program (i.e. use cc -O3 option */ /* or similar option) when you compile. */ /* */ /*>gcc -O3 -o program program.c -lm */ /* */ /* or (for those programs using FFT) */ /* */ /*>gcc -O3 -o program program.c fftsg_h.c -lm */ /* */ /* The time wraps around after 4295 seconds on some computers. */ /* */ /* */ /* Hexadecimal numbers: */ /* -------------------- */ /* */ /* Hexadecimal numbers that can be converted to decimal with the */ /* programs hex2dec.c and fhex2dec.c must be on the following */ /* form (0 <= X <= 9): */ /* */ /* X.A1B678F3EC95D6C4... */ /* */ /* */ /* Algorithm and accuracy: */ /* ----------------------- */ /* */ /* The programs ssehex.c, ssehex2.c and ssehex3.c calculates */ /* SUM 1/k! using the algorithm (C-like notation): */ /* */ /* b=1; */ /* a=0; */ /* */ /* for (i=n; i>=1; i--) { */ /* b=b*i; */ /* a=a+b; */ /* } */ /* */ /* e=a/b; */ /* */ /* Multiplications and divisions with the base (2^32) are made using */ /* the shift operators (<< and >>) and the AND operator (&&) to be */ /* fast. This is also the reason why the calculations are performed */ /* with a hexadecimal (binary) base. */ /* */ /* The final computation of e=a/b is made with the separate program */ /* hexediv.c using Newtons Iteration and FFT-multiplication */ /* for division. */ /* */ /* The programs use long long integer (64 bit) calculations and gives */ /* probably correct result up to 2*10^9 digits of e=exp(1) (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^7 digits. */ /* */ /* Note! If the final decimal result is expected to have n digits */ /* The actual calculation is performed using slightly more digits. */ /* All calculated digits are output from these programs, but only */ /* the n desired decimal digits can be expected to be correct. */ /* */ /* WARNING! */ /* -------- */ /* */ /* This is valid for hexediv.c and fhex2dec.c: */ /* */ /* If you want to make calculations with much more than 4 million */ /* (2^22) digits the convolutions in the FFT multiplications can be */ /* saturated which will give erroneous results or no convergence */ /* at all. If this occur then set mulversion=2 which */ /* is 2 times slower than mulversion=1 but use the base 16 */ /* instead of 256. This means that the convolution will be */ /* saturated for a much larger number of digits (> 100 millions). */ /* */ /* */ /* Timing, benchmark and memory: */ /* ----------------------------- */ /* */ /* The CPU-times are given for a 1400 MHz Intel Celeron computer */ /* (100 MHz memory buss) and 512 Mbyte memory. The cache memory is */ /* 256 kbyte at full CPU speed. The operating system was Red Hat */ /* Linux 9 with the gcc 3.22 C-compiler. The CPU load was on */ /* average 1.00 (no other programs running). */ /* */ /* The programs were compiled using the following lines: */ /* */ /* For ssehex.c, ssehex2.c, ssehex3.c and hex2dec.c: */ /* */ /*>gcc -O3 -o program program.c -lm */ /* */ /* For hexediv.c and fhex2dec.c: */ /* */ /*>gcc -O3 -o program program.c fftsg_h.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) */ /*---------------------------------------------------------------------*/ /* ssehex 3.46 s 791.00 s O(n^2) 1.2 Mbyte */ /* ssehex2 2.86 s 607.09 s O(n^2) 1.2 Mbyte */ /* ssehex3 3.78 s 707.71 s O(n^2) 1.2 Mbyte */ /* hexediv 5.32 s 62.24 s O(n*log(n)^2) 21.0 Mbyte */ /* hex2dec 1.35 s 266.02 s O(n^2) 0.9 Mbyte */ /* fhex2dec 12.48 s 157.51 s O(n*log(n)^3) 36.0 Mbyte */ /*---------------------------------------------------------------------*/ /* */ /***********************************************************************/ /********************************************************************************/ /* */ /* 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 #include double timetic=(double)(CLOCKS_PER_SEC); void printdec(ss,outfile) unsigned long ss; FILE *outfile; { if (ss > 999) { fprintf(outfile,"%4d",ss); } else if (ss > 99) { fprintf(outfile,"0%3d",ss); } else if (ss > 9) { fprintf(outfile,"00%2d",ss); } else { fprintf(outfile,"000%1d",ss); } /* if */ } /* printdec */ unsigned char ppp(s) unsigned char s; { unsigned char ss; if ((s >= '0') && (s <='9')) { ss=s-48; } else { if (s == 'A') { ss=0x0a; } else if (s == 'B') { ss=0x0b; } else if (s == 'C') { ss=0x0c; } else if (s == 'D') { ss=0x0d; } else if (s == 'E') { ss=0x0e; } else if (s == 'F') { ss=0x0f; } /* if */ } /* if */ return(ss); } /* ppp */ int main(argc,argv) int argc; char *argv[]; { FILE *infile,*outfile; int n,n2,m,i,j,ss,n22,ii; unsigned long aa; double sum,lt; unsigned long long ireg; unsigned char cc,hh,cc0; unsigned long *a; unsigned short *s; clock_t t1,td1; if (argc != 3) { printf("Usage: %s hexfile decfile\n",argv[0]); exit(0); } /* if */ if ((infile = fopen(argv[1],"r")) == NULL) { printf("File not found.\n"); exit(0); } /* if */ rewind(infile); fseek(infile,0,SEEK_END); n=ftell(infile)-3; n2=n/8; a=(unsigned long *)calloc(n2+10,sizeof(unsigned long)); m=(int)(2.408*(double)(n2))+1; s=(unsigned short *)calloc(m+10,sizeof(unsigned short)); printf("hex= %9d dec= %9d n2= %9d\n",n,4*m,n2); rewind(infile); cc0=getc(infile); cc=getc(infile); for (i=n2; i>=1; i--) { aa=0; for (j=7; j>=0; j--) { cc=getc(infile); hh=ppp(cc); aa=aa+((unsigned long)(hh) << (4*j)); } /* for j */ a[i]=aa; } /* for i */ lt=log(10000.0)/(32.0*log(2.0)); n22=n2+1; t1=clock(); sum=0.0; ii=1; for (j=1; j<=m; j++) { ireg=0; for (i=ii; i<=n22; i++) { ireg=(unsigned long long)(a[i])*10000+ireg; a[i]=(unsigned long)(ireg & 0x00000000ffffffff); ireg=ireg >> 32; } /* for i */ s[j]=a[n22]; a[n22]=0; sum=sum+lt; ii=(int)(sum)-6; if (ii < 1) { ii=1; } /* if */ } /* for j */ td1=clock()-t1; printf("%10.2f\n",(double)(td1)/timetic); outfile=fopen(argv[2],"w"); putc(cc0,outfile); putc('.',outfile); fprintf(outfile,"\n"); j=1; for (i=1; i<=m; i++) { printdec(s[i],outfile); if (j == 18) { fprintf(outfile,"\n"); j=0; } /* if */ j++; } /* for i */ fprintf(outfile,"\n"); fclose(infile); fclose(outfile); } /* End */