Merge pull request #11888 from mreineck/add_pocketfft

WIP: Add pocketfft sources to numpy for testing, benchmarks, etc.

12 files changed, 2514 insertions, 2179 deletions

diff --git a/doc/release/1.17.0-notes.rst b/doc/release/1.17.0-notes.rst
index ec6f63945..5c58632c0 100644
--- a/doc/release/1.17.0-notes.rst
+++ b/doc/release/1.17.0-notes.rst
@@ -6,6 +6,7 @@ NumPy 1.17.0 Release Notes
 Highlights
 ==========
 
+* NumPy's FFT implementation has switched to pocketfft
 
 New functions
 =============
@@ -38,6 +39,15 @@ New Features
 Improvements
 ============
 
+replacement of the `fftpack`-based FFT module by the `pocketfft` library
+------------------------------------------------------------------------
+Both implementations have the same ancestor (Fortran77 `FFTPACK` by Paul N.
+Swarztrauber), but `pocketfft` contains additional modifications which
+improve both accuracy and performance in some circumstances. For FFT lengths
+containing large prime factors, `pocketfft` uses Bluestein's algorithm, which
+maintains `O(N log N)` run time complexity instead of deteriorating towards
+`O(N*N)` for prime lengths. Also, accuracy for real-valued FFTs with near-prime
+lengths has improved and is on par with complex-valued FFTs.
 
 Changes
 =======
diff --git a/numpy/fft/README.md b/numpy/fft/README.md
new file mode 100644
index 000000000..7040a2e9b
--- /dev/null
+++ b/numpy/fft/README.md
@@ -0,0 +1,53 @@
+PocketFFT
+---------
+
+This is a heavily modified implementation of FFTPack [1,2], with the following
+advantages:
+
+- strictly C99 compliant
+- more accurate twiddle factor computation
+- very fast plan generation
+- worst case complexity for transform sizes with large prime factors is
+  `N*log(N)`, because Bluestein's algorithm [3] is used for these cases.
+
+License
+-------
+
+3-clause BSD (see LICENSE.md)
+
+
+Some code details
+-----------------
+
+Twiddle factor computation:
+
+- making use of symmetries to reduce number of sin/cos evaluations
+- all angles are reduced to the range `[0; pi/4]` for higher accuracy
+- an adapted implementation of `sincospi()` is used, which actually computes
+  `sin(x)` and `(cos(x)-1)`.
+- if `n` sin/cos pairs are required, the adjusted `sincospi()` is only called
+  `2*sqrt(n)` times; the remaining values are obtained by evaluating the
+  angle addition theorems in a numerically accurate way.
+
+Parallel invocation:
+
+- Plans only contain read-only data; all temporary arrays are allocated and
+  deallocated during an individual FFT execution. This means that a single plan
+  can be used in several threads at the same time.
+
+Efficient codelets are available for the factors:
+
+- 2, 3, 4, 5, 7, 11 for complex-valued FFTs
+- 2, 3, 4, 5 for real-valued FFTs
+
+Larger prime factors are handled by somewhat less efficient, generic routines.
+
+For lengths with very large prime factors, Bluestein's algorithm is used, and
+instead of an FFT of length `n`, a convolution of length `n2 >= 2*n-1`
+is performed, where `n2` is chosen to be highly composite.
+
+
+[1] Swarztrauber, P. 1982, Vectorizing the Fast Fourier Transforms
+    (New York: Academic Press), 51
+[2] https://www.netlib.org/fftpack/
+[3] https://en.wikipedia.org/wiki/Chirp_Z-transform
diff --git a/numpy/fft/__init__.py b/numpy/fft/__init__.py
index 44243b483..64b35bc19 100644
--- a/numpy/fft/__init__.py
+++ b/numpy/fft/__init__.py
@@ -3,7 +3,7 @@ from __future__ import division, absolute_import, print_function
 # To get sub-modules
 from .info import __doc__
 
-from .fftpack import *
+from .pocketfft import *
 from .helper import *
 
 from numpy._pytesttester import PytestTester
