summaryrefslogtreecommitdiff
path: root/appl/math/ffts.b
diff options
context:
space:
mode:
Diffstat (limited to 'appl/math/ffts.b')
-rw-r--r--appl/math/ffts.b639
1 files changed, 639 insertions, 0 deletions
diff --git a/appl/math/ffts.b b/appl/math/ffts.b
new file mode 100644
index 00000000..e6132289
--- /dev/null
+++ b/appl/math/ffts.b
@@ -0,0 +1,639 @@
+implement FFTs;
+include "sys.m";
+ sys: Sys;
+ print: import sys;
+include "math.m";
+ math: Math;
+ cos, sin, Degree, Pi: import math;
+include "ffts.m";
+
+# by r. c. singleton, stanford research institute, sept. 1968
+# translated to limbo by eric grosse, jan 1997
+# arrays at(maxf), ck(maxf), bt(maxf), sk(maxf), and np(maxp)
+# are used for temporary storage. if the available storage
+# is insufficient, the program exits.
+# maxf must be >= the maximum prime factor of n.
+# maxp must be > the number of prime factors of n.
+# in addition, if the square-free portion k of n has two or
+# more prime factors, then maxp must be >= k-1.
+# array storage in nfac for a maximum of 15 prime factors of n.
+# if n has more than one square-free factor, the product of the
+# square-free factors must be <= 210
+
+ffts(a,b:array of real, ntot,n,nspan,isn:int){
+ maxp: con 209;
+ i,ii,inc,j,jc,jf,jj,k,k1,k2,k3,k4,kk:int;
+ ks,kspan,kspnn,kt,m,maxf,nn,nt:int;
+ aa,aj,ajm,ajp,ak,akm,akp,bb,bj,bjm,bjp,bk,bkm,bkp:real;
+ c1,c2,c3,c72,cd,rad,radf,s1,s2,s3,s72,s120,sd:real;
+ maxf = 23;
+ if(math == nil){
+ sys = load Sys Sys->PATH;
+ math = load Math Math->PATH;
+ }
+ nfac := array[12] of int;
+ np := array[maxp] of int;
+ at := array[23] of real;
+ ck := array[23] of real;
+ bt := array[23] of real;
+ sk := array[23] of real;
+
+ if(n<2) return;
+ inc = isn;
+ c72 = cos(72.*Degree);
+ s72 = sin(72.*Degree);
+ s120 = sin(120.*Degree);
+ rad = 2.*Pi;
+ if(isn<0){
+ s72 = -s72;
+ s120 = -s120;
+ rad = -rad;
+ inc = -inc;
+ }
+ nt = inc*ntot;
+ ks = inc*nspan;
+ kspan = ks;
+ nn = nt-inc;
+ jc = ks/n;
+ radf = rad*real(jc)*0.5;
+ i = 0;
+ jf = 0;
+
+ # determine the factors of n
+ m = 0;
+ k = n;
+ while(k==k/16*16){
+ m = m+1;
+ nfac[m] = 4;
+ k = k/16;
+ }
+ j = 3;
+ jj = 9;
+ for(;;)
+ if(k%jj==0){
+ m = m+1;
+ nfac[m] = j;
+ k = k/jj;
+ }else{
+ j = j+2;
+ jj = j*j;
+ if(jj>k)
+ break;
+ }
+ if(k<=4){
+ kt = m;
+ nfac[m+1] = k;
+ if(k!=1)
+ m = m+1;
+ }else{
+ if(k==k/4*4){
+ m = m+1;
+ nfac[m] = 2;
+ k = k/4;
+ }
+ kt = m;
+ j = 2;
+ do{
+ if(k%j==0){
+ m = m+1;
+ nfac[m] = j;
+ k = k/j;
+ }
+ j = ((j+1)/2)*2+1;
+ }while(j<=k);
+ }
+ if(kt!=0){
+ j = kt;
+ do{
+ m = m+1;
+ nfac[m] = nfac[j];
+ j = j-1;
+ }while(j!=0);
+ }
+
+ for(;;){ # compute fourier transform
+ sd = radf/real(kspan);
+ cd = sin(sd);
+ cd = 2.0*cd*cd;
+ sd = sin(sd+sd);
+ kk = 1;
+ i = i+1;
+ if(nfac[i]==2){ # transform for factor of 2 (including rotation factor)
+ kspan = kspan/2;
+ k1 = kspan+2;
+ for(;;){
+ k2 = kk+kspan;
+ ak = a[k2-1];
+ bk = b[k2-1];
+ a[k2-1] = a[kk-1]-ak;
+ b[k2-1] = b[kk-1]-bk;
+ a[kk-1] = a[kk-1]+ak;
+ b[kk-1] = b[kk-1]+bk;
+ kk = k2+kspan;
+ if(kk>nn){
+ kk = kk-nn;
+ if(kk>jc)
+ break;
+ }
+ }
+ if(kk>kspan)
+ break;
+ do{
+ c1 = 1.0-cd;
+ s1 = sd;
+ for(;;){
+ k2 = kk+kspan;
+ ak = a[kk-1]-a[k2-1];
+ bk = b[kk-1]-b[k2-1];
+ a[kk-1] = a[kk-1]+a[k2-1];
+ b[kk-1] = b[kk-1]+b[k2-1];
+ a[k2-1] = c1*ak-s1*bk;
+ b[k2-1] = s1*ak+c1*bk;
+ kk = k2+kspan;
+ if(kk>=nt){
+ k2 = kk-nt;
+ c1 = -c1;
+ kk = k1-k2;
+ if(kk<=k2){
+ ak = c1-(cd*c1+sd*s1);
+ s1 = (sd*c1-cd*s1)+s1;
+ c1 = 2.0-(ak*ak+s1*s1);
+ s1 = c1*s1;
+ c1 = c1*ak;
+ kk = kk+jc;
+ if(kk>=k2)
+ break;
+ }
+ }
+ }
+ k1 = k1+inc+inc;
+ kk = (k1-kspan)/2+jc;
+ }while(kk<=jc+jc);
+ }else{ # transform for factor of 4
+ if(nfac[i]!=4){
+ # transform for odd factors
+ k = nfac[i];
+ kspnn = kspan;
+ kspan = kspan/k;
+ if(k==3)
+ for(;;){
+ # transform for factor of 3 (optional code)
+ k1 = kk+kspan;
+ k2 = k1+kspan;
+ ak = a[kk-1];
+ bk = b[kk-1];
+ aj = a[k1-1]+a[k2-1];
+ bj = b[k1-1]+b[k2-1];
+ a[kk-1] = ak+aj;
+ b[kk-1] = bk+bj;
+ ak = -0.5*aj+ak;
+ bk = -0.5*bj+bk;
+ aj = (a[k1-1]-a[k2-1])*s120;
+ bj = (b[k1-1]-b[k2-1])*s120;
+ a[k1-1] = ak-bj;
+ b[k1-1] = bk+aj;
+ a[k2-1] = ak+bj;
+ b[k2-1] = bk-aj;
+ kk = k2+kspan;
+ if(kk>=nn){
+ kk = kk-nn;
+ if(kk>kspan)
+ break;
+ }
+ }
+ else if(k==5){
+ # transform for factor of 5 (optional code)
+ c2 = c72*c72-s72*s72;
+ s2 = 2.0*c72*s72;
+ for(;;){
+ k1 = kk+kspan;
+ k2 = k1+kspan;
+ k3 = k2+kspan;
+ k4 = k3+kspan;
+ akp = a[k1-1]+a[k4-1];
+ akm = a[k1-1]-a[k4-1];
+ bkp = b[k1-1]+b[k4-1];
+ bkm = b[k1-1]-b[k4-1];
+ ajp = a[k2-1]+a[k3-1];
+ ajm = a[k2-1]-a[k3-1];
+ bjp = b[k2-1]+b[k3-1];
+ bjm = b[k2-1]-b[k3-1];
+ aa = a[kk-1];
+ bb = b[kk-1];
+ a[kk-1] = aa+akp+ajp;
+ b[kk-1] = bb+bkp+bjp;
+ ak = akp*c72+ajp*c2+aa;
+ bk = bkp*c72+bjp*c2+bb;
+ aj = akm*s72+ajm*s2;
+ bj = bkm*s72+bjm*s2;
+ a[k1-1] = ak-bj;
+ a[k4-1] = ak+bj;
+ b[k1-1] = bk+aj;
+ b[k4-1] = bk-aj;
+ ak = akp*c2+ajp*c72+aa;
+ bk = bkp*c2+bjp*c72+bb;
+ aj = akm*s2-ajm*s72;
+ bj = bkm*s2-bjm*s72;
+ a[k2-1] = ak-bj;
+ a[k3-1] = ak+bj;
+ b[k2-1] = bk+aj;
+ b[k3-1] = bk-aj;
+ kk = k4+kspan;
+ if(kk>=nn){
+ kk = kk-nn;
+ if(kk>kspan)
+ break;
+ }
+ }
+ }else{
+ if(k!=jf){
+ jf = k;
+ s1 = rad/real(k);
+ c1 = cos(s1);
+ s1 = sin(s1);
+ if(jf>maxf){
+ sys->fprint(sys->fildes(2),"too many primes for fft");
+ exit;
+ }
+ ck[jf-1] = 1.0;
+ sk[jf-1] = 0.0;
+ j = 1;
+ do{
+ ck[j-1] = ck[k-1]*c1+sk[k-1]*s1;
+ sk[j-1] = ck[k-1]*s1-sk[k-1]*c1;
+ k = k-1;
+ ck[k-1] = ck[j-1];
+ sk[k-1] = -sk[j-1];
+ j = j+1;
+ }while(j<k);
+ }
+ for(;;){
+ k1 = kk;
+ k2 = kk+kspnn;
+ aa = a[kk-1];
+ bb = b[kk-1];
+ ak = aa;
+ bk = bb;
+ j = 1;
+ k1 = k1+kspan;
+ do{
+ k2 = k2-kspan;
+ j = j+1;
+ at[j-1] = a[k1-1]+a[k2-1];
+ ak = at[j-1]+ak;
+ bt[j-1] = b[k1-1]+b[k2-1];
+ bk = bt[j-1]+bk;
+ j = j+1;
+ at[j-1] = a[k1-1]-a[k2-1];
+ bt[j-1] = b[k1-1]-b[k2-1];
+ k1 = k1+kspan;
+ }while(k1<k2);
+ a[kk-1] = ak;
+ b[kk-1] = bk;
+ k1 = kk;
+ k2 = kk+kspnn;
+ j = 1;
+ do{
+ k1 = k1+kspan;
+ k2 = k2-kspan;
+ jj = j;
+ ak = aa;
+ bk = bb;
+ aj = 0.0;
+ bj = 0.0;
+ k = 1;
+ do{
+ k = k+1;
+ ak = at[k-1]*ck[jj-1]+ak;
+ bk = bt[k-1]*ck[jj-1]+bk;
+ k = k+1;
+ aj = at[k-1]*sk[jj-1]+aj;
+ bj = bt[k-1]*sk[jj-1]+bj;
+ jj = jj+j;
+ if(jj>jf)
+ jj = jj-jf;
+ }while(k<jf);
+ k = jf-j;
+ a[k1-1] = ak-bj;
+ b[k1-1] = bk+aj;
+ a[k2-1] = ak+bj;
+ b[k2-1] = bk-aj;
+ j = j+1;
+ }while(j<k);
+ kk = kk+kspnn;
+ if(kk>nn){
+ kk = kk-nn;
+ if(kk>kspan)
+ break;
+ }
+ }
+ }
+ # multiply by rotation factor (except for factors of 2 and 4)
+ if(i==m)
+ break;
+ kk = jc+1;
+ do{
+ c2 = 1.0-cd;
+ s1 = sd;
+ do{
+ c1 = c2;
+ s2 = s1;
+ kk = kk+kspan;
+ for(;;){
+ ak = a[kk-1];
+ a[kk-1] = c2*ak-s2*b[kk-1];
+ b[kk-1] = s2*ak+c2*b[kk-1];
+ kk = kk+kspnn;
+ if(kk>nt){
+ ak = s1*s2;
+ s2 = s1*c2+c1*s2;
+ c2 = c1*c2-ak;
+ kk = kk-nt+kspan;
+ if(kk>kspnn)
+ break;
+ }
+ }
+ c2 = c1-(cd*c1+sd*s1);
+ s1 = s1+(sd*c1-cd*s1);
+ c1 = 2.0-(c2*c2+s1*s1);
+ s1 = c1*s1;
+ c2 = c1*c2;
+ kk = kk-kspnn+jc;
+ }while(kk<=kspan);
+ kk = kk-kspan+jc+inc;
+ }while(kk<=jc+jc);
+ }else{
+ kspnn = kspan;
+ kspan = kspan/4;
+ do{
+ c1 = 1.;
+ s1 = 0.;
+ for(;;){
+ k1 = kk+kspan;
+ k2 = k1+kspan;
+ k3 = k2+kspan;
+ akp = a[kk-1]+a[k2-1];
+ akm = a[kk-1]-a[k2-1];
+ ajp = a[k1-1]+a[k3-1];
+ ajm = a[k1-1]-a[k3-1];
+ a[kk-1] = akp+ajp;
+ ajp = akp-ajp;
+ bkp = b[kk-1]+b[k2-1];
+ bkm = b[kk-1]-b[k2-1];
+ bjp = b[k1-1]+b[k3-1];
+ bjm = b[k1-1]-b[k3-1];
+ b[kk-1] = bkp+bjp;
+ bjp = bkp-bjp;
+ do10 := 0;
+ if(isn<0){
+ akp = akm+bjm;
+ akm = akm-bjm;
+ bkp = bkm-ajm;
+ bkm = bkm+ajm;
+ if(s1!=0.) do10 = 1;
+ }else{
+ akp = akm-bjm;
+ akm = akm+bjm;
+ bkp = bkm+ajm;
+ bkm = bkm-ajm;
+ if(s1!=0.) do10 = 1;
+ }
+ if(do10){
+ a[k1-1] = akp*c1-bkp*s1;
+ b[k1-1] = akp*s1+bkp*c1;
+ a[k2-1] = ajp*c2-bjp*s2;
+ b[k2-1] = ajp*s2+bjp*c2;
+ a[k3-1] = akm*c3-bkm*s3;
+ b[k3-1] = akm*s3+bkm*c3;
+ kk = k3+kspan;
+ if(kk<=nt)
+ continue;
+ }else{
+ a[k1-1] = akp;
+ b[k1-1] = bkp;
+ a[k2-1] = ajp;
+ b[k2-1] = bjp;
+ a[k3-1] = akm;
+ b[k3-1] = bkm;
+ kk = k3+kspan;
+ if(kk<=nt)
+ continue;
+ }
+ c2 = c1-(cd*c1+sd*s1);
+ s1 = (sd*c1-cd*s1)+s1;
+ c1 = 2.0-(c2*c2+s1*s1);
+ s1 = c1*s1;
+ c1 = c1*c2;
+ c2 = c1*c1-s1*s1;
+ s2 = 2.0*c1*s1;
+ c3 = c2*c1-s2*s1;
+ s3 = c2*s1+s2*c1;
+ kk = kk-nt+jc;
+ if(kk>kspan)
+ break;
+ }
+ kk = kk-kspan+inc;
+ }while(kk<=jc);
+ if(kspan==jc)
+ break;
+ }
+ }
+ } # end "compute fourier transform"
+
+ # permute the results to normal order---done in two stages
+ # permutation for square factors of n
+ np[0] = ks;
+ if(kt!=0){
+ k = kt+kt+1;
+ if(m<k)
+ k = k-1;
+ j = 1;
+ np[k] = jc;
+ do{
+ np[j] = np[j-1]/nfac[j];
+ np[k-1] = np[k]*nfac[j];
+ j = j+1;
+ k = k-1;
+ }while(j<k);
+ k3 = np[k];
+ kspan = np[1];
+ kk = jc+1;
+ k2 = kspan+1;
+ j = 1;
+ if(n!=ntot){
+ for(;;){
+ # permutation for multivariate transform
+ k = kk+jc;
+ do{
+ ak = a[kk-1];
+ a[kk-1] = a[k2-1];
+ a[k2-1] = ak;
+ bk = b[kk-1];
+ b[kk-1] = b[k2-1];
+ b[k2-1] = bk;
+ kk = kk+inc;
+ k2 = k2+inc;
+ }while(kk<k);
+ kk = kk+ks-jc;
+ k2 = k2+ks-jc;
+ if(kk>=nt){
+ k2 = k2-nt+kspan;
+ kk = kk-nt+jc;
+ if(k2>=ks)
+ permm: for(;;){
+ k2 = k2-np[j-1];
+ j = j+1;
+ k2 = np[j]+k2;
+ if(k2<=np[j-1]){
+ j = 1;
+ do{
+ if(kk<k2)
+ break permm;
+ kk = kk+jc;
+ k2 = kspan+k2;
+ }while(k2<ks);
+ if(kk>=ks)
+ break permm;
+ }
+ }
+ }
+ }
+ jc = k3;
+ }else{
+ for(;;){
+ # permutation for single-variate transform (optional code)
+ ak = a[kk-1];
+ a[kk-1] = a[k2-1];
+ a[k2-1] = ak;
+ bk = b[kk-1];
+ b[kk-1] = b[k2-1];
+ b[k2-1] = bk;
+ kk = kk+inc;
+ k2 = kspan+k2;
+ if(k2>=ks)
+ perms: for(;;){
+ k2 = k2-np[j-1];
+ j = j+1;
+ k2 = np[j]+k2;
+ if(k2<=np[j-1]){
+ j = 1;
+ do{
+ if(kk<k2)
+ break perms;
+ kk = kk+inc;
+ k2 = kspan+k2;
+ }while(k2<ks);
+ if(kk>=ks)
+ break perms;
+ }
+ }
+ }
+ jc = k3;
+ }
+ }
+ if(2*kt+1>=m)
+ return;
+ kspnn = np[kt];
+ # permutation for square-free factors of n
+ j = m-kt;
+ nfac[j+1] = 1;
+ do{
+ nfac[j] = nfac[j]*nfac[j+1];
+ j = j-1;
+ }while(j!=kt);
+ kt = kt+1;
+ nn = nfac[kt]-1;
+ if(nn<=maxp){
+ jj = 0;
+ j = 0;
+ for(;;){
+ k2 = nfac[kt];
+ k = kt+1;
+ kk = nfac[k];
+ j = j+1;
+ if(j>nn)
+ break;
+ for(;;){
+ jj = kk+jj;
+ if(jj<k2)
+ break;
+ jj = jj-k2;
+ k2 = kk;
+ k = k+1;
+ kk = nfac[k];
+ }
+ np[j-1] = jj;
+ }
+ # determine the permutation cycles of length greater than 1
+ j = 0;
+ for(;;){
+ j = j+1;
+ kk = np[j-1];
+ if(kk>=0)
+ if(kk==j){
+ np[j-1] = -j;
+ if(j==nn)
+ break;
+ }else{
+ do{
+ k = kk;
+ kk = np[k-1];
+ np[k-1] = -kk;
+ }while(kk!=j);
+ k3 = kk;
+ }
+ }
+ maxf = inc*maxf;
+ for(;;){
+ j = k3+1;
+ nt = nt-kspnn;
+ ii = nt-inc+1;
+ if(nt<0)
+ break;
+ for(;;){
+ j = j-1;
+ if(np[j-1]>=0){
+ jj = jc;
+ do{
+ kspan = jj;
+ if(jj>maxf)
+ kspan = maxf;
+ jj = jj-kspan;
+ k = np[j-1];
+ kk = jc*k+ii+jj;
+ k1 = kk+kspan;
+ k2 = 0;
+ do{
+ k2 = k2+1;
+ at[k2-1] = a[k1-1];
+ bt[k2-1] = b[k1-1];
+ k1 = k1-inc;
+ }while(k1!=kk);
+ do{
+ k1 = kk+kspan;
+ k2 = k1-jc*(k+np[k-1]);
+ k = -np[k-1];
+ do{
+ a[k1-1] = a[k2-1];
+ b[k1-1] = b[k2-1];
+ k1 = k1-inc;
+ k2 = k2-inc;
+ }while(k1!=kk);
+ kk = k2;
+ }while(k!=j);
+ k1 = kk+kspan;
+ k2 = 0;
+ do{
+ k2 = k2+1;
+ a[k1-1] = at[k2-1];
+ b[k1-1] = bt[k2-1];
+ k1 = k1-inc;
+ }while(k1!=kk);
+ }while(jj!=0);
+ if(j==1)
+ break;
+ }
+ }
+ }
+ }
+}
generated by cgit v1.2.3 (git 2.46.0) at 2025-12-10 16:09:34 +0000