#include <stdio.h>
#include <math.h>

/*fast fourier transform                                                   */
/*     THE COOLEY-TUKEY FAST FOURIER TRANSFORM                             */
/*     EVALUATES COMPLEX FOURIER SERIES FOR COMPLEX OR REAL FUNCTIONS.     */
/*     THAT IS, IT COMPUTES                                                */
/*     FTRAN(J1,J2,...)=SUM(DATA(I1,I2,...)*W1**(I1-1)*(J1-1)              */
/*                                         *W2**(I2-1)*(J2-1)*...),        */
/*     WHERE W1=EXP(-2*PI*SQRT(-1)/NN(1)), W2=EXP(-2*PI*SQRT(-1)/NN(2)),   */
/*     ETC. AND I1 AND J1 RUN FROM 1 TO NN(1), I2 AND J2 RUN FROM 1 TO     */
/*     NN(2), ETC.  THERE IS NO LIMIT ON THE DIMENSIONALITY (NUMBER OF     */
/*     SUBSCRIPTS) OF THE ARRAY OF DATA.  THE PROGRAM WILL PERFORM         */
/*     A THREE-DIMENSIONAL FOURIER TRANSFORM AS EASILY AS A ONE-DIMEN-     */
/*     SIONAL ONE, THO IN A PROPORTIONATELY GREATER TIME.  AN INVERSE      */
/*     TRANSFORM CAN BE PERFORMED, IN WHICH THE SIGN IN THE EXPONENTIALS   */
/*     IS +, INSTEAD OF -.  IF AN INVERSE TRANSFORM IS PERFORMED UPON      */
/*     AN ARRAY OF TRANSFORMED DATA, THE ORIGINAL DATA WILL REAPPEAR,      */
/*     MULTIPLIED BY NN(1)*NN(2)*...  THE ARRAY OF INPUT DATA MAY BE       */
/*     REAL OR COMPLEX, AT THE PROGRAMMERS OPTION, WITH A SAVING OF        */
/*     ABOUT THIRTY PER CENT IN RUNNING TIME FOR REAL OVER COMPLEX.        */
/*     (FOR FASTEST TRANSFORM OF REAL DATA, NN(1) SHOULD BE EVEN.)         */
/*     THE TRANSFORM VALUES ARE ALWAYS COMPLEX, AND ARE RETURNED IN THE    */
/*     ORIGINAL ARRAY OF DATA, REPLACING THE INPUT DATA.  THE LENGTH       */
/*     OF EACH DIMENSION OF THE DATA ARRAY MAY BE ANY INTEGER.  THE        */
/*     PROGRAM RUNS FASTER ON COMPOSITE INTEGERS THAN ON PRIMES, AND IS    */
/*     PARTICULARLY FAST ON NUMBERS RICH IN FACTORS OF TWO.                */
/*     TIMING IS IN FACT GIVEN BY THE FOLLOWING FORMULA.  LET NTOT BE THE  */
/*     TOTAL NUMBER OF POINTS (REAL OR COMPLEX) IN THE DATA ARRAY, THAT    */
/*     IS, NTOT=NN(1)*NN(2)*...  DECOMPOSE NTOT INTO ITS PRIME FACTORS,    */
/*     SUCH AS 2**K2 * 3**K3 * 5**K5 * ...  LET SUM2 BE THE SUM OF ALL     */
/*     THE FACTORS OF TWO IN NTOT, THAT IS, SUM2 = 2*K2.  LET SUMF BE      */
/*     THE SUM OF ALL OTHER FACTORS OF NTOT, THAT IS, SUMF = 3*K3+5*K5+..  */
/*     THE TIME TAKEN BY A MULTIDIMENSIONAL TRANSFORM ON THESE NTOT DATA   */
/*     IS T = T0 + T1*NTOT + T2*NTOT*SUM2 + T3*NTOT*SUMF.  FOR THE PAR-    */
/*     TICULAR IMPLEMENTATION FORTRAN 32 ON THE CDC 3300 (FLOATING POINT   */
/*     ADD TIME = SIX MICROSECONDS),                                       */
/*     T = 3000 + 600*NTOT + 50*NTOT*SUM2 + 175*NTOT*SUMF MICROSECONDS     */
/*     ON COMPLEX DATA.                                                    */
/*     IMPLEMENTATION OF THE DEFINITION BY SUMMATION WILL RUN IN A TIME    */
/*     PROPORTIONAL TO NTOT**2.  FOR HIGHLY COMPOSITE NTOT, THE SAVINGS    */
/*     OFFERED BY COOLEY-TUKEY CAN BE DRAMATIC.  A MATRIX 100 BY 100 WILL  */
/*     BE TRANSFORMED IN TIME PROPORTIONAL TO 10000*(2+2+2+2+5+5+5+5) =    */
/*     280,000 (ASSUMING T2 AND T3 TO BE ROUGHLY COMPARABLE) VERSUS        */
/*     10000**2 = 100,000,000 FOR THE STRAIGHTFORWARD TECHNIQUE.           */
/*     THE COOLEY-TUKEY ALGORITHM PLACES TWO RESTRICTIONS UPON THE         */
/*     NATURE OF THE DATA BEYOND THE USUAL RESTRICTION THAT                */
/*     THE DATA FROM ONE CYCLE OF A PERIODIC FUNCTION.  THEY ARE--         */
/*     1.  THE NUMBER OF INPUT DATA AND THE NUMBER OF TRANSFORM VALUES     */
/*     MUST BE THE SAME.                                                   */
/*     2. CONSIDERING THE DATA TO BE IN THE TIME DOMAIN,                   */
/*     THEY MUST BE EQUI-SPACED AT INTERVALS OF DT.  FURTHER, THE TRANS-   */
/*     FORM VALUES, CONSIDERED TO BE IN FREQUENCY SPACE, WILL BE EQUI-     */
/*     SPACED FROM 0 TO 2*PI*(NN(I)-1)/(NN(I)*DT) AT INTERVALS OF          */
/*     2*PI/(NN(I)*DT) FOR EACH DIMENSION OF LENGTH NN(I).  OF COURSE,     */
/*     DT NEED NOT BE THE SAME FOR EVERY DIMENSION.                        */

