Index: /trunk/Ohana/src/opihi/cmd.data/fft1d.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/fft1d.c	(revision 16093)
+++ /trunk/Ohana/src/opihi/cmd.data/fft1d.c	(revision 16094)
@@ -3,6 +3,5 @@
 int fft1d (int argc, char **argv) {
   
-  int i, Npix, ZeroImaginary;
-  float *t1, *t2, *temp;
+  int Npix, Nbit, ZeroImaginary;
   Vector *Ire, *Iim, *Ore, *Oim;
 
@@ -19,7 +18,6 @@
   }    
   if ((Ire = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((Ore = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((Oim = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE);
-
+  
+  // check the input data (match lengths? binary length?)
   Npix = Ire[0].Nelements;
   if (!ZeroImaginary && (Npix != Iim[0].Nelements)) {
@@ -27,28 +25,14 @@
     return (FALSE);
   }
-
-  if (!IsBinary (Npix)) {
+  if (!IsBinary (Npix, &Nbit)) {
     gprint (GP_ERR, "Npix is not a binary number!\n");
     return (FALSE);
   }
-  
-  ALLOCATE (temp, float, 2*Npix);
-  if (ZeroImaginary) {
-    t1 = Ire[0].elements;
-    for (i = 0; i < Npix; i++, t1++) {
-      temp[2*i  ] = *t1;
-      temp[2*i+1] = 0;
-    }
-  } else {
-    t1 = Ire[0].elements;
-    t2 = Iim[0].elements;
-    for (i = 0; i < Npix; i++, t1++, t2++) {
-      temp[2*i  ] = *t1;
-      temp[2*i+1] = *t2;
-    }
-  }    
-    
-  fft (temp, Npix, 1); 
 
+  // select or create the output vectors
+  if ((Ore = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Oim = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  // allocate sufficient output space
   Ore[0].Nelements = Npix;
   Oim[0].Nelements = Npix;
@@ -56,15 +40,16 @@
   REALLOCATE (Oim[0].elements, float, Npix);
  
-  t1 = Ore[0].elements;
-  t2 = Oim[0].elements;
-  for (i = 0; i < Npix; i++, t1++, t2++) {
-    *t1 = temp[2*i  ] / Npix;
-    *t2 = temp[2*i+1] / Npix;
+  // copy data to output vectors (fft is done in place)
+  memcpy (Ore[0].elements, Ire[0].elements, Npix*sizeof(float));
+
+  // copy imaginary vector or create a zero vector
+  if (ZeroImaginary) {
+    memset (Oim[0].elements, 0, Npix*sizeof(float));
+  } else {
+    memcpy (Oim[0].elements, Iim[0].elements, Npix*sizeof(float));
   }    
-  
-  free (temp);
+
+  fft (Ore[0].elements, Oim[0].elements, Npix, Nbit); 
   
   return (TRUE);
 }
-
-  
Index: /trunk/Ohana/src/opihi/include/data.h
===================================================================
--- /trunk/Ohana/src/opihi/include/data.h	(revision 16093)
+++ /trunk/Ohana/src/opihi/include/data.h	(revision 16094)
@@ -84,6 +84,5 @@
 
 /* in fft.c */
-void fft (float *Data, int N, int isign);
-void fftold (float *Data, int N, int isign);
+void fft (float *dataRe, float *dataIm, int N, int Nbit);
 void fftN (float *data, int *nn, int ndim, int isign);
 int IsBinary (int N);
Index: /trunk/Ohana/src/opihi/lib.data/fft.c
===================================================================
--- /trunk/Ohana/src/opihi/lib.data/fft.c	(revision 16093)
+++ /trunk/Ohana/src/opihi/lib.data/fft.c	(revision 16094)
@@ -1,270 +1,101 @@
 # include "data.h"
 
-#define FSWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
-
-void fft (float *Data, int N, int isign) {
-
-  int n,mmax,m,j,istep,i;
-  double wtemp,wr,wpr,wpi,wi,theta;
-  float tempr,tempi, *data;
-
-  data = Data;
-  n = N << 1;
-  j = 0;
-
-  for (i = 0; i < n; i+=2) {
-    if (j > i) {
-      FSWAP (data[j], data[i]);
-      FSWAP (data[j+1], data[i+1]);
-    }
-    m = n >> 1;
-    while (m >= 2 && j >= m) {
-      j -= m;
-      m >>= 1;
-    }
-    j += m;
-  }
-  mmax = 2;
-  while (n > mmax) {
-    istep = 2*mmax;
-    theta = 6.28318530717959 / (isign*mmax);
-    wtemp = sin(0.5*theta);
-    wpr = -2.0*wtemp*wtemp;
-    wpi = sin(theta);
-    wr = 1.0;
-    wi = 0.0;
-    for (m = 0; m < mmax; m+=2) {
-      for (i = m; i < n; i+=istep) {
-	j = i + mmax;
-	tempr = wr*data[j] - wi*data[j+1];
-	tempi = wr*data[j+1] + wi*data[j];
-	data[j] = data[i] - tempr;
-	data[j+1] = data[i+1] - tempi;
-	data[i] += tempr;
-	data[i+1] += tempi;
-      }
-      wr = (wtemp = wr)*wpr - wi*wpi+wr;
-      wi = wi*wpr + wtemp*wpi + wi;
-    }
-    mmax = istep;
-  }
+// fft based on code by Douglas L. Jones (see note at EOF). modified for Ohana C style
+void fftN (float *data, int *nn, int ndim, int isign) {
+  
+  fprintf (stderr, "broken\n");
+  return;
 }
 
-void fftold (float *Data, int N, int isign) {
+void fft (double *x, double *y, int n, int Nbit) {
 
-  int n,mmax,m,j,istep,i;
-  double wtemp,wr,wpr,wpi,wi,theta;
-  float tempr,tempi, *data;
+  int i,j,k,n1,n2;
+  double c,s,e,a,t1,t2;        
+         
+  // bit-reverse
+  j = 0; 
+  n2 = n/2;
+  for (i = 1; i < n - 1; i++) {
+    n1 = n2;
+    while ( j >= n1 ) {
+      j = j - n1;
+      n1 = n1/2;
+    }
+    j = j + n1;
+               
+    if (i < j) {
+      t1 = x[i];
+      x[i] = x[j];
+      x[j] = t1;
+      t1 = y[i];
+      y[i] = y[j];
+      y[j] = t1;
+    }
+  }
+  
+                                          
+  n1 = 0; /* FFT */
+  n2 = 1;
+                                             
+  for (i=0; i < Nbit; i++) {
+    n1 = n2;
+    n2 = n2 + n2;
+    e = -6.283185307179586/n2;
+    a = 0.0;
+                                             
+    for (j=0; j < n1; j++) {
+      c = cos(a);
+      s = sin(a);
+      a = a + e;
+                                            
+      for (k=j; k < n; k=k+n2) {
+	t1 = c*x[k+n1] - s*y[k+n1];
+	t2 = s*x[k+n1] + c*y[k+n1];
+	x[k+n1] = x[k] - t1;
+	y[k+n1] = y[k] - t2;
+	x[k] = x[k] + t1;
+	y[k] = y[k] + t2;
+      }
+    }
+  }
+                                      
+  return;
+}                          
 
-  data = Data - 1;
-  n = N << 1;
-  j = 1;
+int IsBinary (int N, int *nbin) {
 
-  for (i = 1; i < n; i+=2) {
-    if (j > i) {
-      FSWAP (data[j], data[i]);
-      FSWAP (data[j+1], data[i+1]);
+  // check that a number is binary (2^nbit)
+
+  Nbit = 0;
+  Nset = 0;
+  for (i = 0; i < 8*sizeof(int); i++) {
+    if (N & 0x01) {
+      Nset ++;
+      Nbin = i;
     }
-    m = n >> 1;
-    while (m >= 2 && j > m) {
-      j -= m;
-      m >>= 1;
-    }
-    j += m;
+    N >>= 1;
   }
-  mmax = 2;
-  while (n > mmax) {
-    istep = 2*mmax;
-    theta = 6.28318530717959 / (isign*mmax);
-    wtemp = sin(0.5*theta);
-    wpr = -2.0*wtemp*wtemp;
-    wpi = sin(theta);
-    wr = 1.0;
-    wi = 0.0;
-    for (m = 1; m < mmax; m+=2) {
-      for (i = m; i <= n; i+=istep) {
-	j = i + mmax;
-	tempr = wr*data[j] - wi*data[j+1];
-	tempi = wr*data[j+1] + wi*data[j];
-	data[j] = data[i] - tempr;
-	data[j+1] = data[i+1] - tempi;
-	data[i] += tempr;
-	data[i+1] += tempi;
-      }
-      wr = (wtemp = wr)*wpr - wi*wpi+wr;
-      wi = wi*wpr + wtemp*wpi + wi;
-    }
-    mmax = istep;
-  }
+  *nbin = Nbit;
+  if (Nset > 1) return (FALSE);
 }
 
-/* convert indices to zero reference */
-void fftN (float *data, int *nn, int ndim, int isign) {
+/**********************************************************/
+/* fft.c                                                  */
+/* (c) Douglas L. Jones                                   */
+/* University of Illinois at Urbana-Champaign             */
+/* January 19, 1992                                       */
+/*                                                        */
+/*   fft: in-place radix-2 DIT DFT of a complex input     */
+/*                                                        */
+/*   input:                                               */
+/* n: length of FFT: must be a power of two               */
+/* m: n = 2**m                                            */
+/*   input/output                                         */
+/* x: double array of length n with real part of data     */
+/* y: double array of length n with imag part of data     */
+/*                                                        */
+/*   Permission to copy and use this program is granted   */
+/*   under a Creative Commons "Attribution" license       */
+/*   http://creativecommons.org/licenses/by/1.0/          */
+/**********************************************************/
 
-  int i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
-  int ibit,idim,k1,k2,n,nprev,nrem,ntot;
-  float tempi,tempr;
-  double theta,wi,wpi,wpr,wr,wtemp;
-
-  ntot = 1;
-  for (idim = 0; idim < ndim; idim++) ntot *= nn[idim];
-
-  nprev = 1;
-  for (idim = ndim - 1; idim >= 0; idim--) {
-    n  =  nn[idim];
-    nrem = ntot / (n*nprev);
-    ip1 = nprev << 1;
-    ip2 = ip1*n;
-    ip3 = ip2*nrem;
-    i2rev = 0;
-    for (i2 = 0; i2 < ip2; i2+=ip1) {
-      if (i2 < i2rev) {
-	for (i1 = i2; i1 <= i2+ip1-2; i1+=2) {
-	  for (i3 = i1; i3 < ip3; i3+=ip2) {	
-	    i3rev = i2rev+i3-i2;
-	    FSWAP(data[i3],data[i3rev]);
-	    FSWAP(data[i3+1],data[i3rev+1]);
-	  }
-	}
-      }
-      ibit = ip2 >> 1;
-      while (ibit >= ip1 && i2rev >= ibit) {
-	i2rev -= ibit;
-	ibit >>= 1;
-      }
-      i2rev += ibit;
-    }
-    ifp1 = ip1;
-    while (ifp1 < ip2) {
-      ifp2 = ifp1 << 1;
-      theta = isign*6.28318530717959/(ifp2/ip1);
-      wtemp = sin(0.5*theta);
-      wpr = -2.0*wtemp*wtemp;
-      wpi = sin(theta);
-      wr = 1.0;
-      wi = 0.0;
-      for (i3 = 0; i3 < ifp1; i3+=ip1) {
-	for (i1 = i3; i1 <= i3+ip1-2; i1+=2) {
-	  for (i2 = i1; i2 < ip3; i2+=ifp2) {
-	    k1 = i2;
-	    k2 = k1+ifp1;
-	    tempr = wr*data[k2]-wi*data[k2+1];
-	    tempi = wr*data[k2+1]+wi*data[k2];
-	    data[k2] = data[k1]-tempr;
-	    data[k2+1] = data[k1+1]-tempi;
-	    data[k1] += tempr;
-	    data[k1+1] += tempi;
-	  }
-	}
-	wr = (wtemp = wr)*wpr-wi*wpi+wr;
-	wi = wi*wpr+wtemp*wpi+wi;
-      }
-      ifp1 = ifp2;
-    }
-    nprev *= n;
-  }
-}
-
-#undef FSWAP
-
-/* based on the PRESS routine, this fft takes an array from data[0] to data[2N-1] */
-/* this function takes Data = h(t) and replaces it in situ with H(F) or vice versa.
-   There are assumed to be 2*N input values with 
-   Data[0,2,4,...] the real and Data[1,3,5,...] the imaginary ones.
-   the output is ordered the same.  
-
-   for h(t), values are in time sequence order.
-   for H(F), values are in order F = 0, 1/N, ... 1/2 - 1/N, +/- 1/2, -1/2 + 1/N, ... -1/N 
-
-   no normalization is performed, so a signal of amplitude A sin (w_k * t) will be 
-      give an H(F) value of 0.5 * A * N, and the DC term will have H(0) = A * N.
-
-*/ 
-
-int IsBinary (int N) {
-
-  int i, nbit;
-
-  /* check if number is a binary number */
-  nbit = 0;
-  for (i = 0; i < 8*sizeof(N); i++) {
-    nbit += (N & 0x01);
-    N = (N >> 1);
-  }
-  if (nbit == 1) { 
-    return (1); 
-  } else { 
-    return (0); 
-  }
-
-}  
-
-#define FSWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
-
-void fourn (float *data, int *nn, int ndim, int isign) {
-
-  int i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
-  int ibit,idim,k1,k2,n,nprev,nrem,ntot;
-  float tempi,tempr;
-  double theta,wi,wpi,wpr,wr,wtemp;
-
-  ntot = 1;
-  for (idim = 1; idim <= ndim; idim++)
-    ntot *= nn[idim];
-  nprev = 1;
-  for (idim = ndim; idim >= 1; idim--) {
-    n = nn[idim];
-    nrem = ntot/(n*nprev);
-    ip1 = nprev << 1;
-    ip2 = ip1*n;
-    ip3 = ip2*nrem;
-    i2rev = 1;
-    for (i2 = 1; i2 <= ip2; i2+=ip1) {
-      if (i2 < i2rev) {
-	for (i1 = i2; i1 <= i2+ip1-2; i1+=2) {
-	  for (i3 = i1; i3 <= ip3;i3+=ip2) {
-	    i3rev = i2rev+i3-i2;
-	    FSWAP(data[i3],data[i3rev]);
-	    FSWAP(data[i3+1],data[i3rev+1]);
-	  }
-	}
-      }
-      ibit = ip2 >> 1;
-      while (ibit >= ip1 && i2rev > ibit) {
-	i2rev -=  ibit;
-	ibit >>= 1;
-      }
-      i2rev += ibit;
-    }
-    ifp1 = ip1;
-    while (ifp1 < ip2) {
-      ifp2 = ifp1 << 1;
-      theta = isign*6.28318530717959/(ifp2/ip1);
-      wtemp = sin(0.5*theta);
-      wpr  =  -2.0*wtemp*wtemp;
-      wpi = sin(theta);
-      wr = 1.0;
-      wi = 0.0;
-      for (i3 = 1;i3<=ifp1;i3+=ip1) {
-	for (i1 = i3;i1<=i3+ip1-2;i1+=2) {
-	  for (i2 = i1;i2<=ip3;i2+=ifp2) {
-	    k1 = i2;
-	    k2 = k1+ifp1;
-	    tempr = wr*data[k2]-wi*data[k2+1];
-	    tempi = wr*data[k2+1]+wi*data[k2];
-	    data[k2] = data[k1]-tempr;
-	    data[k2+1] = data[k1+1]-tempi;
-	    data[k1] += tempr;
-	    data[k1+1] += tempi;
-	  }
-	}
-	wr = (wtemp = wr)*wpr-wi*wpi+wr;
-	wi = wi*wpr+wtemp*wpi+wi;
-      }
-      ifp1 = ifp2;
-    }
-    nprev *= n;
-  }
-}
-
-#undef FSWAP
Index: /trunk/Ohana/src/opihi/lib.data/fftold.c
===================================================================
--- /trunk/Ohana/src/opihi/lib.data/fftold.c	(revision 16094)
+++ /trunk/Ohana/src/opihi/lib.data/fftold.c	(revision 16094)
@@ -0,0 +1,270 @@
+# include "data.h"
+
+#define FSWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
+
+void fft (float *Data, int N, int isign) {
+
+  int n,mmax,m,j,istep,i;
+  double wtemp,wr,wpr,wpi,wi,theta;
+  float tempr,tempi, *data;
+
+  data = Data;
+  n = N << 1;
+  j = 0;
+
+  for (i = 0; i < n; i+=2) {
+    if (j > i) {
+      FSWAP (data[j], data[i]);
+      FSWAP (data[j+1], data[i+1]);
+    }
+    m = n >> 1;
+    while (m >= 2 && j >= m) {
+      j -= m;
+      m >>= 1;
+    }
+    j += m;
+  }
+  mmax = 2;
+  while (n > mmax) {
+    istep = 2*mmax;
+    theta = 6.28318530717959 / (isign*mmax);
+    wtemp = sin(0.5*theta);
+    wpr = -2.0*wtemp*wtemp;
+    wpi = sin(theta);
+    wr = 1.0;
+    wi = 0.0;
+    for (m = 0; m < mmax; m+=2) {
+      for (i = m; i < n; i+=istep) {
+	j = i + mmax;
+	tempr = wr*data[j] - wi*data[j+1];
+	tempi = wr*data[j+1] + wi*data[j];
+	data[j] = data[i] - tempr;
+	data[j+1] = data[i+1] - tempi;
+	data[i] += tempr;
+	data[i+1] += tempi;
+      }
+      wr = (wtemp = wr)*wpr - wi*wpi+wr;
+      wi = wi*wpr + wtemp*wpi + wi;
+    }
+    mmax = istep;
+  }
+}
+
+void fftold (float *Data, int N, int isign) {
+
+  int n,mmax,m,j,istep,i;
+  double wtemp,wr,wpr,wpi,wi,theta;
+  float tempr,tempi, *data;
+
+  data = Data - 1;
+  n = N << 1;
+  j = 1;
+
+  for (i = 1; i < n; i+=2) {
+    if (j > i) {
+      FSWAP (data[j], data[i]);
+      FSWAP (data[j+1], data[i+1]);
+    }
+    m = n >> 1;
+    while (m >= 2 && j > m) {
+      j -= m;
+      m >>= 1;
+    }
+    j += m;
+  }
+  mmax = 2;
+  while (n > mmax) {
+    istep = 2*mmax;
+    theta = 6.28318530717959 / (isign*mmax);
+    wtemp = sin(0.5*theta);
+    wpr = -2.0*wtemp*wtemp;
+    wpi = sin(theta);
+    wr = 1.0;
+    wi = 0.0;
+    for (m = 1; m < mmax; m+=2) {
+      for (i = m; i <= n; i+=istep) {
+	j = i + mmax;
+	tempr = wr*data[j] - wi*data[j+1];
+	tempi = wr*data[j+1] + wi*data[j];
+	data[j] = data[i] - tempr;
+	data[j+1] = data[i+1] - tempi;
+	data[i] += tempr;
+	data[i+1] += tempi;
+      }
+      wr = (wtemp = wr)*wpr - wi*wpi+wr;
+      wi = wi*wpr + wtemp*wpi + wi;
+    }
+    mmax = istep;
+  }
+}
+
+/* convert indices to zero reference */
+void fftN (float *data, int *nn, int ndim, int isign) {
+
+  int i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
+  int ibit,idim,k1,k2,n,nprev,nrem,ntot;
+  float tempi,tempr;
+  double theta,wi,wpi,wpr,wr,wtemp;
+
+  ntot = 1;
+  for (idim = 0; idim < ndim; idim++) ntot *= nn[idim];
+
+  nprev = 1;
+  for (idim = ndim - 1; idim >= 0; idim--) {
+    n  =  nn[idim];
+    nrem = ntot / (n*nprev);
+    ip1 = nprev << 1;
+    ip2 = ip1*n;
+    ip3 = ip2*nrem;
+    i2rev = 0;
+    for (i2 = 0; i2 < ip2; i2+=ip1) {
+      if (i2 < i2rev) {
+	for (i1 = i2; i1 <= i2+ip1-2; i1+=2) {
+	  for (i3 = i1; i3 < ip3; i3+=ip2) {	
+	    i3rev = i2rev+i3-i2;
+	    FSWAP(data[i3],data[i3rev]);
+	    FSWAP(data[i3+1],data[i3rev+1]);
+	  }
+	}
+      }
+      ibit = ip2 >> 1;
+      while (ibit >= ip1 && i2rev >= ibit) {
+	i2rev -= ibit;
+	ibit >>= 1;
+      }
+      i2rev += ibit;
+    }
+    ifp1 = ip1;
+    while (ifp1 < ip2) {
+      ifp2 = ifp1 << 1;
+      theta = isign*6.28318530717959/(ifp2/ip1);
+      wtemp = sin(0.5*theta);
+      wpr = -2.0*wtemp*wtemp;
+      wpi = sin(theta);
+      wr = 1.0;
+      wi = 0.0;
+      for (i3 = 0; i3 < ifp1; i3+=ip1) {
+	for (i1 = i3; i1 <= i3+ip1-2; i1+=2) {
+	  for (i2 = i1; i2 < ip3; i2+=ifp2) {
+	    k1 = i2;
+	    k2 = k1+ifp1;
+	    tempr = wr*data[k2]-wi*data[k2+1];
+	    tempi = wr*data[k2+1]+wi*data[k2];
+	    data[k2] = data[k1]-tempr;
+	    data[k2+1] = data[k1+1]-tempi;
+	    data[k1] += tempr;
+	    data[k1+1] += tempi;
+	  }
+	}
+	wr = (wtemp = wr)*wpr-wi*wpi+wr;
+	wi = wi*wpr+wtemp*wpi+wi;
+      }
+      ifp1 = ifp2;
+    }
+    nprev *= n;
+  }
+}
+
+#undef FSWAP
+
+/* based on the PRESS routine, this fft takes an array from data[0] to data[2N-1] */
+/* this function takes Data = h(t) and replaces it in situ with H(F) or vice versa.
+   There are assumed to be 2*N input values with 
+   Data[0,2,4,...] the real and Data[1,3,5,...] the imaginary ones.
+   the output is ordered the same.  
+
+   for h(t), values are in time sequence order.
+   for H(F), values are in order F = 0, 1/N, ... 1/2 - 1/N, +/- 1/2, -1/2 + 1/N, ... -1/N 
+
+   no normalization is performed, so a signal of amplitude A sin (w_k * t) will be 
+      give an H(F) value of 0.5 * A * N, and the DC term will have H(0) = A * N.
+
+*/ 
+
+int IsBinary (int N) {
+
+  int i, nbit;
+
+  /* check if number is a binary number */
+  nbit = 0;
+  for (i = 0; i < 8*sizeof(N); i++) {
+    nbit += (N & 0x01);
+    N = (N >> 1);
+  }
+  if (nbit == 1) { 
+    return (1); 
+  } else { 
+    return (0); 
+  }
+
+}  
+
+#define FSWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
+
+void fourn (float *data, int *nn, int ndim, int isign) {
+
+  int i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
+  int ibit,idim,k1,k2,n,nprev,nrem,ntot;
+  float tempi,tempr;
+  double theta,wi,wpi,wpr,wr,wtemp;
+
+  ntot = 1;
+  for (idim = 1; idim <= ndim; idim++)
+    ntot *= nn[idim];
+  nprev = 1;
+  for (idim = ndim; idim >= 1; idim--) {
+    n = nn[idim];
+    nrem = ntot/(n*nprev);
+    ip1 = nprev << 1;
+    ip2 = ip1*n;
+    ip3 = ip2*nrem;
+    i2rev = 1;
+    for (i2 = 1; i2 <= ip2; i2+=ip1) {
+      if (i2 < i2rev) {
+	for (i1 = i2; i1 <= i2+ip1-2; i1+=2) {
+	  for (i3 = i1; i3 <= ip3;i3+=ip2) {
+	    i3rev = i2rev+i3-i2;
+	    FSWAP(data[i3],data[i3rev]);
+	    FSWAP(data[i3+1],data[i3rev+1]);
+	  }
+	}
+      }
+      ibit = ip2 >> 1;
+      while (ibit >= ip1 && i2rev > ibit) {
+	i2rev -=  ibit;
+	ibit >>= 1;
+      }
+      i2rev += ibit;
+    }
+    ifp1 = ip1;
+    while (ifp1 < ip2) {
+      ifp2 = ifp1 << 1;
+      theta = isign*6.28318530717959/(ifp2/ip1);
+      wtemp = sin(0.5*theta);
+      wpr  =  -2.0*wtemp*wtemp;
+      wpi = sin(theta);
+      wr = 1.0;
+      wi = 0.0;
+      for (i3 = 1;i3<=ifp1;i3+=ip1) {
+	for (i1 = i3;i1<=i3+ip1-2;i1+=2) {
+	  for (i2 = i1;i2<=ip3;i2+=ifp2) {
+	    k1 = i2;
+	    k2 = k1+ifp1;
+	    tempr = wr*data[k2]-wi*data[k2+1];
+	    tempi = wr*data[k2+1]+wi*data[k2];
+	    data[k2] = data[k1]-tempr;
+	    data[k2+1] = data[k1+1]-tempi;
+	    data[k1] += tempr;
+	    data[k1+1] += tempi;
+	  }
+	}
+	wr = (wtemp = wr)*wpr-wi*wpi+wr;
+	wi = wi*wpr+wtemp*wpi+wi;
+      }
+      ifp1 = ifp2;
+    }
+    nprev *= n;
+  }
+}
+
+#undef FSWAP