diff --git a/numpy/fft/fftpack.c b/numpy/fft/fftpack.c
deleted file mode 100644
index 07fa2bf4c..000000000
--- a/numpy/fft/fftpack.c
+++ /dev/null
@@ -1,1536 +0,0 @@
-/*
- * fftpack.c : A set of FFT routines in C.
- * Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber (Version 4, 1985).
-*/
-#define NPY_NO_DEPRECATED_API NPY_API_VERSION
-
-#include <Python.h>
-#include <math.h>
-#include <stdio.h>
-#include <numpy/ndarraytypes.h>
-
-#define DOUBLE
-#ifdef DOUBLE
-#define Treal double
-#else
-#define Treal float
-#endif
-
-#define ref(u,a) u[a]
-
-/* Macros for accurate calculation of the twiddle factors. */
-#define TWOPI 6.283185307179586476925286766559005768391
-#define cos2pi(m, n) cos((TWOPI * (m)) / (n))
-#define sin2pi(m, n) sin((TWOPI * (m)) / (n))
-
-#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
-
-static void sincos2pi(int m, int n, Treal* si, Treal* co)
-/* Calculates sin(2pi * m/n) and cos(2pi * m/n). It is more accurate
- * than the naive calculation as the fraction m/n is reduced to [0, 1/8) first.
- * Due to the symmetry of sin(x) and cos(x) the values for all x can be
- * determined from the function values of the reduced argument in the first
- * octant.
- */
-    {
-        int n8, m8, octant;
-        n8 = 8 * n;
-        m8 = (8 * m) % n8;
-        octant = m8 / n;
-        m8 = m8 % n;
-        switch(octant) {
-            case 0:
-                *co = cos2pi(m8, n8);
-                *si = sin2pi(m8, n8);
-                break;
-            case 1:
-                *co = sin2pi(n-m8, n8);
-                *si = cos2pi(n-m8, n8);
-                break;
-            case 2:
-                *co = -sin2pi(m8, n8);
-                *si = cos2pi(m8, n8);
-                break;
-            case 3:
-                *co = -cos2pi(n-m8, n8);
-                *si = sin2pi(n-m8, n8);
-                break;
-            case 4:
-                *co = -cos2pi(m8, n8);
-                *si = -sin2pi(m8, n8);
-                break;
-            case 5:
-                *co = -sin2pi(n-m8, n8);
-                *si = -cos2pi(n-m8, n8);
-                break;
-            case 6:
-                *co = sin2pi(m8, n8);
-                *si = -cos2pi(m8, n8);
-                break;
-            case 7:
-                *co = cos2pi(n-m8, n8);
-                *si = -sin2pi(n-m8, n8);
-                break;
-        }
-    }
-
-/* ----------------------------------------------------------------------
-   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.86602540378443864676;
-    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.3090169943749474241;
-    static const Treal ti11 = 0.95105651629515357212;
-    static const Treal tr12 = -0.8090169943749474241;
-    static const Treal ti12 = 0.58778525229247312917;
-    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, 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.86602540378443864676;
-    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.86602540378443864676;
-    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.70710678118654752440;
-    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.41421356237309504880;
-    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.3090169943749474241;
-    static const Treal ti11 = 0.95105651629515357212;
-    static const Treal tr12 = -0.8090169943749474241;
-    static const Treal ti12 = 0.58778525229247312917;
-    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.3090169943749474241;
-    static const Treal ti11 = 0.95105651629515357212;
-    static const Treal tr12 = -0.8090169943749474241;
-    static const Treal ti12 = 0.58778525229247312917;
-    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[])
-  {
-    int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;    
-    Treal dc2, ai1, ai2, ar1, ar2, ds2, dcp, dsp, ar1h, ar2h;
-    sincos2pi(1, ip, &dsp, &dcp);
-    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[])
-  {
-    int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
-    Treal dc2, ai1, ai2, ar1, ar2, ds2;
-    int nbd;
-    Treal dcp, dsp, ar1h, ar2h;
-    sincos2pi(1, ip, &dsp, &dcp);
-    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, npy_cfftf, npy_cfftb, cffti1, npy_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 */
-
-
-NPY_VISIBILITY_HIDDEN void npy_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);
-  } /* npy_cfftf */
-
-
-NPY_VISIBILITY_HIDDEN void npy_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);
-  } /* npy_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])
-  {
-    int fi, 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];
-    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;
-        for (ii=4; ii<=idot; ii+=2) {
-          i+= 2;
-          fi+= 1;
-          sincos2pi(fi*ld, n, wa+i, wa+i-1);
-        }
-        if (ip > 5) {
-          wa[i1-1] = wa[i-1];
-          wa[i1] = wa[i];
-        }
-      }
-      l1 = l2;
-    }
-  } /* cffti1 */
-
-
-NPY_VISIBILITY_HIDDEN void npy_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));
-  } /* npy_cffti */
-
-  /* -------------------------------------------------------------------
-rfftf1, rfftb1, npy_rfftf, npy_rfftb, rffti1, npy_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 */
-
-
-static 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 */
-
-
-NPY_VISIBILITY_HIDDEN void npy_rfftf(int n, Treal r[], Treal wsave[])
-  {
-    if (n == 1) return;
-    rfftf1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
-  } /* npy_rfftf */
-
-
-NPY_VISIBILITY_HIDDEN void npy_rfftb(int n, Treal r[], Treal wsave[])
-  {
-    if (n == 1) return;
-    rfftb1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
-  } /* npy_rfftb */
-
-
-static void rffti1(int n, Treal wa[], int ifac[MAXFAC+2])
-  {
-    int fi, 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];
-    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;
-        fi = 0;
-        for (ii = 3; ii <= ido; ii += 2) {
-          i += 2;
-          fi += 1;
-          sincos2pi(fi*ld, n, wa+i-1, wa+i-2);
-        }
-        is += ido;
-      }
-      l1 = l2;
-    }
-  } /* rffti1 */
-
-
-NPY_VISIBILITY_HIDDEN void npy_rffti(int n, Treal wsave[])
-  {
-    if (n == 1) return;
-    rffti1(n, wsave+n, (int*)(wsave+2*n));
-  } /* npy_rffti */
-
-#ifdef __cplusplus
-}
-#endif
diff --git a/numpy/fft/fftpack.h b/numpy/fft/fftpack.h
deleted file mode 100644
index 5e8f4631c..000000000
--- a/numpy/fft/fftpack.h
+++ /dev/null
@@ -1,28 +0,0 @@
-/*
- * 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 NPY_VISIBILITY_HIDDEN void npy_cfftf(int N, Treal data[], const Treal wrk[]);
-extern NPY_VISIBILITY_HIDDEN void npy_cfftb(int N, Treal data[], const Treal wrk[]);
-extern NPY_VISIBILITY_HIDDEN void npy_cffti(int N, Treal wrk[]);
-
-extern NPY_VISIBILITY_HIDDEN void npy_rfftf(int N, Treal data[], const Treal wrk[]);
-extern NPY_VISIBILITY_HIDDEN void npy_rfftb(int N, Treal data[], const Treal wrk[]);
-extern NPY_VISIBILITY_HIDDEN void npy_rffti(int N, Treal wrk[]);
-
-#ifdef __cplusplus
-}
-#endif
diff --git a/numpy/fft/fftpack_litemodule.c b/numpy/fft/fftpack_litemodule.c
deleted file mode 100644
index bd6cfc120..000000000
--- a/numpy/fft/fftpack_litemodule.c
+++ /dev/null
@@ -1,366 +0,0 @@
-#define NPY_NO_DEPRECATED_API NPY_API_VERSION
-
-#include "Python.h"
-#include "numpy/arrayobject.h"
-#include "fftpack.h"
-
-static PyObject *ErrorObject;
-
-static const char fftpack_cfftf__doc__[] = "";
-
-static PyObject *
-fftpack_cfftf(PyObject *NPY_UNUSED(self), PyObject *args)
-{
-    PyObject *op1, *op2;
-    PyArrayObject *data;
-    PyArray_Descr *descr;
-    double *wsave, *dptr;
-    npy_intp nsave;
-    int npts, nrepeats, i;
-
-    if(!PyArg_ParseTuple(args, "OO:cfftf", &op1, &op2)) {
-        return NULL;
-    }
-    data = (PyArrayObject *)PyArray_CopyFromObject(op1,
-            NPY_CDOUBLE, 1, 0);
-    if (data == NULL) {
-        return NULL;
-    }
-    descr = PyArray_DescrFromType(NPY_DOUBLE);
-    if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) {
-        goto fail;
-    }
-    if (data == NULL) {
-        goto fail;
-    }
-
-    npts = PyArray_DIM(data, PyArray_NDIM(data) - 1);
-    if (nsave != npts*4 + 15) {
-        PyErr_SetString(ErrorObject, "invalid work array for fft size");
-        goto fail;
-    }
-
-    nrepeats = PyArray_SIZE(data)/npts;
-    dptr = (double *)PyArray_DATA(data);
-    Py_BEGIN_ALLOW_THREADS;
-    NPY_SIGINT_ON;
-    for (i = 0; i < nrepeats; i++) {
-        npy_cfftf(npts, dptr, wsave);
-        dptr += npts*2;
-    }
-    NPY_SIGINT_OFF;
-    Py_END_ALLOW_THREADS;
-    PyArray_Free(op2, (char *)wsave);
-    return (PyObject *)data;
-
-fail:
-    PyArray_Free(op2, (char *)wsave);
-    Py_DECREF(data);
-    return NULL;
-}
-
-static const char fftpack_cfftb__doc__[] = "";
-
-static PyObject *
-fftpack_cfftb(PyObject *NPY_UNUSED(self), PyObject *args)
-{
-    PyObject *op1, *op2;
-    PyArrayObject *data;
-    PyArray_Descr *descr;
-    double *wsave, *dptr;
-    npy_intp nsave;
-    int npts, nrepeats, i;
-
-    if(!PyArg_ParseTuple(args, "OO:cfftb", &op1, &op2)) {
-        return NULL;
-    }
-    data = (PyArrayObject *)PyArray_CopyFromObject(op1,
-            NPY_CDOUBLE, 1, 0);
-    if (data == NULL) {
-        return NULL;
-    }
-    descr = PyArray_DescrFromType(NPY_DOUBLE);
-    if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) {
-        goto fail;
-    }
-    if (data == NULL) {
-        goto fail;
-    }
-
-    npts = PyArray_DIM(data, PyArray_NDIM(data) - 1);
-    if (nsave != npts*4 + 15) {
-        PyErr_SetString(ErrorObject, "invalid work array for fft size");
-        goto fail;
-    }
-
-    nrepeats = PyArray_SIZE(data)/npts;
-    dptr = (double *)PyArray_DATA(data);
-    Py_BEGIN_ALLOW_THREADS;
-    NPY_SIGINT_ON;
-    for (i = 0; i < nrepeats; i++) {
-        npy_cfftb(npts, dptr, wsave);
-        dptr += npts*2;
-    }
-    NPY_SIGINT_OFF;
-    Py_END_ALLOW_THREADS;
-    PyArray_Free(op2, (char *)wsave);
-    return (PyObject *)data;
-
-fail:
-    PyArray_Free(op2, (char *)wsave);
-    Py_DECREF(data);
-    return NULL;
-}
-
-static const char fftpack_cffti__doc__[] = "";
-
-static PyObject *
-fftpack_cffti(PyObject *NPY_UNUSED(self), PyObject *args)
-{
-    PyArrayObject *op;
-    npy_intp dim;
-    long n;
-
-    if (!PyArg_ParseTuple(args, "l:cffti", &n)) {
-        return NULL;
-    }
-    /*Magic size needed by npy_cffti*/
-    dim = 4*n + 15;
-    /*Create a 1 dimensional array of dimensions of type double*/
-    op = (PyArrayObject *)PyArray_SimpleNew(1, &dim, NPY_DOUBLE);
-    if (op == NULL) {
-        return NULL;
-    }
-
-    Py_BEGIN_ALLOW_THREADS;
-    NPY_SIGINT_ON;
-    npy_cffti(n, (double *)PyArray_DATA((PyArrayObject*)op));
-    NPY_SIGINT_OFF;
-    Py_END_ALLOW_THREADS;
-
-    return (PyObject *)op;
-}
-
-static const char fftpack_rfftf__doc__[] = "";
-
-static PyObject *
-fftpack_rfftf(PyObject *NPY_UNUSED(self), PyObject *args)
-{
-    PyObject *op1, *op2;
-    PyArrayObject *data, *ret;
-    PyArray_Descr *descr;
-    double *wsave = NULL, *dptr, *rptr;
-    npy_intp nsave;
-    int npts, nrepeats, i, rstep;
-
-    if(!PyArg_ParseTuple(args, "OO:rfftf", &op1, &op2)) {
-        return NULL;
-    }
-    data = (PyArrayObject *)PyArray_ContiguousFromObject(op1,
-            NPY_DOUBLE, 1, 0);
-    if (data == NULL) {
-        return NULL;
-    }
-    /* FIXME, direct access changing contents of data->dimensions */
-    npts = PyArray_DIM(data, PyArray_NDIM(data) - 1);
-    PyArray_DIMS(data)[PyArray_NDIM(data) - 1] = npts/2 + 1;
-    ret = (PyArrayObject *)PyArray_Zeros(PyArray_NDIM(data),
-            PyArray_DIMS(data), PyArray_DescrFromType(NPY_CDOUBLE), 0);
-    if (ret == NULL) {
-        goto fail;
-    }
-    PyArray_DIMS(data)[PyArray_NDIM(data) - 1] = npts;
-    rstep = PyArray_DIM(ret, PyArray_NDIM(ret) - 1)*2;
-
-    descr = PyArray_DescrFromType(NPY_DOUBLE);
-    if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -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 *)PyArray_DATA(ret);
-    dptr = (double *)PyArray_DATA(data);
-
-    Py_BEGIN_ALLOW_THREADS;
-    NPY_SIGINT_ON;
-    for (i = 0; i < nrepeats; i++) {
-        memcpy((char *)(rptr+1), dptr, npts*sizeof(double));
-        npy_rfftf(npts, rptr+1, wsave);
-        rptr[0] = rptr[1];
-        rptr[1] = 0.0;
-        rptr += rstep;
-        dptr += npts;
-    }
-    NPY_SIGINT_OFF;
-    Py_END_ALLOW_THREADS;
-    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 const char fftpack_rfftb__doc__[] = "";
-
-static PyObject *
-fftpack_rfftb(PyObject *NPY_UNUSED(self), PyObject *args)
-{
-    PyObject *op1, *op2;
-    PyArrayObject *data, *ret;
-    PyArray_Descr *descr;
-    double *wsave, *dptr, *rptr;
-    npy_intp nsave;
-    int npts, nrepeats, i;
-
-    if(!PyArg_ParseTuple(args, "OO:rfftb", &op1, &op2)) {
-        return NULL;
-    }
-    data = (PyArrayObject *)PyArray_ContiguousFromObject(op1,
-            NPY_CDOUBLE, 1, 0);
-    if (data == NULL) {
-        return NULL;
-    }
-    npts = PyArray_DIM(data, PyArray_NDIM(data) - 1);
-    ret = (PyArrayObject *)PyArray_Zeros(PyArray_NDIM(data), PyArray_DIMS(data),
-            PyArray_DescrFromType(NPY_DOUBLE), 0);
-
-    descr = PyArray_DescrFromType(NPY_DOUBLE);
-    if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -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 *)PyArray_DATA(ret);
-    dptr = (double *)PyArray_DATA(data);
-
-    Py_BEGIN_ALLOW_THREADS;
-    NPY_SIGINT_ON;
-    for (i = 0; i < nrepeats; i++) {
-        memcpy((char *)(rptr + 1), (dptr + 2), (npts - 1)*sizeof(double));
-        rptr[0] = dptr[0];
-        npy_rfftb(npts, rptr, wsave);
-        rptr += npts;
-        dptr += npts*2;
-    }
-    NPY_SIGINT_OFF;
-    Py_END_ALLOW_THREADS;
-    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 const char fftpack_rffti__doc__[] = "";
-
-static PyObject *
-fftpack_rffti(PyObject *NPY_UNUSED(self), PyObject *args)
-{
-  PyArrayObject *op;
-  npy_intp dim;
-  long n;
-
-  if (!PyArg_ParseTuple(args, "l:rffti", &n)) {
-      return NULL;
-  }
-  /*Magic size needed by npy_rffti*/
-  dim = 2*n + 15;
-  /*Create a 1 dimensional array of dimensions of type double*/
-  op = (PyArrayObject *)PyArray_SimpleNew(1, &dim, NPY_DOUBLE);
-  if (op == NULL) {
-      return NULL;
-  }
-  Py_BEGIN_ALLOW_THREADS;
-  NPY_SIGINT_ON;
-  npy_rffti(n, (double *)PyArray_DATA((PyArrayObject*)op));
-  NPY_SIGINT_OFF;
-  Py_END_ALLOW_THREADS;
-
-  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, 0, NULL}          /* sentinel */
-};
-
-#if PY_MAJOR_VERSION >= 3
-static struct PyModuleDef moduledef = {
-        PyModuleDef_HEAD_INIT,
-        "fftpack_lite",
-        NULL,
-        -1,
-        fftpack_methods,
-        NULL,
-        NULL,
-        NULL,
-        NULL
-};
-#endif
-
-/* Initialization function for the module */
-#if PY_MAJOR_VERSION >= 3
-#define RETVAL(x) x
-PyMODINIT_FUNC PyInit_fftpack_lite(void)
-#else
-#define RETVAL(x)
-PyMODINIT_FUNC
-initfftpack_lite(void)
-#endif
-{
-    PyObject *m,*d;
-#if PY_MAJOR_VERSION >= 3
-    m = PyModule_Create(&moduledef);
-#else
-    static const char fftpack_module_documentation[] = "";
-
-    m = Py_InitModule4("fftpack_lite", fftpack_methods,
-            fftpack_module_documentation,
-            (PyObject*)NULL,PYTHON_API_VERSION);
-#endif
-    if (m == NULL) {
-        return RETVAL(NULL);
-    }
-
-    /* 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 */
-
-    return RETVAL(m);
-}
diff --git a/numpy/fft/helper.py b/numpy/fft/helper.py
index 9985f6d4c..a920a4ac0 100644
--- a/numpy/fft/helper.py
+++ b/numpy/fft/helper.py
@@ -4,11 +4,6 @@ Discrete Fourier Transforms - helper.py
 """
 from __future__ import division, absolute_import, print_function
 
-import collections
-try:
-    import threading
-except ImportError:
-    import dummy_threading as threading
 from numpy.compat import integer_types
 from numpy.core import integer, empty, arange, asarray, roll
 from numpy.core.overrides import array_function_dispatch, set_module
@@ -227,99 +222,3 @@ def rfftfreq(n, d=1.0):
     N = n//2 + 1
     results = arange(0, N, dtype=int)
     return results * val
-
-
-class _FFTCache(object):
-    """
-    Cache for the FFT twiddle factors as an LRU (least recently used) cache.
-
-    Parameters
-    ----------
-    max_size_in_mb : int
-        Maximum memory usage of the cache before items are being evicted.
-    max_item_count : int
-        Maximum item count of the cache before items are being evicted.
-
-    Notes
-    -----
-    Items will be evicted if either limit has been reached upon getting and
-    setting. The maximum memory usages is not strictly the given
-    ``max_size_in_mb`` but rather
-    ``max(max_size_in_mb, 1.5 * size_of_largest_item)``. Thus the cache will
-    never be completely cleared - at least one item will remain and a single
-    large item can cause the cache to retain several smaller items even if the
-    given maximum cache size has been exceeded.
-    """
-    def __init__(self, max_size_in_mb, max_item_count):
-        self._max_size_in_bytes = max_size_in_mb * 1024 ** 2
-        self._max_item_count = max_item_count
-        self._dict = collections.OrderedDict()
-        self._lock = threading.Lock()
-
-    def put_twiddle_factors(self, n, factors):
-        """
-        Store twiddle factors for an FFT of length n in the cache.
-
-        Putting multiple twiddle factors for a certain n will store it multiple
-        times.
-
-        Parameters
-        ----------
-        n : int
-            Data length for the FFT.
-        factors : ndarray
-            The actual twiddle values.
-        """
-        with self._lock:
-            # Pop + later add to move it to the end for LRU behavior.
-            # Internally everything is stored in a dictionary whose values are
-            # lists.
-            try:
-                value = self._dict.pop(n)
-            except KeyError:
-                value = []
-            value.append(factors)
-            self._dict[n] = value
-            self._prune_cache()
-
-    def pop_twiddle_factors(self, n):
-        """
-        Pop twiddle factors for an FFT of length n from the cache.
-
-        Will return None if the requested twiddle factors are not available in
-        the cache.
-
-        Parameters
-        ----------
-        n : int
-            Data length for the FFT.
-
-        Returns
-        -------
-        out : ndarray or None
-            The retrieved twiddle factors if available, else None.
-        """
-        with self._lock:
-            if n not in self._dict or not self._dict[n]:
-                return None
-            # Pop + later add to move it to the end for LRU behavior.
-            all_values = self._dict.pop(n)
-            value = all_values.pop()
-            # Only put pack if there are still some arrays left in the list.
-            if all_values:
-                self._dict[n] = all_values
-            return value
-
-    def _prune_cache(self):
-        # Always keep at least one item.
-        while len(self._dict) > 1 and (
-                len(self._dict) > self._max_item_count or self._check_size()):
-            self._dict.popitem(last=False)
-
-    def _check_size(self):
-        item_sizes = [sum(_j.nbytes for _j in _i)
-                      for _i in self._dict.values() if _i]
-        if not item_sizes:
-            return False
-        max_size = max(self._max_size_in_bytes, 1.5 * max(item_sizes))
-        return sum(item_sizes) > max_size
diff --git a/numpy/fft/pocketfft.c b/numpy/fft/pocketfft.c
new file mode 100644
index 000000000..10a741b6f
--- /dev/null
+++ b/numpy/fft/pocketfft.c
@@ -0,0 +1,2398 @@
+/*
+ * This file is part of pocketfft.
+ * Licensed under a 3-clause BSD style license - see LICENSE.md
+ */
+
+/*
+ *  Main implementation file.
+ *
+ *  Copyright (C) 2004-2018 Max-Planck-Society
+ *  \author Martin Reinecke
+ */
+
+#include <math.h>
+#include <string.h>
+#include <stdlib.h>
+
+#include "npy_config.h"
+#define restrict NPY_RESTRICT
+
+#define RALLOC(type,num) \
+  ((type *)malloc((num)*sizeof(type)))
+#define DEALLOC(ptr) \
+  do { free(ptr); (ptr)=NULL; } while(0)
+
+#define SWAP(a,b,type) \
+  do { type tmp_=(a); (a)=(b); (b)=tmp_; } while(0)
+
+#ifdef __GNUC__
+#define NOINLINE __attribute__((noinline))
+#define WARN_UNUSED_RESULT __attribute__ ((warn_unused_result))
+#else
+#define NOINLINE
+#define WARN_UNUSED_RESULT
+#endif
+
+struct cfft_plan_i;
+typedef struct cfft_plan_i * cfft_plan;
+struct rfft_plan_i;
+typedef struct rfft_plan_i * rfft_plan;
+
+// adapted from https://stackoverflow.com/questions/42792939/
+// CAUTION: this function only works for arguments in the range [-0.25; 0.25]!
+static void my_sincosm1pi (double a, double *restrict res)
+  {
+  double s = a * a;
+  /* Approximate cos(pi*x)-1 for x in [-0.25,0.25] */
+  double r =     -1.0369917389758117e-4;
+  r = fma (r, s,  1.9294935641298806e-3);
+  r = fma (r, s, -2.5806887942825395e-2);
+  r = fma (r, s,  2.3533063028328211e-1);
+  r = fma (r, s, -1.3352627688538006e+0);
+  r = fma (r, s,  4.0587121264167623e+0);
+  r = fma (r, s, -4.9348022005446790e+0);
+  double c = r*s;
+  /* Approximate sin(pi*x) for x in [-0.25,0.25] */
+  r =             4.6151442520157035e-4;
+  r = fma (r, s, -7.3700183130883555e-3);
+  r = fma (r, s,  8.2145868949323936e-2);
+  r = fma (r, s, -5.9926452893214921e-1);
+  r = fma (r, s,  2.5501640398732688e+0);
+  r = fma (r, s, -5.1677127800499516e+0);
+  s = s * a;
+  r = r * s;
+  s = fma (a, 3.1415926535897931e+0, r);
+  res[0] = c;
+  res[1] = s;
+  }
+
+NOINLINE static void calc_first_octant(size_t den, double * restrict res)
+  {
+  size_t n = (den+4)>>3;
+  if (n==0) return;
+  res[0]=1.; res[1]=0.;
+  if (n==1) return;
+  size_t l1=(size_t)sqrt(n);
+  for (size_t i=1; i<l1; ++i)
+    my_sincosm1pi((2.*i)/den,&res[2*i]);
+  size_t start=l1;
+  while(start<n)
+    {
+    double cs[2];
+    my_sincosm1pi((2.*start)/den,cs);
+    res[2*start] = cs[0]+1.;
+    res[2*start+1] = cs[1];
+    size_t end = l1;
+    if (start+end>n) end = n-start;
+    for (size_t i=1; i<end; ++i)
+      {
+      double csx[2]={res[2*i], res[2*i+1]};
+      res[2*(start+i)] = ((cs[0]*csx[0] - cs[1]*csx[1] + cs[0]) + csx[0]) + 1.;
+      res[2*(start+i)+1] = (cs[0]*csx[1] + cs[1]*csx[0]) + cs[1] + csx[1];
+      }
+    start += l1;
+    }
+  for (size_t i=1; i<l1; ++i)
+    res[2*i] += 1.;
+  }
+
+NOINLINE static void calc_first_quadrant(size_t n, double * restrict res)
+  {
+  double * restrict p = res+n;
+  calc_first_octant(n<<1, p);
+  size_t ndone=(n+2)>>2;
+  size_t i=0, idx1=0, idx2=2*ndone-2;
+  for (; i+1<ndone; i+=2, idx1+=2, idx2-=2)
+    {
+    res[idx1]   = p[2*i];
+    res[idx1+1] = p[2*i+1];
+    res[idx2]   = p[2*i+3];
+    res[idx2+1] = p[2*i+2];
+    }
+  if (i!=ndone)
+    {
+    res[idx1  ] = p[2*i];
+    res[idx1+1] = p[2*i+1];
+    }
+  }
+
+NOINLINE static void calc_first_half(size_t n, double * restrict res)
+  {
+  int ndone=(n+1)>>1;
+  double * p = res+n-1;
+  calc_first_octant(n<<2, p);
+  int i4=0, in=n, i=0;
+  for (; i4<=in-i4; ++i, i4+=4) // octant 0
+    {
+    res[2*i] = p[2*i4]; res[2*i+1] = p[2*i4+1];
+    }
+  for (; i4-in <= 0; ++i, i4+=4) // octant 1
+    {
+    int xm = in-i4;
+    res[2*i] = p[2*xm+1]; res[2*i+1] = p[2*xm];
+    }
+  for (; i4<=3*in-i4; ++i, i4+=4) // octant 2
+    {
+    int xm = i4-in;
+    res[2*i] = -p[2*xm+1]; res[2*i+1] = p[2*xm];
+    }
+  for (; i<ndone; ++i, i4+=4) // octant 3
+    {
+    int xm = 2*in-i4;
+    res[2*i] = -p[2*xm]; res[2*i+1] = p[2*xm+1];
+    }
+  }
+
+NOINLINE static void fill_first_quadrant(size_t n, double * restrict res)
+  {
+  const double hsqt2 = 0.707106781186547524400844362104849;
+  size_t quart = n>>2;
+  if ((n&7)==0)
+    res[quart] = res[quart+1] = hsqt2;
+  for (size_t i=2, j=2*quart-2; i<quart; i+=2, j-=2)
+    {
+    res[j  ] = res[i+1];
+    res[j+1] = res[i  ];
+    }
+  }
+
+NOINLINE static void fill_first_half(size_t n, double * restrict res)
+  {
+  size_t half = n>>1;
+  if ((n&3)==0)
+    for (size_t i=0; i<half; i+=2)
+      {
+      res[i+half]   = -res[i+1];
+      res[i+half+1] =  res[i  ];
+      }
+  else
+    for (size_t i=2, j=2*half-2; i<half; i+=2, j-=2)
+      {
+      res[j  ] = -res[i  ];
+      res[j+1] =  res[i+1];
+      }
+  }
+
+NOINLINE static void fill_second_half(size_t n, double * restrict res)
+  {
+  if ((n&1)==0)
+    for (size_t i=0; i<n; ++i)
+      res[i+n] = -res[i];
+  else
+    for (size_t i=2, j=2*n-2; i<n; i+=2, j-=2)
+      {
+      res[j  ] =  res[i  ];
+      res[j+1] = -res[i+1];
+      }
+  }
+
+NOINLINE static void sincos_2pibyn_half(size_t n, double * restrict res)
+  {
+  if ((n&3)==0)
+    {
+    calc_first_octant(n, res);
+    fill_first_quadrant(n, res);
+    fill_first_half(n, res);
+    }
+  else if ((n&1)==0)
+    {
+    calc_first_quadrant(n, res);
+    fill_first_half(n, res);
+    }
+  else
+    calc_first_half(n, res);
+  }
+
+NOINLINE static void sincos_2pibyn(size_t n, double * restrict res)
+  {
+  sincos_2pibyn_half(n, res);
+  fill_second_half(n, res);
+  }
+
+NOINLINE static size_t largest_prime_factor (size_t n)
+  {
+  size_t res=1;
+  size_t tmp;
+  while (((tmp=(n>>1))<<1)==n)
+    { res=2; n=tmp; }
+
+  size_t limit=(size_t)sqrt(n+0.01);
+  for (size_t x=3; x<=limit; x+=2)
+  while (((tmp=(n/x))*x)==n)
+    {
+    res=x;
+    n=tmp;
+    limit=(size_t)sqrt(n+0.01);
+    }
+  if (n>1) res=n;
+
+  return res;
+  }
+
+NOINLINE static double cost_guess (size_t n)
+  {
+  const double lfp=1.1; // penalty for non-hardcoded larger factors
+  size_t ni=n;
+  double result=0.;
+  size_t tmp;
+  while (((tmp=(n>>1))<<1)==n)
+    { result+=2; n=tmp; }
+
+  size_t limit=(size_t)sqrt(n+0.01);
+  for (size_t x=3; x<=limit; x+=2)
+  while ((tmp=(n/x))*x==n)
+    {
+    result+= (x<=5) ? x : lfp*x; // penalize larger prime factors
+    n=tmp;
+    limit=(size_t)sqrt(n+0.01);
+    }
+  if (n>1) result+=(n<=5) ? n : lfp*n;
+
+  return result*ni;
+  }
+
+/* returns the smallest composite of 2, 3, 5, 7 and 11 which is >= n */
+NOINLINE static size_t good_size(size_t n)
+  {
+  if (n<=6) return n;
+
+  size_t bestfac=2*n;
+  for (size_t f2=1; f2<bestfac; f2*=2)
+    for (size_t f23=f2; f23<bestfac; f23*=3)
+      for (size_t f235=f23; f235<bestfac; f235*=5)
+        for (size_t f2357=f235; f2357<bestfac; f2357*=7)
+          for (size_t f235711=f2357; f235711<bestfac; f235711*=11)
+            if (f235711>=n) bestfac=f235711;
+  return bestfac;
+  }
+
+typedef struct cmplx {
+  double r,i;
+} cmplx;
+
+#define NFCT 25
+typedef struct cfftp_fctdata
+  {
+  size_t fct;
+  cmplx *tw, *tws;
+  } cfftp_fctdata;
+
+typedef struct cfftp_plan_i
+  {
+  size_t length, nfct;
+  cmplx *mem;
+  cfftp_fctdata fct[NFCT];
+  } cfftp_plan_i;
+typedef struct cfftp_plan_i * cfftp_plan;
+
+#define PMC(a,b,c,d) { a.r=c.r+d.r; a.i=c.i+d.i; b.r=c.r-d.r; b.i=c.i-d.i; }
+#define ADDC(a,b,c) { a.r=b.r+c.r; a.i=b.i+c.i; }
+#define SCALEC(a,b) { a.r*=b; a.i*=b; }
+#define ROT90(a) { double tmp_=a.r; a.r=-a.i; a.i=tmp_; }
+#define ROTM90(a) { double tmp_=-a.r; a.r=a.i; a.i=tmp_; }
+#define CH(a,b,c) ch[(a)+ido*((b)+l1*(c))]
+#define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))]
+#define WA(x,i) wa[(i)-1+(x)*(ido-1)]
+/* a = b*c */
+#define A_EQ_B_MUL_C(a,b,c) { a.r=b.r*c.r-b.i*c.i; a.i=b.r*c.i+b.i*c.r; }
+/* a = conj(b)*c*/
+#define A_EQ_CB_MUL_C(a,b,c) { a.r=b.r*c.r+b.i*c.i; a.i=b.r*c.i-b.i*c.r; }
+
+#define PMSIGNC(a,b,c,d) { a.r=c.r+sign*d.r; a.i=c.i+sign*d.i; b.r=c.r-sign*d.r; b.i=c.i-sign*d.i; }
+/* a = b*c */
+#define MULPMSIGNC(a,b,c) { a.r=b.r*c.r-sign*b.i*c.i; a.i=b.r*c.i+sign*b.i*c.r; }
+/* a *= b */
+#define MULPMSIGNCEQ(a,b) { double xtmp=a.r; a.r=b.r*a.r-sign*b.i*a.i; a.i=b.r*a.i+sign*b.i*xtmp; }
+
+NOINLINE static void pass2b (size_t ido, size_t l1, const cmplx * restrict cc,
+  cmplx * restrict ch, const cmplx * restrict wa)
+  {
+  const size_t cdim=2;
+
+  if (ido==1)
+    for (size_t k=0; k<l1; ++k)
+      PMC (CH(0,k,0),CH(0,k,1),CC(0,0,k),CC(0,1,k))
+  else
+    for (size_t k=0; k<l1; ++k)
+      {
+      PMC (CH(0,k,0),CH(0,k,1),CC(0,0,k),CC(0,1,k))
+      for (size_t i=1; i<ido; ++i)
+        {
+        cmplx t;
+        PMC (CH(i,k,0),t,CC(i,0,k),CC(i,1,k))
+        A_EQ_B_MUL_C (CH(i,k,1),WA(0,i),t)
+        }
+      }
+  }
+
+NOINLINE static void pass2f (size_t ido, size_t l1, const cmplx * restrict cc,
+  cmplx * restrict ch, const cmplx * restrict wa)
+  {
+  const size_t cdim=2;
+
+  if (ido==1)
+    for (size_t k=0; k<l1; ++k)
+      PMC (CH(0,k,0),CH(0,k,1),CC(0,0,k),CC(0,1,k))
+  else
+    for (size_t k=0; k<l1; ++k)
+      {
+      PMC (CH(0,k,0),CH(0,k,1),CC(0,0,k),CC(0,1,k))
+      for (size_t i=1; i<ido; ++i)
+        {
+        cmplx t;
+        PMC (CH(i,k,0),t,CC(i,0,k),CC(i,1,k))
+        A_EQ_CB_MUL_C (CH(i,k,1),WA(0,i),t)
+        }
+      }
+  }
+
+#define PREP3(idx) \
+        cmplx t0 = CC(idx,0,k), t1, t2; \
+        PMC (t1,t2,CC(idx,1,k),CC(idx,2,k)) \
+        CH(idx,k,0).r=t0.r+t1.r; \
+        CH(idx,k,0).i=t0.i+t1.i;
+#define PARTSTEP3a(u1,u2,twr,twi) \
+        { \
+        cmplx ca,cb; \
+        ca.r=t0.r+twr*t1.r; \
+        ca.i=t0.i+twr*t1.i; \
+        cb.i=twi*t2.r; \
+        cb.r=-(twi*t2.i); \
+        PMC(CH(0,k,u1),CH(0,k,u2),ca,cb) \
+        }
+
+#define PARTSTEP3b(u1,u2,twr,twi) \
+        { \
+        cmplx ca,cb,da,db; \
+        ca.r=t0.r+twr*t1.r; \
+        ca.i=t0.i+twr*t1.i; \
+        cb.i=twi*t2.r; \
+        cb.r=-(twi*t2.i); \
+        PMC(da,db,ca,cb) \
+        A_EQ_B_MUL_C (CH(i,k,u1),WA(u1-1,i),da) \
+        A_EQ_B_MUL_C (CH(i,k,u2),WA(u2-1,i),db) \
+        }
+NOINLINE static void pass3b (size_t ido, size_t l1, const cmplx * restrict cc,
+  cmplx * restrict ch, const cmplx * restrict wa)
+  {
+  const size_t cdim=3;
+  const double tw1r=-0.5, tw1i= 0.86602540378443864676;
+
+  if (ido==1)
+    for (size_t k=0; k<l1; ++k)
+      {
+      PREP3(0)
+      PARTSTEP3a(1,2,tw1r,tw1i)
+      }
+  else
+    for (size_t k=0; k<l1; ++k)
+      {
+      {
+      PREP3(0)
+      PARTSTEP3a(1,2,tw1r,tw1i)
+      }
+      for (size_t i=1; i<ido; ++i)
+        {
+        PREP3(i)
+        PARTSTEP3b(1,2,tw1r,tw1i)
+        }
+      }
+  }
+#define PARTSTEP3f(u1,u2,twr,twi) \
+        { \
+        cmplx ca,cb,da,db; \
+        ca.r=t0.r+twr*t1.r; \
+        ca.i=t0.i+twr*t1.i; \
+        cb.i=twi*t2.r; \
+        cb.r=-(twi*t2.i); \
+        PMC(da,db,ca,cb) \
+        A_EQ_CB_MUL_C (CH(i,k,u1),WA(u1-1,i),da) \
+        A_EQ_CB_MUL_C (CH(i,k,u2),WA(u2-1,i),db) \
+        }
+NOINLINE static void pass3f (size_t ido, size_t l1, const cmplx * restrict cc,
+  cmplx * restrict ch, const cmplx * restrict wa)
+  {
+  const size_t cdim=3;
+  const double tw1r=-0.5, tw1i= -0.86602540378443864676;
+
+  if (ido==1)
+    for (size_t k=0; k<l1; ++k)
+      {
+      PREP3(0)
+      PARTSTEP3a(1,2,tw1r,tw1i)
+      }
+  else
+    for (size_t k=0; k<l1; ++k)
+      {
+      {
+      PREP3(0)
+      PARTSTEP3a(1,2,tw1r,tw1i)
+      }
+      for (size_t i=1; i<ido; ++i)
+        {
+        PREP3(i)
+        PARTSTEP3f(1,2,tw1r,tw1i)
+        }
+      }
+  }
+
+NOINLINE static void pass4b (size_t ido, size_t l1, const cmplx * restrict cc,
+  cmplx * restrict ch, const cmplx * restrict wa)
+  {
+  const size_t cdim=4;
+
+  if (ido==1)
+    for (size_t k=0; k<l1; ++k)
+      {
+      cmplx t1, t2, t3, t4;
+      PMC(t2,t1,CC(0,0,k),CC(0,2,k))
+      PMC(t3,t4,CC(0,1,k),CC(0,3,k))
+      ROT90(t4)
+      PMC(CH(0,k,0),CH(0,k,2),t2,t3)
+      PMC(CH(0,k,1),CH(0,k,3),t1,t4)
+      }
+  else
+    for (size_t k=0; k<l1; ++k)
+      {
+      {
+      cmplx t1, t2, t3, t4;
+      PMC(t2,t1,CC(0,0,k),CC(0,2,k))
+      PMC(t3,t4,CC(0,1,k),CC(0,3,k))
+      ROT90(t4)
+      PMC(CH(0,k,0),CH(0,k,2),t2,t3)
+      PMC(CH(0,k,1),CH(0,k,3),t1,t4)
+      }
+      for (size_t i=1; i<ido; ++i)
+        {
+        cmplx c2, c3, c4, t1, t2, t3, t4;
+        cmplx cc0=CC(i,0,k), cc1=CC(i,1,k),cc2=CC(i,2,k),cc3=CC(i,3,k);
+        PMC(t2,t1,cc0,cc2)
+        PMC(t3,t4,cc1,cc3)
+        ROT90(t4)
+        cmplx wa0=WA(0,i), wa1=WA(1,i),wa2=WA(2,i);
+        PMC(CH(i,k,0),c3,t2,t3)
+        PMC(c2,c4,t1,t4)
+        A_EQ_B_MUL_C (CH(i,k,1),wa0,c2)
+        A_EQ_B_MUL_C (CH(i,k,2),wa1,c3)
+        A_EQ_B_MUL_C (CH(i,k,3),wa2,c4)
+        }
+      }
+  }
+NOINLINE static void pass4f (size_t ido, size_t l1, const cmplx * restrict cc,
+  cmplx * restrict ch, const cmplx * restrict wa)
+  {
+  const size_t cdim=4;
+
+  if (ido==1)
+    for (size_t k=0; k<l1; ++k)
+      {
+      cmplx t1, t2, t3, t4;
+      PMC(t2,t1,CC(0,0,k),CC(0,2,k))
+      PMC(t3,t4,CC(0,1,k),CC(0,3,k))
+      ROTM90(t4)
+      PMC(CH(0,k,0),CH(0,k,2),t2,t3)
+      PMC(CH(0,k,1),CH(0,k,3),t1,t4)
+      }
+  else
+    for (size_t k=0; k<l1; ++k)
+      {
+      {
+      cmplx t1, t2, t3, t4;
+      PMC(t2,t1,CC(0,0,k),CC(0,2,k))
+      PMC(t3,t4,CC(0,1,k),CC(0,3,k))
+      ROTM90(t4)
+      PMC(CH(0,k,0),CH(0,k,2),t2,t3)
+      PMC (CH(0,k,1),CH(0,k,3),t1,t4)
+      }
+      for (size_t i=1; i<ido; ++i)
+        {
+        cmplx c2, c3, c4, t1, t2, t3, t4;
+        cmplx cc0=CC(i,0,k), cc1=CC(i,1,k),cc2=CC(i,2,k),cc3=CC(i,3,k);
+        PMC(t2,t1,cc0,cc2)
+        PMC(t3,t4,cc1,cc3)
+        ROTM90(t4)
+        cmplx wa0=WA(0,i), wa1=WA(1,i),wa2=WA(2,i);
+        PMC(CH(i,k,0),c3,t2,t3)
+        PMC(c2,c4,t1,t4)
+        A_EQ_CB_MUL_C (CH(i,k,1),wa0,c2)
+        A_EQ_CB_MUL_C (CH(i,k,2),wa1,c3)
+        A_EQ_CB_MUL_C (CH(i,k,3),wa2,c4)
+        }
+      }
+  }
+
+#define PREP5(idx) \
+        cmplx t0 = CC(idx,0,k), t1, t2, t3, t4; \
+        PMC (t1,t4,CC(idx,1,k),CC(idx,4,k)) \
+        PMC (t2,t3,CC(idx,2,k),CC(idx,3,k)) \
+        CH(idx,k,0).r=t0.r+t1.r+t2.r; \
+        CH(idx,k,0).i=t0.i+t1.i+t2.i;
+
+#define PARTSTEP5a(u1,u2,twar,twbr,twai,twbi) \
+        { \
+        cmplx ca,cb; \
+        ca.r=t0.r+twar*t1.r+twbr*t2.r; \
+        ca.i=t0.i+twar*t1.i+twbr*t2.i; \
+        cb.i=twai*t4.r twbi*t3.r; \
+        cb.r=-(twai*t4.i twbi*t3.i); \
+        PMC(CH(0,k,u1),CH(0,k,u2),ca,cb) \
+        }
+
+#define PARTSTEP5b(u1,u2,twar,twbr,twai,twbi) \
+        { \
+        cmplx ca,cb,da,db; \
+        ca.r=t0.r+twar*t1.r+twbr*t2.r; \
+        ca.i=t0.i+twar*t1.i+twbr*t2.i; \
+        cb.i=twai*t4.r twbi*t3.r; \
+        cb.r=-(twai*t4.i twbi*t3.i); \
+        PMC(da,db,ca,cb) \
+        A_EQ_B_MUL_C (CH(i,k,u1),WA(u1-1,i),da) \
+        A_EQ_B_MUL_C (CH(i,k,u2),WA(u2-1,i),db) \
+        }
+NOINLINE static void pass5b (size_t ido, size_t l1, const cmplx * restrict cc,
+  cmplx * restrict ch, const cmplx * restrict wa)
+  {
+  const size_t cdim=5;
+  const double tw1r= 0.3090169943749474241,
+               tw1i= 0.95105651629515357212,
+               tw2r= -0.8090169943749474241,
+               tw2i= 0.58778525229247312917;
+
+  if (ido==1)
+    for (size_t k=0; k<l1; ++k)
+      {
+      PREP5(0)
+      PARTSTEP5a(1,4,tw1r,tw2r,+tw1i,+tw2i)
+      PARTSTEP5a(2,3,tw2r,tw1r,+tw2i,-tw1i)
+      }
+  else
+    for (size_t k=0; k<l1; ++k)
+      {
+      {
+      PREP5(0)
+      PARTSTEP5a(1,4,tw1r,tw2r,+tw1i,+tw2i)
+      PARTSTEP5a(2,3,tw2r,tw1r,+tw2i,-tw1i)
+      }
+      for (size_t i=1; i<ido; ++i)
+        {
+        PREP5(i)
+        PARTSTEP5b(1,4,tw1r,tw2r,+tw1i,+tw2i)
+        PARTSTEP5b(2,3,tw2r,tw1r,+tw2i,-tw1i)
+        }
+      }
+  }
+#define PARTSTEP5f(u1,u2,twar,twbr,twai,twbi) \
+        { \
+        cmplx ca,cb,da,db; \
+        ca.r=t0.r+twar*t1.r+twbr*t2.r; \
+        ca.i=t0.i+twar*t1.i+twbr*t2.i; \
+        cb.i=twai*t4.r twbi*t3.r; \
+        cb.r=-(twai*t4.i twbi*t3.i); \
+        PMC(da,db,ca,cb) \
+        A_EQ_CB_MUL_C (CH(i,k,u1),WA(u1-1,i),da) \
+        A_EQ_CB_MUL_C (CH(i,k,u2),WA(u2-1,i),db) \
+        }
+NOINLINE static void pass5f (size_t ido, size_t l1, const cmplx * restrict cc,
+  cmplx * restrict ch, const cmplx * restrict wa)
+  {
+  const size_t cdim=5;
+  const double tw1r= 0.3090169943749474241,
+               tw1i= -0.95105651629515357212,
+               tw2r= -0.8090169943749474241,
+               tw2i= -0.58778525229247312917;
+
+  if (ido==1)
+    for (size_t k=0; k<l1; ++k)
+      {
+      PREP5(0)
+      PARTSTEP5a(1,4,tw1r,tw2r,+tw1i,+tw2i)
+      PARTSTEP5a(2,3,tw2r,tw1r,+tw2i,-tw1i)
+      }
+  else
+    for (size_t k=0; k<l1; ++k)
+      {
+      {
+      PREP5(0)
+      PARTSTEP5a(1,4,tw1r,tw2r,+tw1i,+tw2i)
+      PARTSTEP5a(2,3,tw2r,tw1r,+tw2i,-tw1i)
+      }
+      for (size_t i=1; i<ido; ++i)
+        {
+        PREP5(i)
+        PARTSTEP5f(1,4,tw1r,tw2r,+tw1i,+tw2i)
+        PARTSTEP5f(2,3,tw2r,tw1r,+tw2i,-tw1i)
+        }
+      }
+  }
+
+#define PREP7(idx) \
+        cmplx t1 = CC(idx,0,k), t2, t3, t4, t5, t6, t7; \
+        PMC (t2,t7,CC(idx,1,k),CC(idx,6,k)) \
+        PMC (t3,t6,CC(idx,2,k),CC(idx,5,k)) \
+        PMC (t4,t5,CC(idx,3,k),CC(idx,4,k)) \
+        CH(idx,k,0).r=t1.r+t2.r+t3.r+t4.r; \
+        CH(idx,k,0).i=t1.i+t2.i+t3.i+t4.i;
+
+#define PARTSTEP7a0(u1,u2,x1,x2,x3,y1,y2,y3,out1,out2) \
+        { \
+        cmplx ca,cb; \
+        ca.r=t1.r+x1*t2.r+x2*t3.r+x3*t4.r; \
+        ca.i=t1.i+x1*t2.i+x2*t3.i+x3*t4.i; \
+        cb.i=y1*t7.r y2*t6.r y3*t5.r; \
+        cb.r=-(y1*t7.i y2*t6.i y3*t5.i); \
+        PMC(out1,out2,ca,cb) \
+        }
+#define PARTSTEP7a(u1,u2,x1,x2,x3,y1,y2,y3) \
+        PARTSTEP7a0(u1,u2,x1,x2,x3,y1,y2,y3,CH(0,k,u1),CH(0,k,u2))
+#define PARTSTEP7(u1,u2,x1,x2,x3,y1,y2,y3) \
+        { \
+        cmplx da,db; \
+        PARTSTEP7a0(u1,u2,x1,x2,x3,y1,y2,y3,da,db) \
+        MULPMSIGNC (CH(i,k,u1),WA(u1-1,i),da) \
+        MULPMSIGNC (CH(i,k,u2),WA(u2-1,i),db) \
+        }
+
+NOINLINE static void pass7(size_t ido, size_t l1, const cmplx * restrict cc,
+  cmplx * restrict ch, const cmplx * restrict wa, const int sign)
+  {
+  const size_t cdim=7;
+  const double tw1r= 0.623489801858733530525,
+               tw1i= sign * 0.7818314824680298087084,
+               tw2r= -0.222520933956314404289,
+               tw2i= sign * 0.9749279121818236070181,
+               tw3r= -0.9009688679024191262361,
+               tw3i= sign * 0.4338837391175581204758;
+
+  if (ido==1)
+    for (size_t k=0; k<l1; ++k)
+      {
+      PREP7(0)
+      PARTSTEP7a(1,6,tw1r,tw2r,tw3r,+tw1i,+tw2i,+tw3i)
+      PARTSTEP7a(2,5,tw2r,tw3r,tw1r,+tw2i,-tw3i,-tw1i)
+      PARTSTEP7a(3,4,tw3r,tw1r,tw2r,+tw3i,-tw1i,+tw2i)
+      }
+  else
+    for (size_t k=0; k<l1; ++k)
+      {
+      {
+      PREP7(0)
+      PARTSTEP7a(1,6,tw1r,tw2r,tw3r,+tw1i,+tw2i,+tw3i)
+      PARTSTEP7a(2,5,tw2r,tw3r,tw1r,+tw2i,-tw3i,-tw1i)
+      PARTSTEP7a(3,4,tw3r,tw1r,tw2r,+tw3i,-tw1i,+tw2i)
+      }
+      for (size_t i=1; i<ido; ++i)
+        {
+        PREP7(i)
+        PARTSTEP7(1,6,tw1r,tw2r,tw3r,+tw1i,+tw2i,+tw3i)
+        PARTSTEP7(2,5,tw2r,tw3r,tw1r,+tw2i,-tw3i,-tw1i)
+        PARTSTEP7(3,4,tw3r,tw1r,tw2r,+tw3i,-tw1i,+tw2i)
+        }
+      }
+  }
+
+#define PREP11(idx) \
+        cmplx t1 = CC(idx,0,k), t2, t3, t4, t5, t6, t7, t8, t9, t10, t11; \
+        PMC (t2,t11,CC(idx,1,k),CC(idx,10,k)) \
+        PMC (t3,t10,CC(idx,2,k),CC(idx, 9,k)) \
+        PMC (t4,t9 ,CC(idx,3,k),CC(idx, 8,k)) \
+        PMC (t5,t8 ,CC(idx,4,k),CC(idx, 7,k)) \
+        PMC (t6,t7 ,CC(idx,5,k),CC(idx, 6,k)) \
+        CH(idx,k,0).r=t1.r+t2.r+t3.r+t4.r+t5.r+t6.r; \
+        CH(idx,k,0).i=t1.i+t2.i+t3.i+t4.i+t5.i+t6.i;
+
+#define PARTSTEP11a0(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,out1,out2) \
+        { \
+        cmplx ca,cb; \
+        ca.r=t1.r+x1*t2.r+x2*t3.r+x3*t4.r+x4*t5.r+x5*t6.r; \
+        ca.i=t1.i+x1*t2.i+x2*t3.i+x3*t4.i+x4*t5.i+x5*t6.i; \
+        cb.i=y1*t11.r y2*t10.r y3*t9.r y4*t8.r y5*t7.r; \
+        cb.r=-(y1*t11.i y2*t10.i y3*t9.i y4*t8.i y5*t7.i ); \
+        PMC(out1,out2,ca,cb) \
+        }
+#define PARTSTEP11a(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5) \
+        PARTSTEP11a0(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,CH(0,k,u1),CH(0,k,u2))
+#define PARTSTEP11(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5) \
+        { \
+        cmplx da,db; \
+        PARTSTEP11a0(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,da,db) \
+        MULPMSIGNC (CH(i,k,u1),WA(u1-1,i),da) \
+        MULPMSIGNC (CH(i,k,u2),WA(u2-1,i),db) \
+        }
+
+NOINLINE static void pass11 (size_t ido, size_t l1, const cmplx * restrict cc,
+  cmplx * restrict ch, const cmplx * restrict wa, const int sign)
+  {
+  const size_t cdim=11;
+  const double tw1r =        0.8412535328311811688618,
+               tw1i = sign * 0.5406408174555975821076,
+               tw2r =        0.4154150130018864255293,
+               tw2i = sign * 0.9096319953545183714117,
+               tw3r =       -0.1423148382732851404438,
+               tw3i = sign * 0.9898214418809327323761,
+               tw4r =       -0.6548607339452850640569,
+               tw4i = sign * 0.755749574354258283774,
+               tw5r =       -0.9594929736144973898904,
+               tw5i = sign * 0.2817325568414296977114;
+
+  if (ido==1)
+    for (size_t k=0; k<l1; ++k)
+      {
+      PREP11(0)
+      PARTSTEP11a(1,10,tw1r,tw2r,tw3r,tw4r,tw5r,+tw1i,+tw2i,+tw3i,+tw4i,+tw5i)
+      PARTSTEP11a(2, 9,tw2r,tw4r,tw5r,tw3r,tw1r,+tw2i,+tw4i,-tw5i,-tw3i,-tw1i)
+      PARTSTEP11a(3, 8,tw3r,tw5r,tw2r,tw1r,tw4r,+tw3i,-tw5i,-tw2i,+tw1i,+tw4i)
+      PARTSTEP11a(4, 7,tw4r,tw3r,tw1r,tw5r,tw2r,+tw4i,-tw3i,+tw1i,+tw5i,-tw2i)
+      PARTSTEP11a(5, 6,tw5r,tw1r,tw4r,tw2r,tw3r,+tw5i,-tw1i,+tw4i,-tw2i,+tw3i)
+      }
+  else
+    for (size_t k=0; k<l1; ++k)
+      {
+      {
+      PREP11(0)
+      PARTSTEP11a(1,10,tw1r,tw2r,tw3r,tw4r,tw5r,+tw1i,+tw2i,+tw3i,+tw4i,+tw5i)
+      PARTSTEP11a(2, 9,tw2r,tw4r,tw5r,tw3r,tw1r,+tw2i,+tw4i,-tw5i,-tw3i,-tw1i)
+      PARTSTEP11a(3, 8,tw3r,tw5r,tw2r,tw1r,tw4r,+tw3i,-tw5i,-tw2i,+tw1i,+tw4i)
+      PARTSTEP11a(4, 7,tw4r,tw3r,tw1r,tw5r,tw2r,+tw4i,-tw3i,+tw1i,+tw5i,-tw2i)
+      PARTSTEP11a(5, 6,tw5r,tw1r,tw4r,tw2r,tw3r,+tw5i,-tw1i,+tw4i,-tw2i,+tw3i)
+      }
+      for (size_t i=1; i<ido; ++i)
+        {
+        PREP11(i)
+        PARTSTEP11(1,10,tw1r,tw2r,tw3r,tw4r,tw5r,+tw1i,+tw2i,+tw3i,+tw4i,+tw5i)
+        PARTSTEP11(2, 9,tw2r,tw4r,tw5r,tw3r,tw1r,+tw2i,+tw4i,-tw5i,-tw3i,-tw1i)
+        PARTSTEP11(3, 8,tw3r,tw5r,tw2r,tw1r,tw4r,+tw3i,-tw5i,-tw2i,+tw1i,+tw4i)
+        PARTSTEP11(4, 7,tw4r,tw3r,tw1r,tw5r,tw2r,+tw4i,-tw3i,+tw1i,+tw5i,-tw2i)
+        PARTSTEP11(5, 6,tw5r,tw1r,tw4r,tw2r,tw3r,+tw5i,-tw1i,+tw4i,-tw2i,+tw3i)
+        }
+      }
+  }
+
+#define CX(a,b,c) cc[(a)+ido*((b)+l1*(c))]
+#define CX2(a,b) cc[(a)+idl1*(b)]
+#define CH2(a,b) ch[(a)+idl1*(b)]
+
+NOINLINE static int passg (size_t ido, size_t ip, size_t l1,
+  cmplx * restrict cc, cmplx * restrict ch, const cmplx * restrict wa,
+  const cmplx * restrict csarr, const int sign)
+  {
+  const size_t cdim=ip;
+  size_t ipph = (ip+1)/2;
+  size_t idl1 = ido*l1;
+
+  cmplx * restrict wal=RALLOC(cmplx,ip);
+  if (!wal) return -1;
+  wal[0]=(cmplx){1.,0.};
+  for (size_t i=1; i<ip; ++i)
+    wal[i]=(cmplx){csarr[i].r,sign*csarr[i].i};
+
+  for (size_t k=0; k<l1; ++k)
+    for (size_t i=0; i<ido; ++i)
+      CH(i,k,0) = CC(i,0,k);
+  for (size_t j=1, jc=ip-1; j<ipph; ++j, --jc)
+    for (size_t k=0; k<l1; ++k)
+      for (size_t i=0; i<ido; ++i)
+        PMC(CH(i,k,j),CH(i,k,jc),CC(i,j,k),CC(i,jc,k))
+  for (size_t k=0; k<l1; ++k)
+    for (size_t i=0; i<ido; ++i)
+      {
+      cmplx tmp = CH(i,k,0);
+      for (size_t j=1; j<ipph; ++j)
+        ADDC(tmp,tmp,CH(i,k,j))
+      CX(i,k,0) = tmp;
+      }
+  for (size_t l=1, lc=ip-1; l<ipph; ++l, --lc)
+    {
+    // j=0
+    for (size_t ik=0; ik<idl1; ++ik)
+      {
+      CX2(ik,l).r = CH2(ik,0).r+wal[l].r*CH2(ik,1).r+wal[2*l].r*CH2(ik,2).r;
+      CX2(ik,l).i = CH2(ik,0).i+wal[l].r*CH2(ik,1).i+wal[2*l].r*CH2(ik,2).i;
+      CX2(ik,lc).r=-wal[l].i*CH2(ik,ip-1).i-wal[2*l].i*CH2(ik,ip-2).i;
+      CX2(ik,lc).i=wal[l].i*CH2(ik,ip-1).r+wal[2*l].i*CH2(ik,ip-2).r;
+      }
+
+    size_t iwal=2*l;
+    size_t j=3, jc=ip-3;
+    for (; j<ipph-1; j+=2, jc-=2)
+      {
+      iwal+=l; if (iwal>ip) iwal-=ip;
+      cmplx xwal=wal[iwal];
+      iwal+=l; if (iwal>ip) iwal-=ip;
+      cmplx xwal2=wal[iwal];
+      for (size_t ik=0; ik<idl1; ++ik)
+        {
+        CX2(ik,l).r += CH2(ik,j).r*xwal.r+CH2(ik,j+1).r*xwal2.r;
+        CX2(ik,l).i += CH2(ik,j).i*xwal.r+CH2(ik,j+1).i*xwal2.r;
+        CX2(ik,lc).r -= CH2(ik,jc).i*xwal.i+CH2(ik,jc-1).i*xwal2.i;
+        CX2(ik,lc).i += CH2(ik,jc).r*xwal.i+CH2(ik,jc-1).r*xwal2.i;
+        }
+      }
+    for (; j<ipph; ++j, --jc)
+      {
+      iwal+=l; if (iwal>ip) iwal-=ip;
+      cmplx xwal=wal[iwal];
+      for (size_t ik=0; ik<idl1; ++ik)
+        {
+        CX2(ik,l).r += CH2(ik,j).r*xwal.r;
+        CX2(ik,l).i += CH2(ik,j).i*xwal.r;
+        CX2(ik,lc).r -= CH2(ik,jc).i*xwal.i;
+        CX2(ik,lc).i += CH2(ik,jc).r*xwal.i;
+        }
+      }
+    }
+  DEALLOC(wal);
+
+  // shuffling and twiddling
+  if (ido==1)
+    for (size_t j=1, jc=ip-1; j<ipph; ++j, --jc)
+      for (size_t ik=0; ik<idl1; ++ik)
+        {
+        cmplx t1=CX2(ik,j), t2=CX2(ik,jc);
+        PMC(CX2(ik,j),CX2(ik,jc),t1,t2)
+        }
+  else
+    {
+    for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc)
+      for (size_t k=0; k<l1; ++k)
+        {
+        cmplx t1=CX(0,k,j), t2=CX(0,k,jc);
+        PMC(CX(0,k,j),CX(0,k,jc),t1,t2)
+        for (size_t i=1; i<ido; ++i)
+          {
+          cmplx x1, x2;
+          PMC(x1,x2,CX(i,k,j),CX(i,k,jc))
+          size_t idij=(j-1)*(ido-1)+i-1;
+          MULPMSIGNC (CX(i,k,j),wa[idij],x1)
+          idij=(jc-1)*(ido-1)+i-1;
+          MULPMSIGNC (CX(i,k,jc),wa[idij],x2)
+          }
+        }
+    }
+  return 0;
+  }
+
+#undef CH2
+#undef CX2
+#undef CX
+
+NOINLINE WARN_UNUSED_RESULT static int pass_all(cfftp_plan plan, cmplx c[], double fct,
+  const int sign)
+  {
+  if (plan->length==1) return 0;
+  size_t len=plan->length;
+  size_t l1=1, nf=plan->nfct;
+  cmplx *ch = RALLOC(cmplx, len);
+  if (!ch) return -1;
+  cmplx *p1=c, *p2=ch;
+
+  for(size_t k1=0; k1<nf; k1++)
+    {
+    size_t ip=plan->fct[k1].fct;
+    size_t l2=ip*l1;
+    size_t ido = len/l2;
+    if     (ip==4)
+      sign>0 ? pass4b (ido, l1, p1, p2, plan->fct[k1].tw)
+             : pass4f (ido, l1, p1, p2, plan->fct[k1].tw);
+    else if(ip==2)
+      sign>0 ? pass2b (ido, l1, p1, p2, plan->fct[k1].tw)
+             : pass2f (ido, l1, p1, p2, plan->fct[k1].tw);
+    else if(ip==3)
+      sign>0 ? pass3b (ido, l1, p1, p2, plan->fct[k1].tw)
+             : pass3f (ido, l1, p1, p2, plan->fct[k1].tw);
+    else if(ip==5)
+      sign>0 ? pass5b (ido, l1, p1, p2, plan->fct[k1].tw)
+             : pass5f (ido, l1, p1, p2, plan->fct[k1].tw);
+    else if(ip==7)  pass7 (ido, l1, p1, p2, plan->fct[k1].tw, sign);
+    else if(ip==11) pass11(ido, l1, p1, p2, plan->fct[k1].tw, sign);
+    else
+      {
+      if (passg(ido, ip, l1, p1, p2, plan->fct[k1].tw, plan->fct[k1].tws, sign))
+        { DEALLOC(ch); return -1; }
+      SWAP(p1,p2,cmplx *);
+      }
+    SWAP(p1,p2,cmplx *);
+    l1=l2;
+    }
+  if (p1!=c)
+    {
+    if (fct!=1.)
+      for (size_t i=0; i<len; ++i)
+        {
+        c[i].r = ch[i].r*fct;
+        c[i].i = ch[i].i*fct;
+        }
+    else
+      memcpy (c,p1,len*sizeof(cmplx));
+    }
+  else
+    if (fct!=1.)
+      for (size_t i=0; i<len; ++i)
+        {
+        c[i].r *= fct;
+        c[i].i *= fct;
+        }
+  DEALLOC(ch);
+  return 0;
+  }
+
+#undef PMSIGNC
+#undef A_EQ_B_MUL_C
+#undef A_EQ_CB_MUL_C
+#undef MULPMSIGNC
+#undef MULPMSIGNCEQ
+
+#undef WA
+#undef CC
+#undef CH
+#undef ROT90
+#undef SCALEC
+#undef ADDC
+#undef PMC
+
+NOINLINE WARN_UNUSED_RESULT
+static int cfftp_forward(cfftp_plan plan, double c[], double fct)
+  { return pass_all(plan,(cmplx *)c, fct, -1); }
+
+NOINLINE WARN_UNUSED_RESULT
+static int cfftp_backward(cfftp_plan plan, double c[], double fct)
+  { return pass_all(plan,(cmplx *)c, fct, 1); }
+
+NOINLINE WARN_UNUSED_RESULT
+static int cfftp_factorize (cfftp_plan plan)
+  {
+  size_t length=plan->length;
+  size_t nfct=0;
+  while ((length%4)==0)
+    { if (nfct>=NFCT) return -1; plan->fct[nfct++].fct=4; length>>=2; }
+  if ((length%2)==0)
+    {
+    length>>=1;
+    // factor 2 should be at the front of the factor list
+    if (nfct>=NFCT) return -1;
+    plan->fct[nfct++].fct=2;
+    SWAP(plan->fct[0].fct, plan->fct[nfct-1].fct,size_t);
+    }
+  size_t maxl=(size_t)(sqrt((double)length))+1;
+  for (size_t divisor=3; (length>1)&&(divisor<maxl); divisor+=2)
+    if ((length%divisor)==0)
+      {
+      while ((length%divisor)==0)
+        {
+        if (nfct>=NFCT) return -1;
+        plan->fct[nfct++].fct=divisor;
+        length/=divisor;
+        }
+      maxl=(size_t)(sqrt((double)length))+1;
+      }
+  if (length>1) plan->fct[nfct++].fct=length;
+  plan->nfct=nfct;
+  return 0;
+  }
+
+NOINLINE static size_t cfftp_twsize (cfftp_plan plan)
+  {
+  size_t twsize=0, l1=1;
+  for (size_t k=0; k<plan->nfct; ++k)
+    {
+    size_t ip=plan->fct[k].fct, ido= plan->length/(l1*ip);
+    twsize+=(ip-1)*(ido-1);
+    if (ip>11)
+      twsize+=ip;
+    l1*=ip;
+    }
+  return twsize;
+  }
+
+NOINLINE WARN_UNUSED_RESULT static int cfftp_comp_twiddle (cfftp_plan plan)
+  {
+  size_t length=plan->length;
+  double *twid = RALLOC(double, 2*length);
+  if (!twid) return -1;
+  sincos_2pibyn(length, twid);
+  size_t l1=1;
+  size_t memofs=0;
+  for (size_t k=0; k<plan->nfct; ++k)
+    {
+    size_t ip=plan->fct[k].fct, ido= length/(l1*ip);
+    plan->fct[k].tw=plan->mem+memofs;
+    memofs+=(ip-1)*(ido-1);
+    for (size_t j=1; j<ip; ++j)
+      for (size_t i=1; i<ido; ++i)
+        {
+        plan->fct[k].tw[(j-1)*(ido-1)+i-1].r = twid[2*j*l1*i];
+        plan->fct[k].tw[(j-1)*(ido-1)+i-1].i = twid[2*j*l1*i+1];
+        }
+    if (ip>11)
+      {
+      plan->fct[k].tws=plan->mem+memofs;
+      memofs+=ip;
+      for (size_t j=0; j<ip; ++j)
+        {
+        plan->fct[k].tws[j].r = twid[2*j*l1*ido];
+        plan->fct[k].tws[j].i = twid[2*j*l1*ido+1];
+        }
+      }
+    l1*=ip;
+    }
+  DEALLOC(twid);
+  return 0;
+  }
+
+static cfftp_plan make_cfftp_plan (size_t length)
+  {
+  if (length==0) return NULL;
+  cfftp_plan plan = RALLOC(cfftp_plan_i,1);
+  if (!plan) return NULL;
+  plan->length=length;
+  plan->nfct=0;
+  for (size_t i=0; i<NFCT; ++i)
+    plan->fct[i]=(cfftp_fctdata){0,0,0};
+  plan->mem=0;
+  if (length==1) return plan;
+  if (cfftp_factorize(plan)!=0) { DEALLOC(plan); return NULL; }
+  size_t tws=cfftp_twsize(plan);
+  plan->mem=RALLOC(cmplx,tws);
+  if (!plan->mem) { DEALLOC(plan); return NULL; }
+  if (cfftp_comp_twiddle(plan)!=0)
+    { DEALLOC(plan->mem); DEALLOC(plan); return NULL; }
+  return plan;
+  }
+
+static void destroy_cfftp_plan (cfftp_plan plan)
+  {
+  DEALLOC(plan->mem);
+  DEALLOC(plan);
+  }
+
+typedef struct rfftp_fctdata
+  {
+  size_t fct;
+  double *tw, *tws;
+  } rfftp_fctdata;
+
+typedef struct rfftp_plan_i
+  {
+  size_t length, nfct;
+  double *mem;
+  rfftp_fctdata fct[NFCT];
+  } rfftp_plan_i;
+typedef struct rfftp_plan_i * rfftp_plan;
+
+#define WA(x,i) wa[(i)+(x)*(ido-1)]
+#define PM(a,b,c,d) { a=c+d; b=c-d; }
+/* (a+ib) = conj(c+id) * (e+if) */
+#define MULPM(a,b,c,d,e,f) { a=c*e+d*f; b=c*f-d*e; }
+
+#define CC(a,b,c) cc[(a)+ido*((b)+l1*(c))]
+#define CH(a,b,c) ch[(a)+ido*((b)+cdim*(c))]
+
+NOINLINE static void radf2 (size_t ido, size_t l1, const double * restrict cc,
+  double * restrict ch, const double * restrict wa)
+  {
+  const size_t cdim=2;
+
+  for (size_t k=0; k<l1; k++)
+    PM (CH(0,0,k),CH(ido-1,1,k),CC(0,k,0),CC(0,k,1))
+  if ((ido&1)==0)
+    for (size_t k=0; k<l1; k++)
+      {
+      CH(    0,1,k) = -CC(ido-1,k,1);
+      CH(ido-1,0,k) =  CC(ido-1,k,0);
+      }
+  if (ido<=2) return;
+  for (size_t k=0; k<l1; k++)
+    for (size_t i=2; i<ido; i+=2)
+      {
+      size_t ic=ido-i;
+      double tr2, ti2;
+      MULPM (tr2,ti2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1))
+      PM (CH(i-1,0,k),CH(ic-1,1,k),CC(i-1,k,0),tr2)
+      PM (CH(i  ,0,k),CH(ic  ,1,k),ti2,CC(i  ,k,0))
+      }
+  }
+
+NOINLINE static void radf3(size_t ido, size_t l1, const double * restrict cc,
+  double * restrict ch, const double * restrict wa)
+  {
+  const size_t cdim=3;
+  static const double taur=-0.5, taui=0.86602540378443864676;
+
+  for (size_t k=0; k<l1; k++)
+    {
+    double cr2=CC(0,k,1)+CC(0,k,2);
+    CH(0,0,k) = CC(0,k,0)+cr2;
+    CH(0,2,k) = taui*(CC(0,k,2)-CC(0,k,1));
+    CH(ido-1,1,k) = CC(0,k,0)+taur*cr2;
+    }
+  if (ido==1) return;
+  for (size_t k=0; k<l1; k++)
+    for (size_t i=2; i<ido; i+=2)
+      {
+      size_t ic=ido-i;
+      double di2, di3, dr2, dr3;
+      MULPM (dr2,di2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1)) // d2=conj(WA0)*CC1
+      MULPM (dr3,di3,WA(1,i-2),WA(1,i-1),CC(i-1,k,2),CC(i,k,2)) // d3=conj(WA1)*CC2
+      double cr2=dr2+dr3; // c add
+      double ci2=di2+di3;
+      CH(i-1,0,k) = CC(i-1,k,0)+cr2; // c add
+      CH(i  ,0,k) = CC(i  ,k,0)+ci2;
+      double tr2 = CC(i-1,k,0)+taur*cr2; // c add
+      double ti2 = CC(i  ,k,0)+taur*ci2;
+      double tr3 = taui*(di2-di3);  // t3 = taui*i*(d3-d2)?
+      double ti3 = taui*(dr3-dr2);
+      PM(CH(i-1,2,k),CH(ic-1,1,k),tr2,tr3) // PM(i) = t2+t3
+      PM(CH(i  ,2,k),CH(ic  ,1,k),ti3,ti2) // PM(ic) = conj(t2-t3)
+      }
+  }
+
+NOINLINE static void radf4(size_t ido, size_t l1, const double * restrict cc,
+  double * restrict ch, const double * restrict wa)
+  {
+  const size_t cdim=4;
+  static const double hsqt2=0.70710678118654752440;
+
+  for (size_t k=0; k<l1; k++)
+    {
+    double tr1,tr2;
+    PM (tr1,CH(0,2,k),CC(0,k,3),CC(0,k,1))
+    PM (tr2,CH(ido-1,1,k),CC(0,k,0),CC(0,k,2))
+    PM (CH(0,0,k),CH(ido-1,3,k),tr2,tr1)
+    }
+  if ((ido&1)==0)
+    for (size_t k=0; k<l1; k++)
+      {
+      double ti1=-hsqt2*(CC(ido-1,k,1)+CC(ido-1,k,3));
+      double tr1= hsqt2*(CC(ido-1,k,1)-CC(ido-1,k,3));
+      PM (CH(ido-1,0,k),CH(ido-1,2,k),CC(ido-1,k,0),tr1)
+      PM (CH(    0,3,k),CH(    0,1,k),ti1,CC(ido-1,k,2))
+      }
+  if (ido<=2) return;
+  for (size_t k=0; k<l1; k++)
+    for (size_t i=2; i<ido; i+=2)
+      {
+      size_t ic=ido-i;
+      double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
+      MULPM(cr2,ci2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1))
+      MULPM(cr3,ci3,WA(1,i-2),WA(1,i-1),CC(i-1,k,2),CC(i,k,2))
+      MULPM(cr4,ci4,WA(2,i-2),WA(2,i-1),CC(i-1,k,3),CC(i,k,3))
+      PM(tr1,tr4,cr4,cr2)
+      PM(ti1,ti4,ci2,ci4)
+      PM(tr2,tr3,CC(i-1,k,0),cr3)
+      PM(ti2,ti3,CC(i  ,k,0),ci3)
+      PM(CH(i-1,0,k),CH(ic-1,3,k),tr2,tr1)
+      PM(CH(i  ,0,k),CH(ic  ,3,k),ti1,ti2)
+      PM(CH(i-1,2,k),CH(ic-1,1,k),tr3,ti4)
+      PM(CH(i  ,2,k),CH(ic  ,1,k),tr4,ti3)
+      }
+  }
+
+NOINLINE static void radf5(size_t ido, size_t l1, const double * restrict cc,
+  double * restrict ch, const double * restrict wa)
+  {
+  const size_t cdim=5;
+  static const double tr11= 0.3090169943749474241, ti11=0.95105651629515357212,
+                      tr12=-0.8090169943749474241, ti12=0.58778525229247312917;
+
+  for (size_t k=0; k<l1; k++)
+    {
+    double cr2, cr3, ci4, ci5;
+    PM (cr2,ci5,CC(0,k,4),CC(0,k,1))
+    PM (cr3,ci4,CC(0,k,3),CC(0,k,2))
+    CH(0,0,k)=CC(0,k,0)+cr2+cr3;
+    CH(ido-1,1,k)=CC(0,k,0)+tr11*cr2+tr12*cr3;
+    CH(0,2,k)=ti11*ci5+ti12*ci4;
+    CH(ido-1,3,k)=CC(0,k,0)+tr12*cr2+tr11*cr3;
+    CH(0,4,k)=ti12*ci5-ti11*ci4;
+    }
+  if (ido==1) return;
+  for (size_t k=0; k<l1;++k)
+    for (size_t i=2; i<ido; i+=2)
+      {
+      double ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3,
+         dr4, dr5, cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
+      size_t ic=ido-i;
+      MULPM (dr2,di2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1))
+      MULPM (dr3,di3,WA(1,i-2),WA(1,i-1),CC(i-1,k,2),CC(i,k,2))
+      MULPM (dr4,di4,WA(2,i-2),WA(2,i-1),CC(i-1,k,3),CC(i,k,3))
+      MULPM (dr5,di5,WA(3,i-2),WA(3,i-1),CC(i-1,k,4),CC(i,k,4))
+      PM(cr2,ci5,dr5,dr2)
+      PM(ci2,cr5,di2,di5)
+      PM(cr3,ci4,dr4,dr3)
+      PM(ci3,cr4,di3,di4)
+      CH(i-1,0,k)=CC(i-1,k,0)+cr2+cr3;
+      CH(i  ,0,k)=CC(i  ,k,0)+ci2+ci3;
+      tr2=CC(i-1,k,0)+tr11*cr2+tr12*cr3;
+      ti2=CC(i  ,k,0)+tr11*ci2+tr12*ci3;
+      tr3=CC(i-1,k,0)+tr12*cr2+tr11*cr3;
+      ti3=CC(i  ,k,0)+tr12*ci2+tr11*ci3;
+      MULPM(tr5,tr4,cr5,cr4,ti11,ti12)
+      MULPM(ti5,ti4,ci5,ci4,ti11,ti12)
+      PM(CH(i-1,2,k),CH(ic-1,1,k),tr2,tr5)
+      PM(CH(i  ,2,k),CH(ic  ,1,k),ti5,ti2)
+      PM(CH(i-1,4,k),CH(ic-1,3,k),tr3,tr4)
+      PM(CH(i  ,4,k),CH(ic  ,3,k),ti4,ti3)
+      }
+  }
+
+#undef CC
+#undef CH
+#define C1(a,b,c) cc[(a)+ido*((b)+l1*(c))]
+#define C2(a,b) cc[(a)+idl1*(b)]
+#define CH2(a,b) ch[(a)+idl1*(b)]
+#define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))]
+#define CH(a,b,c) ch[(a)+ido*((b)+l1*(c))]
+NOINLINE static void radfg(size_t ido, size_t ip, size_t l1,
+  double * restrict cc, double * restrict ch, const double * restrict wa,
+  const double * restrict csarr)
+  {
+  const size_t cdim=ip;
+  size_t ipph=(ip+1)/2;
+  size_t idl1 = ido*l1;
+
+  if (ido>1)
+    {
+    for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc)              // 114
+      {
+      size_t is=(j-1)*(ido-1),
+             is2=(jc-1)*(ido-1);
+      for (size_t k=0; k<l1; ++k)                            // 113
+        {
+        size_t idij=is;
+        size_t idij2=is2;
+        for (size_t i=1; i<=ido-2; i+=2)                      // 112
+          {
+          double t1=C1(i,k,j ), t2=C1(i+1,k,j ),
+                 t3=C1(i,k,jc), t4=C1(i+1,k,jc);
+          double x1=wa[idij]*t1 + wa[idij+1]*t2,
+                 x2=wa[idij]*t2 - wa[idij+1]*t1,
+                 x3=wa[idij2]*t3 + wa[idij2+1]*t4,
+                 x4=wa[idij2]*t4 - wa[idij2+1]*t3;
+          C1(i  ,k,j ) = x1+x3;
+          C1(i  ,k,jc) = x2-x4;
+          C1(i+1,k,j ) = x2+x4;
+          C1(i+1,k,jc) = x3-x1;
+          idij+=2;
+          idij2+=2;
+          }
+        }
+      }
+    }
+
+  for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc)                // 123
+    for (size_t k=0; k<l1; ++k)                              // 122
+      {
+      double t1=C1(0,k,j), t2=C1(0,k,jc);
+      C1(0,k,j ) = t1+t2;
+      C1(0,k,jc) = t2-t1;
+      }
+
+//everything in C
+//memset(ch,0,ip*l1*ido*sizeof(double));
+
+  for (size_t l=1,lc=ip-1; l<ipph; ++l,--lc)                 // 127
+    {
+    for (size_t ik=0; ik<idl1; ++ik)                         // 124
+      {
+      CH2(ik,l ) = C2(ik,0)+csarr[2*l]*C2(ik,1)+csarr[4*l]*C2(ik,2);
+      CH2(ik,lc) = csarr[2*l+1]*C2(ik,ip-1)+csarr[4*l+1]*C2(ik,ip-2);
+      }
+    size_t iang = 2*l;
+    size_t j=3, jc=ip-3;
+    for (; j<ipph-3; j+=4,jc-=4)              // 126
+      {
+      iang+=l; if (iang>=ip) iang-=ip;
+      double ar1=csarr[2*iang], ai1=csarr[2*iang+1];
+      iang+=l; if (iang>=ip) iang-=ip;
+      double ar2=csarr[2*iang], ai2=csarr[2*iang+1];
+      iang+=l; if (iang>=ip) iang-=ip;
+      double ar3=csarr[2*iang], ai3=csarr[2*iang+1];
+      iang+=l; if (iang>=ip) iang-=ip;
+      double ar4=csarr[2*iang], ai4=csarr[2*iang+1];
+      for (size_t ik=0; ik<idl1; ++ik)                       // 125
+        {
+        CH2(ik,l ) += ar1*C2(ik,j )+ar2*C2(ik,j +1)
+                     +ar3*C2(ik,j +2)+ar4*C2(ik,j +3);
+        CH2(ik,lc) += ai1*C2(ik,jc)+ai2*C2(ik,jc-1)
+                     +ai3*C2(ik,jc-2)+ai4*C2(ik,jc-3);
+        }
+      }
+    for (; j<ipph-1; j+=2,jc-=2)              // 126
+      {
+      iang+=l; if (iang>=ip) iang-=ip;
+      double ar1=csarr[2*iang], ai1=csarr[2*iang+1];
+      iang+=l; if (iang>=ip) iang-=ip;
+      double ar2=csarr[2*iang], ai2=csarr[2*iang+1];
+      for (size_t ik=0; ik<idl1; ++ik)                       // 125
+        {
+        CH2(ik,l ) += ar1*C2(ik,j )+ar2*C2(ik,j +1);
+        CH2(ik,lc) += ai1*C2(ik,jc)+ai2*C2(ik,jc-1);
+        }
+      }
+    for (; j<ipph; ++j,--jc)              // 126
+      {
+      iang+=l; if (iang>=ip) iang-=ip;
+      double ar=csarr[2*iang], ai=csarr[2*iang+1];
+      for (size_t ik=0; ik<idl1; ++ik)                       // 125
+        {
+        CH2(ik,l ) += ar*C2(ik,j );
+        CH2(ik,lc) += ai*C2(ik,jc);
+        }
+      }
+    }
+  for (size_t ik=0; ik<idl1; ++ik)                         // 101
+    CH2(ik,0) = C2(ik,0);
+  for (size_t j=1; j<ipph; ++j)                              // 129
+    for (size_t ik=0; ik<idl1; ++ik)                         // 128
+      CH2(ik,0) += C2(ik,j);
+
+// everything in CH at this point!
+//memset(cc,0,ip*l1*ido*sizeof(double));
+
+  for (size_t k=0; k<l1; ++k)                                // 131
+    for (size_t i=0; i<ido; ++i)                             // 130
+      CC(i,0,k) = CH(i,k,0);
+
+  for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc)                // 137
+    {
+    size_t j2=2*j-1;
+    for (size_t k=0; k<l1; ++k)                              // 136
+      {
+      CC(ido-1,j2,k) = CH(0,k,j);
+      CC(0,j2+1,k) = CH(0,k,jc);
+      }
+    }
+
+  if (ido==1) return;
+
+  for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc)                // 140
+    {
+    size_t j2=2*j-1;
+    for(size_t k=0; k<l1; ++k)                               // 139
+      for(size_t i=1, ic=ido-i-2; i<=ido-2; i+=2, ic-=2)      // 138
+        {
+        CC(i   ,j2+1,k) = CH(i  ,k,j )+CH(i  ,k,jc);
+        CC(ic  ,j2  ,k) = CH(i  ,k,j )-CH(i  ,k,jc);
+        CC(i+1 ,j2+1,k) = CH(i+1,k,j )+CH(i+1,k,jc);
+        CC(ic+1,j2  ,k) = CH(i+1,k,jc)-CH(i+1,k,j );
+        }
+    }
+  }
+#undef C1
+#undef C2
+#undef CH2
+
+#undef CH
+#undef CC
+#define CH(a,b,c) ch[(a)+ido*((b)+l1*(c))]
+#define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))]
+
+NOINLINE static void radb2(size_t ido, size_t l1, const double * restrict cc,
+  double * restrict ch, const double * restrict wa)
+  {
+  const size_t cdim=2;
+
+  for (size_t k=0; k<l1; k++)
+    PM (CH(0,k,0),CH(0,k,1),CC(0,0,k),CC(ido-1,1,k))
+  if ((ido&1)==0)
+    for (size_t k=0; k<l1; k++)
+      {
+      CH(ido-1,k,0) = 2.*CC(ido-1,0,k);
+      CH(ido-1,k,1) =-2.*CC(0    ,1,k);
+      }
+  if (ido<=2) return;
+  for (size_t k=0; k<l1;++k)
+    for (size_t i=2; i<ido; i+=2)
+      {
+      size_t ic=ido-i;
+      double ti2, tr2;
+      PM (CH(i-1,k,0),tr2,CC(i-1,0,k),CC(ic-1,1,k))
+      PM (ti2,CH(i  ,k,0),CC(i  ,0,k),CC(ic  ,1,k))
+      MULPM (CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),ti2,tr2)
+      }
+  }
+
+NOINLINE static void radb3(size_t ido, size_t l1, const double * restrict cc,
+  double * restrict ch, const double * restrict wa)
+  {
+  const size_t cdim=3;
+  static const double taur=-0.5, taui=0.86602540378443864676;
+
+  for (size_t k=0; k<l1; k++)
+    {
+    double tr2=2.*CC(ido-1,1,k);
+    double cr2=CC(0,0,k)+taur*tr2;
+    CH(0,k,0)=CC(0,0,k)+tr2;
+    double ci3=2.*taui*CC(0,2,k);
+    PM (CH(0,k,2),CH(0,k,1),cr2,ci3);
+    }
+  if (ido==1) return;
+  for (size_t k=0; k<l1; k++)
+    for (size_t i=2; i<ido; i+=2)
+      {
+      size_t ic=ido-i;
+      double tr2=CC(i-1,2,k)+CC(ic-1,1,k); // t2=CC(I) + conj(CC(ic))
+      double ti2=CC(i  ,2,k)-CC(ic  ,1,k);
+      double cr2=CC(i-1,0,k)+taur*tr2;     // c2=CC +taur*t2
+      double ci2=CC(i  ,0,k)+taur*ti2;
+      CH(i-1,k,0)=CC(i-1,0,k)+tr2;         // CH=CC+t2
+      CH(i  ,k,0)=CC(i  ,0,k)+ti2;
+      double cr3=taui*(CC(i-1,2,k)-CC(ic-1,1,k));// c3=taui*(CC(i)-conj(CC(ic)))
+      double ci3=taui*(CC(i  ,2,k)+CC(ic  ,1,k));
+      double di2, di3, dr2, dr3;
+      PM(dr3,dr2,cr2,ci3) // d2= (cr2-ci3, ci2+cr3) = c2+i*c3
+      PM(di2,di3,ci2,cr3) // d3= (cr2+ci3, ci2-cr3) = c2-i*c3
+      MULPM(CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),di2,dr2) // ch = WA*d2
+      MULPM(CH(i,k,2),CH(i-1,k,2),WA(1,i-2),WA(1,i-1),di3,dr3)
+      }
+  }
+
+NOINLINE static void radb4(size_t ido, size_t l1, const double * restrict cc,
+  double * restrict ch, const double * restrict wa)
+  {
+  const size_t cdim=4;
+  static const double sqrt2=1.41421356237309504880;
+
+  for (size_t k=0; k<l1; k++)
+    {
+    double tr1, tr2;
+    PM (tr2,tr1,CC(0,0,k),CC(ido-1,3,k))
+    double tr3=2.*CC(ido-1,1,k);
+    double tr4=2.*CC(0,2,k);
+    PM (CH(0,k,0),CH(0,k,2),tr2,tr3)
+    PM (CH(0,k,3),CH(0,k,1),tr1,tr4)
+    }
+  if ((ido&1)==0)
+    for (size_t k=0; k<l1; k++)
+      {
+      double tr1,tr2,ti1,ti2;
+      PM (ti1,ti2,CC(0    ,3,k),CC(0    ,1,k))
+      PM (tr2,tr1,CC(ido-1,0,k),CC(ido-1,2,k))
+      CH(ido-1,k,0)=tr2+tr2;
+      CH(ido-1,k,1)=sqrt2*(tr1-ti1);
+      CH(ido-1,k,2)=ti2+ti2;
+      CH(ido-1,k,3)=-sqrt2*(tr1+ti1);
+      }
+  if (ido<=2) return;
+  for (size_t k=0; k<l1;++k)
+    for (size_t i=2; i<ido; i+=2)
+      {
+      double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
+      size_t ic=ido-i;
+      PM (tr2,tr1,CC(i-1,0,k),CC(ic-1,3,k))
+      PM (ti1,ti2,CC(i  ,0,k),CC(ic  ,3,k))
+      PM (tr4,ti3,CC(i  ,2,k),CC(ic  ,1,k))
+      PM (tr3,ti4,CC(i-1,2,k),CC(ic-1,1,k))
+      PM (CH(i-1,k,0),cr3,tr2,tr3)
+      PM (CH(i  ,k,0),ci3,ti2,ti3)
+      PM (cr4,cr2,tr1,tr4)
+      PM (ci2,ci4,ti1,ti4)
+      MULPM (CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),ci2,cr2)
+      MULPM (CH(i,k,2),CH(i-1,k,2),WA(1,i-2),WA(1,i-1),ci3,cr3)
+      MULPM (CH(i,k,3),CH(i-1,k,3),WA(2,i-2),WA(2,i-1),ci4,cr4)
+      }
+  }
+
+NOINLINE static void radb5(size_t ido, size_t l1, const double * restrict cc,
+  double * restrict ch, const double * restrict wa)
+  {
+  const size_t cdim=5;
+  static const double tr11= 0.3090169943749474241, ti11=0.95105651629515357212,
+                      tr12=-0.8090169943749474241, ti12=0.58778525229247312917;
+
+  for (size_t k=0; k<l1; k++)
+    {
+    double ti5=CC(0,2,k)+CC(0,2,k);
+    double ti4=CC(0,4,k)+CC(0,4,k);
+    double tr2=CC(ido-1,1,k)+CC(ido-1,1,k);
+    double tr3=CC(ido-1,3,k)+CC(ido-1,3,k);
+    CH(0,k,0)=CC(0,0,k)+tr2+tr3;
+    double cr2=CC(0,0,k)+tr11*tr2+tr12*tr3;
+    double cr3=CC(0,0,k)+tr12*tr2+tr11*tr3;
+    double ci4, ci5;
+    MULPM(ci5,ci4,ti5,ti4,ti11,ti12)
+    PM(CH(0,k,4),CH(0,k,1),cr2,ci5)
+    PM(CH(0,k,3),CH(0,k,2),cr3,ci4)
+    }
+  if (ido==1) return;
+  for (size_t k=0; k<l1;++k)
+    for (size_t i=2; i<ido; i+=2)
+      {
+      size_t ic=ido-i;
+      double tr2, tr3, tr4, tr5, ti2, ti3, ti4, ti5;
+      PM(tr2,tr5,CC(i-1,2,k),CC(ic-1,1,k))
+      PM(ti5,ti2,CC(i  ,2,k),CC(ic  ,1,k))
+      PM(tr3,tr4,CC(i-1,4,k),CC(ic-1,3,k))
+      PM(ti4,ti3,CC(i  ,4,k),CC(ic  ,3,k))
+      CH(i-1,k,0)=CC(i-1,0,k)+tr2+tr3;
+      CH(i  ,k,0)=CC(i  ,0,k)+ti2+ti3;
+      double cr2=CC(i-1,0,k)+tr11*tr2+tr12*tr3;
+      double ci2=CC(i  ,0,k)+tr11*ti2+tr12*ti3;
+      double cr3=CC(i-1,0,k)+tr12*tr2+tr11*tr3;
+      double ci3=CC(i  ,0,k)+tr12*ti2+tr11*ti3;
+      double ci4, ci5, cr5, cr4;
+      MULPM(cr5,cr4,tr5,tr4,ti11,ti12)
+      MULPM(ci5,ci4,ti5,ti4,ti11,ti12)
+      double dr2, dr3, dr4, dr5, di2, di3, di4, di5;
+      PM(dr4,dr3,cr3,ci4)
+      PM(di3,di4,ci3,cr4)
+      PM(dr5,dr2,cr2,ci5)
+      PM(di2,di5,ci2,cr5)
+      MULPM(CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),di2,dr2)
+      MULPM(CH(i,k,2),CH(i-1,k,2),WA(1,i-2),WA(1,i-1),di3,dr3)
+      MULPM(CH(i,k,3),CH(i-1,k,3),WA(2,i-2),WA(2,i-1),di4,dr4)
+      MULPM(CH(i,k,4),CH(i-1,k,4),WA(3,i-2),WA(3,i-1),di5,dr5)
+      }
+  }
+
+#undef CC
+#undef CH
+#define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))]
+#define CH(a,b,c) ch[(a)+ido*((b)+l1*(c))]
+#define C1(a,b,c) cc[(a)+ido*((b)+l1*(c))]
+#define C2(a,b) cc[(a)+idl1*(b)]
+#define CH2(a,b) ch[(a)+idl1*(b)]
+
+NOINLINE static void radbg(size_t ido, size_t ip, size_t l1,
+  double * restrict cc, double * restrict ch, const double * restrict wa,
+  const double * restrict csarr)
+  {
+  const size_t cdim=ip;
+  size_t ipph=(ip+1)/ 2;
+  size_t idl1 = ido*l1;
+
+  for (size_t k=0; k<l1; ++k)        // 102
+    for (size_t i=0; i<ido; ++i)     // 101
+      CH(i,k,0) = CC(i,0,k);
+  for (size_t j=1, jc=ip-1; j<ipph; ++j, --jc)   // 108
+    {
+    size_t j2=2*j-1;
+    for (size_t k=0; k<l1; ++k)
+      {
+      CH(0,k,j ) = 2*CC(ido-1,j2,k);
+      CH(0,k,jc) = 2*CC(0,j2+1,k);
+      }
+    }
+
+  if (ido!=1)
+    {
+    for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc)   // 111
+      {
+      size_t j2=2*j-1;
+      for (size_t k=0; k<l1; ++k)
+        for (size_t i=1, ic=ido-i-2; i<=ido-2; i+=2, ic-=2)      // 109
+          {
+          CH(i  ,k,j ) = CC(i  ,j2+1,k)+CC(ic  ,j2,k);
+          CH(i  ,k,jc) = CC(i  ,j2+1,k)-CC(ic  ,j2,k);
+          CH(i+1,k,j ) = CC(i+1,j2+1,k)-CC(ic+1,j2,k);
+          CH(i+1,k,jc) = CC(i+1,j2+1,k)+CC(ic+1,j2,k);
+          }
+      }
+    }
+  for (size_t l=1,lc=ip-1; l<ipph; ++l,--lc)
+    {
+    for (size_t ik=0; ik<idl1; ++ik)
+      {
+      C2(ik,l ) = CH2(ik,0)+csarr[2*l]*CH2(ik,1)+csarr[4*l]*CH2(ik,2);
+      C2(ik,lc) = csarr[2*l+1]*CH2(ik,ip-1)+csarr[4*l+1]*CH2(ik,ip-2);
+      }
+    size_t iang=2*l;
+    size_t j=3,jc=ip-3;
+    for(; j<ipph-3; j+=4,jc-=4)
+      {
+      iang+=l; if(iang>ip) iang-=ip;
+      double ar1=csarr[2*iang], ai1=csarr[2*iang+1];
+      iang+=l; if(iang>ip) iang-=ip;
+      double ar2=csarr[2*iang], ai2=csarr[2*iang+1];
+      iang+=l; if(iang>ip) iang-=ip;
+      double ar3=csarr[2*iang], ai3=csarr[2*iang+1];
+      iang+=l; if(iang>ip) iang-=ip;
+      double ar4=csarr[2*iang], ai4=csarr[2*iang+1];
+      for (size_t ik=0; ik<idl1; ++ik)
+        {
+        C2(ik,l ) += ar1*CH2(ik,j )+ar2*CH2(ik,j +1)
+                    +ar3*CH2(ik,j +2)+ar4*CH2(ik,j +3);
+        C2(ik,lc) += ai1*CH2(ik,jc)+ai2*CH2(ik,jc-1)
+                    +ai3*CH2(ik,jc-2)+ai4*CH2(ik,jc-3);
+        }
+      }
+    for(; j<ipph-1; j+=2,jc-=2)
+      {
+      iang+=l; if(iang>ip) iang-=ip;
+      double ar1=csarr[2*iang], ai1=csarr[2*iang+1];
+      iang+=l; if(iang>ip) iang-=ip;
+      double ar2=csarr[2*iang], ai2=csarr[2*iang+1];
+      for (size_t ik=0; ik<idl1; ++ik)
+        {
+        C2(ik,l ) += ar1*CH2(ik,j )+ar2*CH2(ik,j +1);
+        C2(ik,lc) += ai1*CH2(ik,jc)+ai2*CH2(ik,jc-1);
+        }
+      }
+    for(; j<ipph; ++j,--jc)
+      {
+      iang+=l; if(iang>ip) iang-=ip;
+      double war=csarr[2*iang], wai=csarr[2*iang+1];
+      for (size_t ik=0; ik<idl1; ++ik)
+        {
+        C2(ik,l ) += war*CH2(ik,j );
+        C2(ik,lc) += wai*CH2(ik,jc);
+        }
+      }
+    }
+  for (size_t j=1; j<ipph; ++j)
+    for (size_t ik=0; ik<idl1; ++ik)
+      CH2(ik,0) += CH2(ik,j);
+  for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc)   // 124
+    for (size_t k=0; k<l1; ++k)
+      {
+      CH(0,k,j ) = C1(0,k,j)-C1(0,k,jc);
+      CH(0,k,jc) = C1(0,k,j)+C1(0,k,jc);
+      }
+
+  if (ido==1) return;
+
+  for (size_t j=1, jc=ip-1; j<ipph; ++j, --jc)  // 127
+    for (size_t k=0; k<l1; ++k)
+      for (size_t i=1; i<=ido-2; i+=2)
+        {
+        CH(i  ,k,j ) = C1(i  ,k,j)-C1(i+1,k,jc);
+        CH(i  ,k,jc) = C1(i  ,k,j)+C1(i+1,k,jc);
+        CH(i+1,k,j ) = C1(i+1,k,j)+C1(i  ,k,jc);
+        CH(i+1,k,jc) = C1(i+1,k,j)-C1(i  ,k,jc);
+        }
+
+// All in CH
+
+  for (size_t j=1; j<ip; ++j)
+    {
+    size_t is = (j-1)*(ido-1);
+    for (size_t k=0; k<l1; ++k)
+      {
+      size_t idij = is;
+      for (size_t i=1; i<=ido-2; i+=2)
+        {
+        double t1=CH(i,k,j), t2=CH(i+1,k,j);
+        CH(i  ,k,j) = wa[idij]*t1-wa[idij+1]*t2;
+        CH(i+1,k,j) = wa[idij]*t2+wa[idij+1]*t1;
+        idij+=2;
+        }
+      }
+    }
+  }
+#undef C1
+#undef C2
+#undef CH2
+
+#undef CC
+#undef CH
+#undef PM
+#undef MULPM
+#undef WA
+
+static void copy_and_norm(double *c, double *p1, size_t n, double fct)
+  {
+  if (p1!=c)
+    {
+    if (fct!=1.)
+      for (size_t i=0; i<n; ++i)
+        c[i] = fct*p1[i];
+    else
+      memcpy (c,p1,n*sizeof(double));
+    }
+  else
+    if (fct!=1.)
+      for (size_t i=0; i<n; ++i)
+        c[i] *= fct;
+  }
+
+WARN_UNUSED_RESULT
+static int rfftp_forward(rfftp_plan plan, double c[], double fct)
+  {
+  if (plan->length==1) return 0;
+  size_t n=plan->length;
+  size_t l1=n, nf=plan->nfct;
+  double *ch = RALLOC(double, n);
+  if (!ch) return -1;
+  double *p1=c, *p2=ch;
+
+  for(size_t k1=0; k1<nf;++k1)
+    {
+    size_t k=nf-k1-1;
+    size_t ip=plan->fct[k].fct;
+    size_t ido=n / l1;
+    l1 /= ip;
+    if(ip==4)
+      radf4(ido, l1, p1, p2, plan->fct[k].tw);
+    else if(ip==2)
+      radf2(ido, l1, p1, p2, plan->fct[k].tw);
+    else if(ip==3)
+      radf3(ido, l1, p1, p2, plan->fct[k].tw);
+    else if(ip==5)
+      radf5(ido, l1, p1, p2, plan->fct[k].tw);
+    else
+      {
+      radfg(ido, ip, l1, p1, p2, plan->fct[k].tw, plan->fct[k].tws);
+      SWAP (p1,p2,double *);
+      }
+    SWAP (p1,p2,double *);
+    }
+  copy_and_norm(c,p1,n,fct);
+  DEALLOC(ch);
+  return 0;
+  }
+
+WARN_UNUSED_RESULT
+static int rfftp_backward(rfftp_plan plan, double c[], double fct)
+  {
+  if (plan->length==1) return 0;
+  size_t n=plan->length;
+  size_t l1=1, nf=plan->nfct;
+  double *ch = RALLOC(double, n);
+  if (!ch) return -1;
+  double *p1=c, *p2=ch;
+
+  for(size_t k=0; k<nf; k++)
+    {
+    size_t ip = plan->fct[k].fct,
+           ido= n/(ip*l1);
+    if(ip==4)
+      radb4(ido, l1, p1, p2, plan->fct[k].tw);
+    else if(ip==2)
+      radb2(ido, l1, p1, p2, plan->fct[k].tw);
+    else if(ip==3)
+      radb3(ido, l1, p1, p2, plan->fct[k].tw);
+    else if(ip==5)
+      radb5(ido, l1, p1, p2, plan->fct[k].tw);
+    else
+      radbg(ido, ip, l1, p1, p2, plan->fct[k].tw, plan->fct[k].tws);
+    SWAP (p1,p2,double *);
+    l1*=ip;
+    }
+  copy_and_norm(c,p1,n,fct);
+  DEALLOC(ch);
+  return 0;
+  }
+
+WARN_UNUSED_RESULT
+static int rfftp_factorize (rfftp_plan plan)
+  {
+  size_t length=plan->length;
+  size_t nfct=0;
+  while ((length%4)==0)
+    { if (nfct>=NFCT) return -1; plan->fct[nfct++].fct=4; length>>=2; }
+  if ((length%2)==0)
+    {
+    length>>=1;
+    // factor 2 should be at the front of the factor list
+    if (nfct>=NFCT) return -1;
+    plan->fct[nfct++].fct=2;
+    SWAP(plan->fct[0].fct, plan->fct[nfct-1].fct,size_t);
+    }
+  size_t maxl=(size_t)(sqrt((double)length))+1;
+  for (size_t divisor=3; (length>1)&&(divisor<maxl); divisor+=2)
+    if ((length%divisor)==0)
+      {
+      while ((length%divisor)==0)
+        {
+        if (nfct>=NFCT) return -1;
+        plan->fct[nfct++].fct=divisor;
+        length/=divisor;
+        }
+      maxl=(size_t)(sqrt((double)length))+1;
+      }
+  if (length>1) plan->fct[nfct++].fct=length;
+  plan->nfct=nfct;
+  return 0;
+  }
+
+static size_t rfftp_twsize(rfftp_plan plan)
+  {
+  size_t twsize=0, l1=1;
+  for (size_t k=0; k<plan->nfct; ++k)
+    {
+    size_t ip=plan->fct[k].fct, ido= plan->length/(l1*ip);
+    twsize+=(ip-1)*(ido-1);
+    if (ip>5) twsize+=2*ip;
+    l1*=ip;
+    }
+  return twsize;
+  return 0;
+  }
+
+WARN_UNUSED_RESULT NOINLINE static int rfftp_comp_twiddle (rfftp_plan plan)
+  {
+  size_t length=plan->length;
+  double *twid = RALLOC(double, 2*length);
+  if (!twid) return -1;
+  sincos_2pibyn_half(length, twid);
+  size_t l1=1;
+  double *ptr=plan->mem;
+  for (size_t k=0; k<plan->nfct; ++k)
+    {
+    size_t ip=plan->fct[k].fct, ido=length/(l1*ip);
+    if (k<plan->nfct-1) // last factor doesn't need twiddles
+      {
+      plan->fct[k].tw=ptr; ptr+=(ip-1)*(ido-1);
+      for (size_t j=1; j<ip; ++j)
+        for (size_t i=1; i<=(ido-1)/2; ++i)
+          {
+          plan->fct[k].tw[(j-1)*(ido-1)+2*i-2] = twid[2*j*l1*i];
+          plan->fct[k].tw[(j-1)*(ido-1)+2*i-1] = twid[2*j*l1*i+1];
+          }
+      }
+    if (ip>5) // special factors required by *g functions
+      {
+      plan->fct[k].tws=ptr; ptr+=2*ip;
+      plan->fct[k].tws[0] = 1.;
+      plan->fct[k].tws[1] = 0.;
+      for (size_t i=1; i<=(ip>>1); ++i)
+        {
+        plan->fct[k].tws[2*i  ] = twid[2*i*(length/ip)];
+        plan->fct[k].tws[2*i+1] = twid[2*i*(length/ip)+1];
+        plan->fct[k].tws[2*(ip-i)  ] = twid[2*i*(length/ip)];
+        plan->fct[k].tws[2*(ip-i)+1] = -twid[2*i*(length/ip)+1];
+        }
+      }
+    l1*=ip;
+    }
+  DEALLOC(twid);
+  return 0;
+  }
+
+NOINLINE static rfftp_plan make_rfftp_plan (size_t length)
+  {
+  if (length==0) return NULL;
+  rfftp_plan plan = RALLOC(rfftp_plan_i,1);
+  if (!plan) return NULL;
+  plan->length=length;
+  plan->nfct=0;
+  plan->mem=NULL;
+  for (size_t i=0; i<NFCT; ++i)
+    plan->fct[i]=(rfftp_fctdata){0,0,0};
+  if (length==1) return plan;
+  if (rfftp_factorize(plan)!=0) { DEALLOC(plan); return NULL; }
+  size_t tws=rfftp_twsize(plan);
+  plan->mem=RALLOC(double,tws);
+  if (!plan->mem) { DEALLOC(plan); return NULL; }
+  if (rfftp_comp_twiddle(plan)!=0)
+    { DEALLOC(plan->mem); DEALLOC(plan); return NULL; }
+  return plan;
+  }
+
+NOINLINE static void destroy_rfftp_plan (rfftp_plan plan)
+  {
+  DEALLOC(plan->mem);
+  DEALLOC(plan);
+  }
+
+typedef struct fftblue_plan_i
+  {
+  size_t n, n2;
+  cfftp_plan plan;
+  double *mem;
+  double *bk, *bkf;
+  } fftblue_plan_i;
+typedef struct fftblue_plan_i * fftblue_plan;
+
+NOINLINE static fftblue_plan make_fftblue_plan (size_t length)
+  {
+  fftblue_plan plan = RALLOC(fftblue_plan_i,1);
+  if (!plan) return NULL;
+  plan->n = length;
+  plan->n2 = good_size(plan->n*2-1);
+  plan->mem = RALLOC(double, 2*plan->n+2*plan->n2);
+  if (!plan->mem) { DEALLOC(plan); return NULL; }
+  plan->bk  = plan->mem;
+  plan->bkf = plan->bk+2*plan->n;
+
+/* initialize b_k */
+  double *tmp = RALLOC(double,4*plan->n);
+  if (!tmp) { DEALLOC(plan->mem); DEALLOC(plan); return NULL; }
+  sincos_2pibyn(2*plan->n,tmp);
+  plan->bk[0] = 1;
+  plan->bk[1] = 0;
+
+  size_t coeff=0;
+  for (size_t m=1; m<plan->n; ++m)
+    {
+    coeff+=2*m-1;
+    if (coeff>=2*plan->n) coeff-=2*plan->n;
+    plan->bk[2*m  ] = tmp[2*coeff  ];
+    plan->bk[2*m+1] = tmp[2*coeff+1];
+    }
+
+  /* initialize the zero-padded, Fourier transformed b_k. Add normalisation. */
+  double xn2 = 1./plan->n2;
+  plan->bkf[0] = plan->bk[0]*xn2;
+  plan->bkf[1] = plan->bk[1]*xn2;
+  for (size_t m=2; m<2*plan->n; m+=2)
+    {
+    plan->bkf[m]   = plan->bkf[2*plan->n2-m]   = plan->bk[m]   *xn2;
+    plan->bkf[m+1] = plan->bkf[2*plan->n2-m+1] = plan->bk[m+1] *xn2;
+    }
+  for (size_t m=2*plan->n;m<=(2*plan->n2-2*plan->n+1);++m)
+    plan->bkf[m]=0.;
+  plan->plan=make_cfftp_plan(plan->n2);
+  if (!plan->plan)
+    { DEALLOC(tmp); DEALLOC(plan->mem); DEALLOC(plan); return NULL; }
+  if (cfftp_forward(plan->plan,plan->bkf,1.)!=0)
+    { DEALLOC(tmp); DEALLOC(plan->mem); DEALLOC(plan); return NULL; }
+  DEALLOC(tmp);
+
+  return plan;
+  }
+
+NOINLINE static void destroy_fftblue_plan (fftblue_plan plan)
+  {
+  DEALLOC(plan->mem);
+  destroy_cfftp_plan(plan->plan);
+  DEALLOC(plan);
+  }
+
+NOINLINE WARN_UNUSED_RESULT
+static int fftblue_fft(fftblue_plan plan, double c[], int isign, double fct)
+  {
+  size_t n=plan->n;
+  size_t n2=plan->n2;
+  double *bk  = plan->bk;
+  double *bkf = plan->bkf;
+  double *akf = RALLOC(double, 2*n2);
+  if (!akf) return -1;
+
+/* initialize a_k and FFT it */
+  if (isign>0)
+    for (size_t m=0; m<2*n; m+=2)
+      {
+      akf[m]   = c[m]*bk[m]   - c[m+1]*bk[m+1];
+      akf[m+1] = c[m]*bk[m+1] + c[m+1]*bk[m];
+      }
+  else
+    for (size_t m=0; m<2*n; m+=2)
+      {
+      akf[m]   = c[m]*bk[m]   + c[m+1]*bk[m+1];
+      akf[m+1] =-c[m]*bk[m+1] + c[m+1]*bk[m];
+      }
+  for (size_t m=2*n; m<2*n2; ++m)
+    akf[m]=0;
+
+  if (cfftp_forward (plan->plan,akf,fct)!=0)
+    { DEALLOC(akf); return -1; }
+
+/* do the convolution */
+  if (isign>0)
+    for (size_t m=0; m<2*n2; m+=2)
+      {
+      double im = -akf[m]*bkf[m+1] + akf[m+1]*bkf[m];
+      akf[m  ]  =  akf[m]*bkf[m]   + akf[m+1]*bkf[m+1];
+      akf[m+1]  = im;
+      }
+  else
+    for (size_t m=0; m<2*n2; m+=2)
+      {
+      double im = akf[m]*bkf[m+1] + akf[m+1]*bkf[m];
+      akf[m  ]  = akf[m]*bkf[m]   - akf[m+1]*bkf[m+1];
+      akf[m+1]  = im;
+      }
+
+/* inverse FFT */
+  if (cfftp_backward (plan->plan,akf,1.)!=0)
+    { DEALLOC(akf); return -1; }
+
+/* multiply by b_k */
+  if (isign>0)
+    for (size_t m=0; m<2*n; m+=2)
+      {
+      c[m]   = bk[m]  *akf[m] - bk[m+1]*akf[m+1];
+      c[m+1] = bk[m+1]*akf[m] + bk[m]  *akf[m+1];
+      }
+  else
+    for (size_t m=0; m<2*n; m+=2)
+      {
+      c[m]   = bk[m]  *akf[m] + bk[m+1]*akf[m+1];
+      c[m+1] =-bk[m+1]*akf[m] + bk[m]  *akf[m+1];
+      }
+  DEALLOC(akf);
+  return 0;
+  }
+
+WARN_UNUSED_RESULT
+static int cfftblue_backward(fftblue_plan plan, double c[], double fct)
+  { return fftblue_fft(plan,c,1,fct); }
+
+WARN_UNUSED_RESULT
+static int cfftblue_forward(fftblue_plan plan, double c[], double fct)
+  { return fftblue_fft(plan,c,-1,fct); }
+
+WARN_UNUSED_RESULT
+static int rfftblue_backward(fftblue_plan plan, double c[], double fct)
+  {
+  size_t n=plan->n;
+  double *tmp = RALLOC(double,2*n);
+  if (!tmp) return -1;
+  tmp[0]=c[0];
+  tmp[1]=0.;
+  memcpy (tmp+2,c+1, (n-1)*sizeof(double));
+  if ((n&1)==0) tmp[n+1]=0.;
+  for (size_t m=2; m<n; m+=2)
+    {
+    tmp[2*n-m]=tmp[m];
+    tmp[2*n-m+1]=-tmp[m+1];
+    }
+  if (fftblue_fft(plan,tmp,1,fct)!=0)
+    { DEALLOC(tmp); return -1; }
+  for (size_t m=0; m<n; ++m)
+    c[m] = tmp[2*m];
+  DEALLOC(tmp);
+  return 0;
+  }
+
+WARN_UNUSED_RESULT
+static int rfftblue_forward(fftblue_plan plan, double c[], double fct)
+  {
+  size_t n=plan->n;
+  double *tmp = RALLOC(double,2*n);
+  if (!tmp) return -1;
+  for (size_t m=0; m<n; ++m)
+    {
+    tmp[2*m] = c[m];
+    tmp[2*m+1] = 0.;
+    }
+  if (fftblue_fft(plan,tmp,-1,fct)!=0)
+    { DEALLOC(tmp); return -1; }
+  c[0] = tmp[0];
+  memcpy (c+1, tmp+2, (n-1)*sizeof(double));
+  DEALLOC(tmp);
+  return 0;
+  }
+
+typedef struct cfft_plan_i
+  {
+  cfftp_plan packplan;
+  fftblue_plan blueplan;
+  } cfft_plan_i;
+
+static cfft_plan make_cfft_plan (size_t length)
+  {
+  if (length==0) return NULL;
+  cfft_plan plan = RALLOC(cfft_plan_i,1);
+  if (!plan) return NULL;
+  plan->blueplan=0;
+  plan->packplan=0;
+  if ((length<50) || (largest_prime_factor(length)<=sqrt(length)))
+    {
+    plan->packplan=make_cfftp_plan(length);
+    if (!plan->packplan) { DEALLOC(plan); return NULL; }
+    return plan;
+    }
+  double comp1 = cost_guess(length);
+  double comp2 = 2*cost_guess(good_size(2*length-1));
+  comp2*=1.5; /* fudge factor that appears to give good overall performance */
+  if (comp2<comp1) // use Bluestein
+    {
+    plan->blueplan=make_fftblue_plan(length);
+    if (!plan->blueplan) { DEALLOC(plan); return NULL; }
+    }
+  else
+    {
+    plan->packplan=make_cfftp_plan(length);
+    if (!plan->packplan) { DEALLOC(plan); return NULL; }
+    }
+  return plan;
+  }
+
+static void destroy_cfft_plan (cfft_plan plan)
+  {
+  if (plan->blueplan)
+    destroy_fftblue_plan(plan->blueplan);
+  if (plan->packplan)
+    destroy_cfftp_plan(plan->packplan);
+  DEALLOC(plan);
+  }
+
+WARN_UNUSED_RESULT static int cfft_backward(cfft_plan plan, double c[], double fct)
+  {
+  if (plan->packplan)
+    return cfftp_backward(plan->packplan,c,fct);
+  // if (plan->blueplan)
+  return cfftblue_backward(plan->blueplan,c,fct);
+  }
+
+WARN_UNUSED_RESULT static int cfft_forward(cfft_plan plan, double c[], double fct)
+  {
+  if (plan->packplan)
+    return cfftp_forward(plan->packplan,c,fct);
+  // if (plan->blueplan)
+  return cfftblue_forward(plan->blueplan,c,fct);
+  }
+
+typedef struct rfft_plan_i
+  {
+  rfftp_plan packplan;
+  fftblue_plan blueplan;
+  } rfft_plan_i;
+
+static rfft_plan make_rfft_plan (size_t length)
+  {
+  if (length==0) return NULL;
+  rfft_plan plan = RALLOC(rfft_plan_i,1);
+  if (!plan) return NULL;
+  plan->blueplan=0;
+  plan->packplan=0;
+  if ((length<50) || (largest_prime_factor(length)<=sqrt(length)))
+    {
+    plan->packplan=make_rfftp_plan(length);
+    if (!plan->packplan) { DEALLOC(plan); return NULL; }
+    return plan;
+    }
+  double comp1 = 0.5*cost_guess(length);
+  double comp2 = 2*cost_guess(good_size(2*length-1));
+  comp2*=1.5; /* fudge factor that appears to give good overall performance */
+  if (comp2<comp1) // use Bluestein
+    {
+    plan->blueplan=make_fftblue_plan(length);
+    if (!plan->blueplan) { DEALLOC(plan); return NULL; }
+    }
+  else
+    {
+    plan->packplan=make_rfftp_plan(length);
+    if (!plan->packplan) { DEALLOC(plan); return NULL; }
+    }
+  return plan;
+  }
+
+static void destroy_rfft_plan (rfft_plan plan)
+  {
+  if (plan->blueplan)
+    destroy_fftblue_plan(plan->blueplan);
+  if (plan->packplan)
+    destroy_rfftp_plan(plan->packplan);
+  DEALLOC(plan);
+  }
+
+WARN_UNUSED_RESULT static int rfft_backward(rfft_plan plan, double c[], double fct)
+  {
+  if (plan->packplan)
+    return rfftp_backward(plan->packplan,c,fct);
+  else // if (plan->blueplan)
+    return rfftblue_backward(plan->blueplan,c,fct);
+  }
+
+WARN_UNUSED_RESULT static int rfft_forward(rfft_plan plan, double c[], double fct)
+  {
+  if (plan->packplan)
+    return rfftp_forward(plan->packplan,c,fct);
+  else // if (plan->blueplan)
+    return rfftblue_forward(plan->blueplan,c,fct);
+  }
+
+#define NPY_NO_DEPRECATED_API NPY_API_VERSION
+
+#include "Python.h"
+#include "numpy/arrayobject.h"
+
+static PyObject *
+execute_complex(PyObject *a1, int is_forward, double fct)
+{
+    PyArrayObject *data = (PyArrayObject *)PyArray_CopyFromObject(a1, NPY_CDOUBLE, 1, 0);
+    if (!data) return NULL;
+
+    int npts = PyArray_DIM(data, PyArray_NDIM(data) - 1);
+    cfft_plan plan=NULL;
+
+    int nrepeats = PyArray_SIZE(data)/npts;
+    double *dptr = (double *)PyArray_DATA(data);
+    int fail=0;
+    Py_BEGIN_ALLOW_THREADS;
+    NPY_SIGINT_ON;
+    plan = make_cfft_plan(npts);
+    if (!plan) fail=1;
+    if (!fail)
+      for (int i = 0; i < nrepeats; i++) {
+          int res = is_forward ?
+            cfft_forward(plan, dptr, fct) : cfft_backward(plan, dptr, fct);
+          if (res!=0) { fail=1; break; }
+          dptr += npts*2;
+      }
+    if (plan) destroy_cfft_plan(plan);
+    NPY_SIGINT_OFF;
+    Py_END_ALLOW_THREADS;
+    if (fail) {
+      Py_XDECREF(data);
+      return PyErr_NoMemory();
+    }
+    return (PyObject *)data;
+}
+
+static PyObject *
+execute_real_forward(PyObject *a1, double fct)
+{
+    rfft_plan plan=NULL;
+    int fail = 0;
+    PyArrayObject *data = (PyArrayObject *)PyArray_ContiguousFromObject(a1,
+            NPY_DOUBLE, 1, 0);
+    if (!data) return NULL;
+
+    int ndim = PyArray_NDIM(data);
+    const npy_intp *odim = PyArray_DIMS(data);
+    int npts = odim[ndim - 1];
+    npy_intp *tdim=(npy_intp *)malloc(ndim*sizeof(npy_intp));
+    if (!tdim)
+      { Py_XDECREF(data); return NULL; }
+    for (int d=0; d<ndim-1; ++d)
+      tdim[d] = odim[d];
+    tdim[ndim-1] = npts/2 + 1;
+    PyArrayObject *ret = (PyArrayObject *)PyArray_Empty(ndim,
+            tdim, PyArray_DescrFromType(NPY_CDOUBLE), 0);
+    free(tdim);
+    if (!ret) fail=1;
+    if (!fail) {
+      int rstep = PyArray_DIM(ret, PyArray_NDIM(ret) - 1)*2;
+
+      int nrepeats = PyArray_SIZE(data)/npts;
+      double *rptr = (double *)PyArray_DATA(ret),
+             *dptr = (double *)PyArray_DATA(data);
+
+      Py_BEGIN_ALLOW_THREADS;
+      NPY_SIGINT_ON;
+      plan = make_rfft_plan(npts);
+      if (!plan) fail=1;
+      if (!fail)
+        for (int i = 0; i < nrepeats; i++) {
+            rptr[rstep-1] = 0.0;
+            memcpy((char *)(rptr+1), dptr, npts*sizeof(double));
+            if (rfft_forward(plan, rptr+1, fct)!=0) {fail=1; break;}
+            rptr[0] = rptr[1];
+            rptr[1] = 0.0;
+            rptr += rstep;
+            dptr += npts;
+      }
+      if (plan) destroy_rfft_plan(plan);
+      NPY_SIGINT_OFF;
+      Py_END_ALLOW_THREADS;
+    }
+    if (fail) {
+      Py_XDECREF(data);
+      Py_XDECREF(ret);
+      return PyErr_NoMemory();
+    }
+    Py_DECREF(data);
+    return (PyObject *)ret;
+}
+static PyObject *
+execute_real_backward(PyObject *a1, double fct)
+{
+    rfft_plan plan=NULL;
+    PyArrayObject *data = (PyArrayObject *)PyArray_ContiguousFromObject(a1,
+            NPY_CDOUBLE, 1, 0);
+    if (!data) return NULL;
+    int npts = PyArray_DIM(data, PyArray_NDIM(data) - 1);
+    PyArrayObject *ret = (PyArrayObject *)PyArray_Empty(PyArray_NDIM(data),
+            PyArray_DIMS(data), PyArray_DescrFromType(NPY_DOUBLE), 0);
+    int fail = 0;
+    if (!ret) fail=1;
+    if (!fail) {
+      int nrepeats = PyArray_SIZE(ret)/npts;
+      double *rptr = (double *)PyArray_DATA(ret),
+             *dptr = (double *)PyArray_DATA(data);
+
+      Py_BEGIN_ALLOW_THREADS;
+      NPY_SIGINT_ON;
+      plan = make_rfft_plan(npts);
+      if (!plan) fail=1;
+      if (!fail) {
+        for (int i = 0; i < nrepeats; i++) {
+          memcpy((char *)(rptr + 1), (dptr + 2), (npts - 1)*sizeof(double));
+          rptr[0] = dptr[0];
+          if (rfft_backward(plan, rptr, fct)!=0) {fail=1; break;}
+          rptr += npts;
+          dptr += npts*2;
+        }
+      }
+      if (plan) destroy_rfft_plan(plan);
+      NPY_SIGINT_OFF;
+      Py_END_ALLOW_THREADS;
+    }
+    if (fail) {
+      Py_XDECREF(data);
+      Py_XDECREF(ret);
+      return PyErr_NoMemory();
+    }
+    Py_DECREF(data);
+    return (PyObject *)ret;
+}
+
+static PyObject *
+execute_real(PyObject *a1, int is_forward, double fct)
+{
+    return is_forward ? execute_real_forward(a1, fct)
+                      : execute_real_backward(a1, fct);
+}
+
+static const char execute__doc__[] = "";
+
+static PyObject *
+execute(PyObject *NPY_UNUSED(self), PyObject *args)
+{
+    PyObject *a1;
+    int is_real, is_forward;
+    double fct;
+
+    if(!PyArg_ParseTuple(args, "Oiid:execute", &a1, &is_real, &is_forward, &fct)) {
+        return NULL;
+    }
+
+    return is_real ? execute_real(a1, is_forward, fct)
+                   : execute_complex(a1, is_forward, fct);
+}
+
+/* List of methods defined in the module */
+
+static struct PyMethodDef methods[] = {
+    {"execute",   execute,   1, execute__doc__},
+    {NULL, NULL, 0, NULL}          /* sentinel */
+};
+
+#if PY_MAJOR_VERSION >= 3
+static struct PyModuleDef moduledef = {
+        PyModuleDef_HEAD_INIT,
+        "pocketfft_internal",
+        NULL,
+        -1,
+        methods,
+        NULL,
+        NULL,
+        NULL,
+        NULL
+};
+#endif
+
+/* Initialization function for the module */
+#if PY_MAJOR_VERSION >= 3
+#define RETVAL(x) x
+PyMODINIT_FUNC PyInit_pocketfft_internal(void)
+#else
+#define RETVAL(x)
+PyMODINIT_FUNC
+initpocketfft_internal(void)
+#endif
+{
+    PyObject *m;
+#if PY_MAJOR_VERSION >= 3
+    m = PyModule_Create(&moduledef);
+#else
+    static const char module_documentation[] = "";
+
+    m = Py_InitModule4("pocketfft_internal", methods,
+            module_documentation,
+            (PyObject*)NULL,PYTHON_API_VERSION);
+#endif
+    if (m == NULL) {
+        return RETVAL(NULL);
+    }
+
+    /* Import the array object */
+    import_array();
+
+    /* XXXX Add constants here */
+
+    return RETVAL(m);
+}
diff --git a/numpy/fft/fftpack.py b/numpy/fft/pocketfft.py
index d0df6fb48..794d13937 100644
--- a/numpy/fft/fftpack.py
+++ b/numpy/fft/pocketfft.py
@@ -26,9 +26,6 @@ n = n-dimensional transform
 (Note: 2D routines are just nD routines with different default
 behavior.)
 
-The underlying code for these functions is an f2c-translated and modified
-version of the FFTPACK routines.
-
 """
 from __future__ import division, absolute_import, print_function
 
