/***********************************************************************/ /* */ /* O */ /* OOOOOOOOOOOOO */ /* OOOOOOOOOOOOO */ /* O OO OO */ /* OO OO */ /* OO OO */ /* OO OO */ /* OO OO */ /* OO OO */ /* */ /* */ /* Program: */ /* -------- */ /* */ /* sspieln2.c */ /* */ /* */ /* Version 1.21: */ /* ------------- */ /* */ /* First version: September 1993 */ /* 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: */ /* Minor improvements to some procedures. */ /* */ /* 2003-07-18: */ /* The multiplication version variable (mulversion) is now set */ /* (to 1 or 2) in the file "mulver.txt". The program is renamed */ /* to "sspieln2.c". Faster methods for calculating e=exp(1) */ /* are now moved to the program "sseln2.c". */ /* */ /* 2003-06-30: */ /* Improved FFT-multiplication (lmul). Inversions, divisions and */ /* square roots up to 25 % faster. Inversions and divisions with */ /* 4-th order convergence. */ /* */ /* 2003-06-24: */ /* Very minor changes. */ /* */ /* 2003-06-15: */ /* Very minor (non critical) changes. The program now use the FFT */ /* fftsg_h.c by Takuya Ooura. */ /* */ /* 2003-02-19: */ /* All unused variables removed from the source code. */ /* */ /* 2002-08-19: */ /* Minor changes to this text header. */ /* */ /* 2000-09-02: */ /* Very minor changes. */ /* */ /* 1999-11-05: */ /* Version for Takuya Ooura's very fast FFT. */ /* */ /* 1999-09-25: */ /* Minor speed improvements. The C-code has also been a little bit */ /* "prettier" and the procedures easier to use (fewer temporary */ /* parameters). */ /* */ /* 1999-09-18: */ /* Significant speed improvements. Improvements to inversions, */ /* divisions and square roots. Borwein's 2-th order algorithm for */ /* pi added. */ /* */ /* 1999-07-28: */ /* AGM algorithm for ln(2) and Newtons Iteration for e=exp(1) with */ /* AGM for ln(x). */ /* */ /* Sometimes during 1995-1996: */ /* 4-th order algorithms added for square roots. */ /* */ /* First version September 1993: */ /* Algorithms for pi (Borwein's 4-th order and Gauss-Legendre AGM). */ /* Simple algorithm for e=exp(1) using SUM 1/k!. Algorithms for */ /* the Golden Section and sqrt(n). FFT-multiplication and Newtons */ /* Iteration for inversion, division and square roots. */ /* */ /* */ /* 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! */ /* */ /* */ /* Function: */ /* --------- */ /* */ /* This C program calculates pi, e=exp(1), ln(2), the Golden */ /* Section = (sqrt(5)+1)/2 and sqrt(n) to about 2^power decimals. */ /* */ /* If for example (power of 2) = 10 the calculations are made with */ /* about 2^10 = 1024 digits. */ /* */ /* The program can quite easily be modified for other calculations. */ /* */ /* */ /* Algorithms: */ /* ----------- */ /* */ /* The program uses (see also the references above): */ /* */ /* * FFT for multiplications. */ /* */ /* * 4-th order methods for inversions and divisions. */ /* */ /* * 4-th order methods for square and cube roots. */ /* */ /* * Borweins 4-th order or 2-th order or Gauss-Legendre */ /* (AGM) algorithm for Pi. */ /* */ /* * Simple sum or Newtons Iteration for e=exp(1) (AGM for ln(x)). */ /* */ /* * AGM for ln(2). */ /* */ /* AGM means Arithmetic Geometric Mean. */ /* */ /* */ /* Accuracy: */ /* --------- */ /* */ /* All printed digits are probably correct. The results for */ /* all the numbers have been verified to 2^20 digits. */ /* */ /* */ /* 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 and Linux and I recommend that you */ /* optimize the program (i.e. use cc -O3 option or similar option) */ /* when you compile (gcc if you use a GNU compiler). */ /* */ /*>cc -O3 -o sspieln2 sspieln2.c fftsg_h.c -lm */ /* */ /* Observe that on some computers (HP) the code will be faster if the */ /* variables in the FFT are declared as static. Then remove the */ /* comments around words "static" in the FFT code. */ /* */ /* The time wraps around after 4295 seconds on some computers. */ /* */ /* */ /* Output: */ /* ------- */ /* */ /* Output is to the file "sspieln2.txt". This file name can easily be */ /* changed. */ /* */ /* */ /* Timing, benchmark and memory: */ /* ----------------------------- */ /* */ /* 1. 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 sspieln2 sspieln2.c fftsg_h.c -lm */ /* */ /* TC below means Time Complexity. NI means Newtons Iteration. */ /* The times have been measured for "power of 2" equals 17 and 20 with */ /* the "time" function in Linux. The variable mulversion = 1. */ /* */ /* Intel Celeron 1400 MHz */ /* ---------------------- */ /* */ /* Algorithm 2^17 2^20 TC Memory/(2^20) */ /*---------------------------------------------------------------------*/ /* Pi, Borwein 4-th 26.89 s 384.24 s O(n*log(n)^3) 13.4 Mbyte */ /* Pi, Borwein 2-th 32.91 s 372.56 s O(n*log(n)^3) 14.9 Mbyte */ /* Pi, AGM 22.92 s 289.99 s O(n*log(n)^3) 13.4 Mbyte */ /* e, SUM 1/k! 24.87 s 1922.22 s O(n^2) 12.8 Mbyte */ /* e, NI, AGM ln(x) 731.08 s --- s O(n*log(n)^4) 33.8 Mbyte */ /* ln(2), AGM 102.31 s 1287.43 s O(n*log(n)^3) 30.8 Mbyte */ /* Golden Section 1.11 s 12.11 s O(n*log(n)^2) 11.8 Mbyte */ /* sqrt(2) 0.99 s 9.45 s O(n*log(n)^2) 11.3 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 typedef struct biginttag { short *coef; int size,maxsize,exp0,sgn0; } bigint; double timetic=(double)(CLOCKS_PER_SEC); char version_string[80]="Program sspieln2.c, version 1.21, 2007-08-13"; char file_name[20]="sspieln2.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 pi=3.14159265358979323; double fftmaxerror; 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 */ int main(argc,argv) int argc; char *argv[]; { FILE *outfile,*mvfile; int ff,il,ib,dl,ll2,ll,lsize; int i,j,k,length,ilength,dec,pot,l,nr,n; int qqi,remi,basint,l2,l22; int dl2,decc,mm,decc2,maxmem; double frac,el,le,dll,lt,r,q2; double jd,bbas,bbasi,iif; int nn[5]; int sss[5]; clock_t time1,timediff1,time2,timediff2,time3,timediff3; short *t1,*t2,*a,*b,*c,*d,*e,*f,*g,*h,*y,*y2,*logy; int *aa,*bb; 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 */ printf("\n\n"); printf(" O\n"); printf(" OOOOOOOOOOOOO\n"); printf(" OOOOOOOOOOOOO\n"); printf(" O OO OO\n"); printf(" OO OO\n"); printf(" OO OO\n"); printf(" OO OO\n"); printf(" OO OO\n"); printf(" OO OO\n"); printf("\n\n"); printf("%s\n",version_string); 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"); printf("0: Quit\n\n"); printf("1: Pi, Borweins 4-th order\n\n"); printf("2: Pi, Borweins 2-th order\n\n"); printf("3: Pi, Gauss-Legendre, AGM\n\n"); printf("4: e = exp(1), e = SUM 1/k!\n\n"); printf("5: e = exp(1), Newtons Iteration (using AGM for ln(x))\n\n"); printf("6: ln(2), AGM\n\n"); printf("7: The Golden Section = (sqrt(5)+1)/2\n\n"); printf("8: Sqrt(n), (n integer)\n\n"); printf("Input nr: "); scanf("%d",&nr); printf("\n\n"); if (nr == 1) { printf("Pi, Borweins 4-th order\n"); printf("=======================\n\n"); printf("Input power of 2: "); scanf("%d",&pot); printf("\n\n"); pot=pot-2; length=1; for (i=1; i<=pot; i++) { length=2*length; } /* for i */ ilength=2; frac=0.995; nn[0]=length; nn[1]=ilength; nn[2]=length-ilength; nn[3]=(int)(frac*(double)(length-ilength))-2; l=length; lsize=length; a = (short *)calloc(lsize+2,sizeof(short)); b = (short *)calloc(lsize+2,sizeof(short)); c = (short *)calloc(lsize+2,sizeof(short)); d = (short *)calloc(lsize+2,sizeof(short)); e = (short *)calloc(lsize+2,sizeof(short)); f = (short *)calloc(lsize+2,sizeof(short)); time2=clock(); printf("Decimals: %10d FFT: %10d\n",4*nn[3],4*mulversion*2*length/4); maxmem=2*(6+4)*length+8*4*mulversion*length; maxmem=maxmem/1024; printf("\nMax memory usage: %10d kbyte\n",maxmem); printf("\nCalculating ...\n"); fftmaxerror=0.0; lzero(a,nn); lzero(b,nn); lzero(c,nn); lzero(d,nn); lzero(e,nn); lzero(f,nn); linttol(d,2,0,nn); time3=clock(); lsqrt(c,d,nn); timediff3=clock()-time3; printf("\nTime sqrt: %10.2f\n",(double)(timediff3)/timetic); linttol(d,1,0,nn); lsub(b,c,d,nn); linttol(e,6,0,nn); lcmul(d,c,4,0,nn); lsub(a,e,d,nn); printf("\n it decimals scpu\n\n"); time1=clock(); k=1; iif=8.0; do { time3=clock(); lmul(c,b,b,nn); lmul(c,c,c,nn); linttol(d,1,0,nn); lsub(e,d,c,nn); lsqrtsqrt(b,e,nn); lsub(e,d,b,nn); ladd(c,d,b,nn); ldivv(b,e,c,nn); ladd(c,b,d,nn); lmul(c,c,c,nn); lmul(c,c,c,nn); lmul(c,a,c,nn); lmul(e,b,b,nn); ladd(f,b,e,nn); ladd(e,f,d,nn); lmul(f,e,b,nn); lcmul2(d,f,iif,0,nn); lsub(a,c,d,nn); dec=lgetdec(d,nn); timediff3=clock()-time3; printf("%6d %12d %10.2f\n",k,dec,(double)(timediff3)/timetic); iif=4.0*iif; k++; } while (lcompeps(d,nn)==1); timediff1=clock()-time1; printf("\nTime main loop: %10.2f\n",(double)(timediff1)/timetic); time3=clock(); linv(c,a,nn); timediff3=clock()-time3; printf("\nTime div: %10.2f\n",(double)(timediff3)/timetic); timediff2=clock()-time2; printf("\nTime total: %10.2f\n",(double)(timediff2)/timetic); printf("\n"); printf("FFT max error (should be < 0.25): %16.14f\n",fftmaxerror); printf("\n"); printf("Calculation finished, output to %s.\n\n",file_name); outfile=fopen(file_name,"w"); lfprint(c,nn,outfile); fclose(outfile); printf("\n"); } /* if */ if (nr == 2) { printf("Pi, Borweins 2-th order\n"); printf("=======================\n\n"); printf("Input power of 2: "); scanf("%d",&pot); printf("\n\n"); pot=pot-2; length=1; for (i=1; i<=pot; i++) { length=2*length; } /* for i */ ilength=2; frac=0.995; nn[0]=length; nn[1]=ilength; nn[2]=length-ilength; nn[3]=(int)(frac*(double)(length-ilength))-2; l=length; lsize=length; a = (short *)calloc(lsize+2,sizeof(short)); b = (short *)calloc(lsize+2,sizeof(short)); c = (short *)calloc(lsize+2,sizeof(short)); d = (short *)calloc(lsize+2,sizeof(short)); e = (short *)calloc(lsize+2,sizeof(short)); f = (short *)calloc(lsize+2,sizeof(short)); g = (short *)calloc(lsize+2,sizeof(short)); h = (short *)calloc(lsize+2,sizeof(short)); t1 = (short *)calloc(lsize+2,sizeof(short)); time2=clock(); printf("Decimals: %10d FFT: %10d\n",4*nn[3],4*mulversion*2*length/4); maxmem=2*(9+4)*length+8*4*mulversion*length; maxmem=maxmem/1024; printf("\nMax memory usage: %10d kbyte\n",maxmem); printf("\nCalculating ...\n"); fftmaxerror=0.0; lzero(a,nn); lzero(b,nn); lzero(c,nn); lzero(d,nn); lzero(e,nn); lzero(f,nn); lzero(g,nn); lzero(h,nn); time3=clock(); linttol(d,2,0,nn); lsqrtinv(g,d,nn); lmul(e,g,d,nn); lass(h,e,nn); linttol(t1,2,0,nn); ladd(c,t1,e,nn); lsqrtinv(f,e,nn); lmul(e,f,f,nn); lmul(e,e,f,nn); lmul(b,e,d,nn); lass(a,b,nn); timediff3=clock()-time3; printf("\nTime init: %10.2f\n",(double)(timediff3)/timetic); printf("\n it decimals scpu\n\n"); time1=clock(); k=1; do { time3=clock(); ladd(t1,a,f,nn); lcmul(h,t1,5000,-1,nn); lsqrtinv(f,h,nn); lmul(a,h,f,nn); lmul(t1,b,a,nn); ladd(e,t1,f,nn); linttol(t1,1,0,nn); ladd(d,b,t1,nn); linv(g,d,nn); lmul(b,e,g,nn); linttol(t1,1,0,nn); ladd(e,h,t1,nn); lmul(g,g,e,nn); lmul(c,c,g,nn); linttol(t1,1,0,nn); lsub(d,g,t1,nn); dec=lgetdec(d,nn); timediff3=clock()-time3; printf("%6d %12d %10.2f\n",k,dec,(double)(timediff3)/timetic); k++; } while (lcompeps(d,nn)==1); timediff1=clock()-time1; printf("\nTime main loop: %10.2f\n",(double)(timediff1)/timetic); timediff2=clock()-time2; printf("\nTime total: %10.2f\n",(double)(timediff2)/timetic); printf("\n"); printf("FFT max error (should be < 0.25): %16.14f\n",fftmaxerror); printf("\n"); printf("Calculation finished, output to %s.\n\n",file_name); outfile=fopen(file_name,"w"); lfprint(c,nn,outfile); fclose(outfile); printf("\n"); } /* if */ if (nr == 3) { printf("Pi, Gauss-Legendre, AGM\n"); printf("=======================\n\n"); printf("Input power of 2: "); scanf("%d",&pot); printf("\n\n"); pot=pot-2; length=1; for (i=1; i<=pot; i++) { length=2*length; } /* for i */ ilength=2; frac=0.995; nn[0]=length; nn[1]=ilength; nn[2]=length-ilength; nn[3]=(int)(frac*(double)(length-ilength))-2; l=length; lsize=length; a = (short *)calloc(lsize+2,sizeof(short)); b = (short *)calloc(lsize+2,sizeof(short)); c = (short *)calloc(lsize+2,sizeof(short)); d = (short *)calloc(lsize+2,sizeof(short)); e = (short *)calloc(lsize+2,sizeof(short)); t1 = (short *)calloc(lsize+2,sizeof(short)); time2=clock(); printf("Decimals: %10d FFT: %10d\n",4*nn[3],4*mulversion*2*length/4); maxmem=2*(6+4)*length+8*4*mulversion*length; maxmem=maxmem/1024; printf("\nMax memory usage: %10d kbyte\n",maxmem); printf("\nCalculating ...\n"); fftmaxerror=0.0; lzero(a,nn); lzero(b,nn); lzero(c,nn); lzero(d,nn); lzero(e,nn); linttol(a,1,0,nn); linttol(d,2,0,nn); time3=clock(); lsqrtinv(b,d,nn); timediff3=clock()-time3; printf("\nTime sqrt: %10.2f\n",(double)(timediff3)/timetic); printf("\n it decimals scpu\n\n"); linttol(e,2500,-1,nn); time1=clock(); k=1; iif=1.0; do { time3=clock(); lass(d,a,nn); lmul(c,a,b,nn); ladd(t1,a,b,nn); lcmul(a,t1,5000,-1,nn); lsqrt(b,c,nn); lsub(c,d,a,nn); lmul(d,c,c,nn); lcmul2(c,d,iif,0,nn); lsub(t1,e,c,nn); lass(e,t1,nn); dec=lgetdec(d,nn); timediff3=clock()-time3; printf("%6d %12d %10.2f\n",k,dec,(double)(timediff3)/timetic); iif=2.0*iif; k++; } while (lcompeps(d,nn)==1); timediff1=clock()-time1; printf("\nTime main loop: %10.2f\n",(double)(timediff1)/timetic); ladd(t1,a,b,nn); lmul(a,t1,t1,nn); lcmul2(b,e,4.0,0,nn); time3=clock(); ldivv(c,a,b,nn); timediff3=clock()-time3; printf("\nTime div: %10.2f\n",(double)(timediff3)/timetic); timediff2=clock()-time2; printf("\nTime total: %10.2f\n",(double)(timediff2)/timetic); printf("\n"); printf("FFT max error (should be < 0.25): %16.14f\n",fftmaxerror); printf("\n"); printf("Calculation finished, output to %s.\n\n",file_name); outfile=fopen(file_name,"w"); lfprint(c,nn,outfile); fclose(outfile); printf("\n"); } /* if */ if (nr == 4) { printf("e = exp(1), e = SUM 1/k!\n"); printf("========================\n\n"); printf("Input power of 2: "); scanf("%d",&pot); printf("\n\n"); pot=pot-2; length=1; for (i=1; i<=pot; i++) { length=2*length; } /* for i */ ilength=2; nn[0]=length; nn[1]=ilength; nn[2]=length-ilength; nn[3]=(int)(0.995*(double)(length-ilength))-2; l=length; l2=l/2; l22=l2+4; lsize=length; a = (short *)calloc(lsize+2,sizeof(short)); b = (short *)calloc(lsize+2,sizeof(short)); c = (short *)calloc(lsize+2,sizeof(short)); aa = (int *)calloc(l22+1,sizeof(int)); bb = (int *)calloc(l22+1,sizeof(int)); printf("Decimals: %10d FFT: %10d\n\n",4*nn[3],4*mulversion*2*length/4); time2=clock(); lzero(a,nn); lzero(b,nn); lzero(c,nn); for (i=1; i<=l22; i++) { aa[i]=0; bb[i]=0; } /* for i */ bbas=100000000.0; bbasi=0.00000001; basint=100000000; lt=1.0/log(10.0); el=0.0; le=(double)(4*nn[2]); i=1; while (el=1; j--) { jd=(double)(j); dll=dll+lt*log(jd); ll=(int)(0.125*dll+0.5)+3; r=0.0; for (i=1; i<=ll; i++) { q2=(double)(aa[i])*jd+r+0.5; r=(double)((int)(bbasi*q2)); aa[i]=(int)(q2-bbas*r); } /* for i */ remi=0; for (i=1; i<=ll; i++) { qqi=bb[i]+aa[i]+remi; if (qqi >= basint) { bb[i]=qqi-basint; remi=1; } else { bb[i]=qqi; remi=0; } /* if */ } /* for i */ } /* for j */ timediff1=clock()-time1; printf("\nTime main loop: %10.2f\n",(double)(timediff1)/timetic); ll2=ib+(int)(0.25*dll+0.5)+2; k=1; for (i=ib; i<=ll2; i=i+2) { a[i]=aa[k] % 10000; a[i+1]=aa[k] / 10000; b[i]=bb[k] % 10000; b[i+1]=bb[k] / 10000; k++; } /* for i */ a[l+1]=1; b[l+1]=1; time3=clock(); ldivv(c,b,a,nn); timediff3=clock()-time3; printf("\nTime div: %10.2f\n",(double)(timediff3)/timetic); timediff2=clock()-time2; printf("\nTime total e = exp(1): %10.2f\n",(double)(timediff2)/timetic); printf("\n"); printf("FFT max error (should be < 0.25): %16.14f\n",fftmaxerror); printf("\n"); printf("Calculation finished, output to %s.\n\n",file_name); outfile=fopen(file_name,"w"); lfprint(c,nn,outfile); fclose(outfile); printf("\n"); } /* if */ if (nr == 5) { printf("e=exp(1), Newtons Iteration (using AGM for ln(x))\n"); printf("=================================================\n\n"); printf("Input power of 2: "); scanf("%d",&pot); printf("\n\n"); pot=pot-1; length=1; for (i=1; i<=pot; i++) { length=2*length; } /* for i */ ilength=4; frac=0.995; nn[0]=length; nn[1]=ilength; nn[2]=length-ilength; nn[3]=(int)(frac*(double)(length-ilength))-2; sss[0]=length; sss[1]=ilength; sss[2]=length/2-ilength; if ((pot+1) < 12) { sss[3]=(int)(frac*(double)(length/2-ilength/2))-10; } else { sss[3]=(int)(frac*(double)(length/2-ilength/2))-2; } /* if */ dl2=nn[2]/2; l=length; lsize=length; a = (short *)calloc(lsize+2,sizeof(short)); b = (short *)calloc(lsize+2,sizeof(short)); c = (short *)calloc(lsize+2,sizeof(short)); d = (short *)calloc(lsize+2,sizeof(short)); e = (short *)calloc(lsize+2,sizeof(short)); f = (short *)calloc(lsize+2,sizeof(short)); g = (short *)calloc(lsize+2,sizeof(short)); h = (short *)calloc(lsize+2,sizeof(short)); y = (short *)calloc(lsize+2,sizeof(short)); y2 = (short *)calloc(lsize+2,sizeof(short)); logy = (short *)calloc(lsize+2,sizeof(short)); t1 = (short *)calloc(lsize+2,sizeof(short)); t2 = (short *)calloc(lsize+2,sizeof(short)); time2=clock(); printf("Decimals: %10d FFT: %10d\n",4*sss[3],4*mulversion*2*length/4); maxmem=2*(13+4)*length+8*4*mulversion*length; maxmem=maxmem/1024; printf("\nMax memory usage: %10d kbyte\n",maxmem); printf("\nCalculating ...\n"); fftmaxerror=0.0; lzero(a,nn); lzero(b,nn); lzero(c,nn); lzero(d,nn); lzero(e,nn); lzero(f,nn); lzero(g,nn); lzero(h,nn); linttol(a,1,0,nn); linttol(b,1,-dl2,nn); time3=clock(); time1=clock(); k=1; iif=1.0; do { lmul(c,a,b,nn); ladd(t1,a,b,nn); lcmul(a,t1,5000,-1,nn); decc=lgetdec(c,nn); decc2=decc; decc=decc/4; if ((decc % 2) == 1) { decc=decc-1; } /* if */ lzero(f,nn); l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; mm=l-il+1; for (i=1; i<=(l-decc); i++) { f[mm]=c[mm-decc]; mm--; } /* for i */ lsqrt(h,f,sss); decc=decc/2; lzero(b,nn); mm=l-il+1; for (i=1; i<=(l-decc); i++) { b[mm-decc]=h[mm]; mm--; } /* for i */ lmul(t1,a,a,nn); lsub(d,t1,c,nn); lcmul2(t2,d,iif,0,nn); ladd(t1,e,t2,nn); lass(e,t1,nn); dec=lgetdec(d,nn); iif=2.0*iif; k++; } while (lcompeps(d,sss)==1); linttol(b,5000,-1,nn); lsub(c,b,e,nn); linv(f,c,sss); timediff1=clock()-time1; printf("\nTime init: %10.2f\n",(double)(timediff1)/timetic); lzero(a,nn); lzero(b,nn); lzero(c,nn); lzero(d,nn); lzero(e,nn); /* lzero(f,nn); */ lzero(g,nn); lzero(h,nn); lzero(y,nn); lzero(y2,nn); lzero(logy,nn); y[l+1]=1; y[l-il+1]=2; y[l-il]=7182; y[l-il-1]=8182; y[l-il-2]=8459; y[l-il-3]=452; y[l-il-4]=3536; y[l-il-5]=287; y[l-il-6]=4713; y[l-il-7]=5266; printf("\n it decimals scpu\n\n"); k=1; do { time1=clock(); lass(y2,y,nn); linv(c,y,sss); lcmul(b,c,1,-dl2,nn); linttol(a,1,0,nn); lzero(e,nn); iif=1.0; do { lmul(c,a,b,nn); ladd(t1,a,b,nn); lcmul(a,t1,5000,-1,nn); decc=lgetdec(c,nn); decc2=decc; decc=decc/4; if ((decc % 2) == 1) { decc=decc-1; } /* if */ lzero(g,nn); l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; mm=l-il+1; for (i=1; i<=(l-decc); i++) { g[mm]=c[mm-decc]; mm--; } /* for i */ lsqrt(h,g,sss); decc=decc/2; lzero(b,nn); mm=l-il+1; for (i=1; i<=(l-decc); i++) { b[mm-decc]=h[mm]; mm--; } /* for i */ lmul(t1,a,a,nn); lsub(d,t1,c,nn); lcmul2(t2,d,iif,0,nn); ladd(t1,e,t2,nn); lass(e,t1,nn); dec=lgetdec(d,nn); iif=2.0*iif; } while (lcompeps(d,sss)==1); linttol(b,5000,-1,nn); lsub(c,b,e,nn); time3=clock(); linv(g,c,sss); timediff3=clock()-time3; lsub(logy,g,f,nn); linttol(c,2,0,nn); lsub(d,c,logy,nn); lmul(y,d,y,nn); lsub(d,y,y2,nn); dec=lgetdec(d,nn); timediff1=clock()-time1; printf("%6d %12d %10.2f\n",k,dec,(double)(timediff1)/timetic); k++; } while (lcompeps(d,sss)==1); lass(c,y,nn); timediff2=clock()-time2; printf("\nTime total: %10.2f\n",(double)(timediff2)/timetic); printf("\n"); printf("FFT max error (should be < 0.25): %16.14f\n",fftmaxerror); printf("\n"); printf("Calculation finished, output to %s.\n\n",file_name); outfile=fopen(file_name,"w"); lfprint(c,sss,outfile); fclose(outfile); printf("\n"); } /* if */ if (nr == 6) { printf("ln(2), AGM\n"); printf("==========\n\n"); printf("Input power of 2: "); scanf("%d",&pot); printf("\n\n"); pot=pot-1; length=1; for (i=1; i<=pot; i++) { length=2*length; } /* for i */ ilength=4; frac=0.995; nn[0]=length; nn[1]=ilength; nn[2]=length-ilength; nn[3]=(int)(frac*(double)(length-ilength))-2; sss[0]=length; sss[1]=ilength; sss[2]=length/2-ilength; if ((pot+1) < 12) { sss[3]=(int)(frac*(double)(length/2-ilength/2))-10; } else { sss[3]=(int)(frac*(double)(length/2-ilength/2))-2; } /* if */ dl2=nn[2]/2; l=length; lsize=length; a = (short *)calloc(lsize+2,sizeof(short)); b = (short *)calloc(lsize+2,sizeof(short)); c = (short *)calloc(lsize+2,sizeof(short)); d = (short *)calloc(lsize+2,sizeof(short)); e = (short *)calloc(lsize+2,sizeof(short)); f = (short *)calloc(lsize+2,sizeof(short)); g = (short *)calloc(lsize+2,sizeof(short)); h = (short *)calloc(lsize+2,sizeof(short)); t1 = (short *)calloc(lsize+2,sizeof(short)); t2 = (short *)calloc(lsize+2,sizeof(short)); time2=clock(); printf("Decimals: %10d FFT: %10d\n",4*sss[3],4*mulversion*2*length/4); maxmem=2*(10+4)*length+8*4*mulversion*length; maxmem=maxmem/1024; printf("\nMax memory usage: %10d kbyte\n",maxmem); printf("\nCalculating ...\n"); fftmaxerror=0.0; lzero(a,nn); lzero(b,nn); lzero(c,nn); lzero(d,nn); lzero(e,nn); lzero(f,nn); lzero(g,nn); lzero(h,nn); linttol(a,1,0,nn); linttol(b,1,-dl2,nn); time3=clock(); printf("\n it decimals scpu\n\n"); time1=clock(); k=1; iif=1.0; do { time3=clock(); lmul(c,a,b,nn); ladd(t1,a,b,nn); lcmul(a,t1,5000,-1,nn); decc=lgetdec(c,nn); decc2=decc; decc=decc/4; if ((decc % 2) == 1) { decc=decc-1; } /* if */ lzero(f,nn); l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; mm=l-il+1; for (i=1; i<=(l-decc); i++) { f[mm]=c[mm-decc]; mm--; } /* for i */ lsqrt(h,f,sss); decc=decc/2; lzero(b,nn); mm=l-il+1; for (i=1; i<=(l-decc); i++) { b[mm-decc]=h[mm]; mm--; } /* for i */ lmul(t1,a,a,nn); lsub(d,t1,c,nn); lcmul2(t2,d,iif,0,nn); ladd(t1,e,t2,nn); lass(e,t1,nn); dec=lgetdec(d,nn); timediff3=clock()-time3; printf("%6d %12d %10.2f\n",k,dec,(double)(timediff3)/timetic); iif=2.0*iif; k++; } while (lcompeps(d,sss)==1); timediff1=clock()-time1; printf("\nTime main loop: %10.2f\n",(double)(timediff1)/timetic); linttol(b,5000,-1,nn); lsub(c,b,e,nn); time3=clock(); linv(f,c,sss); timediff3=clock()-time3; printf("\nTime div: %10.2f\n",(double)(timediff3)/timetic); lzero(a,nn); lzero(b,nn); lzero(c,nn); lzero(d,nn); lzero(e,nn); /* lzero(f,nn); */ lzero(g,nn); lzero(h,nn); linttol(a,1,0,nn); linttol(b,5000,-dl2-1,nn); time3=clock(); printf("\n it decimals scpu\n\n"); time1=clock(); k=1; iif=1.0; do { time3=clock(); lmul(c,a,b,nn); ladd(t1,a,b,nn); lcmul(a,t1,5000,-1,nn); decc=lgetdec(c,nn); decc2=decc; decc=decc/4; if ((decc % 2) == 1) { decc=decc-1; } /* if */ lzero(g,nn); l=nn[0]; il=nn[1]; dl=nn[2]; ib=l-il-dl+1; mm=l-il+1; for (i=1; i<=(l-decc); i++) { g[mm]=c[mm-decc]; mm--; } /* for i */ lsqrt(h,g,sss); decc=decc/2; lzero(b,nn); mm=l-il+1; for (i=1; i<=(l-decc); i++) { b[mm-decc]=h[mm]; mm--; } /* for i */ lmul(t1,a,a,nn); lsub(d,t1,c,nn); lcmul2(t2,d,iif,0,nn); ladd(t1,e,t2,nn); lass(e,t1,nn); dec=lgetdec(d,nn); timediff3=clock()-time3; printf("%6d %12d %10.2f\n",k,dec,(double)(timediff3)/timetic); iif=2.0*iif; k++; } while (lcompeps(d,sss)==1); timediff1=clock()-time1; printf("\nTime main loop: %10.2f\n",(double)(timediff1)/timetic); linttol(b,5000,-1,nn); lsub(c,b,e,nn); time3=clock(); linv(g,c,sss); timediff3=clock()-time3; printf("\nTime div: %10.2f\n",(double)(timediff3)/timetic); lsub(c,g,f,nn); timediff2=clock()-time2; printf("\nTime total: %10.2f\n",(double)(timediff2)/timetic); printf("\n"); printf("FFT max error (should be < 0.25): %16.14f\n",fftmaxerror); printf("\n"); printf("Calculation finished, output to %s.\n\n",file_name); outfile=fopen(file_name,"w"); lfprint(c,sss,outfile); fclose(outfile); printf("\n"); } /* if */ if (nr == 7) { printf("The Golden Section = (sqrt(5)+1)/2\n"); printf("==================================\n\n"); printf("Input power of 2: "); scanf("%d",&pot); printf("\n\n"); pot=pot-2; length=1; for (i=1; i<=pot; i++) { length=2*length; } /* for i */ ilength=2; frac=0.995; nn[0]=length; nn[1]=ilength; nn[2]=length-ilength; nn[3]=(int)(frac*(double)(length-ilength))-2; l=length; lsize=length; a = (short *)calloc(lsize+2,sizeof(short)); b = (short *)calloc(lsize+2,sizeof(short)); c = (short *)calloc(lsize+2,sizeof(short)); time2=clock(); printf("Decimals: %10d FFT: %10d\n",4*nn[3],4*mulversion*2*length/4); maxmem=2*(3+4)*length+8*4*mulversion*length; maxmem=maxmem/1024; printf("\nMax memory usage: %10d kbyte\n",maxmem); printf("\nCalculating ...\n"); fftmaxerror=0.0; lzero(a,nn); lzero(b,nn); lzero(c,nn); time3=clock(); linttol(a,5,0,nn); lsqrt(b,a,nn); linttol(c,1,0,nn); ladd(a,b,c,nn); lcmul(b,a,5000,-1,nn); timediff3=clock()-time3; printf("\nTime GS: %10.2f\n",(double)(timediff3)/timetic); printf("\n"); printf("FFT max error (should be < 0.25): %16.14f\n",fftmaxerror); printf("\n"); printf("Calculation finished, output to %s.\n\n",file_name); outfile=fopen(file_name,"w"); lfprint(b,nn,outfile); fclose(outfile); printf("\n"); } /* if */ if (nr == 8) { printf("Sqrt(n), n integer\n"); printf("==================\n\n"); printf("Input n (1 <= n <= 9999): "); scanf("%d",&n); if ((n < 1) || (n > 9999)) { exit(0); } /* if */ printf("\n\n"); printf("Input power of 2: "); scanf("%d",&pot); printf("\n\n"); pot=pot-2; length=1; for (i=1; i<=pot; i++) { length=2*length; } /* for i */ ilength=2; frac=0.995; nn[0]=length; nn[1]=ilength; nn[2]=length-ilength; nn[3]=(int)(frac*(double)(length-ilength))-2; l=length; lsize=length; a = (short *)calloc(lsize+2,sizeof(short)); b = (short *)calloc(lsize+2,sizeof(short)); time2=clock(); printf("Decimals: %10d FFT: %10d\n",4*nn[3],4*mulversion*2*length/4); maxmem=2*(2+4)*length+8*4*mulversion*length; maxmem=maxmem/1024; printf("\nMax memory usage: %10d kbyte\n",maxmem); printf("\nCalculating ...\n"); fftmaxerror=0.0; lzero(a,nn); lzero(b,nn); time3=clock(); linttol(a,n,0,nn); lsqrt(b,a,nn); timediff3=clock()-time3; printf("\nTime Sqrt: %10.2f\n",(double)(timediff3)/timetic); printf("\n"); printf("FFT max error (should be < 0.25): %16.14f\n",fftmaxerror); printf("\n"); printf("Calculation finished, output to %s.\n\n",file_name); outfile=fopen(file_name,"w"); lfprint(b,nn,outfile); fclose(outfile); printf("\n"); } /* if */ } /* End */