Many name-changes in oldnumeric. This may break some numpy code that was using the oldnumeric interface.

10 files changed, 2302 insertions, 0 deletions

diff --git a/numpy/fft/__init__.py b/numpy/fft/__init__.py
new file mode 100644
index 000000000..0c5a7f2ed
--- /dev/null
+++ b/numpy/fft/__init__.py
@@ -0,0 +1,9 @@
+# To get sub-modules
+from info import __doc__
+
+from fftpack import *
+from helper import *
+
+def test(level=1, verbosity=1):
+    from numpy.testing import NumpyTest
+    return NumpyTest().test(level, verbosity)
diff --git a/numpy/fft/fftpack.c b/numpy/fft/fftpack.c
new file mode 100644
index 000000000..3e5e7d2ed
--- /dev/null
+++ b/numpy/fft/fftpack.c
@@ -0,0 +1,1501 @@
+/*
+fftpack.c : A set of FFT routines in C.
+Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber (Version 4, 1985).
+
+*/
+
+/* isign is +1 for backward and -1 for forward transforms */
+
+#include <math.h>
+#include <stdio.h>
+#define DOUBLE
+
+#ifdef DOUBLE
+#define Treal double
+#else
+#define Treal float
+#endif
+
+
+#define ref(u,a) u[a]
+
+#define MAXFAC 13    /* maximum number of factors in factorization of n */
+#define NSPECIAL 4   /* number of factors for which we have special-case routines */
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+/* ----------------------------------------------------------------------
+   passf2, passf3, passf4, passf5, passf. Complex FFT passes fwd and bwd.
+---------------------------------------------------------------------- */
+
+static void passf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[], int isign)
+  /* isign==+1 for backward transform */
+  {
+    int i, k, ah, ac;
+    Treal ti2, tr2;
+    if (ido <= 2) {
+      for (k=0; k<l1; k++) {
+        ah = k*ido;
+        ac = 2*k*ido;
+        ch[ah]              = ref(cc,ac) + ref(cc,ac + ido);
+        ch[ah + ido*l1]     = ref(cc,ac) - ref(cc,ac + ido);
+        ch[ah+1]            = ref(cc,ac+1) + ref(cc,ac + ido + 1);
+        ch[ah + ido*l1 + 1] = ref(cc,ac+1) - ref(cc,ac + ido + 1);
+      }
+    } else {
+      for (k=0; k<l1; k++) {
+        for (i=0; i<ido-1; i+=2) {
+          ah = i + k*ido;
+          ac = i + 2*k*ido;
+          ch[ah]   = ref(cc,ac) + ref(cc,ac + ido);
+          tr2      = ref(cc,ac) - ref(cc,ac + ido);
+          ch[ah+1] = ref(cc,ac+1) + ref(cc,ac + 1 + ido);
+          ti2      = ref(cc,ac+1) - ref(cc,ac + 1 + ido);
+          ch[ah+l1*ido+1] = wa1[i]*ti2 + isign*wa1[i+1]*tr2;
+          ch[ah+l1*ido]   = wa1[i]*tr2 - isign*wa1[i+1]*ti2;
+        }
+      }
+    }
+  } /* passf2 */
+
+
+static void passf3(int ido, int l1, const Treal cc[], Treal ch[],
+      const Treal wa1[], const Treal wa2[], int isign)
+  /* isign==+1 for backward transform */
+  {
+    static const Treal taur = -0.5;
+    static const Treal taui = 0.866025403784439;
+    int i, k, ac, ah;
+    Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
+    if (ido == 2) {
+      for (k=1; k<=l1; k++) {
+        ac = (3*k - 2)*ido;
+        tr2 = ref(cc,ac) + ref(cc,ac + ido);
+        cr2 = ref(cc,ac - ido) + taur*tr2;
+        ah = (k - 1)*ido;
+        ch[ah] = ref(cc,ac - ido) + tr2;
+
+        ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);
+        ci2 = ref(cc,ac - ido + 1) + taur*ti2;
+        ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;
+
+        cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));
+        ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));
+        ch[ah + l1*ido] = cr2 - ci3;
+        ch[ah + 2*l1*ido] = cr2 + ci3;
+        ch[ah + l1*ido + 1] = ci2 + cr3;
+        ch[ah + 2*l1*ido + 1] = ci2 - cr3;
+      }
+    } else {
+      for (k=1; k<=l1; k++) {
+        for (i=0; i<ido-1; i+=2) {
+          ac = i + (3*k - 2)*ido;
+          tr2 = ref(cc,ac) + ref(cc,ac + ido);
+          cr2 = ref(cc,ac - ido) + taur*tr2;
+          ah = i + (k-1)*ido;
+          ch[ah] = ref(cc,ac - ido) + tr2;
+          ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);
+          ci2 = ref(cc,ac - ido + 1) + taur*ti2;
+          ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;
+          cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));
+          ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));
+          dr2 = cr2 - ci3;
+          dr3 = cr2 + ci3;
+          di2 = ci2 + cr3;
+          di3 = ci2 - cr3;
+          ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
+          ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
+          ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
+          ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
+        }
+      }
+    }
+  } /* passf3 */
+
+
+static void passf4(int ido, int l1, const Treal cc[], Treal ch[],
+      const Treal wa1[], const Treal wa2[], const Treal wa3[], int isign)
+  /* isign == -1 for forward transform and +1 for backward transform */
+  {
+    int i, k, ac, ah;
+    Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
+    if (ido == 2) {
+      for (k=0; k<l1; k++) {
+        ac = 4*k*ido + 1;
+        ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);
+        ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);
+        tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);
+        ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);
+        tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);
+        tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);
+        ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);
+        tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);
+        ah = k*ido;
+        ch[ah] = tr2 + tr3;
+        ch[ah + 2*l1*ido] = tr2 - tr3;
+        ch[ah + 1] = ti2 + ti3;
+        ch[ah + 2*l1*ido + 1] = ti2 - ti3;
+        ch[ah + l1*ido] = tr1 + isign*tr4;
+        ch[ah + 3*l1*ido] = tr1 - isign*tr4;
+        ch[ah + l1*ido + 1] = ti1 + isign*ti4;
+        ch[ah + 3*l1*ido + 1] = ti1 - isign*ti4;
+      }
+    } else {
+      for (k=0; k<l1; k++) {
+        for (i=0; i<ido-1; i+=2) {
+          ac = i + 1 + 4*k*ido;
+          ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);
+          ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);
+          ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);
+          tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);
+          tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);
+          tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);
+          ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);
+          tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);
+          ah = i + k*ido;
+          ch[ah] = tr2 + tr3;
+          cr3 = tr2 - tr3;
+          ch[ah + 1] = ti2 + ti3;
+          ci3 = ti2 - ti3;
+          cr2 = tr1 + isign*tr4;
+          cr4 = tr1 - isign*tr4;
+          ci2 = ti1 + isign*ti4;
+          ci4 = ti1 - isign*ti4;
+          ch[ah + l1*ido] = wa1[i]*cr2 - isign*wa1[i + 1]*ci2;
+          ch[ah + l1*ido + 1] = wa1[i]*ci2 + isign*wa1[i + 1]*cr2;
+          ch[ah + 2*l1*ido] = wa2[i]*cr3 - isign*wa2[i + 1]*ci3;
+          ch[ah + 2*l1*ido + 1] = wa2[i]*ci3 + isign*wa2[i + 1]*cr3;
+          ch[ah + 3*l1*ido] = wa3[i]*cr4 -isign*wa3[i + 1]*ci4;
+          ch[ah + 3*l1*ido + 1] = wa3[i]*ci4 + isign*wa3[i + 1]*cr4;
+        }
+      }
+    }
+  } /* passf4 */
+
+
+static void passf5(int ido, int l1, const Treal cc[], Treal ch[],
+      const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[], int isign)
+  /* isign == -1 for forward transform and +1 for backward transform */
+  {
+    static const Treal tr11 = 0.309016994374947;
+    static const Treal ti11 = 0.951056516295154;
+    static const Treal tr12 = -0.809016994374947;
+    static const Treal ti12 = 0.587785252292473;
+    int i, k, ac, ah;
+    Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
+        ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
+    if (ido == 2) {
+      for (k = 1; k <= l1; ++k) {
+        ac = (5*k - 4)*ido + 1;
+        ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);
+        ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);
+        ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);
+        ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);
+        tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);
+        tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);
+        tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);
+        tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);
+        ah = (k - 1)*ido;
+        ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;
+        ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;
+        cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;
+        ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;
+        cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;
+        ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;
+        cr5 = isign*(ti11*tr5 + ti12*tr4);
+        ci5 = isign*(ti11*ti5 + ti12*ti4);
+        cr4 = isign*(ti12*tr5 - ti11*tr4);
+        ci4 = isign*(ti12*ti5 - ti11*ti4);
+        ch[ah + l1*ido] = cr2 - ci5;
+        ch[ah + 4*l1*ido] = cr2 + ci5;
+        ch[ah + l1*ido + 1] = ci2 + cr5;
+        ch[ah + 2*l1*ido + 1] = ci3 + cr4;
+        ch[ah + 2*l1*ido] = cr3 - ci4;
+        ch[ah + 3*l1*ido] = cr3 + ci4;
+        ch[ah + 3*l1*ido + 1] = ci3 - cr4;
+        ch[ah + 4*l1*ido + 1] = ci2 - cr5;
+      }
+    } else {
+      for (k=1; k<=l1; k++) {
+        for (i=0; i<ido-1; i+=2) {
+          ac = i + 1 + (k*5 - 4)*ido;
+          ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);
+          ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);
+          ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);
+          ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);
+          tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);
+          tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);
+          tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);
+          tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);
+          ah = i + (k - 1)*ido;
+          ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;
+          ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;
+          cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;
+
+          ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;
+          cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;
+
+          ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;
+          cr5 = isign*(ti11*tr5 + ti12*tr4);
+          ci5 = isign*(ti11*ti5 + ti12*ti4);
+          cr4 = isign*(ti12*tr5 - ti11*tr4);
+          ci4 = isign*(ti12*ti5 - ti11*ti4);
+          dr3 = cr3 - ci4;
+          dr4 = cr3 + ci4;
+          di3 = ci3 + cr4;
+          di4 = ci3 - cr4;
+          dr5 = cr2 + ci5;
+          dr2 = cr2 - ci5;
+          di5 = ci2 - cr5;
+          di2 = ci2 + cr5;
+          ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
+          ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
+          ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
+          ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
+          ch[ah + 3*l1*ido] = wa3[i]*dr4 - isign*wa3[i+1]*di4;
+          ch[ah + 3*l1*ido + 1] = wa3[i]*di4 + isign*wa3[i+1]*dr4;
+          ch[ah + 4*l1*ido] = wa4[i]*dr5 - isign*wa4[i+1]*di5;
+          ch[ah + 4*l1*ido + 1] = wa4[i]*di5 + isign*wa4[i+1]*dr5;
+        }
+      }
+    }
+  } /* passf5 */
+
+
+static void passf(int *nac, int ido, int ip, int l1, int idl1,
+      Treal cc[], Treal ch[],
+      const Treal wa[], int isign)
+  /* isign is -1 for forward transform and +1 for backward transform */
+  {
+    int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, nt, idj, idl, inc,idp;
+    Treal wai, war;
+
+    idot = ido / 2;
+    nt = ip*idl1;
+    ipph = (ip + 1) / 2;
+    idp = ip*ido;
+    if (ido >= l1) {
+      for (j=1; j<ipph; j++) {
+        jc = ip - j;
+        for (k=0; k<l1; k++) {
+          for (i=0; i<ido; i++) {
+            ch[i + (k + j*l1)*ido] =
+                ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*ip)*ido);
+            ch[i + (k + jc*l1)*ido] =
+                ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*ip)*ido);
+          }
+        }
+      }
+      for (k=0; k<l1; k++)
+        for (i=0; i<ido; i++)
+          ch[i + k*ido] = ref(cc,i + k*ip*ido);
+    } else {
+      for (j=1; j<ipph; j++) {
+        jc = ip - j;
+        for (i=0; i<ido; i++) {
+          for (k=0; k<l1; k++) {
+            ch[i + (k + j*l1)*ido] = ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*
+                ip)*ido);
+            ch[i + (k + jc*l1)*ido] = ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*
+                ip)*ido);
+          }
+        }
+      }
+      for (i=0; i<ido; i++)
+        for (k=0; k<l1; k++)
+          ch[i + k*ido] = ref(cc,i + k*ip*ido);
+    }
+
+    idl = 2 - ido;
+    inc = 0;
+    for (l=1; l<ipph; l++) {
+      lc = ip - l;
+      idl += ido;
+      for (ik=0; ik<idl1; ik++) {
+        cc[ik + l*idl1] = ch[ik] + wa[idl - 2]*ch[ik + idl1];
+        cc[ik + lc*idl1] = isign*wa[idl-1]*ch[ik + (ip-1)*idl1];
+      }
+      idlj = idl;
+      inc += ido;
+      for (j=2; j<ipph; j++) {
+        jc = ip - j;
+        idlj += inc;
+        if (idlj > idp) idlj -= idp;
+        war = wa[idlj - 2];
+        wai = wa[idlj-1];
+        for (ik=0; ik<idl1; ik++) {
+          cc[ik + l*idl1] += war*ch[ik + j*idl1];
+          cc[ik + lc*idl1] += isign*wai*ch[ik + jc*idl1];
+        }
+      }
+    }
+    for (j=1; j<ipph; j++)
+      for (ik=0; ik<idl1; ik++)
+        ch[ik] += ch[ik + j*idl1];
+    for (j=1; j<ipph; j++) {
+      jc = ip - j;
+      for (ik=1; ik<idl1; ik+=2) {
+        ch[ik - 1 + j*idl1] = cc[ik - 1 + j*idl1] - cc[ik + jc*idl1];
+        ch[ik - 1 + jc*idl1] = cc[ik - 1 + j*idl1] + cc[ik + jc*idl1];
+        ch[ik + j*idl1] = cc[ik + j*idl1] + cc[ik - 1 + jc*idl1];
+        ch[ik + jc*idl1] = cc[ik + j*idl1] - cc[ik - 1 + jc*idl1];
+      }
+    }
+    *nac = 1;
+    if (ido == 2) return;
+    *nac = 0;
+    for (ik=0; ik<idl1; ik++)
+      cc[ik] = ch[ik];
+    for (j=1; j<ip; j++) {
+      for (k=0; k<l1; k++) {
+        cc[(k + j*l1)*ido + 0] = ch[(k + j*l1)*ido + 0];
+        cc[(k + j*l1)*ido + 1] = ch[(k + j*l1)*ido + 1];
+      }
+    }
+    if (idot <= l1) {
+      idij = 0;
+      for (j=1; j<ip; j++) {
+        idij += 2;
+        for (i=3; i<ido; i+=2) {
+          idij += 2;
+          for (k=0; k<l1; k++) {
+            cc[i - 1 + (k + j*l1)*ido] =
+                wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
+                isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
+            cc[i + (k + j*l1)*ido] =
+                wa[idij - 2]*ch[i + (k + j*l1)*ido] +
+                isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
+          }
+        }
+      }
+    } else {
+      idj = 2 - ido;
+      for (j=1; j<ip; j++) {
+        idj += ido;
+        for (k = 0; k < l1; k++) {
+          idij = idj;
+          for (i=3; i<ido; i+=2) {
+            idij += 2;
+            cc[i - 1 + (k + j*l1)*ido] =
+                wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
+                isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
+            cc[i + (k + j*l1)*ido] =
+                wa[idij - 2]*ch[i + (k + j*l1)*ido] +
+                isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
+          }
+        }
+      }
+    }
+  } /* passf */
+
+
+  /* ----------------------------------------------------------------------
+radf2,radb2, radf3,radb3, radf4,radb4, radf5,radb5, radfg,radbg.
+Treal FFT passes fwd and bwd.
+---------------------------------------------------------------------- */
+
+static void radf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])
+  {
+    int i, k, ic;
+    Treal ti2, tr2;
+    for (k=0; k<l1; k++) {
+      ch[2*k*ido] =
+          ref(cc,k*ido) + ref(cc,(k + l1)*ido);
+      ch[(2*k+1)*ido + ido-1] =
+          ref(cc,k*ido) - ref(cc,(k + l1)*ido);
+    }
+    if (ido < 2) return;
+    if (ido != 2) {
+      for (k=0; k<l1; k++) {
+        for (i=2; i<ido; i+=2) {
+          ic = ido - i;
+          tr2 = wa1[i - 2]*ref(cc, i-1 + (k + l1)*ido) + wa1[i - 1]*ref(cc, i + (k + l1)*ido);
+          ti2 = wa1[i - 2]*ref(cc, i + (k + l1)*ido) - wa1[i - 1]*ref(cc, i-1 + (k + l1)*ido);
+          ch[i + 2*k*ido] = ref(cc,i + k*ido) + ti2;
+          ch[ic + (2*k+1)*ido] = ti2 - ref(cc,i + k*ido);
+          ch[i - 1 + 2*k*ido] = ref(cc,i - 1 + k*ido) + tr2;
+          ch[ic - 1 + (2*k+1)*ido] = ref(cc,i - 1 + k*ido) - tr2;
+        }
+      }
+      if (ido % 2 == 1) return;
+    }
+    for (k=0; k<l1; k++) {
+      ch[(2*k+1)*ido] = -ref(cc,ido-1 + (k + l1)*ido);
+      ch[ido-1 + 2*k*ido] = ref(cc,ido-1 + k*ido);
+    }
+  } /* radf2 */
+
+
+static void radb2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])
+  {
+    int i, k, ic;
+    Treal ti2, tr2;
+    for (k=0; k<l1; k++) {
+      ch[k*ido] =
+          ref(cc,2*k*ido) + ref(cc,ido-1 + (2*k+1)*ido);
+      ch[(k + l1)*ido] =
+          ref(cc,2*k*ido) - ref(cc,ido-1 + (2*k+1)*ido);
+    }
+    if (ido < 2) return;
+    if (ido != 2) {
+      for (k = 0; k < l1; ++k) {
+        for (i = 2; i < ido; i += 2) {
+          ic = ido - i;
+          ch[i-1 + k*ido] =
+              ref(cc,i-1 + 2*k*ido) + ref(cc,ic-1 + (2*k+1)*ido);
+          tr2 = ref(cc,i-1 + 2*k*ido) - ref(cc,ic-1 + (2*k+1)*ido);
+          ch[i + k*ido] =
+              ref(cc,i + 2*k*ido) - ref(cc,ic + (2*k+1)*ido);
+          ti2 = ref(cc,i + (2*k)*ido) + ref(cc,ic + (2*k+1)*ido);
+          ch[i-1 + (k + l1)*ido] =
+              wa1[i - 2]*tr2 - wa1[i - 1]*ti2;
+          ch[i + (k + l1)*ido] =
+              wa1[i - 2]*ti2 + wa1[i - 1]*tr2;
+        }
+      }
+      if (ido % 2 == 1) return;
+    }
+    for (k = 0; k < l1; k++) {
+      ch[ido-1 + k*ido] = 2*ref(cc,ido-1 + 2*k*ido);
+      ch[ido-1 + (k + l1)*ido] = -2*ref(cc,(2*k+1)*ido);
+    }
+  } /* radb2 */
+
+
+static void radf3(int ido, int l1, const Treal cc[], Treal ch[],
+      const Treal wa1[], const Treal wa2[])
+  {
+    static const Treal taur = -0.5;
+    static const Treal taui = 0.866025403784439;
+    int i, k, ic;
+    Treal ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
+    for (k=0; k<l1; k++) {
+      cr2 = ref(cc,(k + l1)*ido) + ref(cc,(k + 2*l1)*ido);
+      ch[3*k*ido] = ref(cc,k*ido) + cr2;
+      ch[(3*k+2)*ido] = taui*(ref(cc,(k + l1*2)*ido) - ref(cc,(k + l1)*ido));
+      ch[ido-1 + (3*k + 1)*ido] = ref(cc,k*ido) + taur*cr2;
+    }
+    if (ido == 1) return;
+    for (k=0; k<l1; k++) {
+      for (i=2; i<ido; i+=2) {
+        ic = ido - i;
+        dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) +
+            wa1[i - 1]*ref(cc,i + (k + l1)*ido);
+        di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
+        dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + l1*2)*ido) + wa2[i - 1]*ref(cc,i + (k + l1*2)*ido);
+        di3 = wa2[i - 2]*ref(cc,i + (k + l1*2)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + l1*2)*ido);
+        cr2 = dr2 + dr3;
+        ci2 = di2 + di3;
+        ch[i - 1 + 3*k*ido] = ref(cc,i - 1 + k*ido) + cr2;
+        ch[i + 3*k*ido] = ref(cc,i + k*ido) + ci2;
+        tr2 = ref(cc,i - 1 + k*ido) + taur*cr2;
+        ti2 = ref(cc,i + k*ido) + taur*ci2;
+        tr3 = taui*(di2 - di3);
+        ti3 = taui*(dr3 - dr2);
+        ch[i - 1 + (3*k + 2)*ido] = tr2 + tr3;
+        ch[ic - 1 + (3*k + 1)*ido] = tr2 - tr3;
+        ch[i + (3*k + 2)*ido] = ti2 + ti3;
+        ch[ic + (3*k + 1)*ido] = ti3 - ti2;
+      }
+    }
+  } /* radf3 */
+
+
+static void radb3(int ido, int l1, const Treal cc[], Treal ch[],
+      const Treal wa1[], const Treal wa2[])
+  {
+    static const Treal taur = -0.5;
+    static const Treal taui = 0.866025403784439;
+    int i, k, ic;
+    Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
+    for (k=0; k<l1; k++) {
+      tr2 = 2*ref(cc,ido-1 + (3*k + 1)*ido);
+      cr2 = ref(cc,3*k*ido) + taur*tr2;
+      ch[k*ido] = ref(cc,3*k*ido) + tr2;
+      ci3 = 2*taui*ref(cc,(3*k + 2)*ido);
+      ch[(k + l1)*ido] = cr2 - ci3;
+      ch[(k + 2*l1)*ido] = cr2 + ci3;
+    }
+    if (ido == 1) return;
+    for (k=0; k<l1; k++) {
+      for (i=2; i<ido; i+=2) {
+        ic = ido - i;
+        tr2 = ref(cc,i - 1 + (3*k + 2)*ido) + ref(cc,ic - 1 + (3*k + 1)*ido);
+        cr2 = ref(cc,i - 1 + 3*k*ido) + taur*tr2;
+        ch[i - 1 + k*ido] = ref(cc,i - 1 + 3*k*ido) + tr2;
+        ti2 = ref(cc,i + (3*k + 2)*ido) - ref(cc,ic + (3*k + 1)*ido);
+        ci2 = ref(cc,i + 3*k*ido) + taur*ti2;
+        ch[i + k*ido] = ref(cc,i + 3*k*ido) + ti2;
+        cr3 = taui*(ref(cc,i - 1 + (3*k + 2)*ido) - ref(cc,ic - 1 + (3*k + 1)*ido));
+        ci3 = taui*(ref(cc,i + (3*k + 2)*ido) + ref(cc,ic + (3*k + 1)*ido));
+        dr2 = cr2 - ci3;
+        dr3 = cr2 + ci3;
+        di2 = ci2 + cr3;
+        di3 = ci2 - cr3;
+        ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
+        ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
+        ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
+        ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
+      }
+    }
+  } /* radb3 */
+
+
+static void radf4(int ido, int l1, const Treal cc[], Treal ch[],
+      const Treal wa1[], const Treal wa2[], const Treal wa3[])
+  {
+    static const Treal hsqt2 = 0.7071067811865475;
+    int i, k, ic;
+    Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
+    for (k=0; k<l1; k++) {
+      tr1 = ref(cc,(k + l1)*ido) + ref(cc,(k + 3*l1)*ido);
+      tr2 = ref(cc,k*ido) + ref(cc,(k + 2*l1)*ido);
+      ch[4*k*ido] = tr1 + tr2;
+      ch[ido-1 + (4*k + 3)*ido] = tr2 - tr1;
+      ch[ido-1 + (4*k + 1)*ido] = ref(cc,k*ido) - ref(cc,(k + 2*l1)*ido);
+      ch[(4*k + 2)*ido] = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + l1)*ido);
+    }
+    if (ido < 2) return;
+    if (ido != 2) {
+      for (k=0; k<l1; k++) {
+        for (i=2; i<ido; i += 2) {
+          ic = ido - i;
+          cr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
+          ci2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
+          cr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*
+              ido);
+          ci3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*
+              ido);
+          cr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*
+              ido);
+          ci4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*
+              ido);
+          tr1 = cr2 + cr4;
+          tr4 = cr4 - cr2;
+          ti1 = ci2 + ci4;
+          ti4 = ci2 - ci4;
+          ti2 = ref(cc,i + k*ido) + ci3;
+          ti3 = ref(cc,i + k*ido) - ci3;
+          tr2 = ref(cc,i - 1 + k*ido) + cr3;
+          tr3 = ref(cc,i - 1 + k*ido) - cr3;
+          ch[i - 1 + 4*k*ido] = tr1 + tr2;
+          ch[ic - 1 + (4*k + 3)*ido] = tr2 - tr1;
+          ch[i + 4*k*ido] = ti1 + ti2;
+          ch[ic + (4*k + 3)*ido] = ti1 - ti2;
+          ch[i - 1 + (4*k + 2)*ido] = ti4 + tr3;
+          ch[ic - 1 + (4*k + 1)*ido] = tr3 - ti4;
+          ch[i + (4*k + 2)*ido] = tr4 + ti3;
+          ch[ic + (4*k + 1)*ido] = tr4 - ti3;
+        }
+      }
+      if (ido % 2 == 1) return;
+    }
+    for (k=0; k<l1; k++) {
+      ti1 = -hsqt2*(ref(cc,ido-1 + (k + l1)*ido) + ref(cc,ido-1 + (k + 3*l1)*ido));
+      tr1 = hsqt2*(ref(cc,ido-1 + (k + l1)*ido) - ref(cc,ido-1 + (k + 3*l1)*ido));
+      ch[ido-1 + 4*k*ido] = tr1 + ref(cc,ido-1 + k*ido);
+      ch[ido-1 + (4*k + 2)*ido] = ref(cc,ido-1 + k*ido) - tr1;
+      ch[(4*k + 1)*ido] = ti1 - ref(cc,ido-1 + (k + 2*l1)*ido);
+      ch[(4*k + 3)*ido] = ti1 + ref(cc,ido-1 + (k + 2*l1)*ido);
+    }
+  } /* radf4 */
+
+
+static void radb4(int ido, int l1, const Treal cc[], Treal ch[],
+      const Treal wa1[], const Treal wa2[], const Treal wa3[])
+  {
+    static const Treal sqrt2 = 1.414213562373095;
+    int i, k, ic;
+    Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
+    for (k = 0; k < l1; k++) {
+      tr1 = ref(cc,4*k*ido) - ref(cc,ido-1 + (4*k + 3)*ido);
+      tr2 = ref(cc,4*k*ido) + ref(cc,ido-1 + (4*k + 3)*ido);
+      tr3 = ref(cc,ido-1 + (4*k + 1)*ido) + ref(cc,ido-1 + (4*k + 1)*ido);
+      tr4 = ref(cc,(4*k + 2)*ido) + ref(cc,(4*k + 2)*ido);
+      ch[k*ido] = tr2 + tr3;
+      ch[(k + l1)*ido] = tr1 - tr4;
+      ch[(k + 2*l1)*ido] = tr2 - tr3;
+      ch[(k + 3*l1)*ido] = tr1 + tr4;
+    }
+    if (ido < 2) return;
+    if (ido != 2) {
+      for (k = 0; k < l1; ++k) {
+        for (i = 2; i < ido; i += 2) {
+          ic = ido - i;
+          ti1 = ref(cc,i + 4*k*ido) + ref(cc,ic + (4*k + 3)*ido);
+          ti2 = ref(cc,i + 4*k*ido) - ref(cc,ic + (4*k + 3)*ido);
+          ti3 = ref(cc,i + (4*k + 2)*ido) - ref(cc,ic + (4*k + 1)*ido);
+          tr4 = ref(cc,i + (4*k + 2)*ido) + ref(cc,ic + (4*k + 1)*ido);
+          tr1 = ref(cc,i - 1 + 4*k*ido) - ref(cc,ic - 1 + (4*k + 3)*ido);
+          tr2 = ref(cc,i - 1 + 4*k*ido) + ref(cc,ic - 1 + (4*k + 3)*ido);
+          ti4 = ref(cc,i - 1 + (4*k + 2)*ido) - ref(cc,ic - 1 + (4*k + 1)*ido);
+          tr3 = ref(cc,i - 1 + (4*k + 2)*ido) + ref(cc,ic - 1 + (4*k + 1)*ido);
+          ch[i - 1 + k*ido] = tr2 + tr3;
+          cr3 = tr2 - tr3;
+          ch[i + k*ido] = ti2 + ti3;
+          ci3 = ti2 - ti3;
+          cr2 = tr1 - tr4;
+          cr4 = tr1 + tr4;
+          ci2 = ti1 + ti4;
+          ci4 = ti1 - ti4;
+          ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*cr2 - wa1[i - 1]*ci2;
+          ch[i + (k + l1)*ido] = wa1[i - 2]*ci2 + wa1[i - 1]*cr2;
+          ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*cr3 - wa2[i - 1]*ci3;
+          ch[i + (k + 2*l1)*ido] = wa2[i - 2]*ci3 + wa2[i - 1]*cr3;
+          ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*cr4 - wa3[i - 1]*ci4;
+          ch[i + (k + 3*l1)*ido] = wa3[i - 2]*ci4 + wa3[i - 1]*cr4;
+        }
+      }
+      if (ido % 2 == 1) return;
+    }
+    for (k = 0; k < l1; k++) {
+      ti1 = ref(cc,(4*k + 1)*ido) + ref(cc,(4*k + 3)*ido);
+      ti2 = ref(cc,(4*k + 3)*ido) - ref(cc,(4*k + 1)*ido);
+      tr1 = ref(cc,ido-1 + 4*k*ido) - ref(cc,ido-1 + (4*k + 2)*ido);
+      tr2 = ref(cc,ido-1 + 4*k*ido) + ref(cc,ido-1 + (4*k + 2)*ido);
+      ch[ido-1 + k*ido] = tr2 + tr2;
+      ch[ido-1 + (k + l1)*ido] = sqrt2*(tr1 - ti1);
+      ch[ido-1 + (k + 2*l1)*ido] = ti2 + ti2;
+      ch[ido-1 + (k + 3*l1)*ido] = -sqrt2*(tr1 + ti1);
+    }
+  } /* radb4 */
+
+
+static void radf5(int ido, int l1, const Treal cc[], Treal ch[],
+      const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])
+  {
+    static const Treal tr11 = 0.309016994374947;
+    static const Treal ti11 = 0.951056516295154;
+    static const Treal tr12 = -0.809016994374947;
+    static const Treal ti12 = 0.587785252292473;
+    int i, k, ic;
+    Treal ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5,
+        cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
+    for (k = 0; k < l1; k++) {
+      cr2 = ref(cc,(k + 4*l1)*ido) + ref(cc,(k + l1)*ido);
+      ci5 = ref(cc,(k + 4*l1)*ido) - ref(cc,(k + l1)*ido);
+      cr3 = ref(cc,(k + 3*l1)*ido) + ref(cc,(k + 2*l1)*ido);
+      ci4 = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + 2*l1)*ido);
+      ch[5*k*ido] = ref(cc,k*ido) + cr2 + cr3;
+      ch[ido-1 + (5*k + 1)*ido] = ref(cc,k*ido) + tr11*cr2 + tr12*cr3;
+      ch[(5*k + 2)*ido] = ti11*ci5 + ti12*ci4;
+      ch[ido-1 + (5*k + 3)*ido] = ref(cc,k*ido) + tr12*cr2 + tr11*cr3;
+      ch[(5*k + 4)*ido] = ti12*ci5 - ti11*ci4;
+    }
+    if (ido == 1) return;
+    for (k = 0; k < l1; ++k) {
+      for (i = 2; i < ido; i += 2) {
+        ic = ido - i;
+        dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
+        di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
+        dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*ido);
+        di3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*ido);
+        dr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*ido);
+        di4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*ido);
+        dr5 = wa4[i - 2]*ref(cc,i - 1 + (k + 4*l1)*ido) + wa4[i - 1]*ref(cc,i + (k + 4*l1)*ido);
+        di5 = wa4[i - 2]*ref(cc,i + (k + 4*l1)*ido) - wa4[i - 1]*ref(cc,i - 1 + (k + 4*l1)*ido);
+        cr2 = dr2 + dr5;
+        ci5 = dr5 - dr2;
+        cr5 = di2 - di5;
+        ci2 = di2 + di5;
+        cr3 = dr3 + dr4;
+        ci4 = dr4 - dr3;
+        cr4 = di3 - di4;
+        ci3 = di3 + di4;
+        ch[i - 1 + 5*k*ido] = ref(cc,i - 1 + k*ido) + cr2 + cr3;
+        ch[i + 5*k*ido] = ref(cc,i + k*ido) + ci2 + ci3;
+        tr2 = ref(cc,i - 1 + k*ido) + tr11*cr2 + tr12*cr3;
+        ti2 = ref(cc,i + k*ido) + tr11*ci2 + tr12*ci3;
+        tr3 = ref(cc,i - 1 + k*ido) + tr12*cr2 + tr11*cr3;
+        ti3 = ref(cc,i + k*ido) + tr12*ci2 + tr11*ci3;
+        tr5 = ti11*cr5 + ti12*cr4;
+        ti5 = ti11*ci5 + ti12*ci4;
+        tr4 = ti12*cr5 - ti11*cr4;
+        ti4 = ti12*ci5 - ti11*ci4;
+        ch[i - 1 + (5*k + 2)*ido] = tr2 + tr5;
+        ch[ic - 1 + (5*k + 1)*ido] = tr2 - tr5;
+        ch[i + (5*k + 2)*ido] = ti2 + ti5;
+        ch[ic + (5*k + 1)*ido] = ti5 - ti2;
+        ch[i - 1 + (5*k + 4)*ido] = tr3 + tr4;
+        ch[ic - 1 + (5*k + 3)*ido] = tr3 - tr4;
+        ch[i + (5*k + 4)*ido] = ti3 + ti4;
+        ch[ic + (5*k + 3)*ido] = ti4 - ti3;
+      }
+    }
+  } /* radf5 */
+
+
+static void radb5(int ido, int l1, const Treal cc[], Treal ch[],
+      const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])
+  {
+    static const Treal tr11 = 0.309016994374947;
+    static const Treal ti11 = 0.951056516295154;
+    static const Treal tr12 = -0.809016994374947;
+    static const Treal ti12 = 0.587785252292473;
+    int i, k, ic;
+    Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
+        ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
+    for (k = 0; k < l1; k++) {
+      ti5 = 2*ref(cc,(5*k + 2)*ido);
+      ti4 = 2*ref(cc,(5*k + 4)*ido);
+      tr2 = 2*ref(cc,ido-1 + (5*k + 1)*ido);
+      tr3 = 2*ref(cc,ido-1 + (5*k + 3)*ido);
+      ch[k*ido] = ref(cc,5*k*ido) + tr2 + tr3;
+      cr2 = ref(cc,5*k*ido) + tr11*tr2 + tr12*tr3;
+      cr3 = ref(cc,5*k*ido) + tr12*tr2 + tr11*tr3;
+      ci5 = ti11*ti5 + ti12*ti4;
+      ci4 = ti12*ti5 - ti11*ti4;
+      ch[(k + l1)*ido] = cr2 - ci5;
+      ch[(k + 2*l1)*ido] = cr3 - ci4;
+      ch[(k + 3*l1)*ido] = cr3 + ci4;
+      ch[(k + 4*l1)*ido] = cr2 + ci5;
+    }
+    if (ido == 1) return;
+    for (k = 0; k < l1; ++k) {
+      for (i = 2; i < ido; i += 2) {
+        ic = ido - i;
+        ti5 = ref(cc,i + (5*k + 2)*ido) + ref(cc,ic + (5*k + 1)*ido);
+        ti2 = ref(cc,i + (5*k + 2)*ido) - ref(cc,ic + (5*k + 1)*ido);
+        ti4 = ref(cc,i + (5*k + 4)*ido) + ref(cc,ic + (5*k + 3)*ido);
+        ti3 = ref(cc,i + (5*k + 4)*ido) - ref(cc,ic + (5*k + 3)*ido);
+        tr5 = ref(cc,i - 1 + (5*k + 2)*ido) - ref(cc,ic - 1 + (5*k + 1)*ido);
+        tr2 = ref(cc,i - 1 + (5*k + 2)*ido) + ref(cc,ic - 1 + (5*k + 1)*ido);
+        tr4 = ref(cc,i - 1 + (5*k + 4)*ido) - ref(cc,ic - 1 + (5*k + 3)*ido);
+        tr3 = ref(cc,i - 1 + (5*k + 4)*ido) + ref(cc,ic - 1 + (5*k + 3)*ido);
+        ch[i - 1 + k*ido] = ref(cc,i - 1 + 5*k*ido) + tr2 + tr3;
+        ch[i + k*ido] = ref(cc,i + 5*k*ido) + ti2 + ti3;
+        cr2 = ref(cc,i - 1 + 5*k*ido) + tr11*tr2 + tr12*tr3;
+
+        ci2 = ref(cc,i + 5*k*ido) + tr11*ti2 + tr12*ti3;
+        cr3 = ref(cc,i - 1 + 5*k*ido) + tr12*tr2 + tr11*tr3;
+
+        ci3 = ref(cc,i + 5*k*ido) + tr12*ti2 + tr11*ti3;
+        cr5 = ti11*tr5 + ti12*tr4;
+        ci5 = ti11*ti5 + ti12*ti4;
+        cr4 = ti12*tr5 - ti11*tr4;
+        ci4 = ti12*ti5 - ti11*ti4;
+        dr3 = cr3 - ci4;
+        dr4 = cr3 + ci4;
+        di3 = ci3 + cr4;
+        di4 = ci3 - cr4;
+        dr5 = cr2 + ci5;
+        dr2 = cr2 - ci5;
+        di5 = ci2 - cr5;
+        di2 = ci2 + cr5;
+        ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
+        ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
+        ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
+        ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
+        ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*dr4 - wa3[i - 1]*di4;
+        ch[i + (k + 3*l1)*ido] = wa3[i - 2]*di4 + wa3[i - 1]*dr4;
+        ch[i - 1 + (k + 4*l1)*ido] = wa4[i - 2]*dr5 - wa4[i - 1]*di5;
+        ch[i + (k + 4*l1)*ido] = wa4[i - 2]*di5 + wa4[i - 1]*dr5;
+      }
+    }
+  } /* radb5 */
+
+
+static void radfg(int ido, int ip, int l1, int idl1,
+      Treal cc[], Treal ch[], const Treal wa[])
+  {
+    static const Treal twopi = 6.28318530717959;
+    int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;
+    Treal dc2, ai1, ai2, ar1, ar2, ds2, dcp, arg, dsp, ar1h, ar2h;
+    arg = twopi / ip;
+    dcp = cos(arg);
+    dsp = sin(arg);
+    ipph = (ip + 1) / 2;
+    nbd = (ido - 1) / 2;
+    if (ido != 1) {
+      for (ik=0; ik<idl1; ik++) ch[ik] = cc[ik];
+      for (j=1; j<ip; j++)
+        for (k=0; k<l1; k++)
+          ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido];
+      if (nbd <= l1) {
+        is = -ido;
+        for (j=1; j<ip; j++) {
+          is += ido;
+          idij = is-1;
+          for (i=2; i<ido; i+=2) {
+            idij += 2;
+            for (k=0; k<l1; k++) {
+              ch[i - 1 + (k + j*l1)*ido] =
+                  wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
+              ch[i + (k + j*l1)*ido] =
+                  wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
+            }
+          }
+        }
+      } else {
+        is = -ido;
+        for (j=1; j<ip; j++) {
+          is += ido;
+          for (k=0; k<l1; k++) {
+            idij = is-1;
+            for (i=2; i<ido; i+=2) {
+              idij += 2;
+              ch[i - 1 + (k + j*l1)*ido] =
+                  wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
+              ch[i + (k + j*l1)*ido] =
+                  wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
+            }
+          }
+        }
+      }
+      if (nbd >= l1) {
+        for (j=1; j<ipph; j++) {
+          jc = ip - j;
+          for (k=0; k<l1; k++) {
+            for (i=2; i<ido; i+=2) {
+              cc[i - 1 + (k + j*l1)*ido] = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
+              cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
+              cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
+              cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
+            }
+          }
+        }
+      } else {
+        for (j=1; j<ipph; j++) {
+          jc = ip - j;
+          for (i=2; i<ido; i+=2) {
+            for (k=0; k<l1; k++) {
+              cc[i - 1 + (k + j*l1)*ido] =
+                  ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
+              cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
+              cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
+              cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
+            }
+          }
+        }
+      }
+    } else {  /* now ido == 1 */
+      for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];
+    }
+    for (j=1; j<ipph; j++) {
+      jc = ip - j;
+      for (k=0; k<l1; k++) {
+        cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido] + ch[(k + jc*l1)*ido];
+        cc[(k + jc*l1)*ido] = ch[(k + jc*l1)*ido] - ch[(k + j*l1)*ido];
+      }
+    }
+
+    ar1 = 1;
+    ai1 = 0;
+    for (l=1; l<ipph; l++) {
+      lc = ip - l;
+      ar1h = dcp*ar1 - dsp*ai1;
+      ai1 = dcp*ai1 + dsp*ar1;
+      ar1 = ar1h;
+      for (ik=0; ik<idl1; ik++) {
+        ch[ik + l*idl1] = cc[ik] + ar1*cc[ik + idl1];
+        ch[ik + lc*idl1] = ai1*cc[ik + (ip-1)*idl1];
+      }
+      dc2 = ar1;
+      ds2 = ai1;
+      ar2 = ar1;
+      ai2 = ai1;
+      for (j=2; j<ipph; j++) {
+        jc = ip - j;
+        ar2h = dc2*ar2 - ds2*ai2;
+        ai2 = dc2*ai2 + ds2*ar2;
+        ar2 = ar2h;
+        for (ik=0; ik<idl1; ik++) {
+          ch[ik + l*idl1] += ar2*cc[ik + j*idl1];
+          ch[ik + lc*idl1] += ai2*cc[ik + jc*idl1];
+        }
+      }
+    }
+    for (j=1; j<ipph; j++)
+      for (ik=0; ik<idl1; ik++)
+        ch[ik] += cc[ik + j*idl1];
+
+    if (ido >= l1) {
+      for (k=0; k<l1; k++) {
+        for (i=0; i<ido; i++) {
+          ref(cc,i + k*ip*ido) = ch[i + k*ido];
+        }
+      }
+    } else {
+      for (i=0; i<ido; i++) {
+        for (k=0; k<l1; k++) {
+          ref(cc,i + k*ip*ido) = ch[i + k*ido];
+        }
+      }
+    }
+    for (j=1; j<ipph; j++) {
+      jc = ip - j;
+      j2 = 2*j;
+      for (k=0; k<l1; k++) {
+        ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) =
+            ch[(k + j*l1)*ido];
+        ref(cc,(j2 + k*ip)*ido) =
+            ch[(k + jc*l1)*ido];
+      }
+    }
+    if (ido == 1) return;
+    if (nbd >= l1) {
+      for (j=1; j<ipph; j++) {
+        jc = ip - j;
+        j2 = 2*j;
+        for (k=0; k<l1; k++) {
+          for (i=2; i<ido; i+=2) {
+            ic = ido - i;
+            ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
+            ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
+            ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
+            ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
+          }
+        }
+      }
+    } else {
+      for (j=1; j<ipph; j++) {
+        jc = ip - j;
+        j2 = 2*j;
+        for (i=2; i<ido; i+=2) {
+          ic = ido - i;
+          for (k=0; k<l1; k++) {
+            ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
+            ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
+            ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
+            ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
+          }
+        }
+      }
+    }
+  } /* radfg */
+
+
+static void radbg(int ido, int ip, int l1, int idl1,
+      Treal cc[], Treal ch[], const Treal wa[])
+  {
+    static const Treal twopi = 6.28318530717959;
+    int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
+    Treal dc2, ai1, ai2, ar1, ar2, ds2;
+    int nbd;
+    Treal dcp, arg, dsp, ar1h, ar2h;
+    arg = twopi / ip;
+    dcp = cos(arg);
+    dsp = sin(arg);
+    nbd = (ido - 1) / 2;
+    ipph = (ip + 1) / 2;
+    if (ido >= l1) {
+      for (k=0; k<l1; k++) {
+        for (i=0; i<ido; i++) {
+          ch[i + k*ido] = ref(cc,i + k*ip*ido);
+        }
+      }
+    } else {
+      for (i=0; i<ido; i++) {
+        for (k=0; k<l1; k++) {
+          ch[i + k*ido] = ref(cc,i + k*ip*ido);
+        }
+      }
+    }
+    for (j=1; j<ipph; j++) {
+      jc = ip - j;
+      j2 = 2*j;
+      for (k=0; k<l1; k++) {
+        ch[(k + j*l1)*ido] = ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) + ref(cc,ido-1 + (j2 - 1 + k*ip)*
+            ido);
+        ch[(k + jc*l1)*ido] = ref(cc,(j2 + k*ip)*ido) + ref(cc,(j2 + k*ip)*ido);
+      }
+    }
+
+    if (ido != 1) {
+      if (nbd >= l1) {
+        for (j=1; j<ipph; j++) {
+          jc = ip - j;
+          for (k=0; k<l1; k++) {
+            for (i=2; i<ido; i+=2) {
+              ic = ido - i;
+              ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,
+                  ic - 1 + (2*j - 1 + k*ip)*ido);
+              ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -
+                  ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
+              ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic
+                  + (2*j - 1 + k*ip)*ido);
+              ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic
+                  + (2*j - 1 + k*ip)*ido);
+            }
+          }
+        }
+      } else {
+        for (j=1; j<ipph; j++) {
+          jc = ip - j;
+          for (i=2; i<ido; i+=2) {
+            ic = ido - i;
+            for (k=0; k<l1; k++) {
+              ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,
+                  ic - 1 + (2*j - 1 + k*ip)*ido);
+              ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -
+                  ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
+              ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic
+                  + (2*j - 1 + k*ip)*ido);
+              ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic
+                  + (2*j - 1 + k*ip)*ido);
+            }
+          }
+        }
+      }
+    }
+
+    ar1 = 1;
+    ai1 = 0;
+    for (l=1; l<ipph; l++) {
+      lc = ip - l;
+      ar1h = dcp*ar1 - dsp*ai1;
+      ai1 = dcp*ai1 + dsp*ar1;
+      ar1 = ar1h;
+      for (ik=0; ik<idl1; ik++) {
+        cc[ik + l*idl1] = ch[ik] + ar1*ch[ik + idl1];
+        cc[ik + lc*idl1] = ai1*ch[ik + (ip-1)*idl1];
+      }
+      dc2 = ar1;
+      ds2 = ai1;
+      ar2 = ar1;
+      ai2 = ai1;
+      for (j=2; j<ipph; j++) {
+        jc = ip - j;
+        ar2h = dc2*ar2 - ds2*ai2;
+        ai2 = dc2*ai2 + ds2*ar2;
+        ar2 = ar2h;
+        for (ik=0; ik<idl1; ik++) {
+          cc[ik + l*idl1] += ar2*ch[ik + j*idl1];
+          cc[ik + lc*idl1] += ai2*ch[ik + jc*idl1];
+        }
+      }
+    }
+    for (j=1; j<ipph; j++) {
+      for (ik=0; ik<idl1; ik++) {
+        ch[ik] += ch[ik + j*idl1];
+      }
+    }
+    for (j=1; j<ipph; j++) {
+      jc = ip - j;
+      for (k=0; k<l1; k++) {
+        ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido] - cc[(k + jc*l1)*ido];
+        ch[(k + jc*l1)*ido] = cc[(k + j*l1)*ido] + cc[(k + jc*l1)*ido];
+      }
+    }
+
+    if (ido == 1) return;
+    if (nbd >= l1) {
+      for (j=1; j<ipph; j++) {
+        jc = ip - j;
+        for (k=0; k<l1; k++) {
+          for (i=2; i<ido; i+=2) {
+            ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
+            ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] + cc[i + (k + jc*l1)*ido];
+            ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
+            ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
+          }
+        }
+      }
+    } else {
+      for (j=1; j<ipph; j++) {
+        jc = ip - j;
+        for (i=2; i<ido; i+=2) {
+          for (k=0; k<l1; k++) {
+            ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
+            ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j *l1)*ido] + cc[i + (k + jc*l1)*ido];
+            ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
+            ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
+          }
+        }
+      }
+    }
+    for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];
+    for (j=1; j<ip; j++)
+      for (k=0; k<l1; k++)
+        cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido];
+    if (nbd <= l1) {
+      is = -ido;
+      for (j=1; j<ip; j++) {
+        is += ido;
+        idij = is-1;
+        for (i=2; i<ido; i+=2) {
+          idij += 2;
+          for (k=0; k<l1; k++) {
+            cc[i - 1 + (k + j*l1)*ido] = wa[idij - 1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*
+                ch[i + (k + j*l1)*ido];
+            cc[i + (k + j*l1)*ido] = wa[idij - 1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
+          }
+        }
+      }
+    } else {
+      is = -ido;
+      for (j=1; j<ip; j++) {
+        is += ido;
+        for (k=0; k<l1; k++) {
+          idij = is - 1;
+          for (i=2; i<ido; i+=2) {
+            idij += 2;
+            cc[i - 1 + (k + j*l1)*ido] = wa[idij-1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*
+                ch[i + (k + j*l1)*ido];
+            cc[i + (k + j*l1)*ido] = wa[idij-1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
+          }
+        }
+      }
+    }
+  } /* radbg */
+
+  /* ----------------------------------------------------------------------
+cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
+---------------------------------------------------------------------- */
+
+static void cfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2], int isign)
+  {
+    int idot, i;
+    int k1, l1, l2;
+    int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;
+    Treal *cinput, *coutput;
+    nf = ifac[1];
+    na = 0;
+    l1 = 1;
+    iw = 0;
+    for (k1=2; k1<=nf+1; k1++) {
+      ip = ifac[k1];
+      l2 = ip*l1;
+      ido = n / l2;
+      idot = ido + ido;
+      idl1 = idot*l1;
+      if (na) {
+        cinput = ch;
+        coutput = c;
+      } else {
+        cinput = c;
+        coutput = ch;
+      }
+      switch (ip) {
+      case 4:
+        ix2 = iw + idot;
+        ix3 = ix2 + idot;
+        passf4(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], isign);
+        na = !na;
+        break;
+      case 2:
+        passf2(idot, l1, cinput, coutput, &wa[iw], isign);
+        na = !na;
+        break;
+      case 3:
+        ix2 = iw + idot;
+        passf3(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], isign);
+        na = !na;
+        break;
+      case 5:
+        ix2 = iw + idot;
+        ix3 = ix2 + idot;
+        ix4 = ix3 + idot;
+        passf5(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
+        na = !na;
+        break;
+      default:
+        passf(&nac, idot, ip, l1, idl1, cinput, coutput, &wa[iw], isign);
+        if (nac != 0) na = !na;
+      }
+      l1 = l2;
+      iw += (ip - 1)*idot;
+    }
+    if (na == 0) return;
+    for (i=0; i<2*n; i++) c[i] = ch[i];
+  } /* cfftf1 */
+
+
+void cfftf(int n, Treal c[], Treal wsave[])
+  {
+    int iw1, iw2;
+    if (n == 1) return;
+    iw1 = 2*n;
+    iw2 = iw1 + 2*n;
+    cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), -1);
+  } /* cfftf */
+
+
+void cfftb(int n, Treal c[], Treal wsave[])
+  {
+    int iw1, iw2;
+    if (n == 1) return;
+    iw1 = 2*n;
+    iw2 = iw1 + 2*n;
+    cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), +1);
+  } /* cfftb */
+
+
+static void factorize(int n, int ifac[MAXFAC+2], const int ntryh[NSPECIAL])
+  /* Factorize n in factors in ntryh and rest. On exit,
+ifac[0] contains n and ifac[1] contains number of factors,
+the factors start from ifac[2]. */
+  {
+    int ntry=3, i, j=0, ib, nf=0, nl=n, nq, nr;
+startloop:
+    if (j < NSPECIAL)
+      ntry = ntryh[j];
+    else
+      ntry+= 2;
+    j++;
+    do {
+      nq = nl / ntry;
+      nr = nl - ntry*nq;
+      if (nr != 0) goto startloop;
+      nf++;
+      ifac[nf + 1] = ntry;
+      nl = nq;
+      if (ntry == 2 && nf != 1) {
+        for (i=2; i<=nf; i++) {
+          ib = nf - i + 2;
+          ifac[ib + 1] = ifac[ib];
+        }
+        ifac[2] = 2;
+      }
+    } while (nl != 1);
+    ifac[0] = n;
+    ifac[1] = nf;
+  }
+
+
+static void cffti1(int n, Treal wa[], int ifac[MAXFAC+2])
+  {
+    static const Treal twopi = 6.28318530717959;
+    Treal arg, argh, argld, fi;
+    int idot, i, j;
+    int i1, k1, l1, l2;
+    int ld, ii, nf, ip;
+    int ido, ipm;
+
+    static const int ntryh[NSPECIAL] = {
+      3,4,2,5    }; /* Do not change the order of these. */
+
+    factorize(n,ifac,ntryh);
+    nf = ifac[1];
+    argh = twopi/(Treal)n;
+    i = 1;
+    l1 = 1;
+    for (k1=1; k1<=nf; k1++) {
+      ip = ifac[k1+1];
+      ld = 0;
+      l2 = l1*ip;
+      ido = n / l2;
+      idot = ido + ido + 2;
+      ipm = ip - 1;
+      for (j=1; j<=ipm; j++) {
+        i1 = i;
+        wa[i-1] = 1;
+        wa[i] = 0;
+        ld += l1;
+        fi = 0;
+        argld = ld*argh;
+        for (ii=4; ii<=idot; ii+=2) {
+          i+= 2;
+          fi+= 1;
+          arg = fi*argld;
+          wa[i-1] = cos(arg);
+          wa[i] = sin(arg);
+        }
+        if (ip > 5) {
+          wa[i1-1] = wa[i-1];
+          wa[i1] = wa[i];
+        }
+      }
+      l1 = l2;
+    }
+  } /* cffti1 */
+
+
+void cffti(int n, Treal wsave[])
+ {
+    int iw1, iw2;
+    if (n == 1) return;
+    iw1 = 2*n;
+    iw2 = iw1 + 2*n;
+    cffti1(n, wsave+iw1, (int*)(wsave+iw2));
+  } /* cffti */
+
+  /* ----------------------------------------------------------------------
+rfftf1, rfftb1, rfftf, rfftb, rffti1, rffti. Treal FFTs.
+---------------------------------------------------------------------- */
+
+static void rfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])
+  {
+    int i;
+    int k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
+    Treal *cinput, *coutput;
+    nf = ifac[1];
+    na = 1;
+    l2 = n;
+    iw = n-1;
+    for (k1 = 1; k1 <= nf; ++k1) {
+      kh = nf - k1;
+      ip = ifac[kh + 2];
+      l1 = l2 / ip;
+      ido = n / l2;
+      idl1 = ido*l1;
+      iw -= (ip - 1)*ido;
+      na = !na;
+      if (na) {
+        cinput = ch;
+        coutput = c;
+      } else {
+        cinput = c;
+        coutput = ch;
+      }
+      switch (ip) {
+      case 4:
+        ix2 = iw + ido;
+        ix3 = ix2 + ido;
+        radf4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
+        break;
+      case 2:
+        radf2(ido, l1, cinput, coutput, &wa[iw]);
+        break;
+      case 3:
+        ix2 = iw + ido;
+        radf3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
+        break;
+      case 5:
+        ix2 = iw + ido;
+        ix3 = ix2 + ido;
+        ix4 = ix3 + ido;
+        radf5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
+        break;
+      default:
+        if (ido == 1)
+          na = !na;
+        if (na == 0) {
+          radfg(ido, ip, l1, idl1, c, ch, &wa[iw]);
+          na = 1;
+        } else {
+          radfg(ido, ip, l1, idl1, ch, c, &wa[iw]);
+          na = 0;
+        }
+      }
+      l2 = l1;
+    }
+    if (na == 1) return;
+    for (i = 0; i < n; i++) c[i] = ch[i];
+  } /* rfftf1 */
+
+
+void rfftb1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])
+  {
+    int i;
+    int k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
+    Treal *cinput, *coutput;
+    nf = ifac[1];
+    na = 0;
+    l1 = 1;
+    iw = 0;
+    for (k1=1; k1<=nf; k1++) {
+      ip = ifac[k1 + 1];
+      l2 = ip*l1;
+      ido = n / l2;
+      idl1 = ido*l1;
+      if (na) {
+        cinput = ch;
+        coutput = c;
+      } else {
+        cinput = c;
+        coutput = ch;
+      }
+      switch (ip) {
+      case 4:
+        ix2 = iw + ido;
+        ix3 = ix2 + ido;
+        radb4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
+        na = !na;
+        break;
+      case 2:
+        radb2(ido, l1, cinput, coutput, &wa[iw]);
+        na = !na;
+        break;
+      case 3:
+        ix2 = iw + ido;
+        radb3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
+        na = !na;
+        break;
+      case 5:
+        ix2 = iw + ido;
+        ix3 = ix2 + ido;
+        ix4 = ix3 + ido;
+        radb5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
+        na = !na;
+        break;
+      default:
+        radbg(ido, ip, l1, idl1, cinput, coutput, &wa[iw]);
+        if (ido == 1) na = !na;
+      }
+      l1 = l2;
+      iw += (ip - 1)*ido;
+    }
+    if (na == 0) return;
+    for (i=0; i<n; i++) c[i] = ch[i];
+  } /* rfftb1 */
+
+
+void rfftf(int n, Treal r[], Treal wsave[])
+  {
+    if (n == 1) return;
+    rfftf1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
+  } /* rfftf */
+
+
+void rfftb(int n, Treal r[], Treal wsave[])
+  {
+    if (n == 1) return;
+    rfftb1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
+  } /* rfftb */
+
+
+static void rffti1(int n, Treal wa[], int ifac[MAXFAC+2])
+  {
+    static const Treal twopi = 6.28318530717959;
+    Treal arg, argh, argld, fi;
+    int i, j;
+    int k1, l1, l2;
+    int ld, ii, nf, ip, is;
+    int ido, ipm, nfm1;
+    static const int ntryh[NSPECIAL] = {
+      4,2,3,5    }; /* Do not change the order of these. */
+    factorize(n,ifac,ntryh);
+    nf = ifac[1];
+    argh = twopi / n;
+    is = 0;
+    nfm1 = nf - 1;
+    l1 = 1;
+    if (nfm1 == 0) return;
+    for (k1 = 1; k1 <= nfm1; k1++) {
+      ip = ifac[k1 + 1];
+      ld = 0;
+      l2 = l1*ip;
+      ido = n / l2;
+      ipm = ip - 1;
+      for (j = 1; j <= ipm; ++j) {
+        ld += l1;
+        i = is;
+        argld = (Treal) ld*argh;
+        fi = 0;
+        for (ii = 3; ii <= ido; ii += 2) {
+          i += 2;
+          fi += 1;
+          arg = fi*argld;
+          wa[i - 2] = cos(arg);
+          wa[i - 1] = sin(arg);
+        }
+        is += ido;
+      }
+      l1 = l2;
+    }
+  } /* rffti1 */
+
+
+void rffti(int n, Treal wsave[])
+  {
+    if (n == 1) return;
+    rffti1(n, wsave+n, (int*)(wsave+2*n));
+  } /* rffti */
+
+#ifdef __cplusplus
+}
+#endif
diff --git a/numpy/fft/fftpack.h b/numpy/fft/fftpack.h
new file mode 100644
index 000000000..d134784a2
--- /dev/null
+++ b/numpy/fft/fftpack.h
@@ -0,0 +1,28 @@
+/*
+ * This file is part of tela the Tensor Language.
+ * Copyright (c) 1994-1995 Pekka Janhunen
+ */
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#define DOUBLE
+
+#ifdef DOUBLE
+#define Treal double
+#else
+#define Treal float
+#endif
+
+extern void cfftf(int N, Treal data[], const Treal wrk[]);
+extern void cfftb(int N, Treal data[], const Treal wrk[]);
+extern void cffti(int N, Treal wrk[]);
+
+extern void rfftf(int N, Treal data[], const Treal wrk[]);
+extern void rfftb(int N, Treal data[], const Treal wrk[]);
+extern void rffti(int N, Treal wrk[]);
+
+#ifdef __cplusplus
+}
+#endif
diff --git a/numpy/fft/fftpack.py b/numpy/fft/fftpack.py
new file mode 100644
index 000000000..1cb24f2b7
--- /dev/null
+++ b/numpy/fft/fftpack.py
@@ -0,0 +1,323 @@
+"""
+Discrete Fourier Transforms - FFT.py
+
+The underlying code for these functions is an f2c translated and modified
+version of the FFTPACK routines.
+
+fft(a, n=None, axis=-1)
+ifft(a, n=None, axis=-1)
+rfft(a, n=None, axis=-1)
+irfft(a, n=None, axis=-1)
+hfft(a, n=None, axis=-1)
+ihfft(a, n=None, axis=-1)
+fftn(a, s=None, axes=None)
+ifftn(a, s=None, axes=None)
+rfftn(a, s=None, axes=None)
+irfftn(a, s=None, axes=None)
+fft2(a, s=None, axes=(-2,-1))
+ifft2(a, s=None, axes=(-2, -1))
+rfft2(a, s=None, axes=(-2,-1))
+irfft2(a, s=None, axes=(-2, -1))
+"""
+__all__ = ['fft','ifft', 'rfft', 'irfft', 'hfft', 'ihfft', 'rfftn',
+           'irfftn', 'rfft2', 'irfft2', 'fft2', 'ifft2', 'fftn', 'ifftn',
+           'refft', 'irefft','refftn','irefftn', 'refft2', 'irefft2']
+
+from numpy.core import asarray, zeros, swapaxes, shape, conjugate, \
+     take
+import fftpack_lite as fftpack
+from helper import *
+
+_fft_cache = {}
+_real_fft_cache = {}
+
+def _raw_fft(a, n=None, axis=-1, init_function=fftpack.cffti,
+             work_function=fftpack.cfftf, fft_cache = _fft_cache ):
+    a = asarray(a)
+
+    if n == None: n = a.shape[axis]
+
+    try:
+        wsave = fft_cache[n]
+    except(KeyError):
+        wsave = init_function(n)
+        fft_cache[n] = wsave
+
+    if a.shape[axis] != n:
+        s = list(a.shape)
+        if s[axis] > n:
+            index = [slice(None)]*len(s)
+            index[axis] = slice(0,n)
+            a = a[index]
+        else:
+            index = [slice(None)]*len(s)
+            index[axis] = slice(0,s[axis])
+            s[axis] = n
+            z = zeros(s, a.dtype.char)
+            z[index] = a
+            a = z
+
+    if axis != -1:
+        a = swapaxes(a, axis, -1)
+    r = work_function(a, wsave)
+    if axis != -1:
+        r = swapaxes(r, axis, -1)
+    return r
+
+
+def fft(a, n=None, axis=-1):
+    """fft(a, n=None, axis=-1)
+
+    Will return the n point discrete Fourier transform of a. n defaults to the
+    length of a. If n is larger than a, then a will be zero-padded to make up
+    the difference. If n is smaller than a, the first n items in a will be
+    used.
+
+    The packing of the result is "standard": If A = fft(a, n), then A[0]
+    contains the zero-frequency term, A[1:n/2+1] contains the
+    positive-frequency terms, and A[n/2+1:] contains the negative-frequency
+    terms, in order of decreasingly negative frequency. So for an 8-point
+    transform, the frequencies of the result are [ 0, 1, 2, 3, 4, -3, -2, -1].
+
+    This is most efficient for n a power of two. This also stores a cache of
+    working memory for different sizes of fft's, so you could theoretically
+    run into memory problems if you call this too many times with too many
+    different n's."""
+
+    return _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftf, _fft_cache)
+
+
+def ifft(a, n=None, axis=-1):
+    """ifft(a, n=None, axis=-1)
+
+    Will return the n point inverse discrete Fourier transform of a.  n
+    defaults to the length of a. If n is larger than a, then a will be
+    zero-padded to make up the difference. If n is smaller than a, then a will
+    be truncated to reduce its size.
+
+    The input array is expected to be packed the same way as the output of
+    fft, as discussed in it's documentation.
+
+    This is the inverse of fft: ifft(fft(a)) == a within numerical
+    accuracy.
+
+    This is most efficient for n a power of two. This also stores a cache of
+    working memory for different sizes of fft's, so you could theoretically
+    run into memory problems if you call this too many times with too many
+    different n's."""
+
+    a = asarray(a).astype(complex)
+    if n == None:
+        n = shape(a)[axis]
+    return _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftb, _fft_cache) / n
+
+
+def rfft(a, n=None, axis=-1):
+    """rfft(a, n=None, axis=-1)
+
+    Will return the n point discrete Fourier transform of the real valued
+    array a. n defaults to the length of a. n is the length of the input, not
+    the output.
+
+    The returned array will be the nonnegative frequency terms of the
+    Hermite-symmetric, complex transform of the real array. So for an 8-point
+    transform, the frequencies in the result are [ 0, 1, 2, 3, 4]. The first
+    term will be real, as will the last if n is even. The negative frequency
+    terms are not needed because they are the complex conjugates of the
+    positive frequency terms. (This is what I mean when I say
+    Hermite-symmetric.)
+
+    This is most efficient for n a power of two."""
+
+    a = asarray(a).astype(float)
+    return _raw_fft(a, n, axis, fftpack.rffti, fftpack.rfftf, _real_fft_cache)
+
+
+def irfft(a, n=None, axis=-1):
+    """irfft(a, n=None, axis=-1)
+
+    Will return the real valued n point inverse discrete Fourier transform of
+    a, where a contains the nonnegative frequency terms of a Hermite-symmetric
+    sequence. n is the length of the result, not the input. If n is not
+    supplied, the default is 2*(len(a)-1). If you want the length of the
+    result to be odd, you have to say so.
+
+    If you specify an n such that a must be zero-padded or truncated, the
+    extra/removed values will be added/removed at high frequencies. One can
+    thus resample a series to m points via Fourier interpolation by: a_resamp
+    = irfft(rfft(a), m).
+
+    This is the inverse of rfft:
+    irfft(rfft(a), len(a)) == a
+    within numerical accuracy."""
+
+    a = asarray(a).astype(complex)
+    if n == None:
+        n = (shape(a)[axis] - 1) * 2
+    return _raw_fft(a, n, axis, fftpack.rffti, fftpack.rfftb,
+                    _real_fft_cache) / n
+
+
+def hfft(a, n=None, axis=-1):
+    """hfft(a, n=None, axis=-1)
+    ihfft(a, n=None, axis=-1)
+
+    These are a pair analogous to rfft/irfft, but for the
+    opposite case: here the signal is real in the frequency domain and has
+    Hermite symmetry in the time domain. So here it's hermite_fft for which
+    you must supply the length of the result if it is to be odd.
+
+    ihfft(hfft(a), len(a)) == a
+    within numerical accuracy."""
+
+    a = asarray(a).astype(complex)
+    if n == None:
+        n = (shape(a)[axis] - 1) * 2
+    return irfft(conjugate(a), n, axis) * n
+
+
+def ihfft(a, n=None, axis=-1):
+    """hfft(a, n=None, axis=-1)
+    ihfft(a, n=None, axis=-1)
+
+    These are a pair analogous to rfft/irfft, but for the
+    opposite case: here the signal is real in the frequency domain and has
+    Hermite symmetry in the time domain. So here it's hfft for which
+    you must supply the length of the result if it is to be odd.
+
+    ihfft(hfft(a), len(a)) == a
+    within numerical accuracy."""
+
+    a = asarray(a).astype(float)
+    if n == None:
+        n = shape(a)[axis]
+    return conjugate(rfft(a, n, axis))/n
+
+
+def _cook_nd_args(a, s=None, axes=None, invreal=0):
+    if s is None:
+        shapeless = 1
+        if axes == None:
+            s = list(a.shape)
+        else:
+            s = take(a.shape, axes)
+    else:
+        shapeless = 0
+    s = list(s)
+    if axes == None:
+        axes = range(-len(s), 0)
+    if len(s) != len(axes):
+        raise ValueError, "Shape and axes have different lengths."
+    if invreal and shapeless:
+        s[axes[-1]] = (s[axes[-1]] - 1) * 2
+    return s, axes
+
+
+def _raw_fftnd(a, s=None, axes=None, function=fft):
+    a = asarray(a)
+    s, axes = _cook_nd_args(a, s, axes)
+    itl = range(len(axes))
+    itl.reverse()
+    for ii in itl:
+        a = function(a, n=s[ii], axis=axes[ii])
+    return a
+
+
+def fftn(a, s=None, axes=None):
+    """fftn(a, s=None, axes=None)
+
+    The n-dimensional fft of a. s is a sequence giving the shape of the input
+    an result along the transformed axes, as n for fft. Results are packed
+    analogously to fft: the term for zero frequency in all axes is in the
+    low-order corner, while the term for the Nyquist frequency in all axes is
+    in the middle.
+
+    If neither s nor axes is specified, the transform is taken along all
+    axes. If s is specified and axes is not, the last len(s) axes are used.
+    If axes are specified and s is not, the input shape along the specified
+    axes is used. If s and axes are both specified and are not the same
+    length, an exception is raised."""
+
+    return _raw_fftnd(a,s,axes,fft)
+
+def ifftn(a, s=None, axes=None):
+    """ifftn(a, s=None, axes=None)
+
+    The inverse of fftn."""
+
+    return _raw_fftnd(a, s, axes, ifft)
+
+
+def fft2(a, s=None, axes=(-2,-1)):
+    """fft2(a, s=None, axes=(-2,-1))
+
+    The 2d fft of a. This is really just fftn with different default
+    behavior."""
+
+    return _raw_fftnd(a,s,axes,fft)
+
+
+def ifft2(a, s=None, axes=(-2,-1)):
+    """ifft2(a, s=None, axes=(-2, -1))
+
+    The inverse of fft2d. This is really just ifftn with different
+    default behavior."""
+
+    return _raw_fftnd(a, s, axes, ifft)
+
+
+def rfftn(a, s=None, axes=None):
+    """rfftn(a, s=None, axes=None)
+
+    The n-dimensional discrete Fourier transform of a real array a. A real
+    transform as rfft is performed along the axis specified by the last
+    element of axes, then complex transforms as fft are performed along the
+    other axes."""
+
+    a = asarray(a).astype(float)
+    s, axes = _cook_nd_args(a, s, axes)
+    a = rfft(a, s[-1], axes[-1])
+    for ii in range(len(axes)-1):
+        a = fft(a, s[ii], axes[ii])
+    return a
+
+def rfft2(a, s=None, axes=(-2,-1)):
+    """rfft2(a, s=None, axes=(-2,-1))
+
+    The 2d fft of the real valued array a. This is really just rfftn with
+    different default behavior."""
+
+    return rfftn(a, s, axes)
+
+def irfftn(a, s=None, axes=None):
+    """irfftn(a, s=None, axes=None)
+
+    The inverse of rfftn. The transform implemented in ifft is
+    applied along all axes but the last, then the transform implemented in
+    irfft is performed along the last axis. As with
+    irfft, the length of the result along that axis must be
+    specified if it is to be odd."""
+
+    a = asarray(a).astype(complex)
+    s, axes = _cook_nd_args(a, s, axes, invreal=1)
+    for ii in range(len(axes)-1):
+        a = ifft(a, s[ii], axes[ii])
+    a = irfft(a, s[-1], axes[-1])
+    return a
+
+def irfft2(a, s=None, axes=(-2,-1)):
+    """irfft2(a, s=None, axes=(-2, -1))
+
+    The inverse of rfft2. This is really just irfftn with
+    different default behavior."""
+
+    return irfftn(a, s, axes)
+
+# Deprecated names
+from numpy import deprecate
+refft = deprecate(rfft, 'refft', 'rfft')
+irefft = deprecate(irfft, 'irefft', 'irfft')
+refft2 = deprecate(rfft2, 'refft2', 'rfft2')
+irefft2 = deprecate(irfft2, 'irefft2', 'irfft2')
+refftn = deprecate(rfftn, 'refftn', 'rfftn')
+irefftn = deprecate(irfftn, 'irefftn', 'irfftn')
diff --git a/numpy/fft/fftpack_litemodule.c b/numpy/fft/fftpack_litemodule.c
new file mode 100644
index 000000000..f6d9aaf6b
--- /dev/null
+++ b/numpy/fft/fftpack_litemodule.c
@@ -0,0 +1,262 @@
+#include "fftpack.h"
+#include "Python.h"
+#include "numpy/arrayobject.h"
+
+static PyObject *ErrorObject;
+
+/* ----------------------------------------------------- */
+
+static char fftpack_cfftf__doc__[] ="";
+
+PyObject *
+fftpack_cfftf(PyObject *self, PyObject *args)
+{
+  PyObject *op1, *op2;
+  PyArrayObject *data;
+  double *wsave, *dptr;
+  int npts, nsave, nrepeats, i;
+
+  if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) return NULL;
+  data = (PyArrayObject *)PyArray_CopyFromObject(op1, PyArray_CDOUBLE, 1, 0);
+	if (data == NULL) return NULL;	
+  if (PyArray_As1D(&op2, (char **)&wsave, &nsave, PyArray_DOUBLE) == -1) 
+    goto fail;
+  if (data == NULL) goto fail;
+
+  npts = data->dimensions[data->nd-1];
+  if (nsave != npts*4+15) {
+    PyErr_SetString(ErrorObject, "invalid work array for fft size");
+    goto fail;
+  }
+
+  nrepeats = PyArray_SIZE(data)/npts;
+  dptr = (double *)data->data;
+  for (i=0; i<nrepeats; i++) {
+    cfftf(npts, dptr, wsave);
+    dptr += npts*2;
+  }
+  PyArray_Free(op2, (char *)wsave);
+  return (PyObject *)data;
+fail:
+  PyArray_Free(op2, (char *)wsave);
+  Py_DECREF(data);
+  return NULL;
+}
+
+static char fftpack_cfftb__doc__[] ="";
+
+PyObject *
+fftpack_cfftb(PyObject *self, PyObject *args)
+{
+  PyObject *op1, *op2;
+  PyArrayObject *data;
+  double *wsave, *dptr;
+  int npts, nsave, nrepeats, i;
+
+  if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) return NULL;
+  data = (PyArrayObject *)PyArray_CopyFromObject(op1, PyArray_CDOUBLE, 1, 0);
+	if (data == NULL) return NULL;	
+  if (PyArray_As1D(&op2, (char **)&wsave, &nsave, PyArray_DOUBLE) == -1) 
+    goto fail;
+  if (data == NULL) goto fail;
+
+  npts = data->dimensions[data->nd-1];
+  if (nsave != npts*4+15) {
+    PyErr_SetString(ErrorObject, "invalid work array for fft size");
+    goto fail;
+  }
+
+  nrepeats = PyArray_SIZE(data)/npts;
+  dptr = (double *)data->data;
+  for (i=0; i<nrepeats; i++) {
+    cfftb(npts, dptr, wsave);
+    dptr += npts*2;
+  }
+  PyArray_Free(op2, (char *)wsave);
+  return (PyObject *)data;
+fail:
+  PyArray_Free(op2, (char *)wsave);
+  Py_DECREF(data);
+  return NULL;
+}
+
+static char fftpack_cffti__doc__[] ="";
+
+static PyObject *
+fftpack_cffti(PyObject *self, PyObject *args)
+{
+  PyArrayObject *op;
+  int dim, n;
+
+  if (!PyArg_ParseTuple(args, "i", &n)) return NULL;
+
+  dim = 4*n+15; /*Magic size needed by cffti*/
+  /*Create a 1 dimensional array of dimensions of type double*/
+  op = (PyArrayObject *)PyArray_FromDims(1, &dim, PyArray_DOUBLE);
+  if (op == NULL) return NULL;
+
+  cffti(n, (double *)((PyArrayObject*)op)->data);
+
+  return (PyObject *)op;
+}
+
+static char fftpack_rfftf__doc__[] ="";
+
+PyObject *
+fftpack_rfftf(PyObject *self, PyObject *args)
+{
+  PyObject *op1, *op2;
+  PyArrayObject *data, *ret;
+  double *wsave, *dptr, *rptr;
+  int npts, nsave, nrepeats, i, rstep;
+
+  if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) return NULL;
+  data = (PyArrayObject *)PyArray_ContiguousFromObject(op1, PyArray_DOUBLE, 1, 0);
+  if (data == NULL) return NULL;
+  npts = data->dimensions[data->nd-1];
+  data->dimensions[data->nd-1] = npts/2+1;
+  ret = (PyArrayObject *)PyArray_Zeros(data->nd, data->dimensions, 
+					   PyArray_DescrFromType(PyArray_CDOUBLE), 0);
+  data->dimensions[data->nd-1] = npts;
+  rstep = (ret->dimensions[ret->nd-1])*2;
+
+  if (PyArray_As1D(&op2, (char **)&wsave, &nsave, PyArray_DOUBLE) == -1) 
+    goto fail;
+  if (data == NULL || ret == NULL) goto fail;
+
+  if (nsave != npts*2+15) {
+    PyErr_SetString(ErrorObject, "invalid work array for fft size");
+    goto fail;
+  }
+
+  nrepeats = PyArray_SIZE(data)/npts;
+  rptr = (double *)ret->data;
+  dptr = (double *)data->data;
+  
+  for (i=0; i<nrepeats; i++) {
+	memcpy((char *)(rptr+1), dptr, npts*sizeof(double));
+    rfftf(npts, rptr+1, wsave);
+	rptr[0] = rptr[1];
+	rptr[1] = 0.0;
+    rptr += rstep;
+	dptr += npts;
+  }
+  PyArray_Free(op2, (char *)wsave);
+  Py_DECREF(data);
+  return (PyObject *)ret;
+fail:
+  PyArray_Free(op2, (char *)wsave);
+  Py_XDECREF(data);
+  Py_XDECREF(ret);
+  return NULL;
+}
+
+static char fftpack_rfftb__doc__[] ="";
+
+
+PyObject *
+fftpack_rfftb(PyObject *self, PyObject *args)
+{
+  PyObject *op1, *op2;
+  PyArrayObject *data, *ret;
+  double *wsave, *dptr, *rptr;
+  int npts, nsave, nrepeats, i;
+
+  if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) return NULL;
+  data = (PyArrayObject *)PyArray_ContiguousFromObject(op1, PyArray_CDOUBLE, 1, 0);
+  if (data == NULL) return NULL;
+  npts = data->dimensions[data->nd-1];
+  ret = (PyArrayObject *)PyArray_Zeros(data->nd, data->dimensions, 
+					   PyArray_DescrFromType(PyArray_DOUBLE), 0);
+
+  if (PyArray_As1D(&op2, (char **)&wsave, &nsave, PyArray_DOUBLE) == -1) 
+    goto fail;
+  if (data == NULL || ret == NULL) goto fail;
+
+  if (nsave != npts*2+15) {
+    PyErr_SetString(ErrorObject, "invalid work array for fft size");
+    goto fail;
+  }
+
+  nrepeats = PyArray_SIZE(ret)/npts;
+  rptr = (double *)ret->data;
+  dptr = (double *)data->data;
+  
+  for (i=0; i<nrepeats; i++) {
+	memcpy((char *)(rptr+1), (dptr+2), (npts-1)*sizeof(double));
+	rptr[0] = dptr[0];
+    rfftb(npts, rptr, wsave);
+    rptr += npts;
+	dptr += npts*2;
+  }
+  PyArray_Free(op2, (char *)wsave);
+  Py_DECREF(data);
+  return (PyObject *)ret;
+fail:
+  PyArray_Free(op2, (char *)wsave);
+  Py_XDECREF(data);
+  Py_XDECREF(ret);
+  return NULL;
+}
+
+
+static char fftpack_rffti__doc__[] ="";
+
+static PyObject *
+fftpack_rffti(PyObject *self, PyObject *args)
+{
+  PyArrayObject *op;
+  int dim, n;
+
+  if (!PyArg_ParseTuple(args, "i", &n)) return NULL;
+
+  dim = 2*n+15; /*Magic size needed by rffti*/
+  /*Create a 1 dimensional array of dimensions of type double*/
+  op = (PyArrayObject *)PyArray_FromDims(1, &dim, PyArray_DOUBLE);
+  if (op == NULL) return NULL;
+
+  rffti(n, (double *)((PyArrayObject*)op)->data);
+
+  return (PyObject *)op;
+}
+
+
+/* List of methods defined in the module */
+
+static struct PyMethodDef fftpack_methods[] = {
+ {"cfftf",	fftpack_cfftf,	1,	fftpack_cfftf__doc__},
+ {"cfftb",	fftpack_cfftb,	1,	fftpack_cfftb__doc__},
+ {"cffti",	fftpack_cffti,	1,	fftpack_cffti__doc__},
+ {"rfftf",	fftpack_rfftf,	1,	fftpack_rfftf__doc__},
+ {"rfftb",	fftpack_rfftb,	1,	fftpack_rfftb__doc__},
+ {"rffti",	fftpack_rffti,	1,	fftpack_rffti__doc__},
+ {NULL,		NULL}		/* sentinel */
+};
+
+
+/* Initialization function for the module (*must* be called initfftpack) */
+
+static char fftpack_module_documentation[] = 
+""
+;
+
+PyMODINIT_FUNC initfftpack_lite(void)
+{
+	PyObject *m, *d;
+
+	/* Create the module and add the functions */
+	m = Py_InitModule4("fftpack_lite", fftpack_methods,
+		fftpack_module_documentation,
+		(PyObject*)NULL,PYTHON_API_VERSION);
+
+	/* Import the array object */
+	import_array();
+
+	/* Add some symbolic constants to the module */
+	d = PyModule_GetDict(m);
+	ErrorObject = PyErr_NewException("fftpack.error", NULL, NULL);
+	PyDict_SetItemString(d, "error", ErrorObject);
+
+	/* XXXX Add constants here */
+	
+}
diff --git a/numpy/fft/helper.py b/numpy/fft/helper.py
new file mode 100644
index 000000000..4e1da2abd
--- /dev/null
+++ b/numpy/fft/helper.py
@@ -0,0 +1,66 @@
+"""
+Discrete Fourier Transforms - helper.py
+"""
+# Created by Pearu Peterson, September 2002
+
+__all__ = ['fftshift','ifftshift','fftfreq']
+
+from numpy.core import asarray, concatenate, arange, take, \
+    array, integer
+import types
+
+def fftshift(x,axes=None):
+    """ fftshift(x, axes=None) -> y
+
+    Shift zero-frequency component to center of spectrum.
+
+    This function swaps half-spaces for all axes listed (defaults to all).
+
+    Notes:
+      If len(x) is even then the Nyquist component is y[0].
+    """
+    tmp = asarray(x)
+    ndim = len(tmp.shape)
+    if axes is None:
+        axes = range(ndim)
+    y = tmp
+    for k in axes:
+        n = tmp.shape[k]
+        p2 = (n+1)/2
+        mylist = concatenate((arange(p2,n),arange(p2)))
+        y = take(y,mylist,k)
+    return y
+
+
+def ifftshift(x,axes=None):
+    """ ifftshift(x,axes=None) - > y
+
+    Inverse of fftshift.
+    """
+    tmp = asarray(x)
+    ndim = len(tmp.shape)
+    if axes is None:
+        axes = range(ndim)
+    y = tmp
+    for k in axes:
+        n = tmp.shape[k]
+        p2 = n-(n+1)/2
+        mylist = concatenate((arange(p2,n),arange(p2)))
+        y = take(y,mylist,k)
+    return y
+
+def fftfreq(n,d=1.0):
+    """ fftfreq(n, d=1.0) -> f
+
+    DFT sample frequencies
+
+    The returned float array contains the frequency bins in
+    cycles/unit (with zero at the start) given a window length n and a
+    sample spacing d:
+
+      f = [0,1,...,n/2-1,-n/2,...,-1]/(d*n)         if n is even
+      f = [0,1,...,(n-1)/2,-(n-1)/2,...,-1]/(d*n)   if n is odd
+    """
+    assert isinstance(n,types.IntType) or isinstance(n, integer)
+    k = range(0,(n-1)/2+1)+range(-(n/2),0)
+    return array(k,'d')/(n*d)
diff --git a/numpy/fft/info.py b/numpy/fft/info.py
new file mode 100644
index 000000000..4aa9b7d2c
--- /dev/null
+++ b/numpy/fft/info.py
@@ -0,0 +1,29 @@
+"""\
+Core FFT routines
+==================
+
+ Standard FFTs
+
+   fft
+   ifft
+   fft2
+   ifft2
+   fftn
+   ifftn
+
+ Real FFTs
+
+   refft
+   irefft
+   refft2
+   irefft2
+   refftn
+   irefftn
+
+ Hermite FFTs
+
+   hfft
+   ihfft
+"""
+
+depends = ['core']
diff --git a/numpy/fft/old.py b/numpy/fft/old.py
new file mode 100644
index 000000000..86723af95
--- /dev/null
+++ b/numpy/fft/old.py
@@ -0,0 +1,19 @@
+
+__all__ = ['fft', 'fft2d', 'fftnd', 'hermite_fft', 'inverse_fft', 'inverse_fft2d',
+           'inverse_fftnd', 'inverse_hermite_fft', 'inverse_real_fft', 'inverse_real_fft2d',
+           'inverse_real_fftnd', 'real_fft', 'real_fft2d', 'real_fftnd']
+
+from fftpack import fft
+from fftpack import fft2 as fft2d
+from fftpack import fftn as fftnd
+from fftpack import hfft as hermite_fft
+from fftpack import ifft as inverse_fft
+from fftpack import ifft2 as inverse_fft2d
+from fftpack import ifftn as inverse_fftnd
+from fftpack import ihfft as inverse_hermite_fft
+from fftpack import irfft as inverse_real_fft
+from fftpack import irfft2 as inverse_real_fft2d
+from fftpack import irfftn as inverse_real_fftnd
+from fftpack import rfft as real_fft
+from fftpack import rfft2 as real_fft2d
+from fftpack import rfftn as real_fftnd
diff --git a/numpy/fft/setup.py b/numpy/fft/setup.py
new file mode 100644
index 000000000..4d1de069d
--- /dev/null
+++ b/numpy/fft/setup.py
@@ -0,0 +1,20 @@
+
+from os.path import join
+
+def configuration(parent_package='',top_path=None):
+    from numpy.distutils.misc_util import Configuration
+    config = Configuration('fft',parent_package,top_path)
+
+    config.add_data_dir('tests')
+
+    # Configure fftpack_lite
+    config.add_extension('fftpack_lite',
+                         sources=['fftpack_litemodule.c', 'fftpack.c']
+                         )
+
+
+    return config
+
+if __name__ == '__main__':
+    from numpy.distutils.core import setup
+    setup(**configuration(top_path='').todict())
diff --git a/numpy/fft/tests/test_helper.py b/numpy/fft/tests/test_helper.py
new file mode 100644
index 000000000..3d02c01df
--- /dev/null
+++ b/numpy/fft/tests/test_helper.py
@@ -0,0 +1,45 @@
+#!/usr/bin/env python
+# Copied from fftpack.helper by Pearu Peterson, October 2005
+""" Test functions for fftpack.helper module
+"""
+
+import sys
+from numpy.testing import *
+set_package_path()
+from numpy.fft import fftshift,ifftshift,fftfreq
+del sys.path[0]
+
+from numpy import pi
+
+def random(size):
+    return rand(*size)
+
+class test_fftshift(NumpyTestCase):
+
+    def check_definition(self):
+        x = [0,1,2,3,4,-4,-3,-2,-1]
+        y = [-4,-3,-2,-1,0,1,2,3,4]
+        assert_array_almost_equal(fftshift(x),y)
+        assert_array_almost_equal(ifftshift(y),x)
+        x = [0,1,2,3,4,-5,-4,-3,-2,-1]
+        y = [-5,-4,-3,-2,-1,0,1,2,3,4]
+        assert_array_almost_equal(fftshift(x),y)
+        assert_array_almost_equal(ifftshift(y),x)
+
+    def check_inverse(self):
+        for n in [1,4,9,100,211]:
+            x = random((n,))
+            assert_array_almost_equal(ifftshift(fftshift(x)),x)
+
+class test_fftfreq(NumpyTestCase):
+
+    def check_definition(self):
+        x = [0,1,2,3,4,-4,-3,-2,-1]
+        assert_array_almost_equal(9*fftfreq(9),x)
+        assert_array_almost_equal(9*pi*fftfreq(9,pi),x)
+        x = [0,1,2,3,4,-5,-4,-3,-2,-1]
+        assert_array_almost_equal(10*fftfreq(10),x)
+        assert_array_almost_equal(10*pi*fftfreq(10,pi),x)
+
+if __name__ == "__main__":
+    NumpyTest().run()