@@ -37,26 +34,18 @@ __all__ = ['fft', 'ifft', 'rfft', 'irfft', 'hfft', 'ihfft', 'rfftn',
 
 import functools
 
-from numpy.core import (array, asarray, zeros, swapaxes, shape, conjugate,
-                        take, sqrt)
+from numpy.core import asarray, zeros, swapaxes, conjugate, take, sqrt
+from . import pocketfft_internal as pfi
 from numpy.core.multiarray import normalize_axis_index
 from numpy.core import overrides
-from . import fftpack_lite as fftpack
-from .helper import _FFTCache
-
-_fft_cache = _FFTCache(max_size_in_mb=100, max_item_count=32)
-_real_fft_cache = _FFTCache(max_size_in_mb=100, max_item_count=32)
 
 
 array_function_dispatch = functools.partial(
     overrides.array_function_dispatch, module='numpy.fft')
 
 
-def _raw_fft(a, n=None, axis=-1, init_function=fftpack.cffti,
-             work_function=fftpack.cfftf, fft_cache=_fft_cache):
-    a = asarray(a)
+def _raw_fft(a, n, axis, is_real, is_forward, fct):
     axis = normalize_axis_index(axis, a.ndim)
-
     if n is None:
         n = a.shape[axis]
 
@@ -64,15 +53,6 @@ def _raw_fft(a, n=None, axis=-1, init_function=fftpack.cffti,
         raise ValueError("Invalid number of FFT data points (%d) specified."
                          % n)
 
-    # We have to ensure that only a single thread can access a wsave array
-    # at any given time. Thus we remove it from the cache and insert it
-    # again after it has been used. Multiple threads might create multiple
-    # copies of the wsave array. This is intentional and a limitation of
-    # the current C code.
-    wsave = fft_cache.pop_twiddle_factors(n)
-    if wsave is None:
-        wsave = init_function(n)
-
     if a.shape[axis] != n:
         s = list(a.shape)
         if s[axis] > n:
@@ -87,25 +67,22 @@ def _raw_fft(a, n=None, axis=-1, init_function=fftpack.cffti,
             z[tuple(index)] = a
             a = z
 
-    if axis != a.ndim - 1:
+    if axis == a.ndim-1:
+        r = pfi.execute(a, is_real, is_forward, fct)
+    else:
         a = swapaxes(a, axis, -1)
-    r = work_function(a, wsave)
-    if axis != a.ndim - 1:
+        r = pfi.execute(a, is_real, is_forward, fct)
         r = swapaxes(r, axis, -1)
-
-    # As soon as we put wsave back into the cache, another thread could pick it
-    # up and start using it, so we must not do this until after we're
-    # completely done using it ourselves.
-    fft_cache.put_twiddle_factors(n, wsave)
-
     return r
 
 
 def _unitary(norm):
-    if norm not in (None, "ortho"):
-        raise ValueError("Invalid norm value %s, should be None or \"ortho\"."
-                         % norm)
-    return norm is not None
+    if norm is None:
+        return False
+    if norm=="ortho":
+        return True
+    raise ValueError("Invalid norm value %s, should be None or \"ortho\"."
+                     % norm)
 
 
 def _fft_dispatcher(a, n=None, axis=None, norm=None):
@@ -198,12 +175,13 @@ def fft(a, n=None, axis=-1, norm=None):
 
     """
 
-    a = asarray(a).astype(complex, copy=False)
+    a = asarray(a)
     if n is None:
         n = a.shape[axis]
-    output = _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftf, _fft_cache)
-    if _unitary(norm):
-        output *= 1 / sqrt(n)
+    fct = 1
+    if norm is not None and _unitary(norm):
+        fct = 1 / sqrt(n)
+    output = _raw_fft(a, n, axis, False, True, fct)
     return output
 
 
@@ -294,13 +272,14 @@ def ifft(a, n=None, axis=-1, norm=None):
     >>> plt.show()
 
     """
-    # The copy may be required for multithreading.
-    a = array(a, copy=True, dtype=complex)
+    a = asarray(a)
     if n is None:
         n = a.shape[axis]
-    unitary = _unitary(norm)
-    output = _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftb, _fft_cache)
-    return output * (1 / (sqrt(n) if unitary else n))
+    fct = 1/n
+    if norm is not None and _unitary(norm):
+        fct = 1/sqrt(n)
+    output = _raw_fft(a, n, axis, False, False, fct)
+    return output
 
 
 
@@ -383,14 +362,13 @@ def rfft(a, n=None, axis=-1, norm=None):
     exploited to compute only the non-negative frequency terms.
 
     """
-    # The copy may be required for multithreading.
-    a = array(a, copy=True, dtype=float)
-    output = _raw_fft(a, n, axis, fftpack.rffti, fftpack.rfftf,
-                      _real_fft_cache)
-    if _unitary(norm):
+    a = asarray(a)
+    fct = 1
+    if norm is not None and _unitary(norm):
         if n is None:
             n = a.shape[axis]
-        output *= 1 / sqrt(n)
+        fct = 1/sqrt(n)
+    output = _raw_fft(a, n, axis, True, True, fct)
     return output
 
 
@@ -475,14 +453,14 @@ def irfft(a, n=None, axis=-1, norm=None):
     specified, and the output array is purely real.
 
     """
-    # The copy may be required for multithreading.
-    a = array(a, copy=True, dtype=complex)
+    a = asarray(a)
     if n is None:
         n = (a.shape[axis] - 1) * 2
-    unitary = _unitary(norm)
-    output = _raw_fft(a, n, axis, fftpack.rffti, fftpack.rfftb,
-                      _real_fft_cache)
-    return output * (1 / (sqrt(n) if unitary else n))
+    fct = 1/n
+    if norm is not None and _unitary(norm):
+        fct = 1/sqrt(n)
+    output = _raw_fft(a, n, axis, True, False, fct)
+    return output
 
 
 @array_function_dispatch(_fft_dispatcher)
@@ -560,8 +538,7 @@ def hfft(a, n=None, axis=-1, norm=None):
            [ 2., -2.]])
 
     """
-    # The copy may be required for multithreading.
-    a = array(a, copy=True, dtype=complex)
+    a = asarray(a)
     if n is None:
         n = (a.shape[axis] - 1) * 2
     unitary = _unitary(norm)
@@ -621,8 +598,7 @@ def ihfft(a, n=None, axis=-1, norm=None):
     array([ 1.-0.j,  2.-0.j,  3.-0.j,  4.-0.j]) # may vary
 
     """
