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;
				}
			}
		}
	}
}