/***********************************************************************/ /* */ /* */ /* Program: */ /* -------- */ /* */ /* sseln2.c */ /* */ /* */ /* Purpose: */ /* -------- */ /* */ /* This program calculates exp(1) and ln(2) to many decimal digits. */ /* */ /* */ /* Version 1.03: */ /* ------------- */ /* */ /* First version: July 2003 */ /* Latest update: 2007-08-13 */ /* */ /* */ /* 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. Some of the algorithms are found on the web-page "Mathematical */ /* constants and computation" by Xavier Gourdon and Pascal Sebah: */ /* */ /* http://numbers.computation.free.fr/Constants/constants.html */ /* */ /* 2. The web-page "Fast Algorithms and the FEE Method", */ /* by Ekatherina A. Karatsuba: */ /* */ /* http://www.ccas.ru/personal/karatsuba/algen.htm */ /* */ /* 3. The web-page "The Karatsuba Method 'Divide and Conquer'", */ /* by Ekatherina A. Karatsuba: */ /* */ /* http://www.ccas.ru/personal/karatsuba/divcen.htm */ /* */ /* 4. The article "Fast evaluation of transcendental functions", */ /* Problems of Information Transmission, vol. 27, (1991), p. 339-360, */ /* by Ekatherina A. Karatsuba. */ /* */ /* 5. The article "Fast multiprecision evaluation of series of */ /* rational numbers", by Bruno Haible and Thomas Papanikolaou, 1997. */ /* */ /* 6. The program now use a very fast FFT (fftsg_h.c) by Takuya Ooura: */ /* */ /* http://momonga.t.u-tokyo.ac.jp/~ooura/ */ /* */ /* 7. "CLN - Class Library for Numbers" by Bruno Haible et. al: */ /* */ /* http://www.ginac.de/CLN/ */ /* */ /* 8. The documentation "Programs to Calculate some Mathematical */ /* Constants to Large Precision", by Stefan Spaennare. */ /* */ /* */ /* Program history: */ /* ---------------- */ /* */ /* 2007-08-13: */ /* The code has been somewhat prettified and shorter (in main). */ /* The program has the same speed as before. */ /* */ /* 2003-10-09: */ /* Very minor changes. */ /* */ /* 2003-08-27: */ /* Three more series for ln(2) added (method 5 - 7). */ /* Minor improvements to some procedures. Updated benchmarks. */ /* */ /* First version July 2003. */ /* */ /* */ /* Copy-right: */ /* ----------- */ /* */ /* This program is free to copy if this information follows the */ /* program. You are also allowed to modify the program for your */ /* own calculations if you give a reference to the author. */ /* See also notice below! */ /* */ /* */ /* IMPORTANT! */ /* ---------- */ /* */ /* Please don't use the result from this program for important */ /* purposes without doing an independent check of the result! */ /* See also WARNING and notice below! */ /* */ /* */ /* Algorithms: */ /* ----------- */ /* */ /* The program uses (see also the references above): */ /* */ /* * Fast series evaluation for exp(1) and ln(2). */ /* */ /* * FFT for multiplications. */ /* */ /* * 4-th order methods for inversions and divisions. */ /* */ /* * 4-th order methods for square and cube roots. */ /* */ /* */ /* Method = 1: */ /* ----------- */ /* */ /* inf */ /* exp(1); e = SUM 1/k! */ /* k=0 */ /* */ /* Method = 2: */ /* ----------- */ /* */ /* inf */ /* exp(1); 1/e = SUM (-1)^k/k! */ /* k=0 */ /* */ /* ********/ /* Method = 3: */ /* ----------- */ /* */ /* ln(2) = 144*atanh(1/251)+54*atanh(1/449)-38*atanh(1/4801)+62*atanh(1/8749) */ /* */ /* Each iteration adds about 4.8 digits to ln(2). */ /* */ /* Method = 4: */ /* ----------- */ /* */ /* ln(2) = 72*atanh(1/127)+54*atanh(1/449)+34*atanh(1/4801)-10*atanh(1/8749) */ /* */ /* Each iteration adds about 4.2 digits to ln(2). */ /* */ /* Method = 5: */ /* ----------- */ /* */ /* ln(2) = 18*atanh(1/26)-2*atanh(1/4801)+8*atanh(1/8749) */ /* */ /* Each iteration adds about 2.8 digits to ln(2). */ /* */ /* Method = 6: */ /* ----------- */ /* */ /* ln(2) = 14*atanh(1/31)+10*atanh(1/49)+6*atanh(1/161) */ /* */ /* Each iteration adds about 3.0 digits to ln(2). */ /* This series was found by a "brute force" computer search by the author. */ /* */ /* Method = 7: */ /* ----------- */ /* */ /* ln(2) = 6*atanh(1/26)+8*atanh(1/31)+10*atanh(1/49) */ /* */ /* Each iteration adds about 2.8 digits to ln(2). */ /* This series was found by a "brute force" computer search by the author. */ /* */ /* ********/ /* atanh(x) is defined by: */ /* ----------------------- */ /* */ /* 1 1+x inf x^(2k+1) */ /* atanh(x) = arctanh(x) = - log(---) = SUM -------- */ /* 2 1-x k=0 2k+1 */ /* */ /* The series above are truncated to the desired precision. */ /* */ /* */ /* Accuracy: */ /* --------- */ /* */ /* The program calculates about 0.995*2^(power of 2) decimal digits. */ /* */ /* All printed digits are probably correct. The results have */ /* been verified to 0.995*2^20 digits (for mulversion = 1). */ /* */ /* */ /* FFT max error */ /* ------------- */ /* */ /* Set the variable mulversion=2 if "FFT max error" > 0.25. */ /* This makes the program two times slower and requires more memory. */ /* */ /* The variable mulversion is present in the file "mulver.txt" */ /* containing the single value 1 or 2. Default value is 1. */ /* */ /* */ /* WARNING! */ /* -------- */ /* */ /* If you want to make calculations with much more than 2^24 */ /* (16 million) 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 100 */ /* instead of 10000. This means that the convolution will be */ /* saturated for a much larger number of digits (> 100 millions). */ /* */ /* */ /* Compilation: */ /* ------------ */ /* */ /* The program works in Unix, Linux and Windows (DOS using DJGPP). */ /* It is recommended that you optimize the program (i.e. use gcc -O3 */ /* option or similar option) when you compile: */ /* */ /*>gcc -O3 -o sseln2 sseln2.c fftsg_h.c -lm */ /* */ /* The time wraps around after 4295 seconds on some computers. */ /* */ /* */ /* Output: */ /* ------- */ /* */ /* Output is to the files "e.txt" and "ln2.txt". */ /* These file names can easily be changed. */ /* */ /* */ /* 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 9 */ /* Linux with the gcc 3.22 C-compiler. The CPU load was on average */ /* 1.00 (i.e. no other programs running). */ /* */ /* The program was compiled using the line: */ /* */ /*>gcc -O3 -o sseln2 sseln2.c fftsg_h.c -lm */ /* */ /* The times have been measured for "power of 2" equals 17 and 20 with */ /* the "time" function in Linux. TC below means "Time Complexity". */ /* The variable mulversion = 1. */ /* */ /* Intel Celeron 1400 MHz */ /* ---------------------- */ /* */ /* Algorithm 2^17 2^20 TC Memory/(2^20) */ /*---------------------------------------------------------------------*/ /* exp(1), method 1 1.99 s 21.04 s O(n*log(n)^3) 12.8 Mbyte */ /* exp(1), method 2 2.38 s 26.45 s O(n*log(n)^3) 12.8 Mbyte */ /* --------------------------------------------------------------------*/ /* ln(2), method 3 19.17 s 253.59 s O(n*log(n)^3) 16.0 Mbyte */ /* ln(2), method 4 21.44 s 276.08 s O(n*log(n)^3) 16.0 Mbyte */ /* ln(2), method 5 23.51 s 309.69 s O(n*log(n)^3) 16.0 Mbyte */ /* ln(2), method 6 19.54 s 244.00 s O(n*log(n)^3) 16.0 Mbyte */ /* ln(2), method 7 19.61 s 259.85 s O(n*log(n)^3) 16.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); double pi=3.14159265358979323; char version_string[80]="Program sseln2.c, version 1.03, 2007-08-13"; char file_name1[20]="e.txt"; /* Output filename */ char file_name2[20]="ln2.txt"; /* Output filename */ char mv_file_name[20]="mulver.txt"; /* Multiplication version file name */ int mulversion; /* Multiplication version (1 or 2) */ double bas=10000.0; double basi=0.0001; double fftmaxerror; typedef struct biginttag { short *coef; int size,maxsize,exp0,sgn0; } bigint; bigint tmp,tmp2; 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]+d[nn3-j-1]); y2=0.5*(d[j]-d[nn3-j-1]); d[j-1]=x1*y1-y2*(-x2); d[j]=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) short *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 */ a[l+1]=1; } /* lzero */ void linttol(a,man0,exp0,nn) short *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); if (man0 >= 0) { a[l+1]=1; } else { a[l+1]=-1; } /* if */ } /* linttol */ void linttol2(a,man0,exp0,nn) short *a; int man0,exp0; int *nn; { int i,l,il,dl,ib; int tt1,tt2,t1,t2,b; 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 */ b=abs(man0); if (b < 10000) { a[l-il+exp0+1]=b; } else { tt1=b % 100000000; tt2=b / 100000000; t1=tt1 % 10000; t2=tt1 / 10000; a[l-il+exp0+1]=t1; a[l-il+exp0+2]=t2; a[l-il+exp0+3]=tt2; } /* if */ if (man0 >= 0) { a[l+1]=1; } else { a[l+1]=-1; } /* if */ } /* linttol2 */ void ppp(t) int t; { if (t > 999) { printf("%4d",t); } else if (t > 99) { printf("0%3d",t); } else if (t > 9) { printf("00%2d",t); } else { printf("000%1d",t); } /* if */ } /* ppp */ void lprint(a,nn) short *a; int *nn; { int i,j,m,l,il,dl,t; l=nn[0]; il=nn[1]; dl=nn[3]; if (a[l+1]==-1) { printf("-"); } /* if */ i=l; while ((a[i]==0) && (i>(l-il+1))) { i--; } /* while */ m=0; if (i==(l-il+1)) { printf("%4d.",a[i]); m++; for (j=(l-il); j>=(l-il-dl+1); j--) { t=a[j]; ppp(t); m++; if ((m % 18)==0) { printf("\n"); } /* if */ } /* for j */ } else { printf("%d",a[i]); m++; for (j=(i-1); j>=(l-il+1); j--) { t=a[j]; ppp(t); m++; if ((m % 18)==0) { printf("\n"); } /* if */ } /* for j */ printf("."); for (j=(l-il); j>=(l-il-dl+1); j--) { t=a[j]; ppp(t); m++; if ((m % 18)==0) { printf("\n"); } /* if */ } /* for j */ } /* else */ printf("\n\n"); } /* lprint */ void fppp(t,outfile) int t; FILE *outfile; { if (t > 999) { fprintf(outfile,"%4d",t); } else if (t > 99) { fprintf(outfile,"0%3d",t); } else if (t > 9) { fprintf(outfile,"00%2d",t); } else { fprintf(outfile,"000%1d",t); } /* if */ } /* fppp */ void lfprint(a,nn,outfile) short *a; int *nn; FILE *outfile; { int i,j,m,l,il,dl,t; l=nn[0]; il=nn[1]; dl=nn[3]; if (a[l+1]==-1) { fprintf(outfile,"-"); } /* if */ i=l; while ((a[i]==0) && (i>(l-il+1))) { i--; } /* while */ m=0; if (i==(l-il+1)) { fprintf(outfile,"%4d.",a[i]); m++; for (j=(l-il); j>=(l-il-dl+1); j--) { t=a[j]; fppp(t,outfile); m++; if ((m % 18)==0) { fprintf(outfile,"\n"); } /* if */ } /* for j */ } else { fprintf(outfile,"%d",a[i]); m++; for (j=(i-1); j>=(l-il+1); j--) { t=a[j]; fppp(t,outfile); m++; if ((m % 18)==0) { fprintf(outfile,"\n"); } /* if */ } /* for j */ fprintf(outfile,"."); for (j=(l-il); j>=(l-il-dl+1); j--) { t=a[j]; fppp(t,outfile); m++; if ((m % 18)==0) { fprintf(outfile,"\n"); } /* if */ } /* for j */ } /* else */ fprintf(outfile,"\n"); } /* lfprint */ void add(a,b,c,nn) short *a; short *b; short *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=10000; r=0; for (i=ib; i<=l; i++) { q=b[i]+c[i]+r; if (q >= basint) { a[i]=q-basint; r=1; } else { a[i]=q; r=0; } /* if */ } /* for i */ } /* add */ void sub(a,b,c,nn) short *a; short *b; short *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=10000; r=0; for (i=ib; i<=l; i++) { q=b[i]-c[i]-r; if (q < 0) { a[i]=basint+q; r=1; } else { a[i]=q; r=0; } /* if */ } /* for i */ } /* sub */ int labscomp(a,b,nn) short *a; short *b; int *nn; { int i,l,il,dl,ib; l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; i=l; while ((a[i]==b[i]) && (i>=ib)) { i--; } /* while */ if (a[i]>=b[i]) { return(1); } else { return(-1); } /* if */ } /* lasbcomp */ void ladd(a,b,c,nn) short *a; short *b; short *c; int *nn; { int b2,c2,co,l; l=nn[0]; co=labscomp(b,c,nn); b2=b[l+1]; c2=c[l+1]; if ((b2>0) && (c2>0)) { add(a,b,c,nn); a[l+1]=1; } /* if */ if ((b2>0) && (c2<0) && (co==1)) { sub(a,b,c,nn); a[l+1]=1; } /* if */ if ((b2>0) && (c2<0) && (co==-1)) { sub(a,c,b,nn); a[l+1]=-1; } /* if */ if ((b2<0) && (c2>0) && (co==1)) { sub(a,b,c,nn); a[l+1]=-1; } /* if */ if ((b2<0) && (c2>0) && (co==-1)) { sub(a,c,b,nn); a[l+1]=1; } /* if */ if ((b2<0) && (c2<0)) { add(a,b,c,nn); a[l+1]=-1; } /* if */ } /* ladd */ void lsub(a,b,c,nn) short *a; short *b; short *c; int *nn; { int b2,c2,co,l; l=nn[0]; co=labscomp(b,c,nn); b2=b[l+1]; c2=c[l+1]; if ((b2>0) && (c2>0) && (co==1)) { sub(a,b,c,nn); a[l+1]=1; } /* if */ if ((b2>0) && (c2>0) && (co==-1)) { sub(a,c,b,nn); a[l+1]=-1; } /* if */ if ((b2>0) && (c2<0)) { add(a,b,c,nn); a[l+1]=1; } /* if */ if ((b2<0) && (c2>0)) { add(a,b,c,nn); a[l+1]=-1; } /* if */ if ((b2<0) && (c2<0) && (co==1)) { sub(a,c,b,nn); a[l+1]=-1; } /* if */ if ((b2<0) && (c2<0) && (co==-1)) { sub(a,b,c,nn); a[l+1]=1; } /* if */ } /* lsub */ void lcmul(a,b,c,exp0,nn) short *a; short *b; int c; int exp0; int *nn; { int i,j,k,l,il,dl,ib; double q,r,cr; l=nn[0]; il=nn[1]; dl=nn[2]; 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 */ if (c>=0) { a[l+1]=b[l+1]; } else { a[l+1]=-b[l+1]; } /* if */ } /* lcmul */ void lcmul2(a,b,c,exp0,nn) short *a; short *b; double c; int exp0; int *nn; { int i,j,k,l,il,dl,ib; double q,r,cr; l=nn[0]; il=nn[1]; dl=nn[2]; 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 */ if (c>=0.0) { a[l+1]=b[l+1]; } else { a[l+1]=-b[l+1]; } /* if */ } /* lcmul2 */ void ltoreal(a,man0,exp0,nn) short *a; double *man0,*exp0; int *nn; { char s0[20]; char s1[10]; char s2[10]; char s3[10]; char s4[10]; char st2[10]; char st3[10]; char st4[10]; int i,l,il,dl,ib,length,le; double man,ex; l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; i=l; while ((a[i]==0) && (i>=ib)) { i--; } /* while */ if ((i==ib) && (a[i]==0)) { man=0.0; ex=0.0; } else { sprintf(s1,"%d",a[i]); length=strlen(s1); ex=(double)(4*(i-(l-il+1))+length); sprintf(s2,"%d",a[i-1]); sprintf(s3,"%d",a[i-2]); sprintf(s4,"%d",a[i-3]); strcpy(st2,""); le=strlen(s2); for (i=1; i<=(4-le); i++) { strcat(st2,"0"); } /* for i */ strcat(st2,s2); strcpy(st3,""); le=strlen(s3); for (i=1; i<=(4-le); i++) { strcat(st3,"0"); } /* for i */ strcat(st3,s3); strcpy(st4,""); le=strlen(s4); for (i=1; i<=(4-le);i++) { strcat(st4,"0"); } /* for i */ strcat(st4,s4); sprintf(s0,"%s%s%s%s%s","0.",s1,st2,st3,st4); man=atof(s0); } /* if */ *man0=(double)(a[l+1])*man; *exp0=ex; } /* ltoreal */ void lrealtol(a,man0,exp0,nn) short *a; double man0,exp0; int *nn; { char s0[18]; char s1[16]; char z1[7]; char z2[7]; char z3[7]; char z4[7]; double man,ex,exp2,man2; int i,k,l,il,dl,ib,l1,l2,ll; 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 */ man2=man0; ex=exp0; man=fabs(man2); while (man>=1.0) { man=man/10.0; ex=ex+1.0; } /* while */ sprintf(s0,"%17.16f",man); for (i=2; i<=17; i++) { s1[i-2]=s0[i]; } /* for i */ exp2=fabs(ex/(double)(4)); l1=(int)(ex/(double)(4)); l2=(int)(4*(exp2-(double)(abs(l1)))); if (ex<0) { l2=4-l2; l1=l1-1; } /* if */ if (l2==0) { l2=4; l1=l1-1; } /* if */ for (i=0; i<=l2-1; i++) { z1[i]=s1[i]; } /* for i */ z1[i]='\0'; for (i=0; i<=4-1; i++) { z2[i]=s1[l2+i]; } /* for i */ z2[i]='\0'; k=0; for (i=0; i<=4-1; i++) { if ((l2+i+4)<=16) { z3[i]=s1[l2+i+4]; k++; } /* if */ } /* for i */ z3[k]='\0'; for (i=1; i<=4-k; i++) { strcat(z3,"0"); } /* for i */ k=0; for (i=0; i<=4-1; i++) { if ((l2+i+2*4)<=16) { z4[i]=s1[l2+i+2*4]; k++; } /* if */ } /* for i */ z4[k]='\0'; for (i=1; i<=4-k; i++) { strcat(z4,"0"); } /* for i */ ll=l-il+l1+1; a[ll]=atoi(z1); a[ll-1]=atoi(z2); a[ll-2]=atoi(z3); a[ll-3]=atoi(z4); if (man2>=0) { a[l+1]=1; } else { a[l+1]=-1; } /* if */ } /* lrealtol */ int lcompeps(a,nn) short *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 */ int lgetdec(a,nn) short *a; int *nn; { int i,l,il,dl,k,ib; l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; k=0; i=l-il; while ((a[i]==0) && (i>ib)) { k++; i--; } /* while */ return(4*k); } /* lgetdec */ int lcheckdec(a,nn) short *a; int *nn; { int i,l,il,dl,k,ib,ss; l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; k=0; i=l-il; ss=a[i]; while ((a[i]==ss) && (i>ib)) { k++; i--; } /* while */ if ((i==ib) && (ss==0)) { k=0; } /* if */ return(k); } /* lcheckdec */ void lass(a,b,nn) short *a; short *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 lmul(a,b,c,nn) short *a; short *b; short *c; int *nn; { int i,k,l,il,dl,ib,n,n2; int bi,b1,c1; int bn,cn,kkk,bnib,cnib; double isize,rr,xx; double br,bri,maxerror,temp,qqq; double *q; if (mulversion == 1) { l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; i=l; while ((b[i] == 0) && (i > ib)) { i--; } /* while */ bnib=i; bn=bnib-ib; i=l; while ((c[i] == 0) && (i > ib)) { i--; } /* while */ cnib=i; cn=cnib-ib; kkk=bn+cn; if (kkk > dl) { n=1; while (n < kkk) { n=2*n; } /* while */ n2=2*n; q=(double *)malloc((n2+1)*sizeof(double)); isize=1.0/sqrt((double)(n)); for (i=0; i<=n2; i++) { q[i]=0.0; } /* for i */ k=0; for (i=ib; i<=bnib; i++) { q[k]=(double)(b[i])*isize; k=k+2; } /* for i */ k=0; for (i=ib; i<=cnib; i++) { q[k+1]=(double)(c[i])*isize; k=k+2; } /* for i */ cdft(n2,1,q); mul(q,n); cdft(n2,-1,q); maxerror=0.0; rr=0.0; for (i=0; i<=n2; i=i+2) { qqq=q[i]; temp=fabs(qqq-floor(qqq+0.5)); if (temp > maxerror) { maxerror=temp; } /* if */ xx=qqq+rr; rr=floor((xx+0.5)*basi); q[i]=xx-bas*rr; } /* for i */ if (maxerror > fftmaxerror) { fftmaxerror=maxerror; } /* if */ k=2*dl; for (i=ib; i<=l; i++) { a[i]=0; if (k <= n2) { a[i]=(int)(q[k]+0.5); } /* if */ k=k+2; } /* for i */ free(q); } else { for (i=ib; i<=l; i++) { a[i]=0; } /* for i */ } /* if */ a[l+1]=b[l+1]*c[l+1]; } else { l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; i=l; while ((b[i] == 0) && (i > ib)) { i--; } /* while */ bnib=i; bn=bnib-ib; i=l; while ((c[i] == 0) && (i > ib)) { i--; } /* while */ cnib=i; cn=cnib-ib; kkk=bn+cn; if (kkk > dl) { n=1; while (n < kkk) { n=2*n; } /* while */ n=2*n; n2=2*n; q=(double *)malloc((n2+1)*sizeof(double)); br=100.0; bri=0.01; bi=100; isize=1.0/sqrt((double)(n)); for (i=0; i<=n2; i++) { q[i]=0.0; } /* for i */ k=0; for (i=ib; i<=bnib; i++) { b1=b[i] % bi; q[k]=(double)(b1)*isize; b1=b[i] / bi; q[k+2]=(double)(b1)*isize; k=k+4; } /* for i */ k=0; for (i=ib; i<=cnib; i++) { c1=c[i] % bi; q[k+1]=(double)(c1)*isize; c1=c[i] / bi; q[k+3]=(double)(c1)*isize; k=k+4; } /* for i */ cdft(n2,1,q); mul(q,n); cdft(n2,-1,q); maxerror=0.0; rr=0.0; for (i=0; i<=n2; i=i+2) { qqq=q[i]; temp=fabs(qqq-floor(qqq+0.5)); if (temp > maxerror) { maxerror=temp; } /* if */ xx=qqq+rr; rr=floor((xx+0.5)*bri); q[i]=xx-br*rr; } /* for i */ if (maxerror > fftmaxerror) { fftmaxerror=maxerror; } /* if */ k=4*dl; for (i=ib; i<=l; i++) { a[i]=0; if (k <= n2) { a[i]=bi*(int)(q[k+2]+0.5)+(int)(q[k]+0.5); } /* if */ k=k+4; } /* for i */ free(q); } else { for (i=ib; i<=l; i++) { a[i]=0; } /* for i */ } /* if */ a[l+1]=b[l+1]*c[l+1]; } /* if */ } /* lmul */ void linv(a,b,nn) short *a; short *b; int *nn; { double man,ex; int s,l,il,dl,dd,kk,bb,bbb,i; short *x1,*x2,*t1,*t2; int ss[4]; l=nn[0]; il=nn[1]; dl=nn[2]; dd=il+dl; x1=(short *)malloc((l+2)*sizeof(short)); x2=(short *)malloc((l+2)*sizeof(short)); t1=(short *)malloc((l+2)*sizeof(short)); t2=(short *)malloc((l+2)*sizeof(short)); lzero(x1,nn); lzero(x2,nn); lzero(t1,nn); lzero(t2,nn); s=b[l+1]; b[l+1]=1; ltoreal(b,&man,&ex,nn); man=1.0/man; ex=-ex; lrealtol(x1,man,ex,nn); bbb=lcheckdec(b,nn); if (bbb < 16) { bb=32; } else { kk=(int)(log((double)(bbb+2))/log(2.0)); bb=2; for (i=1; i<=kk; i++) { bb=2*bb; } /* for i */ } /* if */ kk=bb; ss[0]=l; ss[1]=il; do { if (kk > dd) { kk=dd; } /* if */ ss[2]=kk-il; lmul(t1,x1,b,ss); linttol(t2,1,0,ss); lsub(x2,t2,t1,ss); lmul(t1,x2,x2,ss); lmul(a,x2,t1,ss); ladd(t2,t1,a,ss); ladd(t1,x2,t2,ss); lmul(t2,t1,x1,ss); ladd(t1,x1,t2,ss); lass(x1,t1,ss); kk=4*kk; } while (lcompeps(x2,nn)==1); lass(a,x1,nn); a[l+1]=s; free(x1); free(x2); free(t1); free(t2); } /* linv */ void ldivv(a,b,c,nn) short *a; short *b; short *c; int *nn; { linv(a,c,nn); lmul(a,a,b,nn); } /* ldivv */ void lsqrtinv(a,b,nn) short *a; short *b; int *nn; { double re0,re1,re2,s10,man,ex,eps; int s,l,il,dl,dd,kk,bb,bbb,i; short *x1,*x2,*t1,*t2; int ss[4]; l=nn[0]; il=nn[1]; dl=nn[2]; dd=il+dl; x1=(short *)malloc((l+2)*sizeof(short)); x2=(short *)malloc((l+2)*sizeof(short)); t1=(short *)malloc((l+2)*sizeof(short)); t2=(short *)malloc((l+2)*sizeof(short)); lzero(x1,nn); lzero(x2,nn); lzero(t1,nn); lzero(t2,nn); eps=0.001; s=b[l+1]; b[l+1]=1; ltoreal(b,&man,&ex,nn); re0=1.0/sqrt(man); re1=-ex/2.0; re2=fabs(re1); s10=sqrt(10.0); if ((re2-(double)((int)(re2)))>eps) { re0=re0/s10; re1=re1+0.5; } /* if */ man=re0; ex=re1; lrealtol(x1,man,ex,nn); bbb=lcheckdec(b,nn); if (bbb < 16) { bb=32; } else { kk=(int)(log((double)(bbb+2))/log(2.0)); bb=2; for (i=1; i<=kk; i++) { bb=2*bb; } /* for i */ } /* if */ kk=bb; ss[0]=l; ss[1]=il; do { if (kk > dd) { kk=dd; } /* if */ ss[2]=kk-il; lmul(t1,x1,b,ss); lmul(t1,x1,t1,ss); linttol(t2,1,0,ss); lsub(x2,t2,t1,ss); lmul(t1,x2,x2,ss); lmul(a,t1,x2,ss); lcmul(t2,a,5,0,ss); lcmul(a,t1,6,0,ss); ladd(t1,a,t2,ss); lcmul(t2,x2,8,0,ss); ladd(a,t1,t2,ss); lmul(t1,a,x1,ss); lcmul(a,t1,625,-1,ss); ladd(t2,x1,a,ss); lass(x1,t2,ss); kk=4*kk; } while (lcompeps(x2,nn)==1); lass(a,x1,nn); free(x1); free(x2); free(t1); free(t2); } /* lsqrtinv */ void lsqrt(a,b,nn) short *a; short *b; int *nn; { lsqrtinv(a,b,nn); lmul(a,a,b,nn); } /* lsqrt */ void lsqrtsqrtinv(a,b,nn) short *a; short *b; int *nn; { double re0,re1,re2,s10,man,ex,eps; int s,l,il,dl,bb,dd,kk,bbb,i; short *x1,*x2,*t1,*t2; int ss[4]; l=nn[0]; il=nn[1]; dl=nn[2]; dd=il+dl; x1=(short *)malloc((l+2)*sizeof(short)); x2=(short *)malloc((l+2)*sizeof(short)); t1=(short *)malloc((l+2)*sizeof(short)); t2=(short *)malloc((l+2)*sizeof(short)); lzero(x1,nn); lzero(x2,nn); lzero(t1,nn); lzero(t2,nn); eps=0.001; s=b[l+1]; b[l+1]=1; ltoreal(b,&man,&ex,nn); re0=1.0/sqrt(sqrt(man)); re1=-ex/4.0; re2=fabs(re1); s10=sqrt(sqrt(10.0)); while ((re2-(double)((int)(re2)))>eps) { re0=re0/s10; re1=re1+0.25; re2=fabs(re1); } /* while */ man=re0; ex=re1; lrealtol(x1,man,ex,nn); bbb=lcheckdec(b,nn); if (bbb < 16) { bb=32; } else { kk=(int)(log((double)(bbb+2))/log(2.0)); bb=2; for (i=1; i<=kk; i++) { bb=2*bb; } /* for i */ } /* if */ kk=bb; ss[0]=l; ss[1]=il; do { if (kk > dd) { kk=dd; } /* if */ ss[2]=kk-il; lmul(t1,x1,x1,ss); lmul(t1,t1,t1,ss); lmul(t1,t1,b,ss); linttol(t2,1,0,ss); lsub(x2,t2,t1,ss); lmul(t1,x2,x2,ss); lmul(a,t1,x2,ss); lcmul(t2,a,15,0,ss); lcmul(a,t1,20,0,ss); ladd(t1,a,t2,ss); lcmul(t2,x2,32,0,ss); ladd(a,t1,t2,ss); lmul(t1,a,x1,ss); lcmul(a,t1,781250,-2,ss); ladd(t2,x1,a,ss); lass(x1,t2,ss); kk=4*kk; } while (lcompeps(x2,nn)==1); lass(a,x1,nn); free(x1); free(x2); free(t1); free(t2); } /* lsqrtsqrtinv */ void lsqrtsqrt(a,b,nn) short *a; short *b; int *nn; { int l,il,dl,dd; short *t1; l=nn[0]; il=nn[1]; dl=nn[2]; dd=il+dl; lsqrtsqrtinv(a,b,nn); t1=(short *)malloc((l+2)*sizeof(short)); lmul(t1,a,a,nn); lmul(t1,t1,a,nn); lmul(a,t1,b,nn); free(t1); } /* lsqrtsqrt */ void linv3(a,nn) short *a; int *nn; { int i,l,il,dl,ib,ll; l=nn[0]; il=nn[1]; dl=nn[2]; ll=l-il+1; ib=l-il-dl+1; for (i=l;i>=ll; i--) { a[i]=0; } /* for i */ for (i=ll-1;i>=ib; i--) { a[i]=3333; } /* for i */ a[l+1]=1; } /* linv3 */ void linv162(a,nn) short *a; int *nn; { int i,l,il,dl,ib,ll,kk; char cc[12]; char ss[12]; cc[0]='7'; cc[1]='2'; cc[2]='8'; cc[3]='3'; cc[4]='9'; cc[5]='5'; cc[6]='0'; cc[7]='6'; cc[8]='1'; cc[9]='\0'; l=nn[0]; il=nn[1]; dl=nn[2]; ll=l-il+1; ib=l-il-dl+1; for (i=l;i>=ll; i--) { a[i]=0; } /* for i */ a[ll-1]=61; kk=0; for (i=ll-2;i>=ib; i--) { ss[0]=cc[(kk % 9)]; ss[1]=cc[((kk+1) % 9)]; ss[2]=cc[((kk+2) % 9)]; ss[3]=cc[((kk+3) % 9)]; ss[4]='\0'; a[i]=atoi(ss); kk=kk+4; } /* for i */ a[l+1]=1; } /* linv162 */ void lcubertinv(a,b,nn) short *a; short *b; int *nn; { double re0,re1,re2,s10,man,ex,eps; int s,l,il,dl,bb,dd,kk,bbb,i; short *x1,*x2,*t1,*t2,*t3; int ss[4]; l=nn[0]; il=nn[1]; dl=nn[2]; dd=il+dl; x1=(short *)malloc((l+2)*sizeof(short)); x2=(short *)malloc((l+2)*sizeof(short)); t1=(short *)malloc((l+2)*sizeof(short)); t2=(short *)malloc((l+2)*sizeof(short)); t3=(short *)malloc((l+2)*sizeof(short)); lzero(x1,nn); lzero(x2,nn); lzero(t1,nn); lzero(t2,nn); linv162(t3,nn); eps=0.001; s=b[l+1]; b[l+1]=1; ltoreal(b,&man,&ex,nn); re0=1.0/exp(log(man)/3.0); re1=-ex/3.0; re2=fabs(re1); s10=exp(log(10.0)/3.0); while ((re2-(double)((int)(re2)))>eps) { re0=re0/s10; re1=re1+1.0/3.0; re2=fabs(re1); } /* while */ man=re0; ex=re1; lrealtol(x1,man,ex,nn); bbb=lcheckdec(b,nn); if (bbb < 16) { bb=32; } else { kk=(int)(log((double)(bbb+2))/log(2.0)); bb=2; for (i=1; i<=kk; i++) { bb=2*bb; } /* for i */ } /* if */ kk=bb; ss[0]=l; ss[1]=il; do { if (kk > dd) { kk=dd; } /* if */ ss[2]=kk-il; lmul(t1,x1,x1,ss); lmul(t1,t1,x1,ss); lmul(t1,t1,b,ss); linttol(t2,1,0,ss); lsub(x2,t2,t1,ss); lmul(t1,x2,x2,ss); lmul(a,t1,x2,ss); lcmul(t2,a,28,0,ss); lcmul(a,t1,36,0,ss); ladd(t1,a,t2,ss); lcmul(t2,x2,54,0,ss); ladd(a,t1,t2,ss); lmul(t1,a,x1,ss); lmul(a,t3,t1,ss); ladd(t2,x1,a,ss); lass(x1,t2,ss); kk=4*kk; } while (lcompeps(x2,nn)==1); lass(a,x1,nn); free(x1); free(x2); free(t1); free(t2); free(t3); } /* lcubertinv */ void lcubert(a,b,nn) short *a; short *b; int *nn; { int l,il,dl,dd; short *t1; l=nn[0]; il=nn[1]; dl=nn[2]; dd=il+dl; lcubertinv(a,b,nn); t1=(short *)malloc((l+2)*sizeof(short)); lmul(t1,a,a,nn); lmul(a,t1,b,nn); free(t1); } /* lcubert */ void biginit(a,maxsize) bigint *a; int maxsize; { a->coef=(short *)malloc((maxsize+1)*sizeof(short)); a->size=1; a->maxsize=maxsize; a->exp0=0; a->sgn0=1; a->coef[1]=0; } /* biginit */ void biginit2(a,maxsize) bigint *a; int maxsize; { int i; a->coef=(short *)malloc((maxsize+1)*sizeof(short)); a->size=1; a->maxsize=maxsize; a->exp0=0; a->sgn0=1; for (i=1; i<=maxsize; i++) { a->coef[i]=0; } /* for i */ } /* biginit2 */ void bigfree(a) bigint *a; { free(a->coef); } /* bigfree */ void inttobig(a,bb) bigint *a; int bb; { int tt1,tt2,t1,t2,b; b=abs(bb); if (b < 10000) { a->coef[1]=b; a->size=1; } else { tt1=b % 100000000; tt2=b / 100000000; t1=tt1 % 10000; t2=tt1 / 10000; a->coef[1]=t1; a->coef[2]=t2; a->coef[3]=tt2; a->size=1; if (t2 > 0) { a->size=2; } /* if */ if (tt2 > 0) { a->size=3; } /* if */ } /* if */ a->exp0=0; if (bb >= 0) { a->sgn0=1; } else { a->sgn0=-1; } /* if */ } /* inttobig */ void bigass(a,b) bigint *a; bigint *b; { int i,n; n=b->size; for (i=1; i<=n; i++) { a->coef[i]=b->coef[i]; } /* for i */ a->size=n; a->exp0=b->exp0; a->sgn0=b->sgn0; } /* bigass */ int bigcomp2(b,c) bigint *b; bigint *c; { int i,n,bn,cn,bb,cc,maxsize,maxs; int bexp,cexp,expmax; int bk,ck,bi,ci,flag; bn=b->size; cn=c->size; maxsize=b->maxsize; maxs=maxsize-2; bexp=b->exp0; cexp=c->exp0; if ((bn+bexp) > (cn+cexp)) { expmax=bexp; } else { expmax=cexp; } /* if */ if (bn > cn) { n=bn; } else { n=cn; } /* if */ bi=expmax-bexp; ci=expmax-cexp; i=n; flag=0; while ((flag == 0) && (i >= 1)) { bk=i+bi; if ((bk <= bn) && (bk >= 1)) { bb=b->coef[bk]; } else { bb=0; } /* if */ ck=i+ci; if ((ck <= cn) && (ck >= 1)) { cc=c->coef[ck]; } else { cc=0; } /* if */ if (bb != cc) { flag=1; } /* if */ i--; } /* while*/ if (bb >= cc) { return(1); } else { return(-1); } /* if */ } /* bigcomp2 */ void bigadd2(a,b,c) bigint *a; bigint *b; bigint *c; { int i,n,bn,cn,bb,cc,maxsize,maxs; int q,r,basint,bexp,cexp,expmax; int bk,ck,bi,ci; bn=b->size; cn=c->size; maxsize=b->maxsize; maxs=maxsize-2; bexp=b->exp0; cexp=c->exp0; if ((bn+bexp) > (cn+cexp)) { expmax=bexp; } else { expmax=cexp; } /* if */ if (bn > cn) { n=bn; } else { n=cn; } /* if */ bi=expmax-bexp; ci=expmax-cexp; basint=10000; r=0; for (i=1; i<=n; i++) { bk=i+bi; if ((bk <= bn) && (bk >= 1)) { bb=b->coef[bk]; } else { bb=0; } /* if */ ck=i+ci; if ((ck <= cn) && (ck >= 1)) { cc=c->coef[ck]; } else { cc=0; } /* if */ q=bb+cc+r; if (q >= basint) { a->coef[i]=q-basint; r=1; } else { a->coef[i]=q; r=0; } /* if */ } /* for i */ if (n >= maxs) { if (r > 0) { for (i=1; i < n; i++) { a->coef[i]=a->coef[i+1]; } /* for i */ a->coef[n]=r; a->size=n; a->exp0=expmax+1; } else { a->size=n; a->exp0=expmax; } /* if */ } else { if (r > 0) { a->coef[n+1]=r; a->size=n+1; a->exp0=expmax; } else { a->size=n; a->exp0=expmax; } /* if */ } /* if */ } /* bigadd2 */ void bigsub2(a,b,c) bigint *a; bigint *b; bigint *c; { int i,n,bn,cn,bb,cc,maxsize,maxs; int q,r,basint,bexp,cexp,expmax; int bk,ck,bi,ci; bn=b->size; cn=c->size; maxsize=b->maxsize; maxs=maxsize-2; bexp=b->exp0; cexp=c->exp0; if ((bn+bexp) > (cn+cexp)) { expmax=bexp; } else { expmax=cexp; } /* if */ if (bn > cn) { n=bn; } else { n=cn; } /* if */ bi=expmax-bexp; ci=expmax-cexp; basint=10000; r=0; for (i=1; i<=n; i++) { bk=i+bi; if ((bk <= bn) && (bk >= 1)) { bb=b->coef[bk]; } else { bb=0; } /* if */ ck=i+ci; if ((ck <= cn) && (ck >= 1)) { cc=c->coef[ck]; } else { cc=0; } /* if */ q=bb-cc-r; if (q < 0) { a->coef[i]=basint+q; r=1; } else { a->coef[i]=q; r=0; } /* if */ } /* for i */ i=n; while ((a->coef[i] == 0) && (i > 1)) { i--; } /* while */ a->size=i; a->exp0=expmax; } /* bigsub2 */ void bigadd(a,b,c) bigint *a; bigint *b; bigint *c; { int b2,c2,co; co=bigcomp2(b,c); b2=b->sgn0; c2=c->sgn0; if ((b2>0) && (c2>0)) { bigadd2(a,b,c); a->sgn0=1; } /* if */ if ((b2>0) && (c2<0) && (co==1)) { bigsub2(a,b,c); a->sgn0=1; } /* if */ if ((b2>0) && (c2<0) && (co==-1)) { bigsub2(a,c,b); a->sgn0=-1; } /* if */ if ((b2<0) && (c2>0) && (co==1)) { bigsub2(a,b,c); a->sgn0=-1; } /* if */ if ((b2<0) && (c2>0) && (co==-1)) { bigsub2(a,c,b); a->sgn0=1; } /* if */ if ((b2<0) && (c2<0)) { bigadd2(a,b,c); a->sgn0=-1; } /* if */ } /* bigadd */ void bigsub(a,b,c,nn) bigint *a; bigint *b; bigint *c; int *nn; { int b2,c2,co; co=bigcomp2(b,c); b2=b->sgn0; c2=c->sgn0; if ((b2>0) && (c2>0) && (co==1)) { bigsub2(a,b,c); a->sgn0=1; } /* if */ if ((b2>0) && (c2>0) && (co==-1)) { bigsub2(a,c,b); a->sgn0=-1; } /* if */ if ((b2>0) && (c2<0)) { bigadd2(a,b,c); a->sgn0=1; } /* if */ if ((b2<0) && (c2>0)) { bigadd2(a,b,c); a->sgn0=-1; } /* if */ if ((b2<0) && (c2<0) && (co==1)) { bigsub2(a,c,b); a->sgn0=-1; } /* if */ if ((b2<0) && (c2<0) && (co==-1)) { bigsub2(a,b,c); a->sgn0=1; } /* if */ } /* bigsub */ void bigmul(a,b,c) bigint *a; bigint *b; bigint *c; { double *q; int i,k,n,n2; int bi,b1,c1; int bn,cn,kkk,maxsize; int bexp,cexp,exptot,gg,maxs; double isize,rr,xx; double br,bri,maxerror,qqq,temp; if (mulversion == 1) { bn=b->size; cn=c->size; bexp=b->exp0; cexp=c->exp0; exptot=bexp+cexp; maxsize=b->maxsize; maxs=maxsize-2; kkk=bn+cn; n=1; while (n < kkk) { n=2*n; } /* while */ n2=2*n; q=(double *)malloc((n2+1)*sizeof(double)); for (i=0; i<=n2; i++) { q[i]=0.0; } /* for i */ isize=1.0/sqrt((double)(n)); k=0; for (i=1; i<=bn; i++) { q[k]=(double)(b->coef[i])*isize; k=k+2; } /* for i */ k=0; for (i=1; i<=cn; i++) { q[k+1]=(double)(c->coef[i])*isize; k=k+2; } /* for i */ cdft(n2,1,q); mul(q,n); cdft(n2,-1,q); maxerror=0.0; rr=0.0; for (i=0; i<=n2; i=i+2) { qqq=q[i]; temp=fabs(qqq-floor(qqq+0.5)); if (temp > maxerror) { maxerror=temp; } /* if */ xx=qqq+rr; rr=floor((xx+0.5)*basi); q[i]=xx-bas*rr; } /* for i */ if (maxerror > fftmaxerror) { fftmaxerror=maxerror; } /* if */ gg=1; if (kkk > maxs) { /* printf("* %8d %8d\n",kkk,maxs); */ gg=kkk-maxs+1; } /* if */ k=2*gg-2; for (i=gg; i<=kkk; i++) { a->coef[i-gg+1]=(int)(q[k]+0.5); k=k+2; } /* for i */ if (a->coef[kkk-gg+1] > 0) { a->size=kkk-gg+1; a->exp0=exptot+gg-1; } else { a->size=kkk-gg; a->exp0=exptot+gg-1; } /* if */ a->sgn0=(b->sgn0)*(c->sgn0); free(q); } else { br=100.0; bri=0.01; bi=100; bn=b->size; cn=c->size; bexp=b->exp0; cexp=c->exp0; exptot=bexp+cexp; maxsize=b->maxsize; maxs=maxsize-2; kkk=bn+cn; n=1; while (n < kkk) { n=2*n; } /* while */ n=2*n; n2=2*n; q=(double *)malloc((n2+1)*sizeof(double)); for (i=0; i<=n2; i++) { q[i]=0.0; } /* for i */ isize=1.0/sqrt((double)(n)); k=0; for (i=1; i<=bn; i++) { b1=b->coef[i] % bi; q[k]=(double)(b1)*isize; b1=b->coef[i] / bi; q[k+2]=(double)(b1)*isize; k=k+4; } /* for i */ k=0; for (i=1; i<=cn; i++) { c1=c->coef[i] % bi; q[k+1]=(double)(c1)*isize; c1=c->coef[i] / bi; q[k+3]=(double)(c1)*isize; k=k+4; } /* for i */ cdft(n2,1,q); mul(q,n); cdft(n2,-1,q); maxerror=0.0; rr=0.0; for (i=0; i<=n2; i=i+2) { qqq=q[i]; temp=fabs(qqq-floor(qqq+0.5)); if (temp > maxerror) { maxerror=temp; } /* if */ xx=qqq+rr; rr=floor((xx+0.5)*bri); q[i]=xx-br*rr; } /* for i */ if (maxerror > fftmaxerror) { fftmaxerror=maxerror; } /* if */ gg=1; if (kkk > maxs) { /* printf("* %8d %8d\n",kkk,maxs); */ gg=kkk-maxs+1; } /* if */ k=4*gg-4; for (i=gg; i<=kkk; i++) { a->coef[i-gg+1]=bi*(int)(q[k+2]+0.5)+(int)(q[k]+0.5); k=k+4; } /* for i */ if (a->coef[kkk-gg+1] > 0) { a->size=kkk-gg+1; a->exp0=exptot+gg-1; } else { a->size=kkk-gg; a->exp0=exptot+gg-1; } /* if */ a->sgn0=(b->sgn0)*(c->sgn0); free(q); } /* if */ } /* bigmul */ void big2l(a_big,b_big,a_l,b_l,nn) bigint *a_big; bigint *b_big; short *a_l; short *b_l; int *nn; { int i,ii,l,il,dl,ib,k,kk,jj,maxexp; l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; if (b_big->exp0 > a_big->exp0) { maxexp=b_big->exp0; } else { maxexp=a_big->exp0; } /* if */ if (b_big->size > a_big->size) { ii=b_big->size; } else { ii=a_big->size; } /* if */ k=l-il; for (i=ii; i >= 1; i--) { if (k >= 1) { jj=i+maxexp-b_big->exp0; if (jj <= b_big->size) { b_l[k]=b_big->coef[jj]; } /* if */ kk=i+maxexp-a_big->exp0; if (kk <= a_big->size) { a_l[k]=a_big->coef[kk]; } /* if */ } /* if */ k--; } /* for i */ } /* big2l */ double julian_day(year,month,day) int year,month,day; { double f,g,a,jdd; if (month>=3) { f=(double)(year); g=(double)(month); } else { f=(double)(year-1); g=(double)(month+12); } /* if */ a=2.0-floor(f/100.0)+floor(f/400.0); jdd=floor(365.25*f)+floor(30.6001*(g+1.0))+ (double)(day)+a+1720994.5; return(jdd); } /* julian_day */ double get_float_days(ttime) time_t *ttime; { int year,month,day,hour,min,sec; double fdays,fhours; struct tm *tt; tt=localtime(ttime); year=tt->tm_year+1900; month=tt->tm_mon+1; day=tt->tm_mday; hour=tt->tm_hour; min=tt->tm_min; sec=tt->tm_sec; fhours=(double)(hour)+(double)(min)/60.0+(double)(sec)/3600.0; fdays=julian_day(year,month,day)+fhours/24.0; return(fdays); } /* get_float_days */ void day_and_time(ttime) time_t *ttime; { int year,month,day,hour,min,sec; struct tm *tt; tt=localtime(ttime); year=tt->tm_year+1900; month=tt->tm_mon+1; day=tt->tm_mday; hour=tt->tm_hour; min=tt->tm_min; sec=tt->tm_sec; printf("%04d-%02d-%02d %02d:%02d:%02d",year,month,day,hour,min,sec); } /* day_and_time */ void fday_and_time(outfile,ttime) FILE *outfile; time_t *ttime; { int year,month,day,hour,min,sec; struct tm *tt; tt=localtime(ttime); year=tt->tm_year+1900; month=tt->tm_mon+1; day=tt->tm_mday; hour=tt->tm_hour; min=tt->tm_min; sec=tt->tm_sec; fprintf(outfile,"%04d-%02d-%02d %02d:%02d:%02d",year,month,day,hour,min,sec); } /* fday_and_time */ void sum_e(rQ,rT,n1,n2,maxsize) bigint *rQ; bigint *rT; int n1; int n2; int maxsize; { int nm,n; bigint LQ,LT; bigint RQ,RT; if ((n2-n1) == 1) { n=n1; if (n == 0) { inttobig(rQ,1); inttobig(rT,1); } else { inttobig(rQ,n); inttobig(rT,1); } /* if */ } else { biginit(&LQ,maxsize); biginit(<,maxsize); biginit(&RQ,maxsize); biginit(&RT,maxsize); nm=(n1+n2)/2; sum_e(&LQ,<,n1,nm,maxsize); sum_e(&RQ,&RT,nm,n2,maxsize); bigmul(rQ,&LQ,&RQ); bigmul(&tmp,&RQ,<); bigadd(rT,&RT,&tmp); bigfree(&LQ); bigfree(<); bigfree(&RQ); bigfree(&RT); } /* if */ } /* sum_e */ void sum_e_inv(rP,rQ,rT,n1,n2,maxsize) bigint *rP; bigint *rQ; bigint *rT; int n1; int n2; int maxsize; { int nm,n; bigint LP,LQ,LT; bigint RP,RQ,RT; if ((n2-n1) == 1) { n=n1; if (n == 0) { inttobig(rP,1); inttobig(rQ,1); inttobig(rT,1); } else { inttobig(rP,1); rP->sgn0=-1; inttobig(rQ,n); inttobig(rT,1); rT->sgn0=-1; } /* if */ } else { biginit(&LP,maxsize); biginit(&LQ,maxsize); biginit(<,maxsize); biginit(&RP,maxsize); biginit(&RQ,maxsize); biginit(&RT,maxsize); nm=(n1+n2)/2; sum_e_inv(&LP,&LQ,<,n1,nm,maxsize); sum_e_inv(&RP,&RQ,&RT,nm,n2,maxsize); bigmul(rP,&LP,&RP); bigmul(rQ,&LQ,&RQ); bigmul(&tmp,&RQ,<); bigmul(&tmp2,&LP,&RT); bigadd(rT,&tmp,&tmp2); bigfree(&LP); bigfree(&LQ); bigfree(<); bigfree(&RP); bigfree(&RQ); bigfree(&RT); } /* if */ } /* sum_e_inv */ void sum_arctanh(rP,rQ,rB,rT,n1,n2,u,v,maxsize) bigint *rP; bigint *rQ; bigint *rB; bigint *rT; int n1; int n2; int u; int v; int maxsize; { int nm,n; bigint LP,LQ,LB,LT; bigint RP,RQ,RB,RT; if ((n2-n1) == 1) { n=n1; if (n == 0) { inttobig(rP,u); inttobig(rQ,v); inttobig(rB,1); inttobig(rT,u); } else { inttobig(rP,u*u); inttobig(rQ,v*v); inttobig(rB,(2*n+1)); inttobig(rT,u*u); } /* if */ } else { biginit(&LP,maxsize); biginit(&LQ,maxsize); biginit(&LB,maxsize); biginit(<,maxsize); biginit(&RP,maxsize); biginit(&RQ,maxsize); biginit(&RB,maxsize); biginit(&RT,maxsize); nm=(n1+n2)/2; sum_arctanh(&LP,&LQ,&LB,<,n1,nm,u,v,maxsize); sum_arctanh(&RP,&RQ,&RB,&RT,nm,n2,u,v,maxsize); bigmul(rP,&LP,&RP); bigmul(rQ,&LQ,&RQ); bigmul(rB,&LB,&RB); bigmul(&tmp,&RQ,<); bigmul(&tmp,&tmp,&RB); bigmul(&tmp2,&LP,&RT); bigmul(&tmp2,&tmp2,&LB); bigadd(rT,&tmp,&tmp2); bigfree(&LP); bigfree(&LQ); bigfree(&LB); bigfree(<); bigfree(&RP); bigfree(&RQ); bigfree(&RB); bigfree(&RT); } /* if */ } /* sum_arctanh */ void ln2_term(sum,u,v,w,sgn0,n1,n2,nn,bigsize) short *sum; int u; int v; int w; int sgn0; int n1; int n2; int *nn; int bigsize; { int lsize; clock_t time1,timediff1; short *a,*b,*ab; bigint rP; bigint rQ; bigint rB; bigint rT; biginit2(&rP,bigsize); biginit2(&rQ,bigsize); biginit2(&rB,bigsize); biginit2(&rT,bigsize); inttobig(&rP,0); inttobig(&rQ,0); inttobig(&rB,0); inttobig(&rT,0); biginit2(&tmp,bigsize); biginit2(&tmp2,bigsize); inttobig(&tmp,0); inttobig(&tmp2,0); time1=clock(); sum_arctanh(&rP,&rQ,&rB,&rT,n1,n2,u,v,bigsize); timediff1=clock()-time1; printf("Time term, fast series: %10.2f\n",(double)(timediff1)/timetic); bigfree(&rP); lsize=nn[0]; a = (short *)calloc(lsize+2,sizeof(short)); b = (short *)calloc(lsize+2,sizeof(short)); ab = (short *)calloc(lsize+2,sizeof(short)); bigass(&tmp,&rT); bigmul(&tmp2,&rB,&rQ); bigfree(&rQ); bigfree(&rB); bigfree(&rT); lzero(a,nn); lzero(b,nn); lzero(ab,nn); big2l(&tmp,&tmp2,a,b,nn); bigfree(&tmp); bigfree(&tmp2); time1=clock(); ldivv(ab,a,b,nn); timediff1=clock()-time1; printf("Time div: %10.2f\n",(double)(timediff1)/timetic); lcmul(a,ab,w,0,nn); lass(b,sum,nn); if (sgn0 == (-1)) { lsub(sum,b,a,nn); } else { ladd(sum,b,a,nn); } /* if */ free(a); free(b); free(ab); } /* ln2_term */ int main(argc,argv) int argc; char *argv[]; { FILE *outfile,*mvfile; int bigsize,lsize,ff; int i,n,length,ilength,pot,pot0; int mem1,mem2,decimals,decimals2; int n1,n2,method,method1,maxmem; double elapseddays,frac,temp; double lt,el,le; int nn[5]; clock_t time2,timediff2; clock_t time3,timediff3; time_t start_time,stop_time; short *a,*b,*e_result,*ln2_result; bigint rP; bigint rQ; bigint rT; if ((mvfile = fopen(mv_file_name,"r")) == NULL) { printf("\n"); printf("Warning, file %s not found, mulversion now set to 1!\n",mv_file_name); printf("\n"); mulversion=1; } else { fscanf(mvfile,"%d",&mulversion); fclose(mvfile); if ((mulversion != 1) && (mulversion != 2)) { printf("\n"); printf("Warning, wrong value of mulversion in %s, mulversion now set to 1!\n",mv_file_name); printf("\n"); mulversion=1; } /* if */ } /* if */ if (argc != 3) { printf("\n"); printf("Usage: %s methodi power-of-2i\n",argv[0]); printf("\n"); printf("%s\n",version_string); printf("\n"); printf("method = 1 => exp(1); e = SUM 1/k!\n"); printf("method = 2 => exp(1); 1/e = SUM (-1)^k/k!\n"); printf("\n"); printf("method = 3 => ln(2) = 144*atanh(1/251)+54*atanh(1/449)-38*atanh(1/4801)+62*atanh(1/8749)\n"); printf("method = 4 => ln(2) = 72*atanh(1/127)+54*atanh(1/449)+34*atanh(1/4801)-10*atanh(1/8749)\n"); printf("method = 5 => ln(2) = 18*atanh(1/26)-2*atanh(1/4801)+8*atanh(1/8749)\n"); printf("method = 6 => ln(2) = 14*atanh(1/31)+10*atanh(1/49)+6*atanh(1/161)\n"); printf("method = 7 => ln(2) = 6*atanh(1/26)+8*atanh(1/31)+10*atanh(1/49)\n"); printf("\n"); printf("The number is calculated with about 2^(power-of-2) decimal digits.\n"); printf("\n"); printf("The variable mulversion is at present set to %1d.\n",mulversion); printf("It is possible to change the value (to %1d) by editing the file \"%s\".\n", 3-mulversion,mv_file_name); printf("Set the variable to 2 for calculations with more than 2^24 (16 million) digits,\n"); printf("or more important at least when \"FFT max error\" > 0.25.\n"); printf("\n"); exit(0); } /* if */ method1=atoi(argv[1]); pot0=atoi(argv[2]); if ((method1 < 1) || (method1 > 7)) { printf("\n"); printf("Wrong method number.\n"); printf("\n"); exit(0); } /* if */ if (pot0 < 8) { pot0=8; } /* if */ method=method1-2; /* exp(1), method1 = 1 ***************************************************************************/ if (method1 == 1) { printf("\n"); printf("e = exp(1), e = SUM 1/k!\n"); printf("========================\n"); printf("\n"); printf("%s\n",version_string); printf("\n"); pot=pot0; pot=pot-2; length=1; for (i=1; i<=pot; i++) { length=2*length; } /* for i */ bigsize=length; ilength=8; frac=0.995; decimals=(int)(frac*(double)(length-ilength))-2; decimals=4*decimals; printf("Decimals: %10d FFT: %10d\n\n",decimals,2*mulversion*length); lt=1.0/log(10.0); el=0.0; le=(double)(4*(length-6)); i=1; while (el= 1) && (method <= 5)) { printf("\n"); printf("ln(2) fast series\n"); printf("================= \n"); printf("\n"); printf("%s\n",version_string); printf("\n"); pot=pot0; n=1; for (i=1; i<=pot; i++) { n=2*n; } /* for i */ decimals=n; decimals2=decimals; frac=0.995; decimals=decimals/4; decimals=(int)(frac*(double)(decimals-8))-2; decimals=4*decimals; bigsize=n/4; length=n/4; ilength=8; frac=0.995; nn[0]=length; nn[1]=ilength; nn[2]=length-ilength; nn[3]=(int)(frac*(double)(length-ilength))-2; lsize=nn[0]; printf("mulversion = %10d\n",mulversion); printf("Power of 2 = %10d\n",pot); printf("n = %10d\n",n); printf("FFT = %10d\n",2*mulversion*bigsize); printf("maxsize = %10d\n",bigsize); printf("decimals = %10d (%d)\n",decimals,decimals2); printf("\n"); printf("Number of calculated decimals: %10d\n",decimals); printf("\n"); temp=(double)(bigsize); mem1=(int)(3.6*1048576.0*log(temp)*temp/(log(1048576.0/8.0)*1048576.0/8.0)+0.5); mem2=8*2*2*mulversion*bigsize; maxmem=mem1+mem2; maxmem=maxmem/1024; printf("Max approx memory usage: %10d kbyte; mulversion: %3d\n",maxmem,mulversion); printf("\n"); time(&start_time); printf("Start time: "); day_and_time(&start_time); time2=clock(); printf("\n"); printf("\n"); printf("Calculating ln(2) ...\n"); printf("\n"); fftmaxerror=0.0; /* method = 1 *****************************/ if (method == 1) { printf("ln(2) = 144*atanh(1/251)+54*atanh(1/449)-38*atanh(1/4801)+62*atanh(1/8749)\n"); printf("\n"); n1=0; n2=(int)((double)(decimals2)/4.70+0.5)+1; printf("n1 = %10d\n",n1); printf("n2 = %10d\n",n2); printf("\n"); ln2_result = (short *)calloc(lsize+2,sizeof(short)); lzero(ln2_result,nn); ln2_term(ln2_result,1,251,144,1,n1,n2,nn,bigsize); ln2_term(ln2_result,1,449,54,1,n1,n2,nn,bigsize); ln2_term(ln2_result,1,4801,38,(-1),n1,n2,nn,bigsize); ln2_term(ln2_result,1,8749,62,1,n1,n2,nn,bigsize); } /* if */ /* method = 2 *****************************/ if (method == 2) { printf("ln(2) = 72*atanh(1/127)+54*atanh(1/449)+34*atanh(1/4801)-10*atanh(1/8749)\n"); printf("\n"); n1=0; n2=(int)((double)(decimals2)/4.10+0.5)+1; printf("n1 = %10d\n",n1); printf("n2 = %10d\n",n2); printf("\n"); ln2_result = (short *)calloc(lsize+2,sizeof(short)); lzero(ln2_result,nn); ln2_term(ln2_result,1,127,72,1,n1,n2,nn,bigsize); ln2_term(ln2_result,1,449,54,1,n1,n2,nn,bigsize); ln2_term(ln2_result,1,4801,34,1,n1,n2,nn,bigsize); ln2_term(ln2_result,1,8749,10,(-1),n1,n2,nn,bigsize); } /* if */ /* method = 3 *****************************/ if (method == 3) { printf("ln(2) = 18*atanh(1/26)-2*atanh(1/4801)+8*atanh(1/8749)\n"); printf("\n"); n1=0; n2=(int)((double)(decimals2)/2.75+0.5)+1; printf("n1 = %10d\n",n1); printf("n2 = %10d\n",n2); printf("\n"); ln2_result = (short *)calloc(lsize+2,sizeof(short)); lzero(ln2_result,nn); ln2_term(ln2_result,1,26,18,1,n1,n2,nn,bigsize); ln2_term(ln2_result,1,4801,2,(-1),n1,n2,nn,bigsize); ln2_term(ln2_result,1,8749,8,1,n1,n2,nn,bigsize); } /* if */ /* method = 4 *****************************/ if (method == 4) { printf("ln(2) = 14*atanh(1/31)+10*atanh(1/49)+6*atanh(1/161)\n"); printf("\n"); n1=0; n2=(int)((double)(decimals2)/2.95+0.5)+1; printf("n1 = %10d\n",n1); printf("n2 = %10d\n",n2); printf("\n"); ln2_result = (short *)calloc(lsize+2,sizeof(short)); lzero(ln2_result,nn); ln2_term(ln2_result,1,31,14,1,n1,n2,nn,bigsize); ln2_term(ln2_result,1,49,10,1,n1,n2,nn,bigsize); ln2_term(ln2_result,1,161,6,1,n1,n2,nn,bigsize); } /* if */ /* method = 5 *****************************/ if (method == 5) { printf("ln(2) = 6*atanh(1/26)+8*atanh(1/31)+10*atanh(1/49)\n"); printf("\n"); n1=0; n2=(int)((double)(decimals2)/2.75+0.5)+1; printf("n1 = %10d\n",n1); printf("n2 = %10d\n",n2); printf("\n"); ln2_result = (short *)calloc(lsize+2,sizeof(short)); lzero(ln2_result,nn); ln2_term(ln2_result,1,26,6,1,n1,n2,nn,bigsize); ln2_term(ln2_result,1,31,8,1,n1,n2,nn,bigsize); ln2_term(ln2_result,1,49,10,1,n1,n2,nn,bigsize); } /* if */ /* Final result ln(2) ********************/ printf("\n"); printf("FFT max error (should be < 0.25): %16.14f\n",fftmaxerror); timediff2=clock()-time2; printf("\n"); printf("Time TOTAL ln(2): %10.2f\n",(double)(timediff2)/timetic); printf("\n"); time(&stop_time); printf("Stop time: "); day_and_time(&stop_time); printf("\n"); printf("\n"); elapseddays=get_float_days(&stop_time)-get_float_days(&start_time); printf("Elapsed clock time: %18.8lf days = %18.6lf hours = %18.2lf seconds\n", elapseddays,24.0*elapseddays,24.0*3600.0*elapseddays); printf("\n"); printf("Calculation finished, output to %s.\n",file_name2); printf("\n"); outfile=fopen(file_name2,"w"); lfprint(ln2_result,nn,outfile); fclose(outfile); } /* if */ } /* End */