/*     THE CALLING SEQUENCE IS--                                           */
/*     CALL FOURT(DATAR,DATAI,NN,NDIM,IFRWD,ICPLX,WORKR,WORKI)             */

/*     DATAR AND DATAI ARE THE ARRAYS USED TO HOLD THE REAL AND IMAGINARY  */
/*     PARTS OF THE INPUT DATA ON INPUT AND THE TRANSFORM VALUES ON        */
/*     OUTPUT.  THEY ARE FLOATING POINT ARRAYS, MULTIDIMENSIONAL WITH      */
/*     IDENTICAL DIMENSIONALITY AND EXTENT.  THE EXTENT OF EACH DIMENSION  */
/*     IS GIVEN IN THE INTEGER ARRAY NN, OF LENGTH NDIM.  THAT IS,         */
/*     NDIM IS THE DIMENSIONALITY OF THE ARRAYS DATAR AND DATAI.           */
/*     IFRWD IS AN INTEGER USED TO INDICATE THE DIRECTION OF THE FOURIER   */
/*     TRANSFORM.  IT IS NON-ZERO TO INDICATE A FORWARD TRANSFORM          */
/*     (EXPONENTIAL SIGN IS -) AND ZERO TO INDICATE AN INVERSE TRANSFORM   */
/*     (SIGN IS +).  ICPLX IS AN INTEGER TO INDICATE WHETHER THE DATA      */
/*     ARE REAL OR COMPLEX.  IT IS NON-ZERO FOR COMPLEX, ZERO FOR REAL.    */
/*     IF IT IS ZERO (REAL) THE CONTENTS OF ARRAY DATAI WILL BE ASSUMED    */
/*     TO BE ZERO, AND NEED NOT BE EXPLICITLY SET TO ZERO.  AS EXPLAINED   */
/*     ABOVE, THE TRANSFORM RESULTS ARE ALWAYS COMPLEX AND ARE STORED      */
/*     IN DATAR AND DATAI ON RETURN.  WORKR AND WORKI ARE ARRAYS USED      */
/*     FOR WORKING STORAGE.  THEY ARE NOT NECESSARY IF ALL THE DIMENSIONS  */
/*     OF THE DATA ARE POWERS OF TWO.  IN THIS CASE, THE ARRAYS MAY BE     */
/*     REPLACED BY THE NUMBER 0 IN THE CALLING SEQUENCE.  THUS, USE OF     */
/*     POWERS OF TWO CAN FREE A GOOD DEAL OF STORAGE.  IF ANY DIMENSION    */
/*     IS NOT A POWER OF TWO, THESE ARRAYS MUST BE SUPPLIED.  THEY ARE     */
/*     FLOATING POINT, ONE DIMENSIONAL OF LENGTH EQUAL TO THE LARGEST      */
/*     ARRAY DIMENSION, THAT IS, TO THE LARGEST VALUE OF NN(I).            */
/*     WORKR AND WORKI, IF SUPPLIED, MUST NOT BE THE SAME ARRAYS AS DATAR  */
/*     OR DATAI.  ALL SUBSCRIPTS OF ALL ARRAYS BEGIN AT 1.                 */

