text: /*
   * 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
   */
  
  /*
   * floating point version is added by LTE Inc. and does not tested well.
   * it may not work well.
   */
  
  #include <math.h>
  #include <string.h>
  #include <stdlib.h>
  #include <stdio.h> // dbg
  
  #ifdef _MSC_VER
  // Suppress warnings on MSVC
  #pragma warning( disable : 4244 )
  #pragma warning( disable : 4305 )
  #pragma warning( disable : 4204 )
  #pragma warning( disable : 4267 )
  #endif
  
  //#include "npy_config.h"
  //#define restrict NPY_RESTRICT
  
  #define 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;
  
  struct cfft_plan_i_f;
  typedef struct cfft_plan_i_f * cfft_plan_f;
  struct rfft_plan_i_f;
  typedef struct rfft_plan_i_f * rfft_plan_f;
  
  // 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;
    }
  
  static void my_sincosm1pi_f (float a, float *restrict res)
    {
    float s = a * a;
    /* Approximate cos(pi*x)-1 for x in [-0.25,0.25] */
    float 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);
    float 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_octant_f(size_t den, float * 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_f((2.*i)/den,&res[2*i]);
    size_t start=l1;
    while(start<n)
      {
      float cs[2];
      my_sincosm1pi_f((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)
        {
        float 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_quadrant_f(size_t n, float * restrict res)
    {
    float * restrict p = res+n;
    calc_first_octant_f(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 calc_first_half_f(size_t n, float * restrict res)
    {
    int ndone=(n+1)>>1;
    float * p = res+n-1;
    calc_first_octant_f(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_quadrant_f(size_t n, float * restrict res)
    {
    const float 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_first_half_f(size_t n, float * 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 fill_second_half_f(size_t n, float * 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_half_f(size_t n, float * restrict res)
    {
    if ((n&3)==0)
      {
      calc_first_octant_f(n, res);
      fill_first_quadrant_f(n, res);
      fill_first_half_f(n, res);
      }
    else if ((n&1)==0)
      {
      calc_first_quadrant_f(n, res);
      fill_first_half_f(n, res);
      }
    else
      calc_first_half_f(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 void sincos_2pibyn_f(size_t n, float * restrict res)
    {
    sincos_2pibyn_half_f(n, res);
    fill_second_half_f(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;
    }
  
  NOINLINE static float cost_guess_f (size_t n)
    {
    const float lfp=1.1; // penalty for non-hardcoded larger factors
    size_t ni=n;
    float result=0.;
    size_t tmp;
    while (((tmp=(n>>1))<<1)==n)
      { result+=2; n=tmp; }
  
    size_t limit=(size_t)sqrtf(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;
  
  typedef struct cmplx_f {
    float r,i;
  } cmplx_f;
  
  #define NFCT 25
  typedef struct cfftp_fctdata
    {
    size_t fct;
    cmplx *tw, *tws;
    } cfftp_fctdata;
  
  typedef struct cfftp_fctdata_f
    {
    size_t fct;
    cmplx_f *tw, *tws;
    } cfftp_fctdata_f;
  
  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;
  
  typedef struct cfftp_plan_i_f
    {
    size_t length, nfct;
    cmplx_f *mem;
    cfftp_fctdata_f fct[NFCT];
    } cfftp_plan_i_f;
  typedef struct cfftp_plan_i_f * cfftp_plan_f;
  
  #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 pass2b_f (size_t ido, size_t l1, const cmplx_f * restrict cc,
    cmplx_f * restrict ch, const cmplx_f * 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_f 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)
          }
        }
    }
  
  NOINLINE static void pass2f_f (size_t ido, size_t l1, const cmplx_f * restrict cc,
    cmplx_f * restrict ch, const cmplx_f * 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_f 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 PREP3_f(idx) \
          cmplx_f 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 PARTSTEP3a_f(u1,u2,twr,twi) \
          { \
          cmplx_f 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) \
          }
  #define PARTSTEP3b_f(u1,u2,twr,twi) \
          { \
          cmplx_f 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)
          }
        }
    }
  
  NOINLINE static void pass3b_f (size_t ido, size_t l1, const cmplx_f * restrict cc,
    cmplx_f * restrict ch, const cmplx_f * restrict wa)
    {
    const size_t cdim=3;
    const float tw1r=-0.5, tw1i= 0.86602540378443864676;
  
    if (ido==1)
      for (size_t k=0; k<l1; ++k)
        {
        PREP3_f(0)
        PARTSTEP3a_f(1,2,tw1r,tw1i)
        }
    else
      for (size_t k=0; k<l1; ++k)
        {
        {
        PREP3_f(0)
        PARTSTEP3a_f(1,2,tw1r,tw1i)
        }
        for (size_t i=1; i<ido; ++i)
          {
          PREP3_f(i)
          PARTSTEP3b_f(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)
          }
        }
    }
  #define PARTSTEP3f_f(u1,u2,twr,twi) \
          { \
          cmplx_f 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_f (size_t ido, size_t l1, const cmplx_f * restrict cc,
    cmplx_f * restrict ch, const cmplx_f * restrict wa)
    {
    const size_t cdim=3;
    const float tw1r=-0.5, tw1i= -0.86602540378443864676;
  
    if (ido==1)
      for (size_t k=0; k<l1; ++k)
        {
        PREP3_f(0)
        PARTSTEP3a_f(1,2,tw1r,tw1i)
        }
    else
      for (size_t k=0; k<l1; ++k)
        {
        {
        PREP3_f(0)
        PARTSTEP3a_f(1,2,tw1r,tw1i)
        }
        for (size_t i=1; i<ido; ++i)
          {
          PREP3_f(i)
          PARTSTEP3f_f(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 pass4b_f (size_t ido, size_t l1, const cmplx_f * restrict cc,
    cmplx_f * restrict ch, const cmplx_f * restrict wa)
    {
    const size_t cdim=4;
  
    if (ido==1)
      for (size_t k=0; k<l1; ++k)
        {
        cmplx_f 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_f 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_f c2, c3, c4, t1, t2, t3, t4;
          cmplx_f 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_f 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)
          }
        }
    }
  
  NOINLINE static void pass4f_f (size_t ido, size_t l1, const cmplx_f * restrict cc,
    cmplx_f * restrict ch, const cmplx_f * restrict wa)
    {
    const size_t cdim=4;
  
    if (ido==1)
      for (size_t k=0; k<l1; ++k)
        {
        cmplx_f 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_f 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_f c2, c3, c4, t1, t2, t3, t4;
          cmplx_f 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_f 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 PREP5_f(idx) \
          cmplx_f 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_f(u1,u2,twar,twbr,twai,twbi) \
          { \
          cmplx_f 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_f(u1,u2,twar,twbr,twai,twbi) \
          { \
          cmplx_f 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_f (size_t ido, size_t l1, const cmplx_f * restrict cc,
    cmplx_f * restrict ch, const cmplx_f * restrict wa)
    {
    const size_t cdim=5;
    const float tw1r= 0.3090169943749474241,
                 tw1i= 0.95105651629515357212,
                 tw2r= -0.8090169943749474241,
                 tw2i= 0.58778525229247312917;
  
    if (ido==1)
      for (size_t k=0; k<l1; ++k)
        {
        PREP5_f(0)
        PARTSTEP5a_f(1,4,tw1r,tw2r,+tw1i,+tw2i)
        PARTSTEP5a_f(2,3,tw2r,tw1r,+tw2i,-tw1i)
        }
    else
      for (size_t k=0; k<l1; ++k)
        {
        {
        PREP5_f(0)
        PARTSTEP5a_f(1,4,tw1r,tw2r,+tw1i,+tw2i)
        PARTSTEP5a_f(2,3,tw2r,tw1r,+tw2i,-tw1i)
        }
        for (size_t i=1; i<ido; ++i)
          {
          PREP5_f(i)
          PARTSTEP5b_f(1,4,tw1r,tw2r,+tw1i,+tw2i)
          PARTSTEP5b_f(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 PARTSTEP5f_f(u1,u2,twar,twbr,twai,twbi) \
          { \
          cmplx_f 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_f (size_t ido, size_t l1, const cmplx_f * restrict cc,
    cmplx_f * restrict ch, const cmplx_f * restrict wa)
    {
    const size_t cdim=5;
    const float tw1r= 0.3090169943749474241,
                 tw1i= -0.95105651629515357212,
                 tw2r= -0.8090169943749474241,
                 tw2i= -0.58778525229247312917;
  
    if (ido==1)
      for (size_t k=0; k<l1; ++k)
        {
        PREP5_f(0)
        PARTSTEP5a_f(1,4,tw1r,tw2r,+tw1i,+tw2i)
        PARTSTEP5a_f(2,3,tw2r,tw1r,+tw2i,-tw1i)
        }
    else
      for (size_t k=0; k<l1; ++k)
        {
        {
        PREP5_f(0)
        PARTSTEP5a_f(1,4,tw1r,tw2r,+tw1i,+tw2i)
        PARTSTEP5a_f(2,3,tw2r,tw1r,+tw2i,-tw1i)
        }
        for (size_t i=1; i<ido; ++i)
          {
          PREP5_f(i)
          PARTSTEP5f_f(1,4,tw1r,tw2r,+tw1i,+tw2i)
          PARTSTEP5f_f(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 PREP7_f(idx) \
          cmplx_f 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_f(u1,u2,x1,x2,x3,y1,y2,y3,out1,out2) \
          { \
          cmplx_f 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_f(u1,u2,x1,x2,x3,y1,y2,y3) \
          PARTSTEP7a0_f(u1,u2,x1,x2,x3,y1,y2,y3,CH(0,k,u1),CH(0,k,u2))
  #define PARTSTEP7_f(u1,u2,x1,x2,x3,y1,y2,y3) \
          { \
          cmplx_f da,db; \
          PARTSTEP7a0_f(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_f(size_t ido, size_t l1, const cmplx_f * restrict cc,
    cmplx_f * restrict ch, const cmplx_f * restrict wa, const int sign)
    {
    const size_t cdim=7;
    const float 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_f(0)
        PARTSTEP7a_f(1,6,tw1r,tw2r,tw3r,+tw1i,+tw2i,+tw3i)
        PARTSTEP7a_f(2,5,tw2r,tw3r,tw1r,+tw2i,-tw3i,-tw1i)
        PARTSTEP7a_f(3,4,tw3r,tw1r,tw2r,+tw3i,-tw1i,+tw2i)
        }
    else
      for (size_t k=0; k<l1; ++k)
        {
        {
        PREP7_f(0)
        PARTSTEP7a_f(1,6,tw1r,tw2r,tw3r,+tw1i,+tw2i,+tw3i)
        PARTSTEP7a_f(2,5,tw2r,tw3r,tw1r,+tw2i,-tw3i,-tw1i)
        PARTSTEP7a_f(3,4,tw3r,tw1r,tw2r,+tw3i,-tw1i,+tw2i)
        }
        for (size_t i=1; i<ido; ++i)
          {
          PREP7_f(i)
          PARTSTEP7_f(1,6,tw1r,tw2r,tw3r,+tw1i,+tw2i,+tw3i)
          PARTSTEP7_f(2,5,tw2r,tw3r,tw1r,+tw2i,-tw3i,-tw1i)
          PARTSTEP7_f(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 PREP11_f(idx) \
          cmplx_f 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_f(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,out1,out2) \
          { \
          cmplx_f 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_f(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5) \
          PARTSTEP11a0_f(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,CH(0,k,u1),CH(0,k,u2))
  #define PARTSTEP11_f(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5) \
          { \
          cmplx_f da,db; \
          PARTSTEP11a0_f(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_f (size_t ido, size_t l1, const cmplx_f * restrict cc,
    cmplx_f * restrict ch, const cmplx_f * restrict wa, const int sign)
    {
    const size_t cdim=11;
    const float 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_f(0)
        PARTSTEP11a_f(1,10,tw1r,tw2r,tw3r,tw4r,tw5r,+tw1i,+tw2i,+tw3i,+tw4i,+tw5i)
        PARTSTEP11a_f(2, 9,tw2r,tw4r,tw5r,tw3r,tw1r,+tw2i,+tw4i,-tw5i,-tw3i,-tw1i)
        PARTSTEP11a_f(3, 8,tw3r,tw5r,tw2r,tw1r,tw4r,+tw3i,-tw5i,-tw2i,+tw1i,+tw4i)
        PARTSTEP11a_f(4, 7,tw4r,tw3r,tw1r,tw5r,tw2r,+tw4i,-tw3i,+tw1i,+tw5i,-tw2i)
        PARTSTEP11a_f(5, 6,tw5r,tw1r,tw4r,tw2r,tw3r,+tw5i,-tw1i,+tw4i,-tw2i,+tw3i)
        }
    else
      for (size_t k=0; k<l1; ++k)
        {
        {
        PREP11_f(0)
        PARTSTEP11a_f(1,10,tw1r,tw2r,tw3r,tw4r,tw5r,+tw1i,+tw2i,+tw3i,+tw4i,+tw5i)
        PARTSTEP11a_f(2, 9,tw2r,tw4r,tw5r,tw3r,tw1r,+tw2i,+tw4i,-tw5i,-tw3i,-tw1i)
        PARTSTEP11a_f(3, 8,tw3r,tw5r,tw2r,tw1r,tw4r,+tw3i,-tw5i,-tw2i,+tw1i,+tw4i)
        PARTSTEP11a_f(4, 7,tw4r,tw3r,tw1r,tw5r,tw2r,+tw4i,-tw3i,+tw1i,+tw5i,-tw2i)
        PARTSTEP11a_f(5, 6,tw5r,tw1r,tw4r,tw2r,tw3r,+tw5i,-tw1i,+tw4i,-tw2i,+tw3i)
        }
        for (size_t i=1; i<ido; ++i)
          {
          PREP11_f(i)
          PARTSTEP11_f(1,10,tw1r,tw2r,tw3r,tw4r,tw5r,+tw1i,+tw2i,+tw3i,+tw4i,+tw5i)
          PARTSTEP11_f(2, 9,tw2r,tw4r,tw5r,tw3r,tw1r,+tw2i,+tw4i,-tw5i,-tw3i,-tw1i)
          PARTSTEP11_f(3, 8,tw3r,tw5r,tw2r,tw1r,tw4r,+tw3i,-tw5i,-tw2i,+tw1i,+tw4i)
          PARTSTEP11_f(4, 7,tw4r,tw3r,tw1r,tw5r,tw2r,+tw4i,-tw3i,+tw1i,+tw5i,-tw2i)
          PARTSTEP11_f(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;
    }
  
  NOINLINE static int passg_f (size_t ido, size_t ip, size_t l1,
    cmplx_f * restrict cc, cmplx_f * restrict ch, const cmplx_f * restrict wa,
    const cmplx_f * restrict csarr, const int sign)
    {
    const size_t cdim=ip;
    size_t ipph = (ip+1)/2;
    size_t idl1 = ido*l1;
  
    cmplx_f * restrict wal=RALLOC(cmplx_f,ip);
    if (!wal) return -1;
    wal[0]=(cmplx_f){1.,0.};
    for (size_t i=1; i<ip; ++i)
      wal[i]=(cmplx_f){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_f 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_f xwal=wal[iwal];
        iwal+=l; if (iwal>ip) iwal-=ip;
        cmplx_f 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_f 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_f 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_f 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_f 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;
    }
  
  NOINLINE WARN_UNUSED_RESULT static int pass_all_f(cfftp_plan_f plan, cmplx_f c[], float fct,
    const int sign)
    {
    if (plan->length==1) return 0;
    size_t len=plan->length;
    size_t l1=1, nf=plan->nfct;
    cmplx_f *ch = RALLOC(cmplx_f, len);
    if (!ch) return -1;
    cmplx_f *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_f (ido, l1, p1, p2, plan->fct[k1].tw)
               : pass4f_f (ido, l1, p1, p2, plan->fct[k1].tw);
      else if(ip==2)
        sign>0 ? pass2b_f (ido, l1, p1, p2, plan->fct[k1].tw)
               : pass2f_f (ido, l1, p1, p2, plan->fct[k1].tw);
      else if(ip==3)
        sign>0 ? pass3b_f (ido, l1, p1, p2, plan->fct[k1].tw)
               : pass3f_f (ido, l1, p1, p2, plan->fct[k1].tw);
      else if(ip==5)
        sign>0 ? pass5b_f (ido, l1, p1, p2, plan->fct[k1].tw)
               : pass5f_f (ido, l1, p1, p2, plan->fct[k1].tw);
      else if(ip==7)  pass7_f (ido, l1, p1, p2, plan->fct[k1].tw, sign);
      else if(ip==11) pass11_f(ido, l1, p1, p2, plan->fct[k1].tw, sign);
      else
        {
        if (passg_f(ido, ip, l1, p1, p2, plan->fct[k1].tw, plan->fct[k1].tws, sign))
          { DEALLOC(ch); return -1; }
        SWAP(p1,p2,cmplx_f *);
        }
      SWAP(p1,p2,cmplx_f *);
      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_f));
      }
    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_forward_f(cfftp_plan_f plan, float c[], float fct)
    { return pass_all_f(plan,(cmplx_f *)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_backward_f(cfftp_plan_f plan, float c[], float fct)
    { return pass_all_f(plan,(cmplx_f *)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 WARN_UNUSED_RESULT
  static int cfftp_factorize_f (cfftp_plan_f 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)(sqrtf((float)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)(sqrtf((float)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 static size_t cfftp_twsize_f (cfftp_plan_f 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;
    }
  
  NOINLINE WARN_UNUSED_RESULT static int cfftp_comp_twiddle_f (cfftp_plan_f plan)
    {
    size_t length=plan->length;
    float *twid = RALLOC(float, 2*length);
    if (!twid) return -1;
    sincos_2pibyn_f(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 cfftp_plan_f make_cfftp_plan_f (size_t length)
    {
    if (length==0) return NULL;
    cfftp_plan_f plan = RALLOC(cfftp_plan_i_f,1);
    if (!plan) return NULL;
    plan->length=length;
    plan->nfct=0;
    for (size_t i=0; i<NFCT; ++i)
      plan->fct[i]=(cfftp_fctdata_f){0,0,0};
    plan->mem=0;
    if (length==1) return plan;
    if (cfftp_factorize_f(plan)!=0) { DEALLOC(plan); return NULL; }
    size_t tws=cfftp_twsize_f(plan);
    plan->mem=RALLOC(cmplx_f,tws);
    if (!plan->mem) { DEALLOC(plan); return NULL; }
    if (cfftp_comp_twiddle_f(plan)!=0)
      { DEALLOC(plan->mem); DEALLOC(plan); return NULL; }
    return plan;
    }
  
  static void destroy_cfftp_plan (cfftp_plan plan)
    {
    DEALLOC(plan->mem);
    DEALLOC(plan);
    }
  
  static void destroy_cfftp_plan_f (cfftp_plan_f plan)
    {
    DEALLOC(plan->mem);
    DEALLOC(plan);
    }
  
  typedef struct rfftp_fctdata
    {
    size_t fct;
    double *tw, *tws;
    } rfftp_fctdata;
  
  typedef struct rfftp_fctdata_f
    {
    size_t fct;
    float *tw, *tws;
    } rfftp_fctdata_f;
  
  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;
  
  typedef struct rfftp_plan_i_f
    {
    size_t length, nfct;
    float *mem;
    rfftp_fctdata_f fct[NFCT];
    } rfftp_plan_i_f;
  typedef struct rfftp_plan_i_f * rfftp_plan_f;
  
  #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 radf2_f (size_t ido, size_t l1, const float * restrict cc,
    float * restrict ch, const float * 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;
        float 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 radf3_f(size_t ido, size_t l1, const float * restrict cc,
    float * restrict ch, const float * restrict wa)
    {
    const size_t cdim=3;
    static const float taur=-0.5, taui=0.86602540378443864676;
  
    for (size_t k=0; k<l1; k++)
      {
      float 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;
        float 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
        float cr2=dr2+dr3; // c add
        float 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;
        float tr2 = CC(i-1,k,0)+taur*cr2; // c add
        float ti2 = CC(i  ,k,0)+taur*ci2;
        float tr3 = taui*(di2-di3);  // t3 = taui*i*(d3-d2)?
        float 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 radf4_f(size_t ido, size_t l1, const float * restrict cc,
    float * restrict ch, const float * restrict wa)
    {
    const size_t cdim=4;
    static const float hsqt2=0.70710678118654752440f;
  
    for (size_t k=0; k<l1; k++)
      {
      float 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++)
        {
        float ti1=-hsqt2*(CC(ido-1,k,1)+CC(ido-1,k,3));
        float 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;
        float 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)
        }
    }
  
  NOINLINE static void radf5_f(size_t ido, size_t l1, const float * restrict cc,
    float * restrict ch, const float * restrict wa)
    {
    const size_t cdim=5;
    static const float tr11= 0.3090169943749474241, ti11=0.95105651629515357212,
                        tr12=-0.8090169943749474241, ti12=0.58778525229247312917;
  
    for (size_t k=0; k<l1; k++)
      {
      float 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)
        {
        float 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 );
          }
      }
    }
  
  NOINLINE static void radfg_f(size_t ido, size_t ip, size_t l1,
    float * restrict cc, float * restrict ch, const float * restrict wa,
    const float * 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
            {
            float t1=C1(i,k,j ), t2=C1(i+1,k,j ),
                   t3=C1(i,k,jc), t4=C1(i+1,k,jc);
            float 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
        {
        float 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;
        float ar1=csarr[2*iang], ai1=csarr[2*iang+1];
        iang+=l; if (iang>=ip) iang-=ip;
        float ar2=csarr[2*iang], ai2=csarr[2*iang+1];
        iang+=l; if (iang>=ip) iang-=ip;
        float ar3=csarr[2*iang], ai3=csarr[2*iang+1];
        iang+=l; if (iang>=ip) iang-=ip;
        float 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;
        float ar1=csarr[2*iang], ai1=csarr[2*iang+1];
        iang+=l; if (iang>=ip) iang-=ip;
        float 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;
        float 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;
    }
  
  static void copy_and_norm_f(float *c, float *p1, size_t n, double fct)
    {
    if (p1!=c)
      {
      if (fct!=1.f)
        for (size_t i=0; i<n; ++i)
          c[i] = fct*p1[i];
      else
        memcpy (c,p1,n*sizeof(float));
      }
    else
      if (fct!=1.f)
        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_forward_f(rfftp_plan_f plan, float c[], float fct)
    {
    if (plan->length==1) return 0;
    size_t n=plan->length;
    size_t l1=n, nf=plan->nfct;
    float *ch = RALLOC(float, n);
    if (!ch) return -1;
    float *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_f(ido, l1, p1, p2, plan->fct[k].tw);
      else if(ip==2)
        radf2_f(ido, l1, p1, p2, plan->fct[k].tw);
      else if(ip==3)
        radf3_f(ido, l1, p1, p2, plan->fct[k].tw);
      else if(ip==5)
        radf5_f(ido, l1, p1, p2, plan->fct[k].tw);
      else
        {
        radfg_f(ido, ip, l1, p1, p2, plan->fct[k].tw, plan->fct[k].tws);
        SWAP (p1,p2,float *);
        }
      SWAP (p1,p2,float *);
      }
    copy_and_norm_f(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;
    }
  
  WARN_UNUSED_RESULT
  static int rfftp_factorize_f (rfftp_plan_f 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;
    }
  
  static size_t rfftp_twsize_f(rfftp_plan_f 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;
    }
  
  
  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;
    }
  
  WARN_UNUSED_RESULT NOINLINE static int rfftp_comp_twiddle_f (rfftp_plan_f plan)
    {
    size_t length=plan->length;
    float *twid = RALLOC(float, 2*length);
    if (!twid) return -1;
    sincos_2pibyn_half_f(length, twid);
    size_t l1=1;
    float *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 rfftp_plan_f make_rfftp_plan_f (size_t length)
    {
    if (length==0) return NULL;
    rfftp_plan_f plan = RALLOC(rfftp_plan_i_f,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_f){0,0,0};
    if (length==1) return plan;
    if (rfftp_factorize_f(plan)!=0) { DEALLOC(plan); return NULL; }
    size_t tws=rfftp_twsize_f(plan);
    plan->mem=RALLOC(float,tws);
    if (!plan->mem) { DEALLOC(plan); return NULL; }
    if (rfftp_comp_twiddle_f(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);
    }
  
  NOINLINE static void destroy_rfftp_plan_f (rfftp_plan_f 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;
  
  typedef struct fftblue_plan_i_f
    {
    size_t n, n2;
    cfftp_plan_f plan;
    float *mem;
    float *bk, *bkf;
    } fftblue_plan_i_f;
  typedef struct fftblue_plan_i_f * fftblue_plan_f;
  
  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 fftblue_plan_f make_fftblue_plan_f (size_t length)
    {
    fftblue_plan_f plan = RALLOC(fftblue_plan_i_f,1);
    if (!plan) return NULL;
    plan->n = length;
    plan->n2 = good_size(plan->n*2-1);
    plan->mem = RALLOC(float, 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 */
    float *tmp = RALLOC(float,4*plan->n);
    if (!tmp) { DEALLOC(plan->mem); DEALLOC(plan); return NULL; }
    sincos_2pibyn_f(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. */
    float 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_f(plan->n2);
    if (!plan->plan)
      { DEALLOC(tmp); DEALLOC(plan->mem); DEALLOC(plan); return NULL; }
    if (cfftp_forward_f(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 static void destroy_fftblue_plan_f (fftblue_plan_f plan)
    {
    DEALLOC(plan->mem);
    destroy_cfftp_plan_f(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;
    }
  
  NOINLINE WARN_UNUSED_RESULT
  static int fftblue_fft_f(fftblue_plan_f plan, float c[], int isign, float fct)
    {
    size_t n=plan->n;
    size_t n2=plan->n2;
    float *bk  = plan->bk;
    float *bkf = plan->bkf;
    float *akf = RALLOC(float, 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_f (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)
        {
        float 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)
        {
        float 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_f (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_backward_f(fftblue_plan_f plan, float c[], float fct)
    { return fftblue_fft_f(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 cfftblue_forward_f(fftblue_plan_f plan, float c[], float fct)
    { return fftblue_fft_f(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_backward_f(fftblue_plan_f plan, float c[], float fct)
    {
    size_t n=plan->n;
    float *tmp = RALLOC(float,2*n);
    if (!tmp) return -1;
    tmp[0]=c[0];
    tmp[1]=0.;
    memcpy (tmp+2,c+1, (n-1)*sizeof(float));
    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_f(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;
    }
  
  WARN_UNUSED_RESULT
  static int rfftblue_forward_f(fftblue_plan_f plan, float c[], float fct)
    {
    size_t n=plan->n;
    float *tmp = RALLOC(float,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_f(plan,tmp,-1,fct)!=0)
      { DEALLOC(tmp); return -1; }
    c[0] = tmp[0];
    memcpy (c+1, tmp+2, (n-1)*sizeof(float));
    DEALLOC(tmp);
    return 0;
    }
  
  typedef struct cfft_plan_i
    {
    cfftp_plan packplan;
    fftblue_plan blueplan;
    } cfft_plan_i;
  
  typedef struct cfft_plan_i_f
    {
    cfftp_plan_f packplan;
    fftblue_plan_f blueplan;
    } cfft_plan_i_f;
  
  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 cfft_plan_f make_cfft_plan_f (size_t length)
    {
    if (length==0) return NULL;
    cfft_plan_f plan = RALLOC(cfft_plan_i_f,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_f(length);
      if (!plan->packplan) { DEALLOC(plan); return NULL; }
      return plan;
      }
    float comp1 = cost_guess(length);
    float 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_f(length);
      if (!plan->blueplan) { DEALLOC(plan); return NULL; }
      }
    else
      {
      plan->packplan=make_cfftp_plan_f(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_backward_f(cfft_plan_f plan, float c[], float fct)
    {
    if (plan->packplan)
      return cfftp_backward_f(plan->packplan,c,fct);
    // if (plan->blueplan)
    return cfftblue_backward_f(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);
    }
  
  WARN_UNUSED_RESULT static int cfft_forward_f(cfft_plan_f plan, float c[], float fct)
    {
    if (plan->packplan)
      return cfftp_forward_f(plan->packplan,c,fct);
    // if (plan->blueplan)
    return cfftblue_forward_f(plan->blueplan,c,fct);
    }
  
  typedef struct rfft_plan_i
    {
    rfftp_plan packplan;
    fftblue_plan blueplan;
    } rfft_plan_i;
  
  typedef struct rfft_plan_i_f
    {
    rfftp_plan_f packplan;
    fftblue_plan_f blueplan;
    } rfft_plan_i_f;
  
  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 rfft_plan_f make_rfft_plan_f (size_t length)
    {
    if (length==0) return NULL;
    rfft_plan_f plan = RALLOC(rfft_plan_i_f,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_f(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_f(length);
      if (!plan->blueplan) { DEALLOC(plan); return NULL; }
      }
    else
      {
      plan->packplan=make_rfftp_plan_f(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);
    }
  
  static void destroy_rfft_plan_f (rfft_plan_f plan)
    {
    if (plan->blueplan)
      destroy_fftblue_plan_f(plan->blueplan);
    if (plan->packplan)
      destroy_rfftp_plan_f(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);
    }
  
  WARN_UNUSED_RESULT static int rfft_forward_f(rfft_plan_f plan, float c[], float fct)
    {
    if (plan->packplan)
      return rfftp_forward_f(plan->packplan,c,fct);
    else // if (plan->blueplan)
      return rfftblue_forward_f(plan->blueplan,c,fct);
    }
  
  
  #if 0
  #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_FromAny(a1,
              PyArray_DescrFromType(NPY_CDOUBLE), 1, 0,
              NPY_ARRAY_ENSURECOPY | NPY_ARRAY_DEFAULT |
              NPY_ARRAY_ENSUREARRAY | NPY_ARRAY_FORCECAST,
              NULL);
      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_FromAny(a1,
              PyArray_DescrFromType(NPY_DOUBLE), 1, 0,
              NPY_ARRAY_DEFAULT | NPY_ARRAY_ENSUREARRAY | NPY_ARRAY_FORCECAST,
              NULL);
      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_FromAny(a1,
              PyArray_DescrFromType(NPY_CDOUBLE), 1, 0,
              NPY_ARRAY_DEFAULT | NPY_ARRAY_ENSUREARRAY | NPY_ARRAY_FORCECAST,
              NULL);
      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);
  }
  #endif
  
  // C API
  int rfft_forward_1d_array(const double *input, const int fft_len, const int nsamples, const int nrows, const float norm_factor, double *output)
  {
      rfft_plan plan=NULL;
      int fail = 0;
  
      if ((fft_len < 1) || (nsamples < 1) || (nrows < 1) || (input == NULL) || (output == NULL)) {
        // Invalid parameter.
        return -1;
      }
  
      //PyArrayObject *data = (PyArrayObject *)PyArray_FromAny(a1,
      //        PyArray_DescrFromType(NPY_DOUBLE), 1, 0,
      //        NPY_ARRAY_DEFAULT | NPY_ARRAY_ENSUREARRAY | NPY_ARRAY_FORCECAST,
      //        NULL);
      //if (!data) return NULL;
  
      int npts = fft_len;
  
      //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;
  
      int num_output_complex_points = npts / 2 + 1;
  
      int rstep = num_output_complex_points * 2;	// output stride
  
      int nrepeats = nrows;
      double *rptr = output;
      const double *dptr = input;
  
      plan = make_rfft_plan(npts);
      if (!plan) fail=1;
  
      double fct = norm_factor;
  
      //printf("npts = %d, norm = %f\n", npts, norm_factor);
      //printf("rstep = %d\n", rstep);
      //printf("nrepeats = %d\n", nrepeats);
  
      int nprocessed = 0;
      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;
            dptr += nsamples;
            nprocessed++;
      }
      if (plan) destroy_rfft_plan(plan);
  
      if (fail) {
        return nprocessed;
      }
  
      return nrows;
  }
  
  int cfft_backward_1d_array(const double *input, const int ncolumns, const int nrows, const float norm_factor, double *output)
  {
      if (!input) {
          return -1;
      }
  
      if (!output) {
          return -1;
      }
  
      int npts = ncolumns;
      cfft_plan plan=NULL;
  
      int nrepeats = nrows;
      const double *dptr = input;
      int fail=0;
      double fct = norm_factor;
      //Py_BEGIN_ALLOW_THREADS;
      //NPY_SIGINT_ON;
      plan = make_cfft_plan(npts);
      if (!plan) fail=1;
  
      double *out_ptr = output;
  
      if (!fail)
        for (int i = 0; i < nrepeats; i++) {
            printf("i = %d(%d), ncols = %d, nrows = %d\n", i, nrepeats, ncolumns, nrows);
            memcpy(out_ptr, dptr, sizeof(double) * 2 * npts);
            int res = cfft_backward(plan, out_ptr, fct); // in-place
            if (res!=0) { fail=1; break; }
            dptr += npts*2;
            out_ptr += npts*2;
        }
      if (plan) destroy_cfft_plan(plan);
      //NPY_SIGINT_OFF;
      //Py_END_ALLOW_THREADS;
      if (fail) {
        //Py_XDECREF(data);
        //return PyErr_NoMemory();
        return -1;
      }
      //return (PyObject *)data;
  
      return nrepeats;
  }
  