-    # The copy may be required for multithreading.
-    a = array(a, copy=True, dtype=float)
+    a = asarray(a)
     if n is None:
         n = a.shape[axis]
     unitary = _unitary(norm)
@@ -1122,8 +1098,7 @@ def rfftn(a, s=None, axes=None, norm=None):
             [0.+0.j,  0.+0.j]]])
 
     """
-    # The copy may be required for multithreading.
-    a = array(a, copy=True, dtype=float)
+    a = asarray(a)
     s, axes = _cook_nd_args(a, s, axes)
     a = rfft(a, s[-1], axes[-1], norm)
     for ii in range(len(axes)-1):
@@ -1255,8 +1230,7 @@ def irfftn(a, s=None, axes=None, norm=None):
             [1.,  1.]]])
 
     """
-    # The copy may be required for multithreading.
-    a = array(a, copy=True, dtype=complex)
+    a = asarray(a)
     s, axes = _cook_nd_args(a, s, axes, invreal=1)
     for ii in range(len(axes)-1):
         a = ifft(a, s[ii], axes[ii], norm)
diff --git a/numpy/fft/setup.py b/numpy/fft/setup.py
index cd99a82d7..6c3548b65 100644
--- a/numpy/fft/setup.py
+++ b/numpy/fft/setup.py
@@ -7,9 +7,9 @@ def configuration(parent_package='',top_path=None):
 
     config.add_data_dir('tests')
 