/*     THERE ARE NO ERROR MESSAGES OR ERROR HALTS IN THIS PROGRAM.  THE    */
/*     PROGRAM RETURNS IMMEDIATELY IF NDIM OR ANY NN(I) IS LESS THAN ONE.  */

/*     PROGRAM MODIFIED FROM A SUBROUTINE OF BRENNER                       */
/*     10-06-2000, MLR                                                     */


void fourt(double *datar,double *datai, int nn[3], int ndim, int ifrwd, int icplx, double *workr,double *worki)
{
  int ifact[21],ntot,idim,np1,n,np2,m,ntwo,iff,idiv,iquot,irem,inon2,non2p,np0,nprev,icase,ifmin,i,j,jmax,np2hf,i2,i1max,i3,j3,i1,ifp1,ifp2,i2max,i1rng,istep,imin,imax,mmax,mmin,mstep,j1,j2max,j2,jmin,j3max,nhalf;
  double theta,wstpr,wstpi,wminr,wmini,wr,wi,wtemp,thetm,wmstr,wmsti,twowr,sr,si,oldsr,oldsi,stmpr,stmpi,tempr,tempi,difi,difr,sumr,sumi,TWOPI = 6.283185307179586476925286766559;

    ntot = 1;
    for (idim = 0; idim < ndim; idim++) {
      ntot *= nn[idim];
    }

    /*main loop for each dimension*/
    np1 = 1;
    for (idim = 1; idim <= ndim; idim++) {
      n = nn[idim-1];
      np2 = np1*n;

      if (n < 1) {
	goto L920;
      } else if (n == 1) {
	goto L900;
      } 

      /*is n a power of 2 and if not, what are its factors*/
      m = n;
      ntwo = np1;
      iff = 1;
      idiv = 2;

    L10:
      iquot = m/idiv;
      irem = m-idiv*iquot;
      if (iquot < idiv) 
	goto L50;
      if (irem == 0) {
	ntwo *= 2;
	ifact[iff] = idiv;
	iff++;
	m= iquot;
	goto L10;
      } 
      idiv = 3;
      inon2 = iff;

    L30:
      iquot = m/idiv;
      irem = m-idiv*iquot;
      if (iquot < idiv) 
	goto L60;
      if (irem == 0) {
	ifact[iff] = idiv;
	iff++;
	m = iquot;
	goto L30;
      }
      
      idiv += 2;
      goto L30;
      
    L50:
      inon2 = iff;
      if (irem != 0) 
	goto L60;
      ntwo *= 2;
      goto L70;
      
    L60:
      ifact[iff] = m;
      
    L70:
      non2p = np2/ntwo;
      

	/*SEPARATE FOUR CASES-- 
	  1. COMPLEX TRANSFORM 
	  2. REAL TRANSFORM FOR THE 2ND, 3RD, ETC. DIMENSION.  METHOD: TRANSFORM HALF THE DATA, SUPPLYING THE OTHER HALF BY CONJUGATE SYMMETRY.
	  3. REAL TRANSFORM FOR THE 1ST DIMENSION, N ODD.  METHOD: SET THE IMAGINARY PARTS TO ZERO.
	  4. REAL TRANSFORM FOR THE 1ST DIMENSION, N EVEN.  METHOD: TRANSFORM A COMPLEX ARRAY OF LENGTH N/2 WHOSE REAL PARTS ARE THE EVEN NUMBERED REAL VALUES AND WHOSE IMAGINARY PARTS ARE THE ODD-NUMBERED REAL VALUES.  UNSCRAMBLE AND SUPPLY THE SECOND HALF BY CONJUGATE SYMMETRY. */
	
      icase = 1;
      ifmin = 1;
      if (icplx != 0) 
	goto L100;
      icase = 2;
      if (idim > 1) 
	goto L100;
      icase = 3;
      if (ntwo <= np1) 
	goto L100;
      icase = 4;
      ifmin = 2;
      ntwo /= 2;
      n /= 2;
      np2 /= 2;
      ntot /= 2;
      i = 1;
      for (j = 1; j <= ntot; j++) {
	datar[j] = datar[i];
	datai[j] = datar[i+1];
	i += 2;
      }

	/*shuffle data by bit reversal, since n = 2^k. As the shuffling can be done by simple interchange, no working array is needed*/
    L100:
      if (non2p > 1) 
	goto L200;
      np2hf = np2/2;
      j = 1;
      for (i2 = 1; i2 <= np2; i2 += np1) {
	if (j >= i2) 
	  goto L130;
	i1max = i2+np1-1;
	for (i1 = i2; i1 <= i1max; i1++) {
	  for (i3 = i1; i3 <= ntot; i3 += np2) {
	    j3 = j+i3-i2;
	    tempr = datar[i3];
	    tempi = datai[i3];
	    datar[i3] = datar[j3];
	    datai[i3] = datai[j3];
	    datar[j3] = tempr;
	    datai[j3] = tempi;
	  }
	}
	
      L130:
	m = np2hf;
	
      L140:
	if (j <= m) {
	  j += m;
	} else {
	  j -= m;
	  m /= 2;
	  if (m >= np1) 
	    goto L140;
	}
	
      }
      goto L300;

	/*shuffle data by digit reversal for general n*/
    L200:
      for (i1 = 1; i1 <= np1; i1++) {
	for (i3 = i1; i3 <= ntot; i3 += np2) {
	  j = i3;
	  for (i = 1; i <= n; i++) {
	    if (icase != 3) {
	      workr[i] = datar[j];
	      worki[i] = datai[j];
	    } else {
	      workr[i] = datar[j];
	      worki[i] = 0.;
	    }
	    ifp2 = np2;
	    iff = ifmin;
	  L250:
	    ifp1 = ifp2/ifact[iff];
	    j += ifp1;
	    if (j >= i3+ifp2) { 
	      j -= ifp2;
	      ifp2 = ifp1;
	      iff += 1;
	      if (ifp2 > np1) 
		goto L250;
	    }
	  }
	  i2max = i3+np2-np1;
	  i = 1;
	  for (i2 = i3; i2 <= i2max; i2 += np1) {
	    datar[i2] = workr[i];
	    datai[i2] = worki[i];
	    i++;
	  }
	  
	}
      }
      
      /*special case-- W=1*/
    L300:
      i1rng = np1;
      if (icase == 2) 
	i1rng = np0*(1+nprev/2);
      if (ntwo <= np1) 
	goto L600;
      for (i1 = 1; i1 <= i1rng; i1++) {
	imin = np1+i1;
	istep = 2*np1;
	goto L330;
	
      L310:
	j = i1;
	for (i = imin; i <= ntot; i += istep) {
	  tempr = datar[i];
	  tempi = datai[i];
	  datar[i] = datar[j]-tempr;
	  datai[i] = datai[j]-tempi;
	  datar[j] = datar[j]+tempr;
	  datai[j] = datai[j]+tempi;
	  j += istep;
	}
	imin = 2*imin-i1;
	istep *= 2;
	
      L330:
	if (istep <= ntwo) 
	  goto L310;
	
	/*special case-- W = -sqrt(-1)*/
	imin = 3*np1+i1;
	istep = 4*np1;
	goto L420;
	
      L400:
	j = imin-istep/2;
	for (i = imin; i <= ntot; i += istep) {
	  if (ifrwd != 0) {
	    tempr = datai[i];
	    tempi = -datar[i];
	  } else {
	    tempr = -datai[i];
	    tempi = datar[i];
	  }
	  datar[i] = datar[j]-tempr;
	  datai[i] = datai[j]-tempi;
	  datar[j] += tempr;
	  datai[j] += tempi;
	  j += istep;
	}
	
	imin = 2*imin-i1;
	istep *= 2;
	
      L420:
	if (istep <= ntwo) 
	  goto L400;
      }
      
      /*main loop for factors of 2. W=EXP(-2*PI*SQRT(-1)*m/mmax) */
      theta = -TWOPI/8.;
      wstpr = 0.;
      wstpi = -1.;
      if (ifrwd == 0) {
	theta = -theta;
	wstpi = 1.;
      }
      mmax = 8*np1;
      goto L540;
      
    L500:
      wminr = cos(theta);
      wmini = sin(theta);
      wr = wminr;
      wi = wmini;
      mmin = mmax/2+np1;
      mstep = np1*2;
      for (m = mmin; m <= mmax; m += mstep) {
	for (i1 = 1; i1 <= i1rng; i1++) {
	  istep = mmax;
	  imin = m+i1;
	L510:
	  j = imin-istep/2;
	  for (i = imin; i <= ntot; i += istep) {
	    tempr = datar[i]*wr-datai[i]*wi;
	    tempi = datar[i]*wi+datai[i]*wr;
	    datar[i] = datar[j]-tempr;
	    datai[i] = datai[j]-tempi;
	    datar[j] += tempr;
	    datai[j] += tempi;
	    j += istep;
	  }
	  imin = 2*imin-i1;
	  istep *= 2;
	  if (istep <= ntwo) 
	    goto L510;
	}
	wtemp = wr*wstpi;
	wr = wr*wstpr-wi*wstpi;
	wi = wi*wstpr+wtemp;
      }
      wstpr = wminr;
      wstpi = wmini;
      theta /= 2.;
      mmax += mmax;
    L540:
      if (mmax <= ntwo) 
	goto L500;
      
      /*main loop for factors not equal to 2-- W=EXP(-2*PI*SQRT(-1)*(j2-i3)/ifp2)*/
    L600:
      if (non2p <= 1) 
	goto L700;
      ifp1 = ntwo;
      iff = inon2;
    L610:
      ifp2 = ifact[iff]*ifp1;
      theta = -TWOPI/ (double)ifact[iff]; 
      if (ifrwd == 0) 
	theta = -theta;
      thetm = theta/ (double)(ifp1/np1);
      wstpr = cos(theta);
      wstpi = sin(theta);
      wmstr = cos(thetm);
      wmsti = sin(thetm);
      wminr = 1.;
      wmini = 0.;
      
      for (j1 = 1; j1 <= ifp1; j1 += np1) {
	i1max = j1+i1rng-1;
	for (i1 = j1; i1 <= i1max; i1++) {
	  for (i3 = i1; i3 <= ntot; i3 += np2) {
	    i = 1;
	    wr = wminr;
	    wi = wmini;
	    j2max = i3+ifp2-ifp1;
	    for (j2 = i3; j2 <= j2max; j2 += ifp1) {
	      twowr = 2.*wr;
	      jmin = i3;
	      j3max = j2+np2-ifp2;
	      for (j3 = j2 ; j3 <= j3max; j3 += ifp2) {
		j = jmin+ifp2-ifp1;
		sr = datar[j];
		si = datai[j];
		oldsr = 0.;
		oldsi = 0.;
		j -= ifp1;
	      L620:
		stmpr = sr;
		stmpi = si;
		sr = twowr*sr-oldsr+datar[j];
		si = twowr*si-oldsi+datai[j];
		oldsr = stmpr;
		oldsi = stmpi;
		j -= ifp1;
		if (j > jmin) 
		  goto L620;
		workr[i] = wr*sr-wi*si-oldsr+datar[j];
		worki[i] = wi*sr+wr*si-oldsi+datai[j];
		jmin += ifp2;
		i++;
	      }
	      wtemp = wr*wstpi;
	      wr = wr*wstpr-wi*wstpi;
	      wi = wi*wstpr+wtemp;
	    }
	    i = 1;
	    for (j2 = i3; j2 <= j2max; j2 += ifp1) {
	      j3max = j2+np2-ifp2;
	      for (j3 = j2; j3 <= j3max; j3 += ifp2) {
		datar[j3] = workr[i];
		datai[j3] = worki[i];
		i++;
	      }
	    }
	  }
	}
	wtemp = wminr*wmsti;
	wminr = wminr*wmstr-wmini*wmsti;
	wmini = wmini*wmstr+wtemp;
      }
      iff++;
      ifp1 = ifp2;
      if (ifp1 < np2) 
	goto L610;
      
      /*complete a real transform in the 1st dimension, n even, by conjugate symmetries*/
    L700:
      switch (icase) {
      case 1:  
	goto L900;
	break;
      case 2:  
	goto L800;
	break;
      case 3:  
	goto L900;
	break;
      }

      nhalf = n;
      n += n;
      theta = -TWOPI/ (double) n;
      if (ifrwd == 0) 
	theta = -theta;
      wstpr = cos(theta);
      wstpi = sin(theta);
      wr = wstpr;
      wi = wstpi;
      imin = 2;
      jmin = nhalf;
      goto L725;
    L710:
      j = jmin;
      for (i = imin; i <= ntot; i += np2) {
	sumr = (datar[i]+datar[j])/2.;
	sumi = (datai[i]+datai[j])/2.;
	difr = (datar[i]-datar[j])/2.;
	difi = (datai[i]-datai[j])/2.;
	tempr = wr*sumi+wi*difr;
	tempi = wi*sumi-wr*difr;
	datar[i] = sumr+tempr;
	datai[i] = difi+tempi;
	datar[j] = sumr-tempr;
	datai[j] = tempi-difi;
	j += np2;
      }
      imin++;
      jmin--;
      wtemp = wr*wstpi;
      wr = wr*wstpr-wi*wstpi;
      wi = wi*wstpr+wtemp;
    L725:
      if (imin < jmin) {
	goto L710;
      } else if (imin > jmin) {
	goto L740;
      }
      if (ifrwd == 0) 
	goto L740;
      for (i = imin; i <= ntot; i += np2) {
	datai[i] = -datai[i];
      }
    L740:
      np2 *= 2;
      ntot *= 2;
      j = ntot+1;
      imax = ntot/2+1;
    L745:
      imin = imax-nhalf;
      i = imin;
      goto L755;
    L750:
      datar[j] = datar[i];
      datai[j] = -datai[i];
    L755:
      i++;
      j--;
      if (i < imax) 
	goto L750;
      datar[j] = datar[imin]-datai[imin];
      datai[j] = 0.;
      if (i >= j) {
	goto L780;
      } else {
	goto L770;
      }
    L765:
      datar[j] = datar[i];
      datai[j] = datai[i];
    L770:
      i--;
      j--;
      if (i > imin) 
	goto L765;
      datar[j] = datar[imin]+datai[imin];
      datai[j] = 0.;
      imax = imin;
      goto L745;
    L780:
      datar[1] += datai[1];
      datai[1] = 0.;
      goto L900;
      
      /*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/
    L800:
      if (nprev <= 2) 
	goto L900;
      for (i3 = 1; i3 <= ntot; i3 += np2) {
	i2max = i3+np2-np1;
	for (i2 = i3; i2 <= i2max; i2 += np1) {
	  imax = i2+np1-1;
	  imin = i2+i1rng;
	  jmax = 2*i3+np1-imin;
	  if (i2 > i3) 
	    jmax += np2;
	  if (idim > 2) {
	    j = jmax+np0;
	    for (i = imin; i <= imax; i++) {
	      datar[i] = datar[j];
	      datai[i] = -datai[j];
	      j--;
	    }
	  }
	  j = jmax;
	  for (i = imin; i <= imax; i += np0) {
	    datar[i] = datar[j];
	    datai[i] = -datai[j];
	    j -= np0;
	  }
	}
      }

      /*end of loop on each dimension*/
    L900:
      np0 = np1;
      np1 = np2;
      nprev = n;
    }
 L920:
    return;
}