/***********************************************************************/ /* */ /* */ /* Program: */ /* -------- */ /* */ /* fhex2dec.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); int mulversion=1; double bas=256.0; double basi=0.00390625; double pi=3.14159265358979323; void mul(d,n) double *d; int n; { int nn,nn3,nn2,j; double x1,x2,y1,y2; nn=n+1; nn2=2*n+2; nn3=nn2+1; x1=d[0]; y1=d[1]; x2=y1; y1=0.0; y2=0.0; d[0]=x1*x2-y1*y2; d[1]=x1*y2+y1*x2; for (j=3; j<=nn; j=j+2) { x1=0.5*(d[j-1]+d[nn2-j-1]); x2=0.5*(d[j-1]-d[nn2-j-1]); y1=0.5*(d[j+1-1]+d[nn3-j-1]); y2=0.5*(d[j+1-1]-d[nn3-j-1]); d[j-1]=x1*y1-y2*(-x2); d[j+1-1]=x1*(-x2)+y2*y1; d[nn2-j-1]=x1*y1-(-y2)*x2; d[nn3-j-1]=x1*x2+(-y2)*y1; } /* for j */ } /* mul */ void lzero(a,nn) unsigned char *a; int *nn; { int i,l,il,dl,ib; l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; for (i=l;i>=ib; i--) { a[i]=0; } /* for i */ } /* lzero */ void linttol(a,man0,exp0,nn) unsigned char *a; int man0,exp0; int *nn; { int i,l,il,dl,ib; l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; for (i=l;i>=ib; i--) { a[i]=0; } /* for i */ a[l-il+exp0+1]=abs(man0); } /* linttol */ void lass(a,b,nn) unsigned char *a; unsigned char *b; int *nn; { int i,l,il,dl,ib; l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; for (i=l+1; i>=ib; i--) { a[i]=b[i]; } /* for i */ } /* lass */ void add(a,b,c,nn) unsigned char *a; unsigned char *b; unsigned char *c; int *nn; { int i,ib,il,dl,l; int q,r,basint; l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; basint=256; r=0; for (i=ib; i<=l; i++) { q=(int)(b[i])+(int)(c[i])+r; if (q >= basint) { a[i]=(unsigned char)(q-basint); r=1; } else { a[i]=(unsigned char)(q); r=0; } /* if */ } /* for i */ } /* add */ void sub(a,b,c,nn) unsigned char *a; unsigned char *b; unsigned char *c; int *nn; { int i,ib,il,dl,l; int q,r,basint; l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; basint=256; r=0; for (i=ib; i<=l; i++) { q=(int)(b[i])-(int)(c[i])-r; if (q < 0) { a[i]=(unsigned char)(basint+q); r=1; } else { a[i]=(unsigned char)(q); r=0; } /* if */ } /* for i */ } /* sub */ void sub2(a,b,c,i1,in,j1,k1) unsigned char *a; unsigned char *b; unsigned char *c; int i1,in,j1,k1; { int i,j,k; int q,r,basint; basint=256; r=0; j=j1; k=k1; for (i=i1; i<=in; i++) { q=(int)(b[j])-(int)(c[k])-r; if (q < 0) { a[i]=(unsigned char)(basint+q); r=1; } else { a[i]=(unsigned char)(q); r=0; } /* if */ j++; k++; } /* for i */ } /* sub2 */ void shift(a,b,i1,in,j1) unsigned char *a; unsigned char *b; int i1,in,j1; { int i,j; j=j1; for (i=i1; i<=in; i++) { a[i]=b[j]; j++; } /* for i */ } /* shift */ void shift2(a,b,i1,in,j1) unsigned char *a; unsigned char *b; int i1,in,j1; { int i,j; j=j1; for (i=in; i>=i1; i--) { a[i]=b[j]; j--; } /* for i */ } /* shift2 */ void shift3(a,b,i1,in,j1) unsigned char *a; unsigned char *b; int i1,in,j1; { int i,j; j=j1; for (i=i1; i<=in; i++) { a[j]=b[i]; j++; } /* for i */ } /* shift3 */ void shift4(a,b,i1,in,j1) unsigned char *a; unsigned char *b; int i1,in,j1; { int i,j; j=j1; for (i=in; i>=i1; i--) { a[j]=b[i]; j--; } /* for i */ } /* shift4 */ int lcompeps(a,nn) unsigned char *a; int *nn; { int i,l,il,dl,k,ib,n,ll; l=nn[0]; il=nn[1]; dl=nn[2]; n=nn[3]; ib=l-il-dl+1; k=0; ll=l-il-n; for (i=l; i>=ll; i--) { if (a[i]==0) { k++; } /* if */ } /* for i */ if (k==(il+n+1)) { return(-1); } else { return(1); } /* if */ } /* lcompeps */ void lcmul(a,b,c,nn,exp0) unsigned char *a; int b; unsigned char *c; int *nn; int exp0; { int i,j,k,l,il,dl,ib,b2; double q,r,b2r; l=nn[0]; il=nn[1]; dl=nn[2]; b2=b; ib=l-il-dl+1; k=ib+exp0; for (i=ib; i=ib) { a[k]=(int)(q-bas*r); } /* if */ k++; } /* for j */ a[k]=(int)(r+0.5); k++; for (i=k; i<=l; i++) { a[i]=0; } /* for i */ } /* lcmul */ void lmul(a,b,c,nn,q) unsigned char *a; unsigned char *b; unsigned char *c; int *nn; double *q; { int i,k,l,il,dl,ib,n,n2; int bi,b1,c1; double isize,rr,xx; double br,bri; if (mulversion == 1) { l=nn[0]; il=nn[1]; dl=nn[2]; n=2*l; n2=2*n; isize=1.0/sqrt((double)(n)); ib=l-il-dl+1; for (i=0; i=1; j--) { i=(int)(fu); x1[j]=(unsigned char)(i); fu=256.0*(fu-(double)(i)); } /* for j */ do { lass(x2,x1,nn); lcmul(a,2,x1,nn,0); lmul(t1,b,x1,nn,q); lmul(t2,x1,t1,nn,q); sub(x1,a,t2,nn); sub(t1,x1,x2,nn); } while (lcompeps(t1,nn)==1); lass(a,x1,nn); a[l+1]=s; } /* linv */ void ldivv(a,b,c,x1,x2,t1,t2,nn,q) unsigned char *a; unsigned char *b; unsigned char *c; unsigned char *x1; unsigned char *x2; unsigned char *t1; unsigned char *t2; int *nn; double *q; { linv(a,c,x1,x2,t1,t2,nn,q); lmul(a,a,b,nn,q); } /* ldivv */ unsigned char ppp1(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); } /* ppp1 */ void fppp(t,outfile) unsigned long 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 */ void vprint(v,i1,in) unsigned char *v; int i1,in; { int i; for (i=i1; i<=in; i++) { printf("%9d\n",(int)(v[i])); } /* for i */ printf("##############################\n"); } /* vprint */ void vprint2(v,i1,in,ss) unsigned char *v; int i1,in,ss; { int i; for (i=i1; i<=in; i++) { printf("%9d\n",(int)(v[i])); if (i == ss) { printf("................\n"); } /* if */ } /* for i */ printf("##############################\n"); } /* vprint2 */ void vzero(v,i1,in) unsigned char *v; int i1,in; { int i; for (i=i1; i<=in; i++) { v[i]=0; } /* for i */ } /* vzero */ int main(argc,argv) int argc; char *argv[]; { FILE *outfile,*infile; int ll,ll2,ll4,ql; int i,j,k,pot,n,n8; int n2,m; int ii,jj,kk,i2; int d1,d0; int s1,s0,i3,i4,iii,jjj; unsigned char cc,h1,h0,cc0; int gg[5]; int nn[31]; clock_t time1,timediff1,time2,timediff2; unsigned char *v1,*v2,*v3,*v4,*v5,*v6; unsigned long ttt; unsigned char *mmm[31]; unsigned char **ddd[3]; unsigned char *mpddd[3]; unsigned char *mpddd2[3]; double *q; ql=4*mulversion; 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/2; n8=n2/8; m=(int)(2.408*(double)(n8))+1; rewind(infile); pot=(int)(log(8.0*(double)(m))/log(2.0)); printf("pot=%3d; hex=%9d; dec=%9d\n",pot,n,8*m); ll=1; for (i=1; i<=pot; i++) { ll=2*ll; } /* for i */ ll2=2*ll; ll4=2*ll; printf(" "); printf("ll=%9d\n",ll2); jj=ll; ii=1; for (i=1; i<=pot+2; i++) { nn[i]=jj; mmm[i]=(unsigned char *)calloc(jj+1,sizeof(unsigned char)); if ((i >= 1) && (i <= 2)) { ddd[i]=(unsigned char **)calloc(ii+1,sizeof(unsigned char*)); mpddd[i]=(unsigned char *)calloc((ii+1)*(jj+1),sizeof(unsigned char)); mpddd2[i]=mpddd[i]; for (j=0; j<=ii; j++) { ddd[i][j]=mpddd[i]; mpddd[i]=mpddd[i]+jj+1; } /* for j */ } /* if */ jj=jj/2; ii=2*ii; } /* for i */ v1 = (unsigned char *)calloc(ll4+2,sizeof(unsigned char)); v2 = (unsigned char *)calloc(ll4+2,sizeof(unsigned char)); v3 = (unsigned char *)calloc(ll4+2,sizeof(unsigned char)); v4 = (unsigned char *)calloc(ll4+2,sizeof(unsigned char)); v5 = (unsigned char *)calloc(ll4+2,sizeof(unsigned char)); v6 = (unsigned char *)calloc(ll4+2,sizeof(unsigned char)); q = (double *)calloc(ql*ll+1,sizeof(double)); d0=ll; d1=ll/2; s0=ll2; s1=ll2/2; cc0=getc(infile); cc=getc(infile); vzero(v3,1,s0); ii=s1-n2+1; for (i=s1; i>=ii; i--) { cc=getc(infile); h1=ppp1(cc); cc=getc(infile); h0=ppp1(cc); v3[i]=(h1 << 4)+h0; } /* for i */ time1=clock(); time2=clock(); ii=pot-1; gg[0]=nn[ii]; gg[1]=nn[ii]; gg[2]=0; linttol(v1,10,0,gg); for (i=1; i<=2; i++) { lmul(v1,v1,v1,gg,q); } /* for i */ jj=pot-1; for (i=3; i<=pot+1; i++) { gg[0]=nn[jj]; gg[1]=nn[jj]; gg[2]=0; lmul(v1,v1,v1,gg,q); shift(mmm[jj],v1,1,nn[jj],1); jj=jj-1; } /* for i */ free(q); vzero(v2,1,s0); shift3(v2,v1,1,d0,s1+1); gg[0]=ll2; gg[1]=ll2/2; gg[2]=ll2/2; q = (double *)calloc(ql*ll2+1,sizeof(double)); lmul(v2,v3,v2,gg,q); free(q); shift(ddd[1][1],v2,1,d0,s1+1); timediff2=clock()-time2; printf("\n"); printf("Init time= %9.2f\n",(double)(timediff2)/timetic); printf("\n"); d0=ll; d1=ll/2; s0=ll2; s1=ll2/2; ii=1; i3=1; jjj=ll/2; for (i=1; i<=pot-2; i++) { time2=clock(); free(mmm[i]); free(v1); free(v2); free(v3); free(v4); free(v5); free(v6); v1 = (unsigned char *)calloc(s0+2,sizeof(unsigned char)); v2 = (unsigned char *)calloc(s0+2,sizeof(unsigned char)); v3 = (unsigned char *)calloc(s0+2,sizeof(unsigned char)); v4 = (unsigned char *)calloc(s0+2,sizeof(unsigned char)); v5 = (unsigned char *)calloc(s0+2,sizeof(unsigned char)); v6 = (unsigned char *)calloc(s0+2,sizeof(unsigned char)); i2=i+1; i4=3-i3; j=d1; while (mmm[i2][j] == 0) { j--; } /* while */ kk=j; printf("pot=%3d; ii=%9d; jjj2=%9d;",i,ii,d0); vzero(v2,1,s0); vzero(v1,1,s0); vzero(v4,1,s0); vzero(v5,1,s0); vzero(v6,1,s0); shift2(v2,mmm[i2],s1-kk-1-2,s1-2,j); gg[0]=s1; gg[1]=2; gg[2]=s1-2; gg[3]=s1-2-10; q = (double *)calloc(ql*s1+1,sizeof(double)); linv(v3,v2,v1,v4,v5,v6,gg,q); free(q); for (k=1; k<=ii; k++) { jj=2*k-1; vzero(v4,1,s0); shift3(v4,ddd[i3][k],1,d0,s1+1); vzero(v5,1,s0); shift2(v5,v3,1+1,s1+1,s1-1); vzero(v1,1,s0); gg[0]=s0; gg[1]=s0/2; gg[2]=s0/2; q = (double *)calloc(ql*s0+1,sizeof(double)); lmul(v1,v4,v5,gg,q); free(q); shift(ddd[i4][jj+1],v1,1,d1,s1+1+kk); vzero(v4,1,s0); shift(v4,v1,1,1+d1+1,s1+1+kk); vzero(v6,1,s0); shift(v6,mmm[i2],1,1+d1+1,1); vzero(v2,1,s0); gg[0]=s1; gg[1]=s1; gg[2]=0; q = (double *)calloc(ql*s1+1,sizeof(double)); lmul(v2,v4,v6,gg,q); free(q); vzero(v5,1,s0); sub2(v5,ddd[i3][k],v2,1,d0+1,1,1); shift(ddd[i4][jj],v5,1,d1+1,1); } /* for k */ d0=d0/2; d1=d1/2; s0=s0/2; s1=s1/2; if (s0 < 64) { s0=64; s1=s0/2; } /* if */ free(ddd[i3]); free(mpddd2[i3]); ii=2*ii; iii=2*ii; jjj=jjj/2; ddd[i3]=(unsigned char **)calloc(iii+1,sizeof(unsigned char*)); mpddd[i3]=(unsigned char *)calloc((iii+1)*(jjj+1),sizeof(unsigned char)); mpddd2[i3]=mpddd[i3]; for (j=0; j<=iii; j++) { ddd[i3][j]=mpddd[i3]; mpddd[i3]=mpddd[i3]+jjj+1; } /* for j */ if (i3 == 1) { i3=2; } else { i3=1; } /* if */ timediff2=clock()-time2; printf(" time= %9.2f\n",(double)(timediff2)/timetic); } /* for i */ timediff1=clock()-time1; printf("\n"); printf("Time total=%10.2f\n",(double)(timediff1)/timetic); printf("\n"); outfile=fopen(argv[2],"w"); putc(cc0,outfile); fprintf(outfile,".\n"); kk=i4; k=1; for (i=ii; i>=(ii-m); i--) { jj=1; ttt=0; for (j=0; j<=3; j++) { ttt=ttt+((unsigned long)(ddd[kk][i][jj]) << (j*8)); jj++; } /* for j */ fppp(ttt,outfile); if (k == 9) { fprintf(outfile,"\n"); k=0; } /* if */ k++; } /* for i */ fprintf(outfile,"\n"); fclose(outfile); fclose(infile); } /* End */