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

Lines 1–500 of 639. Load more