-    # Configure fftpack_lite
-    config.add_extension('fftpack_lite',
-                         sources=['fftpack_litemodule.c', 'fftpack.c']
+    # Configure pocketfft_internal
+    config.add_extension('pocketfft_internal',
+                         sources=['pocketfft.c']
                          )
 
     return config
diff --git a/numpy/fft/tests/test_helper.py b/numpy/fft/tests/test_helper.py
index 8d315fa02..6613c8002 100644
--- a/numpy/fft/tests/test_helper.py
+++ b/numpy/fft/tests/test_helper.py
@@ -7,7 +7,6 @@ from __future__ import division, absolute_import, print_function
 import numpy as np
 from numpy.testing import assert_array_almost_equal, assert_equal
 from numpy import fft, pi
-from numpy.fft.helper import _FFTCache
 
 
 class TestFFTShift(object):
@@ -168,81 +167,3 @@ class TestIRFFTN(object):
 
         # Should not raise error
         fft.irfftn(a, axes=axes)
-
-
-class TestFFTCache(object):
-
-    def test_basic_behaviour(self):
-        c = _FFTCache(max_size_in_mb=1, max_item_count=4)
-
-        # Put
-        c.put_twiddle_factors(1, np.ones(2, dtype=np.float32))
-        c.put_twiddle_factors(2, np.zeros(2, dtype=np.float32))
-
-        # Get
-        assert_array_almost_equal(c.pop_twiddle_factors(1),
-                                  np.ones(2, dtype=np.float32))
-        assert_array_almost_equal(c.pop_twiddle_factors(2),
-                                  np.zeros(2, dtype=np.float32))
-
-        # Nothing should be left.
-        assert_equal(len(c._dict), 0)
-
-        # Now put everything in twice so it can be retrieved once and each will
-        # still have one item left.
-        for _ in range(2):
-            c.put_twiddle_factors(1, np.ones(2, dtype=np.float32))
-            c.put_twiddle_factors(2, np.zeros(2, dtype=np.float32))
-        assert_array_almost_equal(c.pop_twiddle_factors(1),
-                                  np.ones(2, dtype=np.float32))
-        assert_array_almost_equal(c.pop_twiddle_factors(2),
-                                  np.zeros(2, dtype=np.float32))
-        assert_equal(len(c._dict), 2)
-
-    def test_automatic_pruning(self):
-        # That's around 2600 single precision samples.
-        c = _FFTCache(max_size_in_mb=0.01, max_item_count=4)
-
-        c.put_twiddle_factors(1, np.ones(200, dtype=np.float32))
-        c.put_twiddle_factors(2, np.ones(200, dtype=np.float32))
-        assert_equal(list(c._dict.keys()), [1, 2])
-
-        # This is larger than the limit but should still be kept.
-        c.put_twiddle_factors(3, np.ones(3000, dtype=np.float32))
-        assert_equal(list(c._dict.keys()), [1, 2, 3])
-        # Add one more.
-        c.put_twiddle_factors(4, np.ones(3000, dtype=np.float32))
-        # The other three should no longer exist.
-        assert_equal(list(c._dict.keys()), [4])
-
-        # Now test the max item count pruning.
-        c = _FFTCache(max_size_in_mb=0.01, max_item_count=2)
-        c.put_twiddle_factors(2, np.empty(2))
-        c.put_twiddle_factors(1, np.empty(2))
-        # Can still be accessed.
-        assert_equal(list(c._dict.keys()), [2, 1])
-
-        c.put_twiddle_factors(3, np.empty(2))
-        # 1 and 3 can still be accessed - c[2] has been touched least recently
-        # and is thus evicted.
-        assert_equal(list(c._dict.keys()), [1, 3])
-
-        # One last test. We will add a single large item that is slightly
-        # bigger then the cache size. Some small items can still be added.
-        c = _FFTCache(max_size_in_mb=0.01, max_item_count=5)
-        c.put_twiddle_factors(1, np.ones(3000, dtype=np.float32))
-        c.put_twiddle_factors(2, np.ones(2, dtype=np.float32))
-        c.put_twiddle_factors(3, np.ones(2, dtype=np.float32))
-        c.put_twiddle_factors(4, np.ones(2, dtype=np.float32))
-        assert_equal(list(c._dict.keys()), [1, 2, 3, 4])
-
-        # One more big item. This time it is 6 smaller ones but they are
-        # counted as one big item.
-        for _ in range(6):
-            c.put_twiddle_factors(5, np.ones(500, dtype=np.float32))
-        # '1' no longer in the cache. Rest still in the cache.
-        assert_equal(list(c._dict.keys()), [2, 3, 4, 5])
-
-        # Another big item - should now be the only item in the cache.
-        c.put_twiddle_factors(6, np.ones(4000, dtype=np.float32))
-        assert_equal(list(c._dict.keys()), [6])
diff --git a/numpy/fft/tests/test_fftpack.py b/numpy/fft/tests/test_pocketfft.py
index 8d6cd8407..0552f6afd 100644
--- a/numpy/fft/tests/test_fftpack.py
+++ b/numpy/fft/tests/test_pocketfft.py
@@ -28,6 +28,16 @@ class TestFFTShift(object):
 
 class TestFFT1D(object):
 
+    def test_identity(self):
+        maxlen = 512
+        x = random(maxlen) + 1j*random(maxlen)
+        xr = random(maxlen)
+        for i in range(1,maxlen):
+            assert_array_almost_equal(np.fft.ifft(np.fft.fft(x[0:i])), x[0:i],
+                                      decimal=12)
+            assert_array_almost_equal(np.fft.irfft(np.fft.rfft(xr[0:i]),i),
+                                      xr[0:i], decimal=12)
+
     def test_fft(self):
         x = random(30) + 1j*random(30)
         assert_array_almost_equal(fft1(x), np.fft.fft(x))