<gh_stars>0
# is it a good idea to put it on a top-level here?.. maybe some singleton would be nicer or something...
from cffi import FFI
ffi = FFI()


ffi.cdef('''
//typedef ... Treal; // damn, stuff like this doesn't work with cffi yet... structs are fine, but not the basic types :(
// so, we stick to 'double' here:
void cffti(int n, double wsave[]);
void cfftf(int n, double c[], double wsave[]);
''')


C = ffi.verify('''
    /*
    fftpack.c : A set of FFT routines in C.
    Algorithmically based on Fortran-77 FFTPACK by <NAME> (Version 4, 1985).

    */

    /* isign is +1 for backward and -1 for forward transforms */

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

    #ifdef DOUBLE
    #define Treal double
    #else
    #define Treal float
    #endif


    #define ref(u,a) u[a]

    #define MAXFAC 13    /* maximum number of factors in factorization of n */
    #define NSPECIAL 4   /* number of factors for which we have special-case routines */

    #ifdef __cplusplus
    extern "C" {
    #endif


    /* ----------------------------------------------------------------------
    passf2, passf3, passf4, passf5, passf. Complex FFT passes fwd and bwd.
    ---------------------------------------------------------------------- */

    static void passf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[], int isign)
    /* isign==+1 for backward transform */
    {
        int i, k, ah, ac;
        Treal ti2, tr2;
        if (ido <= 2) {
        for (k=0; k<l1; k++) {
            ah = k*ido;
            ac = 2*k*ido;
            ch[ah]              = ref(cc,ac) + ref(cc,ac + ido);
            ch[ah + ido*l1]     = ref(cc,ac) - ref(cc,ac + ido);
            ch[ah+1]            = ref(cc,ac+1) + ref(cc,ac + ido + 1);
            ch[ah + ido*l1 + 1] = ref(cc,ac+1) - ref(cc,ac + ido + 1);
        }
        } else {
        for (k=0; k<l1; k++) {
            for (i=0; i<ido-1; i+=2) {
            ah = i + k*ido;
            ac = i + 2*k*ido;
            ch[ah]   = ref(cc,ac) + ref(cc,ac + ido);
            tr2      = ref(cc,ac) - ref(cc,ac + ido);
            ch[ah+1] = ref(cc,ac+1) + ref(cc,ac + 1 + ido);
            ti2      = ref(cc,ac+1) - ref(cc,ac + 1 + ido);
            ch[ah+l1*ido+1] = wa1[i]*ti2 + isign*wa1[i+1]*tr2;
            ch[ah+l1*ido]   = wa1[i]*tr2 - isign*wa1[i+1]*ti2;
            }
        }
        }
    } /* passf2 */


    static void passf3(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[], int isign)
    /* isign==+1 for backward transform */
    {
        static const Treal taur = -0.5;
        static const Treal taui = 0.866025403784439;
        int i, k, ac, ah;
        Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
        if (ido == 2) {
        for (k=1; k<=l1; k++) {
            ac = (3*k - 2)*ido;
            tr2 = ref(cc,ac) + ref(cc,ac + ido);
            cr2 = ref(cc,ac - ido) + taur*tr2;
            ah = (k - 1)*ido;
            ch[ah] = ref(cc,ac - ido) + tr2;

            ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);
            ci2 = ref(cc,ac - ido + 1) + taur*ti2;
            ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;

            cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));
            ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));
            ch[ah + l1*ido] = cr2 - ci3;
            ch[ah + 2*l1*ido] = cr2 + ci3;
            ch[ah + l1*ido + 1] = ci2 + cr3;
            ch[ah + 2*l1*ido + 1] = ci2 - cr3;
        }
        } else {
        for (k=1; k<=l1; k++) {
            for (i=0; i<ido-1; i+=2) {
            ac = i + (3*k - 2)*ido;
            tr2 = ref(cc,ac) + ref(cc,ac + ido);
            cr2 = ref(cc,ac - ido) + taur*tr2;
            ah = i + (k-1)*ido;
            ch[ah] = ref(cc,ac - ido) + tr2;
            ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);
            ci2 = ref(cc,ac - ido + 1) + taur*ti2;
            ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;
            cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));
            ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));
            dr2 = cr2 - ci3;
            dr3 = cr2 + ci3;
            di2 = ci2 + cr3;
            di3 = ci2 - cr3;
            ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
            ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
            ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
            ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
            }
        }
        }
    } /* passf3 */


    static void passf4(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[], const Treal wa3[], int isign)
    /* isign == -1 for forward transform and +1 for backward transform */
    {
        int i, k, ac, ah;
        Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
        if (ido == 2) {
        for (k=0; k<l1; k++) {
            ac = 4*k*ido + 1;
            ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);
            ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);
            tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);
            ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);
            tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);
            tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);
            ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);
            tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);
            ah = k*ido;
            ch[ah] = tr2 + tr3;
            ch[ah + 2*l1*ido] = tr2 - tr3;
            ch[ah + 1] = ti2 + ti3;
            ch[ah + 2*l1*ido + 1] = ti2 - ti3;
            ch[ah + l1*ido] = tr1 + isign*tr4;
            ch[ah + 3*l1*ido] = tr1 - isign*tr4;
            ch[ah + l1*ido + 1] = ti1 + isign*ti4;
            ch[ah + 3*l1*ido + 1] = ti1 - isign*ti4;
        }
        } else {
        for (k=0; k<l1; k++) {
            for (i=0; i<ido-1; i+=2) {
            ac = i + 1 + 4*k*ido;
            ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);
            ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);
            ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);
            tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);
            tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);
            tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);
            ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);
            tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);
            ah = i + k*ido;
            ch[ah] = tr2 + tr3;
            cr3 = tr2 - tr3;
            ch[ah + 1] = ti2 + ti3;
            ci3 = ti2 - ti3;
            cr2 = tr1 + isign*tr4;
            cr4 = tr1 - isign*tr4;
            ci2 = ti1 + isign*ti4;
            ci4 = ti1 - isign*ti4;
            ch[ah + l1*ido] = wa1[i]*cr2 - isign*wa1[i + 1]*ci2;
            ch[ah + l1*ido + 1] = wa1[i]*ci2 + isign*wa1[i + 1]*cr2;
            ch[ah + 2*l1*ido] = wa2[i]*cr3 - isign*wa2[i + 1]*ci3;
            ch[ah + 2*l1*ido + 1] = wa2[i]*ci3 + isign*wa2[i + 1]*cr3;
            ch[ah + 3*l1*ido] = wa3[i]*cr4 -isign*wa3[i + 1]*ci4;
            ch[ah + 3*l1*ido + 1] = wa3[i]*ci4 + isign*wa3[i + 1]*cr4;
            }
        }
        }
    } /* passf4 */


    static void passf5(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[], int isign)
    /* isign == -1 for forward transform and +1 for backward transform */
    {
        static const Treal tr11 = 0.309016994374947;
        static const Treal ti11 = 0.951056516295154;
        static const Treal tr12 = -0.809016994374947;
        static const Treal ti12 = 0.587785252292473;
        int i, k, ac, ah;
        Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
            ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
        if (ido == 2) {
        for (k = 1; k <= l1; ++k) {
            ac = (5*k - 4)*ido + 1;
            ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);
            ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);
            ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);
            ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);
            tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);
            tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);
            tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);
            tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);
            ah = (k - 1)*ido;
            ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;
            ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;
            cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;
            ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;
            cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;
            ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;
            cr5 = isign*(ti11*tr5 + ti12*tr4);
            ci5 = isign*(ti11*ti5 + ti12*ti4);
            cr4 = isign*(ti12*tr5 - ti11*tr4);
            ci4 = isign*(ti12*ti5 - ti11*ti4);
            ch[ah + l1*ido] = cr2 - ci5;
            ch[ah + 4*l1*ido] = cr2 + ci5;
            ch[ah + l1*ido + 1] = ci2 + cr5;
            ch[ah + 2*l1*ido + 1] = ci3 + cr4;
            ch[ah + 2*l1*ido] = cr3 - ci4;
            ch[ah + 3*l1*ido] = cr3 + ci4;
            ch[ah + 3*l1*ido + 1] = ci3 - cr4;
            ch[ah + 4*l1*ido + 1] = ci2 - cr5;
        }
        } else {
        for (k=1; k<=l1; k++) {
            for (i=0; i<ido-1; i+=2) {
            ac = i + 1 + (k*5 - 4)*ido;
            ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);
            ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);
            ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);
            ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);
            tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);
            tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);
            tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);
            tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);
            ah = i + (k - 1)*ido;
            ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;
            ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;
            cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;

            ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;
            cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;

            ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;
            cr5 = isign*(ti11*tr5 + ti12*tr4);
            ci5 = isign*(ti11*ti5 + ti12*ti4);
            cr4 = isign*(ti12*tr5 - ti11*tr4);
            ci4 = isign*(ti12*ti5 - ti11*ti4);
            dr3 = cr3 - ci4;
            dr4 = cr3 + ci4;
            di3 = ci3 + cr4;
            di4 = ci3 - cr4;
            dr5 = cr2 + ci5;
            dr2 = cr2 - ci5;
            di5 = ci2 - cr5;
            di2 = ci2 + cr5;
            ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
            ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
            ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
            ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
            ch[ah + 3*l1*ido] = wa3[i]*dr4 - isign*wa3[i+1]*di4;
            ch[ah + 3*l1*ido + 1] = wa3[i]*di4 + isign*wa3[i+1]*dr4;
            ch[ah + 4*l1*ido] = wa4[i]*dr5 - isign*wa4[i+1]*di5;
            ch[ah + 4*l1*ido + 1] = wa4[i]*di5 + isign*wa4[i+1]*dr5;
            }
        }
        }
    } /* passf5 */


    static void passf(int *nac, int ido, int ip, int l1, int idl1,
        Treal cc[], Treal ch[],
        const Treal wa[], int isign)
    /* isign is -1 for forward transform and +1 for backward transform */
    {
        int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, idj, idl, inc,idp;
        Treal wai, war;

        idot = ido / 2;
        /* nt = ip*idl1;*/
        ipph = (ip + 1) / 2;
        idp = ip*ido;
        if (ido >= l1) {
        for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (k=0; k<l1; k++) {
            for (i=0; i<ido; i++) {
                ch[i + (k + j*l1)*ido] =
                    ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*ip)*ido);
                ch[i + (k + jc*l1)*ido] =
                    ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*ip)*ido);
            }
            }
        }
        for (k=0; k<l1; k++)
            for (i=0; i<ido; i++)
            ch[i + k*ido] = ref(cc,i + k*ip*ido);
        } else {
        for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (i=0; i<ido; i++) {
            for (k=0; k<l1; k++) {
                ch[i + (k + j*l1)*ido] = ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*
                    ip)*ido);
                ch[i + (k + jc*l1)*ido] = ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*
                    ip)*ido);
            }
            }
        }
        for (i=0; i<ido; i++)
            for (k=0; k<l1; k++)
            ch[i + k*ido] = ref(cc,i + k*ip*ido);
        }

        idl = 2 - ido;
        inc = 0;
        for (l=1; l<ipph; l++) {
        lc = ip - l;
        idl += ido;
        for (ik=0; ik<idl1; ik++) {
            cc[ik + l*idl1] = ch[ik] + wa[idl - 2]*ch[ik + idl1];
            cc[ik + lc*idl1] = isign*wa[idl-1]*ch[ik + (ip-1)*idl1];
        }
        idlj = idl;
        inc += ido;
        for (j=2; j<ipph; j++) {
            jc = ip - j;
            idlj += inc;
            if (idlj > idp) idlj -= idp;
            war = wa[idlj - 2];
            wai = wa[idlj-1];
            for (ik=0; ik<idl1; ik++) {
            cc[ik + l*idl1] += war*ch[ik + j*idl1];
            cc[ik + lc*idl1] += isign*wai*ch[ik + jc*idl1];
            }
        }
        }
        for (j=1; j<ipph; j++)
        for (ik=0; ik<idl1; ik++)
            ch[ik] += ch[ik + j*idl1];
        for (j=1; j<ipph; j++) {
        jc = ip - j;
        for (ik=1; ik<idl1; ik+=2) {
            ch[ik - 1 + j*idl1] = cc[ik - 1 + j*idl1] - cc[ik + jc*idl1];
            ch[ik - 1 + jc*idl1] = cc[ik - 1 + j*idl1] + cc[ik + jc*idl1];
            ch[ik + j*idl1] = cc[ik + j*idl1] + cc[ik - 1 + jc*idl1];
            ch[ik + jc*idl1] = cc[ik + j*idl1] - cc[ik - 1 + jc*idl1];
        }
        }
        *nac = 1;
        if (ido == 2) return;
        *nac = 0;
        for (ik=0; ik<idl1; ik++)
        cc[ik] = ch[ik];
        for (j=1; j<ip; j++) {
        for (k=0; k<l1; k++) {
            cc[(k + j*l1)*ido + 0] = ch[(k + j*l1)*ido + 0];
            cc[(k + j*l1)*ido + 1] = ch[(k + j*l1)*ido + 1];
        }
        }
        if (idot <= l1) {
        idij = 0;
        for (j=1; j<ip; j++) {
            idij += 2;
            for (i=3; i<ido; i+=2) {
            idij += 2;
            for (k=0; k<l1; k++) {
                cc[i - 1 + (k + j*l1)*ido] =
                    wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
                    isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
                cc[i + (k + j*l1)*ido] =
                    wa[idij - 2]*ch[i + (k + j*l1)*ido] +
                    isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
            }
            }
        }
        } else {
        idj = 2 - ido;
        for (j=1; j<ip; j++) {
            idj += ido;
            for (k = 0; k < l1; k++) {
            idij = idj;
            for (i=3; i<ido; i+=2) {
                idij += 2;
                cc[i - 1 + (k + j*l1)*ido] =
                    wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
                    isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
                cc[i + (k + j*l1)*ido] =
                    wa[idij - 2]*ch[i + (k + j*l1)*ido] +
                    isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
            }
            }
        }
        }
    } /* passf */


    /* ----------------------------------------------------------------------
    radf2,radb2, radf3,radb3, radf4,radb4, radf5,radb5, radfg,radbg.
    Treal FFT passes fwd and bwd.
    ---------------------------------------------------------------------- */

    static void radf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])
    {
        int i, k, ic;
        Treal ti2, tr2;
        for (k=0; k<l1; k++) {
        ch[2*k*ido] =
            ref(cc,k*ido) + ref(cc,(k + l1)*ido);
        ch[(2*k+1)*ido + ido-1] =
            ref(cc,k*ido) - ref(cc,(k + l1)*ido);
        }
        if (ido < 2) return;
        if (ido != 2) {
        for (k=0; k<l1; k++) {
            for (i=2; i<ido; i+=2) {
            ic = ido - i;
            tr2 = wa1[i - 2]*ref(cc, i-1 + (k + l1)*ido) + wa1[i - 1]*ref(cc, i + (k + l1)*ido);
            ti2 = wa1[i - 2]*ref(cc, i + (k + l1)*ido) - wa1[i - 1]*ref(cc, i-1 + (k + l1)*ido);
            ch[i + 2*k*ido] = ref(cc,i + k*ido) + ti2;
            ch[ic + (2*k+1)*ido] = ti2 - ref(cc,i + k*ido);
            ch[i - 1 + 2*k*ido] = ref(cc,i - 1 + k*ido) + tr2;
            ch[ic - 1 + (2*k+1)*ido] = ref(cc,i - 1 + k*ido) - tr2;
            }
        }
        if (ido % 2 == 1) return;
        }
        for (k=0; k<l1; k++) {
        ch[(2*k+1)*ido] = -ref(cc,ido-1 + (k + l1)*ido);
        ch[ido-1 + 2*k*ido] = ref(cc,ido-1 + k*ido);
        }
    } /* radf2 */


    static void radb2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])
    {
        int i, k, ic;
        Treal ti2, tr2;
        for (k=0; k<l1; k++) {
        ch[k*ido] =
            ref(cc,2*k*ido) + ref(cc,ido-1 + (2*k+1)*ido);
        ch[(k + l1)*ido] =
            ref(cc,2*k*ido) - ref(cc,ido-1 + (2*k+1)*ido);
        }
        if (ido < 2) return;
        if (ido != 2) {
        for (k = 0; k < l1; ++k) {
            for (i = 2; i < ido; i += 2) {
            ic = ido - i;
            ch[i-1 + k*ido] =
                ref(cc,i-1 + 2*k*ido) + ref(cc,ic-1 + (2*k+1)*ido);
            tr2 = ref(cc,i-1 + 2*k*ido) - ref(cc,ic-1 + (2*k+1)*ido);
            ch[i + k*ido] =
                ref(cc,i + 2*k*ido) - ref(cc,ic + (2*k+1)*ido);
            ti2 = ref(cc,i + (2*k)*ido) + ref(cc,ic + (2*k+1)*ido);
            ch[i-1 + (k + l1)*ido] =
                wa1[i - 2]*tr2 - wa1[i - 1]*ti2;
            ch[i + (k + l1)*ido] =
                wa1[i - 2]*ti2 + wa1[i - 1]*tr2;
            }
        }
        if (ido % 2 == 1) return;
        }
        for (k = 0; k < l1; k++) {
        ch[ido-1 + k*ido] = 2*ref(cc,ido-1 + 2*k*ido);
        ch[ido-1 + (k + l1)*ido] = -2*ref(cc,(2*k+1)*ido);
        }
    } /* radb2 */


    static void radf3(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[])
    {
        static const Treal taur = -0.5;
        static const Treal taui = 0.866025403784439;
        int i, k, ic;
        Treal ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
        for (k=0; k<l1; k++) {
        cr2 = ref(cc,(k + l1)*ido) + ref(cc,(k + 2*l1)*ido);
        ch[3*k*ido] = ref(cc,k*ido) + cr2;
        ch[(3*k+2)*ido] = taui*(ref(cc,(k + l1*2)*ido) - ref(cc,(k + l1)*ido));
        ch[ido-1 + (3*k + 1)*ido] = ref(cc,k*ido) + taur*cr2;
        }
        if (ido == 1) return;
        for (k=0; k<l1; k++) {
        for (i=2; i<ido; i+=2) {
            ic = ido - i;
            dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) +
                wa1[i - 1]*ref(cc,i + (k + l1)*ido);
            di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
            dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + l1*2)*ido) + wa2[i - 1]*ref(cc,i + (k + l1*2)*ido);
            di3 = wa2[i - 2]*ref(cc,i + (k + l1*2)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + l1*2)*ido);
            cr2 = dr2 + dr3;
            ci2 = di2 + di3;
            ch[i - 1 + 3*k*ido] = ref(cc,i - 1 + k*ido) + cr2;
            ch[i + 3*k*ido] = ref(cc,i + k*ido) + ci2;
            tr2 = ref(cc,i - 1 + k*ido) + taur*cr2;
            ti2 = ref(cc,i + k*ido) + taur*ci2;
            tr3 = taui*(di2 - di3);
            ti3 = taui*(dr3 - dr2);
            ch[i - 1 + (3*k + 2)*ido] = tr2 + tr3;
            ch[ic - 1 + (3*k + 1)*ido] = tr2 - tr3;
            ch[i + (3*k + 2)*ido] = ti2 + ti3;
            ch[ic + (3*k + 1)*ido] = ti3 - ti2;
        }
        }
    } /* radf3 */


    static void radb3(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[])
    {
        static const Treal taur = -0.5;
        static const Treal taui = 0.866025403784439;
        int i, k, ic;
        Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
        for (k=0; k<l1; k++) {
        tr2 = 2*ref(cc,ido-1 + (3*k + 1)*ido);
        cr2 = ref(cc,3*k*ido) + taur*tr2;
        ch[k*ido] = ref(cc,3*k*ido) + tr2;
        ci3 = 2*taui*ref(cc,(3*k + 2)*ido);
        ch[(k + l1)*ido] = cr2 - ci3;
        ch[(k + 2*l1)*ido] = cr2 + ci3;
        }
        if (ido == 1) return;
        for (k=0; k<l1; k++) {
        for (i=2; i<ido; i+=2) {
            ic = ido - i;
            tr2 = ref(cc,i - 1 + (3*k + 2)*ido) + ref(cc,ic - 1 + (3*k + 1)*ido);
            cr2 = ref(cc,i - 1 + 3*k*ido) + taur*tr2;
            ch[i - 1 + k*ido] = ref(cc,i - 1 + 3*k*ido) + tr2;
            ti2 = ref(cc,i + (3*k + 2)*ido) - ref(cc,ic + (3*k + 1)*ido);
            ci2 = ref(cc,i + 3*k*ido) + taur*ti2;
            ch[i + k*ido] = ref(cc,i + 3*k*ido) + ti2;
            cr3 = taui*(ref(cc,i - 1 + (3*k + 2)*ido) - ref(cc,ic - 1 + (3*k + 1)*ido));
            ci3 = taui*(ref(cc,i + (3*k + 2)*ido) + ref(cc,ic + (3*k + 1)*ido));
            dr2 = cr2 - ci3;
            dr3 = cr2 + ci3;
            di2 = ci2 + cr3;
            di3 = ci2 - cr3;
            ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
            ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
            ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
            ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
        }
        }
    } /* radb3 */


    static void radf4(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[], const Treal wa3[])
    {
        static const Treal hsqt2 = 0.7071067811865475;
        int i, k, ic;
        Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
        for (k=0; k<l1; k++) {
        tr1 = ref(cc,(k + l1)*ido) + ref(cc,(k + 3*l1)*ido);
        tr2 = ref(cc,k*ido) + ref(cc,(k + 2*l1)*ido);
        ch[4*k*ido] = tr1 + tr2;
        ch[ido-1 + (4*k + 3)*ido] = tr2 - tr1;
        ch[ido-1 + (4*k + 1)*ido] = ref(cc,k*ido) - ref(cc,(k + 2*l1)*ido);
        ch[(4*k + 2)*ido] = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + l1)*ido);
        }
        if (ido < 2) return;
        if (ido != 2) {
        for (k=0; k<l1; k++) {
            for (i=2; i<ido; i += 2) {
            ic = ido - i;
            cr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
            ci2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
            cr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*
                ido);
            ci3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*
                ido);
            cr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*
                ido);
            ci4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*
                ido);
            tr1 = cr2 + cr4;
            tr4 = cr4 - cr2;
            ti1 = ci2 + ci4;
            ti4 = ci2 - ci4;
            ti2 = ref(cc,i + k*ido) + ci3;
            ti3 = ref(cc,i + k*ido) - ci3;
            tr2 = ref(cc,i - 1 + k*ido) + cr3;
            tr3 = ref(cc,i - 1 + k*ido) - cr3;
            ch[i - 1 + 4*k*ido] = tr1 + tr2;
            ch[ic - 1 + (4*k + 3)*ido] = tr2 - tr1;
            ch[i + 4*k*ido] = ti1 + ti2;
            ch[ic + (4*k + 3)*ido] = ti1 - ti2;
            ch[i - 1 + (4*k + 2)*ido] = ti4 + tr3;
            ch[ic - 1 + (4*k + 1)*ido] = tr3 - ti4;
            ch[i + (4*k + 2)*ido] = tr4 + ti3;
            ch[ic + (4*k + 1)*ido] = tr4 - ti3;
            }
        }
        if (ido % 2 == 1) return;
        }
        for (k=0; k<l1; k++) {
        ti1 = -hsqt2*(ref(cc,ido-1 + (k + l1)*ido) + ref(cc,ido-1 + (k + 3*l1)*ido));
        tr1 = hsqt2*(ref(cc,ido-1 + (k + l1)*ido) - ref(cc,ido-1 + (k + 3*l1)*ido));
        ch[ido-1 + 4*k*ido] = tr1 + ref(cc,ido-1 + k*ido);
        ch[ido-1 + (4*k + 2)*ido] = ref(cc,ido-1 + k*ido) - tr1;
        ch[(4*k + 1)*ido] = ti1 - ref(cc,ido-1 + (k + 2*l1)*ido);
        ch[(4*k + 3)*ido] = ti1 + ref(cc,ido-1 + (k + 2*l1)*ido);
        }
    } /* radf4 */


    static void radb4(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[], const Treal wa3[])
    {
        static const Treal sqrt2 = 1.414213562373095;
        int i, k, ic;
        Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
        for (k = 0; k < l1; k++) {
        tr1 = ref(cc,4*k*ido) - ref(cc,ido-1 + (4*k + 3)*ido);
        tr2 = ref(cc,4*k*ido) + ref(cc,ido-1 + (4*k + 3)*ido);
        tr3 = ref(cc,ido-1 + (4*k + 1)*ido) + ref(cc,ido-1 + (4*k + 1)*ido);
        tr4 = ref(cc,(4*k + 2)*ido) + ref(cc,(4*k + 2)*ido);
        ch[k*ido] = tr2 + tr3;
        ch[(k + l1)*ido] = tr1 - tr4;
        ch[(k + 2*l1)*ido] = tr2 - tr3;
        ch[(k + 3*l1)*ido] = tr1 + tr4;
        }
        if (ido < 2) return;
        if (ido != 2) {
        for (k = 0; k < l1; ++k) {
            for (i = 2; i < ido; i += 2) {
            ic = ido - i;
            ti1 = ref(cc,i + 4*k*ido) + ref(cc,ic + (4*k + 3)*ido);
            ti2 = ref(cc,i + 4*k*ido) - ref(cc,ic + (4*k + 3)*ido);
            ti3 = ref(cc,i + (4*k + 2)*ido) - ref(cc,ic + (4*k + 1)*ido);
            tr4 = ref(cc,i + (4*k + 2)*ido) + ref(cc,ic + (4*k + 1)*ido);
            tr1 = ref(cc,i - 1 + 4*k*ido) - ref(cc,ic - 1 + (4*k + 3)*ido);
            tr2 = ref(cc,i - 1 + 4*k*ido) + ref(cc,ic - 1 + (4*k + 3)*ido);
            ti4 = ref(cc,i - 1 + (4*k + 2)*ido) - ref(cc,ic - 1 + (4*k + 1)*ido);
            tr3 = ref(cc,i - 1 + (4*k + 2)*ido) + ref(cc,ic - 1 + (4*k + 1)*ido);
            ch[i - 1 + k*ido] = tr2 + tr3;
            cr3 = tr2 - tr3;
            ch[i + k*ido] = ti2 + ti3;
            ci3 = ti2 - ti3;
            cr2 = tr1 - tr4;
            cr4 = tr1 + tr4;
            ci2 = ti1 + ti4;
            ci4 = ti1 - ti4;
            ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*cr2 - wa1[i - 1]*ci2;
            ch[i + (k + l1)*ido] = wa1[i - 2]*ci2 + wa1[i - 1]*cr2;
            ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*cr3 - wa2[i - 1]*ci3;
            ch[i + (k + 2*l1)*ido] = wa2[i - 2]*ci3 + wa2[i - 1]*cr3;
            ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*cr4 - wa3[i - 1]*ci4;
            ch[i + (k + 3*l1)*ido] = wa3[i - 2]*ci4 + wa3[i - 1]*cr4;
            }
        }
        if (ido % 2 == 1) return;
        }
        for (k = 0; k < l1; k++) {
        ti1 = ref(cc,(4*k + 1)*ido) + ref(cc,(4*k + 3)*ido);
        ti2 = ref(cc,(4*k + 3)*ido) - ref(cc,(4*k + 1)*ido);
        tr1 = ref(cc,ido-1 + 4*k*ido) - ref(cc,ido-1 + (4*k + 2)*ido);
        tr2 = ref(cc,ido-1 + 4*k*ido) + ref(cc,ido-1 + (4*k + 2)*ido);
        ch[ido-1 + k*ido] = tr2 + tr2;
        ch[ido-1 + (k + l1)*ido] = sqrt2*(tr1 - ti1);
        ch[ido-1 + (k + 2*l1)*ido] = ti2 + ti2;
        ch[ido-1 + (k + 3*l1)*ido] = -sqrt2*(tr1 + ti1);
        }
    } /* radb4 */


    static void radf5(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])
    {
        static const Treal tr11 = 0.309016994374947;
        static const Treal ti11 = 0.951056516295154;
        static const Treal tr12 = -0.809016994374947;
        static const Treal ti12 = 0.587785252292473;
        int i, k, ic;
        Treal ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5,
            cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
        for (k = 0; k < l1; k++) {
        cr2 = ref(cc,(k + 4*l1)*ido) + ref(cc,(k + l1)*ido);
        ci5 = ref(cc,(k + 4*l1)*ido) - ref(cc,(k + l1)*ido);
        cr3 = ref(cc,(k + 3*l1)*ido) + ref(cc,(k + 2*l1)*ido);
        ci4 = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + 2*l1)*ido);
        ch[5*k*ido] = ref(cc,k*ido) + cr2 + cr3;
        ch[ido-1 + (5*k + 1)*ido] = ref(cc,k*ido) + tr11*cr2 + tr12*cr3;
        ch[(5*k + 2)*ido] = ti11*ci5 + ti12*ci4;
        ch[ido-1 + (5*k + 3)*ido] = ref(cc,k*ido) + tr12*cr2 + tr11*cr3;
        ch[(5*k + 4)*ido] = ti12*ci5 - ti11*ci4;
        }
        if (ido == 1) return;
        for (k = 0; k < l1; ++k) {
        for (i = 2; i < ido; i += 2) {
            ic = ido - i;
            dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
            di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
            dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*ido);
            di3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*ido);
            dr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*ido);
            di4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*ido);
            dr5 = wa4[i - 2]*ref(cc,i - 1 + (k + 4*l1)*ido) + wa4[i - 1]*ref(cc,i + (k + 4*l1)*ido);
            di5 = wa4[i - 2]*ref(cc,i + (k + 4*l1)*ido) - wa4[i - 1]*ref(cc,i - 1 + (k + 4*l1)*ido);
            cr2 = dr2 + dr5;
            ci5 = dr5 - dr2;
            cr5 = di2 - di5;
            ci2 = di2 + di5;
            cr3 = dr3 + dr4;
            ci4 = dr4 - dr3;
            cr4 = di3 - di4;
            ci3 = di3 + di4;
            ch[i - 1 + 5*k*ido] = ref(cc,i - 1 + k*ido) + cr2 + cr3;
            ch[i + 5*k*ido] = ref(cc,i + k*ido) + ci2 + ci3;
            tr2 = ref(cc,i - 1 + k*ido) + tr11*cr2 + tr12*cr3;
            ti2 = ref(cc,i + k*ido) + tr11*ci2 + tr12*ci3;
            tr3 = ref(cc,i - 1 + k*ido) + tr12*cr2 + tr11*cr3;
            ti3 = ref(cc,i + k*ido) + tr12*ci2 + tr11*ci3;
            tr5 = ti11*cr5 + ti12*cr4;
            ti5 = ti11*ci5 + ti12*ci4;
            tr4 = ti12*cr5 - ti11*cr4;
            ti4 = ti12*ci5 - ti11*ci4;
            ch[i - 1 + (5*k + 2)*ido] = tr2 + tr5;
            ch[ic - 1 + (5*k + 1)*ido] = tr2 - tr5;
            ch[i + (5*k + 2)*ido] = ti2 + ti5;
            ch[ic + (5*k + 1)*ido] = ti5 - ti2;
            ch[i - 1 + (5*k + 4)*ido] = tr3 + tr4;
            ch[ic - 1 + (5*k + 3)*ido] = tr3 - tr4;
            ch[i + (5*k + 4)*ido] = ti3 + ti4;
            ch[ic + (5*k + 3)*ido] = ti4 - ti3;
        }
        }
    } /* radf5 */


    static void radb5(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])
    {
        static const Treal tr11 = 0.309016994374947;
        static const Treal ti11 = 0.951056516295154;
        static const Treal tr12 = -0.809016994374947;
        static const Treal ti12 = 0.587785252292473;
        int i, k, ic;
        Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
            ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
        for (k = 0; k < l1; k++) {
        ti5 = 2*ref(cc,(5*k + 2)*ido);
        ti4 = 2*ref(cc,(5*k + 4)*ido);
        tr2 = 2*ref(cc,ido-1 + (5*k + 1)*ido);
        tr3 = 2*ref(cc,ido-1 + (5*k + 3)*ido);
        ch[k*ido] = ref(cc,5*k*ido) + tr2 + tr3;
        cr2 = ref(cc,5*k*ido) + tr11*tr2 + tr12*tr3;
        cr3 = ref(cc,5*k*ido) + tr12*tr2 + tr11*tr3;
        ci5 = ti11*ti5 + ti12*ti4;
        ci4 = ti12*ti5 - ti11*ti4;
        ch[(k + l1)*ido] = cr2 - ci5;
        ch[(k + 2*l1)*ido] = cr3 - ci4;
        ch[(k + 3*l1)*ido] = cr3 + ci4;
        ch[(k + 4*l1)*ido] = cr2 + ci5;
        }
        if (ido == 1) return;
        for (k = 0; k < l1; ++k) {
        for (i = 2; i < ido; i += 2) {
            ic = ido - i;
            ti5 = ref(cc,i + (5*k + 2)*ido) + ref(cc,ic + (5*k + 1)*ido);
            ti2 = ref(cc,i + (5*k + 2)*ido) - ref(cc,ic + (5*k + 1)*ido);
            ti4 = ref(cc,i + (5*k + 4)*ido) + ref(cc,ic + (5*k + 3)*ido);
            ti3 = ref(cc,i + (5*k + 4)*ido) - ref(cc,ic + (5*k + 3)*ido);
            tr5 = ref(cc,i - 1 + (5*k + 2)*ido) - ref(cc,ic - 1 + (5*k + 1)*ido);
            tr2 = ref(cc,i - 1 + (5*k + 2)*ido) + ref(cc,ic - 1 + (5*k + 1)*ido);
            tr4 = ref(cc,i - 1 + (5*k + 4)*ido) - ref(cc,ic - 1 + (5*k + 3)*ido);
            tr3 = ref(cc,i - 1 + (5*k + 4)*ido) + ref(cc,ic - 1 + (5*k + 3)*ido);
            ch[i - 1 + k*ido] = ref(cc,i - 1 + 5*k*ido) + tr2 + tr3;
            ch[i + k*ido] = ref(cc,i + 5*k*ido) + ti2 + ti3;
            cr2 = ref(cc,i - 1 + 5*k*ido) + tr11*tr2 + tr12*tr3;

            ci2 = ref(cc,i + 5*k*ido) + tr11*ti2 + tr12*ti3;
            cr3 = ref(cc,i - 1 + 5*k*ido) + tr12*tr2 + tr11*tr3;

            ci3 = ref(cc,i + 5*k*ido) + tr12*ti2 + tr11*ti3;
            cr5 = ti11*tr5 + ti12*tr4;
            ci5 = ti11*ti5 + ti12*ti4;
            cr4 = ti12*tr5 - ti11*tr4;
            ci4 = ti12*ti5 - ti11*ti4;
            dr3 = cr3 - ci4;
            dr4 = cr3 + ci4;
            di3 = ci3 + cr4;
            di4 = ci3 - cr4;
            dr5 = cr2 + ci5;
            dr2 = cr2 - ci5;
            di5 = ci2 - cr5;
            di2 = ci2 + cr5;
            ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
            ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
            ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
            ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
            ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*dr4 - wa3[i - 1]*di4;
            ch[i + (k + 3*l1)*ido] = wa3[i - 2]*di4 + wa3[i - 1]*dr4;
            ch[i - 1 + (k + 4*l1)*ido] = wa4[i - 2]*dr5 - wa4[i - 1]*di5;
            ch[i + (k + 4*l1)*ido] = wa4[i - 2]*di5 + wa4[i - 1]*dr5;
        }
        }
    } /* radb5 */


    static void radfg(int ido, int ip, int l1, int idl1,
        Treal cc[], Treal ch[], const Treal wa[])
    {
        static const Treal twopi = 6.28318530717959;
        int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;
        Treal dc2, ai1, ai2, ar1, ar2, ds2, dcp, arg, dsp, ar1h, ar2h;
        arg = twopi / ip;
        dcp = cos(arg);
        dsp = sin(arg);
        ipph = (ip + 1) / 2;
        nbd = (ido - 1) / 2;
        if (ido != 1) {
        for (ik=0; ik<idl1; ik++) ch[ik] = cc[ik];
        for (j=1; j<ip; j++)
            for (k=0; k<l1; k++)
            ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido];
        if (nbd <= l1) {
            is = -ido;
            for (j=1; j<ip; j++) {
            is += ido;
            idij = is-1;
            for (i=2; i<ido; i+=2) {
                idij += 2;
                for (k=0; k<l1; k++) {
                ch[i - 1 + (k + j*l1)*ido] =
                    wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
                ch[i + (k + j*l1)*ido] =
                    wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
                }
            }
            }
        } else {
            is = -ido;
            for (j=1; j<ip; j++) {
            is += ido;
            for (k=0; k<l1; k++) {
                idij = is-1;
                for (i=2; i<ido; i+=2) {
                idij += 2;
                ch[i - 1 + (k + j*l1)*ido] =
                    wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
                ch[i + (k + j*l1)*ido] =
                    wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
                }
            }
            }
        }
        if (nbd >= l1) {
            for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (k=0; k<l1; k++) {
                for (i=2; i<ido; i+=2) {
                cc[i - 1 + (k + j*l1)*ido] = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
                cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
                cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
                cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
                }
            }
            }
        } else {
            for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (i=2; i<ido; i+=2) {
                for (k=0; k<l1; k++) {
                cc[i - 1 + (k + j*l1)*ido] =
                    ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
                cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
                cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
                cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
                }
            }
            }
        }
        } else {  /* now ido == 1 */
        for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];
        }
        for (j=1; j<ipph; j++) {
        jc = ip - j;
        for (k=0; k<l1; k++) {
            cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido] + ch[(k + jc*l1)*ido];
            cc[(k + jc*l1)*ido] = ch[(k + jc*l1)*ido] - ch[(k + j*l1)*ido];
        }
        }

        ar1 = 1;
        ai1 = 0;
        for (l=1; l<ipph; l++) {
        lc = ip - l;
        ar1h = dcp*ar1 - dsp*ai1;
        ai1 = dcp*ai1 + dsp*ar1;
        ar1 = ar1h;
        for (ik=0; ik<idl1; ik++) {
            ch[ik + l*idl1] = cc[ik] + ar1*cc[ik + idl1];
            ch[ik + lc*idl1] = ai1*cc[ik + (ip-1)*idl1];
        }
        dc2 = ar1;
        ds2 = ai1;
        ar2 = ar1;
        ai2 = ai1;
        for (j=2; j<ipph; j++) {
            jc = ip - j;
            ar2h = dc2*ar2 - ds2*ai2;
            ai2 = dc2*ai2 + ds2*ar2;
            ar2 = ar2h;
            for (ik=0; ik<idl1; ik++) {
            ch[ik + l*idl1] += ar2*cc[ik + j*idl1];
            ch[ik + lc*idl1] += ai2*cc[ik + jc*idl1];
            }
        }
        }
        for (j=1; j<ipph; j++)
        for (ik=0; ik<idl1; ik++)
            ch[ik] += cc[ik + j*idl1];

        if (ido >= l1) {
        for (k=0; k<l1; k++) {
            for (i=0; i<ido; i++) {
            ref(cc,i + k*ip*ido) = ch[i + k*ido];
            }
        }
        } else {
        for (i=0; i<ido; i++) {
            for (k=0; k<l1; k++) {
            ref(cc,i + k*ip*ido) = ch[i + k*ido];
            }
        }
        }
        for (j=1; j<ipph; j++) {
        jc = ip - j;
        j2 = 2*j;
        for (k=0; k<l1; k++) {
            ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) =
                ch[(k + j*l1)*ido];
            ref(cc,(j2 + k*ip)*ido) =
                ch[(k + jc*l1)*ido];
        }
        }
        if (ido == 1) return;
        if (nbd >= l1) {
        for (j=1; j<ipph; j++) {
            jc = ip - j;
            j2 = 2*j;
            for (k=0; k<l1; k++) {
            for (i=2; i<ido; i+=2) {
                ic = ido - i;
                ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
                ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
                ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
                ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
            }
            }
        }
        } else {
        for (j=1; j<ipph; j++) {
            jc = ip - j;
            j2 = 2*j;
            for (i=2; i<ido; i+=2) {
            ic = ido - i;
            for (k=0; k<l1; k++) {
                ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
                ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
                ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
                ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
            }
            }
        }
        }
    } /* radfg */


    static void radbg(int ido, int ip, int l1, int idl1,
        Treal cc[], Treal ch[], const Treal wa[])
    {
        static const Treal twopi = 6.28318530717959;
        int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
        Treal dc2, ai1, ai2, ar1, ar2, ds2;
        int nbd;
        Treal dcp, arg, dsp, ar1h, ar2h;
        arg = twopi / ip;
        dcp = cos(arg);
        dsp = sin(arg);
        nbd = (ido - 1) / 2;
        ipph = (ip + 1) / 2;
        if (ido >= l1) {
        for (k=0; k<l1; k++) {
            for (i=0; i<ido; i++) {
            ch[i + k*ido] = ref(cc,i + k*ip*ido);
            }
        }
        } else {
        for (i=0; i<ido; i++) {
            for (k=0; k<l1; k++) {
            ch[i + k*ido] = ref(cc,i + k*ip*ido);
            }
        }
        }
        for (j=1; j<ipph; j++) {
        jc = ip - j;
        j2 = 2*j;
        for (k=0; k<l1; k++) {
            ch[(k + j*l1)*ido] = ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) + ref(cc,ido-1 + (j2 - 1 + k*ip)*
                ido);
            ch[(k + jc*l1)*ido] = ref(cc,(j2 + k*ip)*ido) + ref(cc,(j2 + k*ip)*ido);
        }
        }

        if (ido != 1) {
        if (nbd >= l1) {
            for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (k=0; k<l1; k++) {
                for (i=2; i<ido; i+=2) {
                ic = ido - i;
                ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,
                    ic - 1 + (2*j - 1 + k*ip)*ido);
                ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -
                    ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
                ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic
                    + (2*j - 1 + k*ip)*ido);
                ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic
                    + (2*j - 1 + k*ip)*ido);
                }
            }
            }
        } else {
            for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (i=2; i<ido; i+=2) {
                ic = ido - i;
                for (k=0; k<l1; k++) {
                ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,
                    ic - 1 + (2*j - 1 + k*ip)*ido);
                ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -
                    ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
                ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic
                    + (2*j - 1 + k*ip)*ido);
                ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic
                    + (2*j - 1 + k*ip)*ido);
                }
            }
            }
        }
        }

        ar1 = 1;
        ai1 = 0;
        for (l=1; l<ipph; l++) {
        lc = ip - l;
        ar1h = dcp*ar1 - dsp*ai1;
        ai1 = dcp*ai1 + dsp*ar1;
        ar1 = ar1h;
        for (ik=0; ik<idl1; ik++) {
            cc[ik + l*idl1] = ch[ik] + ar1*ch[ik + idl1];
            cc[ik + lc*idl1] = ai1*ch[ik + (ip-1)*idl1];
        }
        dc2 = ar1;
        ds2 = ai1;
        ar2 = ar1;
        ai2 = ai1;
        for (j=2; j<ipph; j++) {
            jc = ip - j;
            ar2h = dc2*ar2 - ds2*ai2;
            ai2 = dc2*ai2 + ds2*ar2;
            ar2 = ar2h;
            for (ik=0; ik<idl1; ik++) {
            cc[ik + l*idl1] += ar2*ch[ik + j*idl1];
            cc[ik + lc*idl1] += ai2*ch[ik + jc*idl1];
            }
        }
        }
        for (j=1; j<ipph; j++) {
        for (ik=0; ik<idl1; ik++) {
            ch[ik] += ch[ik + j*idl1];
        }
        }
        for (j=1; j<ipph; j++) {
        jc = ip - j;
        for (k=0; k<l1; k++) {
            ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido] - cc[(k + jc*l1)*ido];
            ch[(k + jc*l1)*ido] = cc[(k + j*l1)*ido] + cc[(k + jc*l1)*ido];
        }
        }

        if (ido == 1) return;
        if (nbd >= l1) {
        for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (k=0; k<l1; k++) {
            for (i=2; i<ido; i+=2) {
                ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
                ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] + cc[i + (k + jc*l1)*ido];
                ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
                ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
            }
            }
        }
        } else {
        for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (i=2; i<ido; i+=2) {
            for (k=0; k<l1; k++) {
                ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
                ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j *l1)*ido] + cc[i + (k + jc*l1)*ido];
                ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
                ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
            }
            }
        }
        }
        for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];
        for (j=1; j<ip; j++)
        for (k=0; k<l1; k++)
            cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido];
        if (nbd <= l1) {
        is = -ido;
        for (j=1; j<ip; j++) {
            is += ido;
            idij = is-1;
            for (i=2; i<ido; i+=2) {
            idij += 2;
            for (k=0; k<l1; k++) {
                cc[i - 1 + (k + j*l1)*ido] = wa[idij - 1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*
                    ch[i + (k + j*l1)*ido];
                cc[i + (k + j*l1)*ido] = wa[idij - 1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
            }
            }
        }
        } else {
        is = -ido;
        for (j=1; j<ip; j++) {
            is += ido;
            for (k=0; k<l1; k++) {
            idij = is - 1;
            for (i=2; i<ido; i+=2) {
                idij += 2;
                cc[i - 1 + (k + j*l1)*ido] = wa[idij-1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*
                    ch[i + (k + j*l1)*ido];
                cc[i + (k + j*l1)*ido] = wa[idij-1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
            }
            }
        }
        }
    } /* radbg */

    /* ----------------------------------------------------------------------
    cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
    ---------------------------------------------------------------------- */

    static void cfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2], int isign)
    {
        int idot, i;
        int k1, l1, l2;
        int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;
        Treal *cinput, *coutput;
        nf = ifac[1];
        na = 0;
        l1 = 1;
        iw = 0;
        for (k1=2; k1<=nf+1; k1++) {
        ip = ifac[k1];
        l2 = ip*l1;
        ido = n / l2;
        idot = ido + ido;
        idl1 = idot*l1;
        if (na) {
            cinput = ch;
            coutput = c;
        } else {
            cinput = c;
            coutput = ch;
        }
        switch (ip) {
        case 4:
            ix2 = iw + idot;
            ix3 = ix2 + idot;
            passf4(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], isign);
            na = !na;
            break;
        case 2:
            passf2(idot, l1, cinput, coutput, &wa[iw], isign);
            na = !na;
            break;
        case 3:
            ix2 = iw + idot;
            passf3(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], isign);
            na = !na;
            break;
        case 5:
            ix2 = iw + idot;
            ix3 = ix2 + idot;
            ix4 = ix3 + idot;
            passf5(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
            na = !na;
            break;
        default:
            passf(&nac, idot, ip, l1, idl1, cinput, coutput, &wa[iw], isign);
            if (nac != 0) na = !na;
        }
        l1 = l2;
        iw += (ip - 1)*idot;
        }
        if (na == 0) return;
        for (i=0; i<2*n; i++) c[i] = ch[i];
    } /* cfftf1 */


    void cfftf(int n, Treal c[], Treal wsave[])
    {
        int iw1, iw2;
        if (n == 1) return;
        iw1 = 2*n;
        iw2 = iw1 + 2*n;
        cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), -1);
    } /* cfftf */


    void cfftb(int n, Treal c[], Treal wsave[])
    {
        int iw1, iw2;
        if (n == 1) return;
        iw1 = 2*n;
        iw2 = iw1 + 2*n;
        cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), +1);
    } /* cfftb */


    static void factorize(int n, int ifac[MAXFAC+2], const int ntryh[NSPECIAL])
    /* Factorize n in factors in ntryh and rest. On exit,
    ifac[0] contains n and ifac[1] contains number of factors,
    the factors start from ifac[2]. */
    {
        int ntry=3, i, j=0, ib, nf=0, nl=n, nq, nr;
    startloop:
        if (j < NSPECIAL)
        ntry = ntryh[j];
        else
        ntry+= 2;
        j++;
        do {
        nq = nl / ntry;
        nr = nl - ntry*nq;
        if (nr != 0) goto startloop;
        nf++;
        ifac[nf + 1] = ntry;
        nl = nq;
        if (ntry == 2 && nf != 1) {
            for (i=2; i<=nf; i++) {
            ib = nf - i + 2;
            ifac[ib + 1] = ifac[ib];
            }
            ifac[2] = 2;
        }
        } while (nl != 1);
        ifac[0] = n;
        ifac[1] = nf;
    }


    static void cffti1(int n, Treal wa[], int ifac[MAXFAC+2])
    {
        static const Treal twopi = 6.28318530717959;
        Treal arg, argh, argld, fi;
        int idot, i, j;
        int i1, k1, l1, l2;
        int ld, ii, nf, ip;
        int ido, ipm;

        static const int ntryh[NSPECIAL] = {
        3,4,2,5    }; /* Do not change the order of these. */

        factorize(n,ifac,ntryh);
        nf = ifac[1];
        argh = twopi/(Treal)n;
        i = 1;
        l1 = 1;
        for (k1=1; k1<=nf; k1++) {
        ip = ifac[k1+1];
        ld = 0;
        l2 = l1*ip;
        ido = n / l2;
        idot = ido + ido + 2;
        ipm = ip - 1;
        for (j=1; j<=ipm; j++) {
            i1 = i;
            wa[i-1] = 1;
            wa[i] = 0;
            ld += l1;
            fi = 0;
            argld = ld*argh;
            for (ii=4; ii<=idot; ii+=2) {
            i+= 2;
            fi+= 1;
            arg = fi*argld;
            wa[i-1] = cos(arg);
            wa[i] = sin(arg);
            }
            if (ip > 5) {
            wa[i1-1] = wa[i-1];
            wa[i1] = wa[i];
            }
        }
        l1 = l2;
        }
    } /* cffti1 */


    void cffti(int n, Treal wsave[])
    {
        int iw1, iw2;
        if (n == 1) return;
        iw1 = 2*n;
        iw2 = iw1 + 2*n;
        cffti1(n, wsave+iw1, (int*)(wsave+iw2));
    } /* cffti */

    /* ----------------------------------------------------------------------
    rfftf1, rfftb1, rfftf, rfftb, rffti1, rffti. Treal FFTs.
    ---------------------------------------------------------------------- */

    static void rfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])
    {
        int i;
        int k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
        Treal *cinput, *coutput;
        nf = ifac[1];
        na = 1;
        l2 = n;
        iw = n-1;
        for (k1 = 1; k1 <= nf; ++k1) {
        kh = nf - k1;
        ip = ifac[kh + 2];
        l1 = l2 / ip;
        ido = n / l2;
        idl1 = ido*l1;
        iw -= (ip - 1)*ido;
        na = !na;
        if (na) {
            cinput = ch;
            coutput = c;
        } else {
            cinput = c;
            coutput = ch;
        }
        switch (ip) {
        case 4:
            ix2 = iw + ido;
            ix3 = ix2 + ido;
            radf4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
            break;
        case 2:
            radf2(ido, l1, cinput, coutput, &wa[iw]);
            break;
        case 3:
            ix2 = iw + ido;
            radf3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
            break;
        case 5:
            ix2 = iw + ido;
            ix3 = ix2 + ido;
            ix4 = ix3 + ido;
            radf5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
            break;
        default:
            if (ido == 1)
            na = !na;
            if (na == 0) {
            radfg(ido, ip, l1, idl1, c, ch, &wa[iw]);
            na = 1;
            } else {
            radfg(ido, ip, l1, idl1, ch, c, &wa[iw]);
            na = 0;
            }
        }
        l2 = l1;
        }
        if (na == 1) return;
        for (i = 0; i < n; i++) c[i] = ch[i];
    } /* rfftf1 */


    void rfftb1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])
    {
        int i;
        int k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
        Treal *cinput, *coutput;
        nf = ifac[1];
        na = 0;
        l1 = 1;
        iw = 0;
        for (k1=1; k1<=nf; k1++) {
        ip = ifac[k1 + 1];
        l2 = ip*l1;
        ido = n / l2;
        idl1 = ido*l1;
        if (na) {
            cinput = ch;
            coutput = c;
        } else {
            cinput = c;
            coutput = ch;
        }
        switch (ip) {
        case 4:
            ix2 = iw + ido;
            ix3 = ix2 + ido;
            radb4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
            na = !na;
            break;
        case 2:
            radb2(ido, l1, cinput, coutput, &wa[iw]);
            na = !na;
            break;
        case 3:
            ix2 = iw + ido;
            radb3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
            na = !na;
            break;
        case 5:
            ix2 = iw + ido;
            ix3 = ix2 + ido;
            ix4 = ix3 + ido;
            radb5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
            na = !na;
            break;
        default:
            radbg(ido, ip, l1, idl1, cinput, coutput, &wa[iw]);
            if (ido == 1) na = !na;
        }
        l1 = l2;
        iw += (ip - 1)*ido;
        }
        if (na == 0) return;
        for (i=0; i<n; i++) c[i] = ch[i];
    } /* rfftb1 */


    void rfftf(int n, Treal r[], Treal wsave[])
    {
        if (n == 1) return;
        rfftf1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
    } /* rfftf */


    void rfftb(int n, Treal r[], Treal wsave[])
    {
        if (n == 1) return;
        rfftb1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
    } /* rfftb */


    static void rffti1(int n, Treal wa[], int ifac[MAXFAC+2])
    {
        static const Treal twopi = 6.28318530717959;
        Treal arg, argh, argld, fi;
        int i, j;
        int k1, l1, l2;
        int ld, ii, nf, ip, is;
        int ido, ipm, nfm1;
        static const int ntryh[NSPECIAL] = {
        4,2,3,5    }; /* Do not change the order of these. */
        factorize(n,ifac,ntryh);
        nf = ifac[1];
        argh = twopi / n;
        is = 0;
        nfm1 = nf - 1;
        l1 = 1;
        if (nfm1 == 0) return;
        for (k1 = 1; k1 <= nfm1; k1++) {
        ip = ifac[k1 + 1];
        ld = 0;
        l2 = l1*ip;
        ido = n / l2;
        ipm = ip - 1;
        for (j = 1; j <= ipm; ++j) {
            ld += l1;
            i = is;
            argld = (Treal) ld*argh;
            fi = 0;
            for (ii = 3; ii <= ido; ii += 2) {
            i += 2;
            fi += 1;
            arg = fi*argld;
            wa[i - 2] = cos(arg);
            wa[i - 1] = sin(arg);
            }
            is += ido;
        }
        l1 = l2;
        }
    } /* rffti1 */


    void rffti(int n, Treal wsave[])
    {
        if (n == 1) return;
        rffti1(n, wsave+n, (int*)(wsave+2*n));
    } /* rffti */

    #ifdef __cplusplus
    }
    #endif
''', libraries=[])

import numpy as np

def cffti(n):
    '''
    void cffti(int n, Treal wsave[])
    '''
    wsave = np.empty(4*n+15, 'double') # don't even ask me why 4 and why 15. Just took them from C-ext wrapper "as is".
    wsave_ptr = wsave.__array_interface__['data'][0]
    C.cffti(n, ffi.cast('double*', wsave_ptr))
    return wsave


def cfftf(a, wsave):
    '''
    void cfftf(int n, Troeal c[], Treal wsave[]);
    '''
    a_copy = a.copy()
    npts = a.shape[-1]
    if wsave.shape[0] != npts * 4 + 15:
        raise ValueError("invalid work array for fft size")
    nrepeats = a_copy.size / npts
    n = ffi.cast('int', npts)
    wsave_ptr = wsave.__array_interface__['data'][0]
    for i in range(nrepeats):
        print 'npts,i,a.shape',npts,i,a.shape,a.dtype
        dptr = a_copy.__array_interface__['data'][0] + npts * i
        a_copy_cdata = ffi.cast('double*', dptr)
        C.cfftf(n, a_copy_cdata, ffi.cast('double*', wsave_ptr))
    return a_copy

# is it a good idea to put it on a top-level here?.. maybe some singleton would be nicer or something...
from cffi import FFI
ffi = FFI()


ffi.cdef('''
//typedef ... Treal; // damn, stuff like this doesn't work with cffi yet... structs are fine, but not the basic types :(
// so, we stick to 'double' here:
void cffti(int n, double wsave[]);
void cfftf(int n, double c[], double wsave[]);
''')


C = ffi.verify('''
    /*
    fftpack.c : A set of FFT routines in C.
    Algorithmically based on Fortran-77 FFTPACK by <NAME> (Version 4, 1985).

    */

    /* isign is +1 for backward and -1 for forward transforms */

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

    #ifdef DOUBLE
    #define Treal double
    #else
    #define Treal float
    #endif


    #define ref(u,a) u[a]

    #define MAXFAC 13    /* maximum number of factors in factorization of n */
    #define NSPECIAL 4   /* number of factors for which we have special-case routines */

    #ifdef __cplusplus
    extern "C" {
    #endif


    /* ----------------------------------------------------------------------
    passf2, passf3, passf4, passf5, passf. Complex FFT passes fwd and bwd.
    ---------------------------------------------------------------------- */

    static void passf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[], int isign)
    /* isign==+1 for backward transform */
    {
        int i, k, ah, ac;
        Treal ti2, tr2;
        if (ido <= 2) {
        for (k=0; k<l1; k++) {
            ah = k*ido;
            ac = 2*k*ido;
            ch[ah]              = ref(cc,ac) + ref(cc,ac + ido);
            ch[ah + ido*l1]     = ref(cc,ac) - ref(cc,ac + ido);
            ch[ah+1]            = ref(cc,ac+1) + ref(cc,ac + ido + 1);
            ch[ah + ido*l1 + 1] = ref(cc,ac+1) - ref(cc,ac + ido + 1);
        }
        } else {
        for (k=0; k<l1; k++) {
            for (i=0; i<ido-1; i+=2) {
            ah = i + k*ido;
            ac = i + 2*k*ido;
            ch[ah]   = ref(cc,ac) + ref(cc,ac + ido);
            tr2      = ref(cc,ac) - ref(cc,ac + ido);
            ch[ah+1] = ref(cc,ac+1) + ref(cc,ac + 1 + ido);
            ti2      = ref(cc,ac+1) - ref(cc,ac + 1 + ido);
            ch[ah+l1*ido+1] = wa1[i]*ti2 + isign*wa1[i+1]*tr2;
            ch[ah+l1*ido]   = wa1[i]*tr2 - isign*wa1[i+1]*ti2;
            }
        }
        }
    } /* passf2 */


    static void passf3(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[], int isign)
    /* isign==+1 for backward transform */
    {
        static const Treal taur = -0.5;
        static const Treal taui = 0.866025403784439;
        int i, k, ac, ah;
        Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
        if (ido == 2) {
        for (k=1; k<=l1; k++) {
            ac = (3*k - 2)*ido;
            tr2 = ref(cc,ac) + ref(cc,ac + ido);
            cr2 = ref(cc,ac - ido) + taur*tr2;
            ah = (k - 1)*ido;
            ch[ah] = ref(cc,ac - ido) + tr2;

            ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);
            ci2 = ref(cc,ac - ido + 1) + taur*ti2;
            ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;

            cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));
            ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));
            ch[ah + l1*ido] = cr2 - ci3;
            ch[ah + 2*l1*ido] = cr2 + ci3;
            ch[ah + l1*ido + 1] = ci2 + cr3;
            ch[ah + 2*l1*ido + 1] = ci2 - cr3;
        }
        } else {
        for (k=1; k<=l1; k++) {
            for (i=0; i<ido-1; i+=2) {
            ac = i + (3*k - 2)*ido;
            tr2 = ref(cc,ac) + ref(cc,ac + ido);
            cr2 = ref(cc,ac - ido) + taur*tr2;
            ah = i + (k-1)*ido;
            ch[ah] = ref(cc,ac - ido) + tr2;
            ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);
            ci2 = ref(cc,ac - ido + 1) + taur*ti2;
            ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;
            cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));
            ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));
            dr2 = cr2 - ci3;
            dr3 = cr2 + ci3;
            di2 = ci2 + cr3;
            di3 = ci2 - cr3;
            ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
            ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
            ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
            ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
            }
        }
        }
    } /* passf3 */


    static void passf4(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[], const Treal wa3[], int isign)
    /* isign == -1 for forward transform and +1 for backward transform */
    {
        int i, k, ac, ah;
        Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
        if (ido == 2) {
        for (k=0; k<l1; k++) {
            ac = 4*k*ido + 1;
            ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);
            ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);
            tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);
            ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);
            tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);
            tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);
            ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);
            tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);
            ah = k*ido;
            ch[ah] = tr2 + tr3;
            ch[ah + 2*l1*ido] = tr2 - tr3;
            ch[ah + 1] = ti2 + ti3;
            ch[ah + 2*l1*ido + 1] = ti2 - ti3;
            ch[ah + l1*ido] = tr1 + isign*tr4;
            ch[ah + 3*l1*ido] = tr1 - isign*tr4;
            ch[ah + l1*ido + 1] = ti1 + isign*ti4;
            ch[ah + 3*l1*ido + 1] = ti1 - isign*ti4;
        }
        } else {
        for (k=0; k<l1; k++) {
            for (i=0; i<ido-1; i+=2) {
            ac = i + 1 + 4*k*ido;
            ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);
            ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);
            ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);
            tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);
            tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);
            tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);
            ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);
            tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);
            ah = i + k*ido;
            ch[ah] = tr2 + tr3;
            cr3 = tr2 - tr3;
            ch[ah + 1] = ti2 + ti3;
            ci3 = ti2 - ti3;
            cr2 = tr1 + isign*tr4;
            cr4 = tr1 - isign*tr4;
            ci2 = ti1 + isign*ti4;
            ci4 = ti1 - isign*ti4;
            ch[ah + l1*ido] = wa1[i]*cr2 - isign*wa1[i + 1]*ci2;
            ch[ah + l1*ido + 1] = wa1[i]*ci2 + isign*wa1[i + 1]*cr2;
            ch[ah + 2*l1*ido] = wa2[i]*cr3 - isign*wa2[i + 1]*ci3;
            ch[ah + 2*l1*ido + 1] = wa2[i]*ci3 + isign*wa2[i + 1]*cr3;
            ch[ah + 3*l1*ido] = wa3[i]*cr4 -isign*wa3[i + 1]*ci4;
            ch[ah + 3*l1*ido + 1] = wa3[i]*ci4 + isign*wa3[i + 1]*cr4;
            }
        }
        }
    } /* passf4 */


    static void passf5(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[], int isign)
    /* isign == -1 for forward transform and +1 for backward transform */
    {
        static const Treal tr11 = 0.309016994374947;
        static const Treal ti11 = 0.951056516295154;
        static const Treal tr12 = -0.809016994374947;
        static const Treal ti12 = 0.587785252292473;
        int i, k, ac, ah;
        Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
            ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
        if (ido == 2) {
        for (k = 1; k <= l1; ++k) {
            ac = (5*k - 4)*ido + 1;
            ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);
            ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);
            ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);
            ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);
            tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);
            tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);
            tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);
            tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);
            ah = (k - 1)*ido;
            ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;
            ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;
            cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;
            ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;
            cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;
            ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;
            cr5 = isign*(ti11*tr5 + ti12*tr4);
            ci5 = isign*(ti11*ti5 + ti12*ti4);
            cr4 = isign*(ti12*tr5 - ti11*tr4);
            ci4 = isign*(ti12*ti5 - ti11*ti4);
            ch[ah + l1*ido] = cr2 - ci5;
            ch[ah + 4*l1*ido] = cr2 + ci5;
            ch[ah + l1*ido + 1] = ci2 + cr5;
            ch[ah + 2*l1*ido + 1] = ci3 + cr4;
            ch[ah + 2*l1*ido] = cr3 - ci4;
            ch[ah + 3*l1*ido] = cr3 + ci4;
            ch[ah + 3*l1*ido + 1] = ci3 - cr4;
            ch[ah + 4*l1*ido + 1] = ci2 - cr5;
        }
        } else {
        for (k=1; k<=l1; k++) {
            for (i=0; i<ido-1; i+=2) {
            ac = i + 1 + (k*5 - 4)*ido;
            ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);
            ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);
            ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);
            ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);
            tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);
            tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);
            tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);
            tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);
            ah = i + (k - 1)*ido;
            ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;
            ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;
            cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;

            ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;
            cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;

            ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;
            cr5 = isign*(ti11*tr5 + ti12*tr4);
            ci5 = isign*(ti11*ti5 + ti12*ti4);
            cr4 = isign*(ti12*tr5 - ti11*tr4);
            ci4 = isign*(ti12*ti5 - ti11*ti4);
            dr3 = cr3 - ci4;
            dr4 = cr3 + ci4;
            di3 = ci3 + cr4;
            di4 = ci3 - cr4;
            dr5 = cr2 + ci5;
            dr2 = cr2 - ci5;
            di5 = ci2 - cr5;
            di2 = ci2 + cr5;
            ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
            ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
            ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
            ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
            ch[ah + 3*l1*ido] = wa3[i]*dr4 - isign*wa3[i+1]*di4;
            ch[ah + 3*l1*ido + 1] = wa3[i]*di4 + isign*wa3[i+1]*dr4;
            ch[ah + 4*l1*ido] = wa4[i]*dr5 - isign*wa4[i+1]*di5;
            ch[ah + 4*l1*ido + 1] = wa4[i]*di5 + isign*wa4[i+1]*dr5;
            }
        }
        }
    } /* passf5 */


    static void passf(int *nac, int ido, int ip, int l1, int idl1,
        Treal cc[], Treal ch[],
        const Treal wa[], int isign)
    /* isign is -1 for forward transform and +1 for backward transform */
    {
        int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, idj, idl, inc,idp;
        Treal wai, war;

        idot = ido / 2;
        /* nt = ip*idl1;*/
        ipph = (ip + 1) / 2;
        idp = ip*ido;
        if (ido >= l1) {
        for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (k=0; k<l1; k++) {
            for (i=0; i<ido; i++) {
                ch[i + (k + j*l1)*ido] =
                    ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*ip)*ido);
                ch[i + (k + jc*l1)*ido] =
                    ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*ip)*ido);
            }
            }
        }
        for (k=0; k<l1; k++)
            for (i=0; i<ido; i++)
            ch[i + k*ido] = ref(cc,i + k*ip*ido);
        } else {
        for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (i=0; i<ido; i++) {
            for (k=0; k<l1; k++) {
                ch[i + (k + j*l1)*ido] = ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*
                    ip)*ido);
                ch[i + (k + jc*l1)*ido] = ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*
                    ip)*ido);
            }
            }
        }
        for (i=0; i<ido; i++)
            for (k=0; k<l1; k++)
            ch[i + k*ido] = ref(cc,i + k*ip*ido);
        }

        idl = 2 - ido;
        inc = 0;
        for (l=1; l<ipph; l++) {
        lc = ip - l;
        idl += ido;
        for (ik=0; ik<idl1; ik++) {
            cc[ik + l*idl1] = ch[ik] + wa[idl - 2]*ch[ik + idl1];
            cc[ik + lc*idl1] = isign*wa[idl-1]*ch[ik + (ip-1)*idl1];
        }
        idlj = idl;
        inc += ido;
        for (j=2; j<ipph; j++) {
            jc = ip - j;
            idlj += inc;
            if (idlj > idp) idlj -= idp;
            war = wa[idlj - 2];
            wai = wa[idlj-1];
            for (ik=0; ik<idl1; ik++) {
            cc[ik + l*idl1] += war*ch[ik + j*idl1];
            cc[ik + lc*idl1] += isign*wai*ch[ik + jc*idl1];
            }
        }
        }
        for (j=1; j<ipph; j++)
        for (ik=0; ik<idl1; ik++)
            ch[ik] += ch[ik + j*idl1];
        for (j=1; j<ipph; j++) {
        jc = ip - j;
        for (ik=1; ik<idl1; ik+=2) {
            ch[ik - 1 + j*idl1] = cc[ik - 1 + j*idl1] - cc[ik + jc*idl1];
            ch[ik - 1 + jc*idl1] = cc[ik - 1 + j*idl1] + cc[ik + jc*idl1];
            ch[ik + j*idl1] = cc[ik + j*idl1] + cc[ik - 1 + jc*idl1];
            ch[ik + jc*idl1] = cc[ik + j*idl1] - cc[ik - 1 + jc*idl1];
        }
        }
        *nac = 1;
        if (ido == 2) return;
        *nac = 0;
        for (ik=0; ik<idl1; ik++)
        cc[ik] = ch[ik];
        for (j=1; j<ip; j++) {
        for (k=0; k<l1; k++) {
            cc[(k + j*l1)*ido + 0] = ch[(k + j*l1)*ido + 0];
            cc[(k + j*l1)*ido + 1] = ch[(k + j*l1)*ido + 1];
        }
        }
        if (idot <= l1) {
        idij = 0;
        for (j=1; j<ip; j++) {
            idij += 2;
            for (i=3; i<ido; i+=2) {
            idij += 2;
            for (k=0; k<l1; k++) {
                cc[i - 1 + (k + j*l1)*ido] =
                    wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
                    isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
                cc[i + (k + j*l1)*ido] =
                    wa[idij - 2]*ch[i + (k + j*l1)*ido] +
                    isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
            }
            }
        }
        } else {
        idj = 2 - ido;
        for (j=1; j<ip; j++) {
            idj += ido;
            for (k = 0; k < l1; k++) {
            idij = idj;
            for (i=3; i<ido; i+=2) {
                idij += 2;
                cc[i - 1 + (k + j*l1)*ido] =
                    wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
                    isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
                cc[i + (k + j*l1)*ido] =
                    wa[idij - 2]*ch[i + (k + j*l1)*ido] +
                    isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
            }
            }
        }
        }
    } /* passf */


    /* ----------------------------------------------------------------------
    radf2,radb2, radf3,radb3, radf4,radb4, radf5,radb5, radfg,radbg.
    Treal FFT passes fwd and bwd.
    ---------------------------------------------------------------------- */

    static void radf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])
    {
        int i, k, ic;
        Treal ti2, tr2;
        for (k=0; k<l1; k++) {
        ch[2*k*ido] =
            ref(cc,k*ido) + ref(cc,(k + l1)*ido);
        ch[(2*k+1)*ido + ido-1] =
            ref(cc,k*ido) - ref(cc,(k + l1)*ido);
        }
        if (ido < 2) return;
        if (ido != 2) {
        for (k=0; k<l1; k++) {
            for (i=2; i<ido; i+=2) {
            ic = ido - i;
            tr2 = wa1[i - 2]*ref(cc, i-1 + (k + l1)*ido) + wa1[i - 1]*ref(cc, i + (k + l1)*ido);
            ti2 = wa1[i - 2]*ref(cc, i + (k + l1)*ido) - wa1[i - 1]*ref(cc, i-1 + (k + l1)*ido);
            ch[i + 2*k*ido] = ref(cc,i + k*ido) + ti2;
            ch[ic + (2*k+1)*ido] = ti2 - ref(cc,i + k*ido);
            ch[i - 1 + 2*k*ido] = ref(cc,i - 1 + k*ido) + tr2;
            ch[ic - 1 + (2*k+1)*ido] = ref(cc,i - 1 + k*ido) - tr2;
            }
        }
        if (ido % 2 == 1) return;
        }
        for (k=0; k<l1; k++) {
        ch[(2*k+1)*ido] = -ref(cc,ido-1 + (k + l1)*ido);
        ch[ido-1 + 2*k*ido] = ref(cc,ido-1 + k*ido);
        }
    } /* radf2 */


    static void radb2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])
    {
        int i, k, ic;
        Treal ti2, tr2;
        for (k=0; k<l1; k++) {
        ch[k*ido] =
            ref(cc,2*k*ido) + ref(cc,ido-1 + (2*k+1)*ido);
        ch[(k + l1)*ido] =
            ref(cc,2*k*ido) - ref(cc,ido-1 + (2*k+1)*ido);
        }
        if (ido < 2) return;
        if (ido != 2) {
        for (k = 0; k < l1; ++k) {
            for (i = 2; i < ido; i += 2) {
            ic = ido - i;
            ch[i-1 + k*ido] =
                ref(cc,i-1 + 2*k*ido) + ref(cc,ic-1 + (2*k+1)*ido);
            tr2 = ref(cc,i-1 + 2*k*ido) - ref(cc,ic-1 + (2*k+1)*ido);
            ch[i + k*ido] =
                ref(cc,i + 2*k*ido) - ref(cc,ic + (2*k+1)*ido);
            ti2 = ref(cc,i + (2*k)*ido) + ref(cc,ic + (2*k+1)*ido);
            ch[i-1 + (k + l1)*ido] =
                wa1[i - 2]*tr2 - wa1[i - 1]*ti2;
            ch[i + (k + l1)*ido] =
                wa1[i - 2]*ti2 + wa1[i - 1]*tr2;
            }
        }
        if (ido % 2 == 1) return;
        }
        for (k = 0; k < l1; k++) {
        ch[ido-1 + k*ido] = 2*ref(cc,ido-1 + 2*k*ido);
        ch[ido-1 + (k + l1)*ido] = -2*ref(cc,(2*k+1)*ido);
        }
    } /* radb2 */


    static void radf3(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[])
    {
        static const Treal taur = -0.5;
        static const Treal taui = 0.866025403784439;
        int i, k, ic;
        Treal ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
        for (k=0; k<l1; k++) {
        cr2 = ref(cc,(k + l1)*ido) + ref(cc,(k + 2*l1)*ido);
        ch[3*k*ido] = ref(cc,k*ido) + cr2;
        ch[(3*k+2)*ido] = taui*(ref(cc,(k + l1*2)*ido) - ref(cc,(k + l1)*ido));
        ch[ido-1 + (3*k + 1)*ido] = ref(cc,k*ido) + taur*cr2;
        }
        if (ido == 1) return;
        for (k=0; k<l1; k++) {
        for (i=2; i<ido; i+=2) {
            ic = ido - i;
            dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) +
                wa1[i - 1]*ref(cc,i + (k + l1)*ido);
            di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
            dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + l1*2)*ido) + wa2[i - 1]*ref(cc,i + (k + l1*2)*ido);
            di3 = wa2[i - 2]*ref(cc,i + (k + l1*2)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + l1*2)*ido);
            cr2 = dr2 + dr3;
            ci2 = di2 + di3;
            ch[i - 1 + 3*k*ido] = ref(cc,i - 1 + k*ido) + cr2;
            ch[i + 3*k*ido] = ref(cc,i + k*ido) + ci2;
            tr2 = ref(cc,i - 1 + k*ido) + taur*cr2;
            ti2 = ref(cc,i + k*ido) + taur*ci2;
            tr3 = taui*(di2 - di3);
            ti3 = taui*(dr3 - dr2);
            ch[i - 1 + (3*k + 2)*ido] = tr2 + tr3;
            ch[ic - 1 + (3*k + 1)*ido] = tr2 - tr3;
            ch[i + (3*k + 2)*ido] = ti2 + ti3;
            ch[ic + (3*k + 1)*ido] = ti3 - ti2;
        }
        }
    } /* radf3 */


    static void radb3(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[])
    {
        static const Treal taur = -0.5;
        static const Treal taui = 0.866025403784439;
        int i, k, ic;
        Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
        for (k=0; k<l1; k++) {
        tr2 = 2*ref(cc,ido-1 + (3*k + 1)*ido);
        cr2 = ref(cc,3*k*ido) + taur*tr2;
        ch[k*ido] = ref(cc,3*k*ido) + tr2;
        ci3 = 2*taui*ref(cc,(3*k + 2)*ido);
        ch[(k + l1)*ido] = cr2 - ci3;
        ch[(k + 2*l1)*ido] = cr2 + ci3;
        }
        if (ido == 1) return;
        for (k=0; k<l1; k++) {
        for (i=2; i<ido; i+=2) {
            ic = ido - i;
            tr2 = ref(cc,i - 1 + (3*k + 2)*ido) + ref(cc,ic - 1 + (3*k + 1)*ido);
            cr2 = ref(cc,i - 1 + 3*k*ido) + taur*tr2;
            ch[i - 1 + k*ido] = ref(cc,i - 1 + 3*k*ido) + tr2;
            ti2 = ref(cc,i + (3*k + 2)*ido) - ref(cc,ic + (3*k + 1)*ido);
            ci2 = ref(cc,i + 3*k*ido) + taur*ti2;
            ch[i + k*ido] = ref(cc,i + 3*k*ido) + ti2;
            cr3 = taui*(ref(cc,i - 1 + (3*k + 2)*ido) - ref(cc,ic - 1 + (3*k + 1)*ido));
            ci3 = taui*(ref(cc,i + (3*k + 2)*ido) + ref(cc,ic + (3*k + 1)*ido));
            dr2 = cr2 - ci3;
            dr3 = cr2 + ci3;
            di2 = ci2 + cr3;
            di3 = ci2 - cr3;
            ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
            ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
            ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
            ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
        }
        }
    } /* radb3 */


    static void radf4(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[], const Treal wa3[])
    {
        static const Treal hsqt2 = 0.7071067811865475;
        int i, k, ic;
        Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
        for (k=0; k<l1; k++) {
        tr1 = ref(cc,(k + l1)*ido) + ref(cc,(k + 3*l1)*ido);
        tr2 = ref(cc,k*ido) + ref(cc,(k + 2*l1)*ido);
        ch[4*k*ido] = tr1 + tr2;
        ch[ido-1 + (4*k + 3)*ido] = tr2 - tr1;
        ch[ido-1 + (4*k + 1)*ido] = ref(cc,k*ido) - ref(cc,(k + 2*l1)*ido);
        ch[(4*k + 2)*ido] = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + l1)*ido);
        }
        if (ido < 2) return;
        if (ido != 2) {
        for (k=0; k<l1; k++) {
            for (i=2; i<ido; i += 2) {
            ic = ido - i;
            cr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
            ci2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
            cr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*
                ido);
            ci3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*
                ido);
            cr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*
                ido);
            ci4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*
                ido);
            tr1 = cr2 + cr4;
            tr4 = cr4 - cr2;
            ti1 = ci2 + ci4;
            ti4 = ci2 - ci4;
            ti2 = ref(cc,i + k*ido) + ci3;
            ti3 = ref(cc,i + k*ido) - ci3;
            tr2 = ref(cc,i - 1 + k*ido) + cr3;
            tr3 = ref(cc,i - 1 + k*ido) - cr3;
            ch[i - 1 + 4*k*ido] = tr1 + tr2;
            ch[ic - 1 + (4*k + 3)*ido] = tr2 - tr1;
            ch[i + 4*k*ido] = ti1 + ti2;
            ch[ic + (4*k + 3)*ido] = ti1 - ti2;
            ch[i - 1 + (4*k + 2)*ido] = ti4 + tr3;
            ch[ic - 1 + (4*k + 1)*ido] = tr3 - ti4;
            ch[i + (4*k + 2)*ido] = tr4 + ti3;
            ch[ic + (4*k + 1)*ido] = tr4 - ti3;
            }
        }
        if (ido % 2 == 1) return;
        }
        for (k=0; k<l1; k++) {
        ti1 = -hsqt2*(ref(cc,ido-1 + (k + l1)*ido) + ref(cc,ido-1 + (k + 3*l1)*ido));
        tr1 = hsqt2*(ref(cc,ido-1 + (k + l1)*ido) - ref(cc,ido-1 + (k + 3*l1)*ido));
        ch[ido-1 + 4*k*ido] = tr1 + ref(cc,ido-1 + k*ido);
        ch[ido-1 + (4*k + 2)*ido] = ref(cc,ido-1 + k*ido) - tr1;
        ch[(4*k + 1)*ido] = ti1 - ref(cc,ido-1 + (k + 2*l1)*ido);
        ch[(4*k + 3)*ido] = ti1 + ref(cc,ido-1 + (k + 2*l1)*ido);
        }
    } /* radf4 */


    static void radb4(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[], const Treal wa3[])
    {
        static const Treal sqrt2 = 1.414213562373095;
        int i, k, ic;
        Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
        for (k = 0; k < l1; k++) {
        tr1 = ref(cc,4*k*ido) - ref(cc,ido-1 + (4*k + 3)*ido);
        tr2 = ref(cc,4*k*ido) + ref(cc,ido-1 + (4*k + 3)*ido);
        tr3 = ref(cc,ido-1 + (4*k + 1)*ido) + ref(cc,ido-1 + (4*k + 1)*ido);
        tr4 = ref(cc,(4*k + 2)*ido) + ref(cc,(4*k + 2)*ido);
        ch[k*ido] = tr2 + tr3;
        ch[(k + l1)*ido] = tr1 - tr4;
        ch[(k + 2*l1)*ido] = tr2 - tr3;
        ch[(k + 3*l1)*ido] = tr1 + tr4;
        }
        if (ido < 2) return;
        if (ido != 2) {
        for (k = 0; k < l1; ++k) {
            for (i = 2; i < ido; i += 2) {
            ic = ido - i;
            ti1 = ref(cc,i + 4*k*ido) + ref(cc,ic + (4*k + 3)*ido);
            ti2 = ref(cc,i + 4*k*ido) - ref(cc,ic + (4*k + 3)*ido);
            ti3 = ref(cc,i + (4*k + 2)*ido) - ref(cc,ic + (4*k + 1)*ido);
            tr4 = ref(cc,i + (4*k + 2)*ido) + ref(cc,ic + (4*k + 1)*ido);
            tr1 = ref(cc,i - 1 + 4*k*ido) - ref(cc,ic - 1 + (4*k + 3)*ido);
            tr2 = ref(cc,i - 1 + 4*k*ido) + ref(cc,ic - 1 + (4*k + 3)*ido);
            ti4 = ref(cc,i - 1 + (4*k + 2)*ido) - ref(cc,ic - 1 + (4*k + 1)*ido);
            tr3 = ref(cc,i - 1 + (4*k + 2)*ido) + ref(cc,ic - 1 + (4*k + 1)*ido);
            ch[i - 1 + k*ido] = tr2 + tr3;
            cr3 = tr2 - tr3;
            ch[i + k*ido] = ti2 + ti3;
            ci3 = ti2 - ti3;
            cr2 = tr1 - tr4;
            cr4 = tr1 + tr4;
            ci2 = ti1 + ti4;
            ci4 = ti1 - ti4;
            ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*cr2 - wa1[i - 1]*ci2;
            ch[i + (k + l1)*ido] = wa1[i - 2]*ci2 + wa1[i - 1]*cr2;
            ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*cr3 - wa2[i - 1]*ci3;
            ch[i + (k + 2*l1)*ido] = wa2[i - 2]*ci3 + wa2[i - 1]*cr3;
            ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*cr4 - wa3[i - 1]*ci4;
            ch[i + (k + 3*l1)*ido] = wa3[i - 2]*ci4 + wa3[i - 1]*cr4;
            }
        }
        if (ido % 2 == 1) return;
        }
        for (k = 0; k < l1; k++) {
        ti1 = ref(cc,(4*k + 1)*ido) + ref(cc,(4*k + 3)*ido);
        ti2 = ref(cc,(4*k + 3)*ido) - ref(cc,(4*k + 1)*ido);
        tr1 = ref(cc,ido-1 + 4*k*ido) - ref(cc,ido-1 + (4*k + 2)*ido);
        tr2 = ref(cc,ido-1 + 4*k*ido) + ref(cc,ido-1 + (4*k + 2)*ido);
        ch[ido-1 + k*ido] = tr2 + tr2;
        ch[ido-1 + (k + l1)*ido] = sqrt2*(tr1 - ti1);
        ch[ido-1 + (k + 2*l1)*ido] = ti2 + ti2;
        ch[ido-1 + (k + 3*l1)*ido] = -sqrt2*(tr1 + ti1);
        }
    } /* radb4 */


    static void radf5(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])
    {
        static const Treal tr11 = 0.309016994374947;
        static const Treal ti11 = 0.951056516295154;
        static const Treal tr12 = -0.809016994374947;
        static const Treal ti12 = 0.587785252292473;
        int i, k, ic;
        Treal ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5,
            cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
        for (k = 0; k < l1; k++) {
        cr2 = ref(cc,(k + 4*l1)*ido) + ref(cc,(k + l1)*ido);
        ci5 = ref(cc,(k + 4*l1)*ido) - ref(cc,(k + l1)*ido);
        cr3 = ref(cc,(k + 3*l1)*ido) + ref(cc,(k + 2*l1)*ido);
        ci4 = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + 2*l1)*ido);
        ch[5*k*ido] = ref(cc,k*ido) + cr2 + cr3;
        ch[ido-1 + (5*k + 1)*ido] = ref(cc,k*ido) + tr11*cr2 + tr12*cr3;
        ch[(5*k + 2)*ido] = ti11*ci5 + ti12*ci4;
        ch[ido-1 + (5*k + 3)*ido] = ref(cc,k*ido) + tr12*cr2 + tr11*cr3;
        ch[(5*k + 4)*ido] = ti12*ci5 - ti11*ci4;
        }
        if (ido == 1) return;
        for (k = 0; k < l1; ++k) {
        for (i = 2; i < ido; i += 2) {
            ic = ido - i;
            dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
            di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
            dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*ido);
            di3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*ido);
            dr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*ido);
            di4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*ido);
            dr5 = wa4[i - 2]*ref(cc,i - 1 + (k + 4*l1)*ido) + wa4[i - 1]*ref(cc,i + (k + 4*l1)*ido);
            di5 = wa4[i - 2]*ref(cc,i + (k + 4*l1)*ido) - wa4[i - 1]*ref(cc,i - 1 + (k + 4*l1)*ido);
            cr2 = dr2 + dr5;
            ci5 = dr5 - dr2;
            cr5 = di2 - di5;
            ci2 = di2 + di5;
            cr3 = dr3 + dr4;
            ci4 = dr4 - dr3;
            cr4 = di3 - di4;
            ci3 = di3 + di4;
            ch[i - 1 + 5*k*ido] = ref(cc,i - 1 + k*ido) + cr2 + cr3;
            ch[i + 5*k*ido] = ref(cc,i + k*ido) + ci2 + ci3;
            tr2 = ref(cc,i - 1 + k*ido) + tr11*cr2 + tr12*cr3;
            ti2 = ref(cc,i + k*ido) + tr11*ci2 + tr12*ci3;
            tr3 = ref(cc,i - 1 + k*ido) + tr12*cr2 + tr11*cr3;
            ti3 = ref(cc,i + k*ido) + tr12*ci2 + tr11*ci3;
            tr5 = ti11*cr5 + ti12*cr4;
            ti5 = ti11*ci5 + ti12*ci4;
            tr4 = ti12*cr5 - ti11*cr4;
            ti4 = ti12*ci5 - ti11*ci4;
            ch[i - 1 + (5*k + 2)*ido] = tr2 + tr5;
            ch[ic - 1 + (5*k + 1)*ido] = tr2 - tr5;
            ch[i + (5*k + 2)*ido] = ti2 + ti5;
            ch[ic + (5*k + 1)*ido] = ti5 - ti2;
            ch[i - 1 + (5*k + 4)*ido] = tr3 + tr4;
            ch[ic - 1 + (5*k + 3)*ido] = tr3 - tr4;
            ch[i + (5*k + 4)*ido] = ti3 + ti4;
            ch[ic + (5*k + 3)*ido] = ti4 - ti3;
        }
        }
    } /* radf5 */


    static void radb5(int ido, int l1, const Treal cc[], Treal ch[],
        const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])
    {
        static const Treal tr11 = 0.309016994374947;
        static const Treal ti11 = 0.951056516295154;
        static const Treal tr12 = -0.809016994374947;
        static const Treal ti12 = 0.587785252292473;
        int i, k, ic;
        Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
            ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
        for (k = 0; k < l1; k++) {
        ti5 = 2*ref(cc,(5*k + 2)*ido);
        ti4 = 2*ref(cc,(5*k + 4)*ido);
        tr2 = 2*ref(cc,ido-1 + (5*k + 1)*ido);
        tr3 = 2*ref(cc,ido-1 + (5*k + 3)*ido);
        ch[k*ido] = ref(cc,5*k*ido) + tr2 + tr3;
        cr2 = ref(cc,5*k*ido) + tr11*tr2 + tr12*tr3;
        cr3 = ref(cc,5*k*ido) + tr12*tr2 + tr11*tr3;
        ci5 = ti11*ti5 + ti12*ti4;
        ci4 = ti12*ti5 - ti11*ti4;
        ch[(k + l1)*ido] = cr2 - ci5;
        ch[(k + 2*l1)*ido] = cr3 - ci4;
        ch[(k + 3*l1)*ido] = cr3 + ci4;
        ch[(k + 4*l1)*ido] = cr2 + ci5;
        }
        if (ido == 1) return;
        for (k = 0; k < l1; ++k) {
        for (i = 2; i < ido; i += 2) {
            ic = ido - i;
            ti5 = ref(cc,i + (5*k + 2)*ido) + ref(cc,ic + (5*k + 1)*ido);
            ti2 = ref(cc,i + (5*k + 2)*ido) - ref(cc,ic + (5*k + 1)*ido);
            ti4 = ref(cc,i + (5*k + 4)*ido) + ref(cc,ic + (5*k + 3)*ido);
            ti3 = ref(cc,i + (5*k + 4)*ido) - ref(cc,ic + (5*k + 3)*ido);
            tr5 = ref(cc,i - 1 + (5*k + 2)*ido) - ref(cc,ic - 1 + (5*k + 1)*ido);
            tr2 = ref(cc,i - 1 + (5*k + 2)*ido) + ref(cc,ic - 1 + (5*k + 1)*ido);
            tr4 = ref(cc,i - 1 + (5*k + 4)*ido) - ref(cc,ic - 1 + (5*k + 3)*ido);
            tr3 = ref(cc,i - 1 + (5*k + 4)*ido) + ref(cc,ic - 1 + (5*k + 3)*ido);
            ch[i - 1 + k*ido] = ref(cc,i - 1 + 5*k*ido) + tr2 + tr3;
            ch[i + k*ido] = ref(cc,i + 5*k*ido) + ti2 + ti3;
            cr2 = ref(cc,i - 1 + 5*k*ido) + tr11*tr2 + tr12*tr3;

            ci2 = ref(cc,i + 5*k*ido) + tr11*ti2 + tr12*ti3;
            cr3 = ref(cc,i - 1 + 5*k*ido) + tr12*tr2 + tr11*tr3;

            ci3 = ref(cc,i + 5*k*ido) + tr12*ti2 + tr11*ti3;
            cr5 = ti11*tr5 + ti12*tr4;
            ci5 = ti11*ti5 + ti12*ti4;
            cr4 = ti12*tr5 - ti11*tr4;
            ci4 = ti12*ti5 - ti11*ti4;
            dr3 = cr3 - ci4;
            dr4 = cr3 + ci4;
            di3 = ci3 + cr4;
            di4 = ci3 - cr4;
            dr5 = cr2 + ci5;
            dr2 = cr2 - ci5;
            di5 = ci2 - cr5;
            di2 = ci2 + cr5;
            ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
            ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
            ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
            ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
            ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*dr4 - wa3[i - 1]*di4;
            ch[i + (k + 3*l1)*ido] = wa3[i - 2]*di4 + wa3[i - 1]*dr4;
            ch[i - 1 + (k + 4*l1)*ido] = wa4[i - 2]*dr5 - wa4[i - 1]*di5;
            ch[i + (k + 4*l1)*ido] = wa4[i - 2]*di5 + wa4[i - 1]*dr5;
        }
        }
    } /* radb5 */


    static void radfg(int ido, int ip, int l1, int idl1,
        Treal cc[], Treal ch[], const Treal wa[])
    {
        static const Treal twopi = 6.28318530717959;
        int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;
        Treal dc2, ai1, ai2, ar1, ar2, ds2, dcp, arg, dsp, ar1h, ar2h;
        arg = twopi / ip;
        dcp = cos(arg);
        dsp = sin(arg);
        ipph = (ip + 1) / 2;
        nbd = (ido - 1) / 2;
        if (ido != 1) {
        for (ik=0; ik<idl1; ik++) ch[ik] = cc[ik];
        for (j=1; j<ip; j++)
            for (k=0; k<l1; k++)
            ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido];
        if (nbd <= l1) {
            is = -ido;
            for (j=1; j<ip; j++) {
            is += ido;
            idij = is-1;
            for (i=2; i<ido; i+=2) {
                idij += 2;
                for (k=0; k<l1; k++) {
                ch[i - 1 + (k + j*l1)*ido] =
                    wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
                ch[i + (k + j*l1)*ido] =
                    wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
                }
            }
            }
        } else {
            is = -ido;
            for (j=1; j<ip; j++) {
            is += ido;
            for (k=0; k<l1; k++) {
                idij = is-1;
                for (i=2; i<ido; i+=2) {
                idij += 2;
                ch[i - 1 + (k + j*l1)*ido] =
                    wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
                ch[i + (k + j*l1)*ido] =
                    wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
                }
            }
            }
        }
        if (nbd >= l1) {
            for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (k=0; k<l1; k++) {
                for (i=2; i<ido; i+=2) {
                cc[i - 1 + (k + j*l1)*ido] = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
                cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
                cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
                cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
                }
            }
            }
        } else {
            for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (i=2; i<ido; i+=2) {
                for (k=0; k<l1; k++) {
                cc[i - 1 + (k + j*l1)*ido] =
                    ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
                cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
                cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
                cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
                }
            }
            }
        }
        } else {  /* now ido == 1 */
        for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];
        }
        for (j=1; j<ipph; j++) {
        jc = ip - j;
        for (k=0; k<l1; k++) {
            cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido] + ch[(k + jc*l1)*ido];
            cc[(k + jc*l1)*ido] = ch[(k + jc*l1)*ido] - ch[(k + j*l1)*ido];
        }
        }

        ar1 = 1;
        ai1 = 0;
        for (l=1; l<ipph; l++) {
        lc = ip - l;
        ar1h = dcp*ar1 - dsp*ai1;
        ai1 = dcp*ai1 + dsp*ar1;
        ar1 = ar1h;
        for (ik=0; ik<idl1; ik++) {
            ch[ik + l*idl1] = cc[ik] + ar1*cc[ik + idl1];
            ch[ik + lc*idl1] = ai1*cc[ik + (ip-1)*idl1];
        }
        dc2 = ar1;
        ds2 = ai1;
        ar2 = ar1;
        ai2 = ai1;
        for (j=2; j<ipph; j++) {
            jc = ip - j;
            ar2h = dc2*ar2 - ds2*ai2;
            ai2 = dc2*ai2 + ds2*ar2;
            ar2 = ar2h;
            for (ik=0; ik<idl1; ik++) {
            ch[ik + l*idl1] += ar2*cc[ik + j*idl1];
            ch[ik + lc*idl1] += ai2*cc[ik + jc*idl1];
            }
        }
        }
        for (j=1; j<ipph; j++)
        for (ik=0; ik<idl1; ik++)
            ch[ik] += cc[ik + j*idl1];

        if (ido >= l1) {
        for (k=0; k<l1; k++) {
            for (i=0; i<ido; i++) {
            ref(cc,i + k*ip*ido) = ch[i + k*ido];
            }
        }
        } else {
        for (i=0; i<ido; i++) {
            for (k=0; k<l1; k++) {
            ref(cc,i + k*ip*ido) = ch[i + k*ido];
            }
        }
        }
        for (j=1; j<ipph; j++) {
        jc = ip - j;
        j2 = 2*j;
        for (k=0; k<l1; k++) {
            ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) =
                ch[(k + j*l1)*ido];
            ref(cc,(j2 + k*ip)*ido) =
                ch[(k + jc*l1)*ido];
        }
        }
        if (ido == 1) return;
        if (nbd >= l1) {
        for (j=1; j<ipph; j++) {
            jc = ip - j;
            j2 = 2*j;
            for (k=0; k<l1; k++) {
            for (i=2; i<ido; i+=2) {
                ic = ido - i;
                ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
                ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
                ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
                ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
            }
            }
        }
        } else {
        for (j=1; j<ipph; j++) {
            jc = ip - j;
            j2 = 2*j;
            for (i=2; i<ido; i+=2) {
            ic = ido - i;
            for (k=0; k<l1; k++) {
                ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
                ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
                ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
                ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
            }
            }
        }
        }
    } /* radfg */


    static void radbg(int ido, int ip, int l1, int idl1,
        Treal cc[], Treal ch[], const Treal wa[])
    {
        static const Treal twopi = 6.28318530717959;
        int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
        Treal dc2, ai1, ai2, ar1, ar2, ds2;
        int nbd;
        Treal dcp, arg, dsp, ar1h, ar2h;
        arg = twopi / ip;
        dcp = cos(arg);
        dsp = sin(arg);
        nbd = (ido - 1) / 2;
        ipph = (ip + 1) / 2;
        if (ido >= l1) {
        for (k=0; k<l1; k++) {
            for (i=0; i<ido; i++) {
            ch[i + k*ido] = ref(cc,i + k*ip*ido);
            }
        }
        } else {
        for (i=0; i<ido; i++) {
            for (k=0; k<l1; k++) {
            ch[i + k*ido] = ref(cc,i + k*ip*ido);
            }
        }
        }
        for (j=1; j<ipph; j++) {
        jc = ip - j;
        j2 = 2*j;
        for (k=0; k<l1; k++) {
            ch[(k + j*l1)*ido] = ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) + ref(cc,ido-1 + (j2 - 1 + k*ip)*
                ido);
            ch[(k + jc*l1)*ido] = ref(cc,(j2 + k*ip)*ido) + ref(cc,(j2 + k*ip)*ido);
        }
        }

        if (ido != 1) {
        if (nbd >= l1) {
            for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (k=0; k<l1; k++) {
                for (i=2; i<ido; i+=2) {
                ic = ido - i;
                ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,
                    ic - 1 + (2*j - 1 + k*ip)*ido);
                ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -
                    ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
                ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic
                    + (2*j - 1 + k*ip)*ido);
                ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic
                    + (2*j - 1 + k*ip)*ido);
                }
            }
            }
        } else {
            for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (i=2; i<ido; i+=2) {
                ic = ido - i;
                for (k=0; k<l1; k++) {
                ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,
                    ic - 1 + (2*j - 1 + k*ip)*ido);
                ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -
                    ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
                ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic
                    + (2*j - 1 + k*ip)*ido);
                ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic
                    + (2*j - 1 + k*ip)*ido);
                }
            }
            }
        }
        }

        ar1 = 1;
        ai1 = 0;
        for (l=1; l<ipph; l++) {
        lc = ip - l;
        ar1h = dcp*ar1 - dsp*ai1;
        ai1 = dcp*ai1 + dsp*ar1;
        ar1 = ar1h;
        for (ik=0; ik<idl1; ik++) {
            cc[ik + l*idl1] = ch[ik] + ar1*ch[ik + idl1];
            cc[ik + lc*idl1] = ai1*ch[ik + (ip-1)*idl1];
        }
        dc2 = ar1;
        ds2 = ai1;
        ar2 = ar1;
        ai2 = ai1;
        for (j=2; j<ipph; j++) {
            jc = ip - j;
            ar2h = dc2*ar2 - ds2*ai2;
            ai2 = dc2*ai2 + ds2*ar2;
            ar2 = ar2h;
            for (ik=0; ik<idl1; ik++) {
            cc[ik + l*idl1] += ar2*ch[ik + j*idl1];
            cc[ik + lc*idl1] += ai2*ch[ik + jc*idl1];
            }
        }
        }
        for (j=1; j<ipph; j++) {
        for (ik=0; ik<idl1; ik++) {
            ch[ik] += ch[ik + j*idl1];
        }
        }
        for (j=1; j<ipph; j++) {
        jc = ip - j;
        for (k=0; k<l1; k++) {
            ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido] - cc[(k + jc*l1)*ido];
            ch[(k + jc*l1)*ido] = cc[(k + j*l1)*ido] + cc[(k + jc*l1)*ido];
        }
        }

        if (ido == 1) return;
        if (nbd >= l1) {
        for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (k=0; k<l1; k++) {
            for (i=2; i<ido; i+=2) {
                ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
                ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] + cc[i + (k + jc*l1)*ido];
                ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
                ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
            }
            }
        }
        } else {
        for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (i=2; i<ido; i+=2) {
            for (k=0; k<l1; k++) {
                ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
                ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j *l1)*ido] + cc[i + (k + jc*l1)*ido];
                ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
                ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
            }
            }
        }
        }
        for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];
        for (j=1; j<ip; j++)
        for (k=0; k<l1; k++)
            cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido];
        if (nbd <= l1) {
        is = -ido;
        for (j=1; j<ip; j++) {
            is += ido;
            idij = is-1;
            for (i=2; i<ido; i+=2) {
            idij += 2;
            for (k=0; k<l1; k++) {
                cc[i - 1 + (k + j*l1)*ido] = wa[idij - 1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*
                    ch[i + (k + j*l1)*ido];
                cc[i + (k + j*l1)*ido] = wa[idij - 1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
            }
            }
        }
        } else {
        is = -ido;
        for (j=1; j<ip; j++) {
            is += ido;
            for (k=0; k<l1; k++) {
            idij = is - 1;
            for (i=2; i<ido; i+=2) {
                idij += 2;
                cc[i - 1 + (k + j*l1)*ido] = wa[idij-1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*
                    ch[i + (k + j*l1)*ido];
                cc[i + (k + j*l1)*ido] = wa[idij-1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
            }
            }
        }
        }
    } /* radbg */

    /* ----------------------------------------------------------------------
    cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
    ---------------------------------------------------------------------- */

    static void cfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2], int isign)
    {
        int idot, i;
        int k1, l1, l2;
        int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;
        Treal *cinput, *coutput;
        nf = ifac[1];
        na = 0;
        l1 = 1;
        iw = 0;
        for (k1=2; k1<=nf+1; k1++) {
        ip = ifac[k1];
        l2 = ip*l1;
        ido = n / l2;
        idot = ido + ido;
        idl1 = idot*l1;
        if (na) {
            cinput = ch;
            coutput = c;
        } else {
            cinput = c;
            coutput = ch;
        }
        switch (ip) {
        case 4:
            ix2 = iw + idot;
            ix3 = ix2 + idot;
            passf4(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], isign);
            na = !na;
            break;
        case 2:
            passf2(idot, l1, cinput, coutput, &wa[iw], isign);
            na = !na;
            break;
        case 3:
            ix2 = iw + idot;
            passf3(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], isign);
            na = !na;
            break;
        case 5:
            ix2 = iw + idot;
            ix3 = ix2 + idot;
            ix4 = ix3 + idot;
            passf5(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
            na = !na;
            break;
        default:
            passf(&nac, idot, ip, l1, idl1, cinput, coutput, &wa[iw], isign);
            if (nac != 0) na = !na;
        }
        l1 = l2;
        iw += (ip - 1)*idot;
        }
        if (na == 0) return;
        for (i=0; i<2*n; i++) c[i] = ch[i];
    } /* cfftf1 */


    void cfftf(int n, Treal c[], Treal wsave[])
    {
        int iw1, iw2;
        if (n == 1) return;
        iw1 = 2*n;
        iw2 = iw1 + 2*n;
        cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), -1);
    } /* cfftf */


    void cfftb(int n, Treal c[], Treal wsave[])
    {
        int iw1, iw2;
        if (n == 1) return;
        iw1 = 2*n;
        iw2 = iw1 + 2*n;
        cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), +1);
    } /* cfftb */


    static void factorize(int n, int ifac[MAXFAC+2], const int ntryh[NSPECIAL])
    /* Factorize n in factors in ntryh and rest. On exit,
    ifac[0] contains n and ifac[1] contains number of factors,
    the factors start from ifac[2]. */
    {
        int ntry=3, i, j=0, ib, nf=0, nl=n, nq, nr;
    startloop:
        if (j < NSPECIAL)
        ntry = ntryh[j];
        else
        ntry+= 2;
        j++;
        do {
        nq = nl / ntry;
        nr = nl - ntry*nq;
        if (nr != 0) goto startloop;
        nf++;
        ifac[nf + 1] = ntry;
        nl = nq;
        if (ntry == 2 && nf != 1) {
            for (i=2; i<=nf; i++) {
            ib = nf - i + 2;
            ifac[ib + 1] = ifac[ib];
            }
            ifac[2] = 2;
        }
        } while (nl != 1);
        ifac[0] = n;
        ifac[1] = nf;
    }


    static void cffti1(int n, Treal wa[], int ifac[MAXFAC+2])
    {
        static const Treal twopi = 6.28318530717959;
        Treal arg, argh, argld, fi;
        int idot, i, j;
        int i1, k1, l1, l2;
        int ld, ii, nf, ip;
        int ido, ipm;

        static const int ntryh[NSPECIAL] = {
        3,4,2,5    }; /* Do not change the order of these. */

        factorize(n,ifac,ntryh);
        nf = ifac[1];
        argh = twopi/(Treal)n;
        i = 1;
        l1 = 1;
        for (k1=1; k1<=nf; k1++) {
        ip = ifac[k1+1];
        ld = 0;
        l2 = l1*ip;
        ido = n / l2;
        idot = ido + ido + 2;
        ipm = ip - 1;
        for (j=1; j<=ipm; j++) {
            i1 = i;
            wa[i-1] = 1;
            wa[i] = 0;
            ld += l1;
            fi = 0;
            argld = ld*argh;
            for (ii=4; ii<=idot; ii+=2) {
            i+= 2;
            fi+= 1;
            arg = fi*argld;
            wa[i-1] = cos(arg);
            wa[i] = sin(arg);
            }
            if (ip > 5) {
            wa[i1-1] = wa[i-1];
            wa[i1] = wa[i];
            }
        }
        l1 = l2;
        }
    } /* cffti1 */


    void cffti(int n, Treal wsave[])
    {
        int iw1, iw2;
        if (n == 1) return;
        iw1 = 2*n;
        iw2 = iw1 + 2*n;
        cffti1(n, wsave+iw1, (int*)(wsave+iw2));
    } /* cffti */

    /* ----------------------------------------------------------------------
    rfftf1, rfftb1, rfftf, rfftb, rffti1, rffti. Treal FFTs.
    ---------------------------------------------------------------------- */

    static void rfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])
    {
        int i;
        int k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
        Treal *cinput, *coutput;
        nf = ifac[1];
        na = 1;
        l2 = n;
        iw = n-1;
        for (k1 = 1; k1 <= nf; ++k1) {
        kh = nf - k1;
        ip = ifac[kh + 2];
        l1 = l2 / ip;
        ido = n / l2;
        idl1 = ido*l1;
        iw -= (ip - 1)*ido;
        na = !na;
        if (na) {
            cinput = ch;
            coutput = c;
        } else {
            cinput = c;
            coutput = ch;
        }
        switch (ip) {
        case 4:
            ix2 = iw + ido;
            ix3 = ix2 + ido;
            radf4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
            break;
        case 2:
            radf2(ido, l1, cinput, coutput, &wa[iw]);
            break;
        case 3:
            ix2 = iw + ido;
            radf3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
            break;
        case 5:
            ix2 = iw + ido;
            ix3 = ix2 + ido;
            ix4 = ix3 + ido;
            radf5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
            break;
        default:
            if (ido == 1)
            na = !na;
            if (na == 0) {
            radfg(ido, ip, l1, idl1, c, ch, &wa[iw]);
            na = 1;
            } else {
            radfg(ido, ip, l1, idl1, ch, c, &wa[iw]);
            na = 0;
            }
        }
        l2 = l1;
        }
        if (na == 1) return;
        for (i = 0; i < n; i++) c[i] = ch[i];
    } /* rfftf1 */


    void rfftb1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])
    {
        int i;
        int k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
        Treal *cinput, *coutput;
        nf = ifac[1];
        na = 0;
        l1 = 1;
        iw = 0;
        for (k1=1; k1<=nf; k1++) {
        ip = ifac[k1 + 1];
        l2 = ip*l1;
        ido = n / l2;
        idl1 = ido*l1;
        if (na) {
            cinput = ch;
            coutput = c;
        } else {
            cinput = c;
            coutput = ch;
        }
        switch (ip) {
        case 4:
            ix2 = iw + ido;
            ix3 = ix2 + ido;
            radb4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
            na = !na;
            break;
        case 2:
            radb2(ido, l1, cinput, coutput, &wa[iw]);
            na = !na;
            break;
        case 3:
            ix2 = iw + ido;
            radb3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
            na = !na;
            break;
        case 5:
            ix2 = iw + ido;
            ix3 = ix2 + ido;
            ix4 = ix3 + ido;
            radb5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
            na = !na;
            break;
        default:
            radbg(ido, ip, l1, idl1, cinput, coutput, &wa[iw]);
            if (ido == 1) na = !na;
        }
        l1 = l2;
        iw += (ip - 1)*ido;
        }
        if (na == 0) return;
        for (i=0; i<n; i++) c[i] = ch[i];
    } /* rfftb1 */


    void rfftf(int n, Treal r[], Treal wsave[])
    {
        if (n == 1) return;
        rfftf1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
    } /* rfftf */


    void rfftb(int n, Treal r[], Treal wsave[])
    {
        if (n == 1) return;
        rfftb1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
    } /* rfftb */


    static void rffti1(int n, Treal wa[], int ifac[MAXFAC+2])
    {
        static const Treal twopi = 6.28318530717959;
        Treal arg, argh, argld, fi;
        int i, j;
        int k1, l1, l2;
        int ld, ii, nf, ip, is;
        int ido, ipm, nfm1;
        static const int ntryh[NSPECIAL] = {
        4,2,3,5    }; /* Do not change the order of these. */
        factorize(n,ifac,ntryh);
        nf = ifac[1];
        argh = twopi / n;
        is = 0;
        nfm1 = nf - 1;
        l1 = 1;
        if (nfm1 == 0) return;
        for (k1 = 1; k1 <= nfm1; k1++) {
        ip = ifac[k1 + 1];
        ld = 0;
        l2 = l1*ip;
        ido = n / l2;
        ipm = ip - 1;
        for (j = 1; j <= ipm; ++j) {
            ld += l1;
            i = is;
            argld = (Treal) ld*argh;
            fi = 0;
            for (ii = 3; ii <= ido; ii += 2) {
            i += 2;
            fi += 1;
            arg = fi*argld;
            wa[i - 2] = cos(arg);
            wa[i - 1] = sin(arg);
            }
            is += ido;
        }
        l1 = l2;
        }
    } /* rffti1 */


    void rffti(int n, Treal wsave[])
    {
        if (n == 1) return;
        rffti1(n, wsave+n, (int*)(wsave+2*n));
    } /* rffti */

    #ifdef __cplusplus
    }
    #endif
''', libraries=[])

import numpy as np

def cffti(n):
    '''
    void cffti(int n, Treal wsave[])
    '''
    wsave = np.empty(4*n+15, 'double') # don't even ask me why 4 and why 15. Just took them from C-ext wrapper "as is".
    wsave_ptr = wsave.__array_interface__['data'][0]
    C.cffti(n, ffi.cast('double*', wsave_ptr))
    return wsave


def cfftf(a, wsave):
    '''
    void cfftf(int n, Troeal c[], Treal wsave[]);
    '''
    a_copy = a.copy()
    npts = a.shape[-1]
    if wsave.shape[0] != npts * 4 + 15:
        raise ValueError("invalid work array for fft size")
    nrepeats = a_copy.size / npts
    n = ffi.cast('int', npts)
    wsave_ptr = wsave.__array_interface__['data'][0]
    for i in range(nrepeats):
        print 'npts,i,a.shape',npts,i,a.shape,a.dtype
        dptr = a_copy.__array_interface__['data'][0] + npts * i
        a_copy_cdata = ffi.cast('double*', dptr)
        C.cfftf(n, a_copy_cdata, ffi.cast('double*', wsave_ptr))
    return a_copy
# is it a good idea to put it on a top-level here?.. maybe some singleton would be nicer or something... //typedef ... Treal; // damn, stuff like this doesn't work with cffi yet... structs are fine, but not the basic types :( // so, we stick to 'double' here: void cffti(int n, double wsave[]); void cfftf(int n, double c[], double wsave[]); /*     fftpack.c : A set of FFT routines in C.     Algorithmically based on Fortran-77 FFTPACK by <NAME> (Version 4, 1985).      */      /* isign is +1 for backward and -1 for forward transforms */      #include <math.h>     #include <stdio.h>     #define DOUBLE      #ifdef DOUBLE     #define Treal double     #else     #define Treal float     #endif       #define ref(u,a) u[a]      #define MAXFAC 13    /* maximum number of factors in factorization of n */     #define NSPECIAL 4   /* number of factors for which we have special-case routines */      #ifdef __cplusplus     extern "C" {     #endif       /* ----------------------------------------------------------------------     passf2, passf3, passf4, passf5, passf. Complex FFT passes fwd and bwd.     ---------------------------------------------------------------------- */      static void passf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[], int isign)     /* isign==+1 for backward transform */     {         int i, k, ah, ac;         Treal ti2, tr2;         if (ido <= 2) {         for (k=0; k<l1; k++) {             ah = k*ido;             ac = 2*k*ido;             ch[ah]              = ref(cc,ac) + ref(cc,ac + ido);             ch[ah + ido*l1]     = ref(cc,ac) - ref(cc,ac + ido);             ch[ah+1]            = ref(cc,ac+1) + ref(cc,ac + ido + 1);             ch[ah + ido*l1 + 1] = ref(cc,ac+1) - ref(cc,ac + ido + 1);         }         } else {         for (k=0; k<l1; k++) {             for (i=0; i<ido-1; i+=2) {             ah = i + k*ido;             ac = i + 2*k*ido;             ch[ah]   = ref(cc,ac) + ref(cc,ac + ido);             tr2      = ref(cc,ac) - ref(cc,ac + ido);             ch[ah+1] = ref(cc,ac+1) + ref(cc,ac + 1 + ido);             ti2      = ref(cc,ac+1) - ref(cc,ac + 1 + ido);             ch[ah+l1*ido+1] = wa1[i]*ti2 + isign*wa1[i+1]*tr2;             ch[ah+l1*ido]   = wa1[i]*tr2 - isign*wa1[i+1]*ti2;             }         }         }     } /* passf2 */       static void passf3(int ido, int l1, const Treal cc[], Treal ch[],         const Treal wa1[], const Treal wa2[], int isign)     /* isign==+1 for backward transform */     {         static const Treal taur = -0.5;         static const Treal taui = 0.866025403784439;         int i, k, ac, ah;         Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;         if (ido == 2) {         for (k=1; k<=l1; k++) {             ac = (3*k - 2)*ido;             tr2 = ref(cc,ac) + ref(cc,ac + ido);             cr2 = ref(cc,ac - ido) + taur*tr2;             ah = (k - 1)*ido;             ch[ah] = ref(cc,ac - ido) + tr2;              ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);             ci2 = ref(cc,ac - ido + 1) + taur*ti2;             ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;              cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));             ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));             ch[ah + l1*ido] = cr2 - ci3;             ch[ah + 2*l1*ido] = cr2 + ci3;             ch[ah + l1*ido + 1] = ci2 + cr3;             ch[ah + 2*l1*ido + 1] = ci2 - cr3;         }         } else {         for (k=1; k<=l1; k++) {             for (i=0; i<ido-1; i+=2) {             ac = i + (3*k - 2)*ido;             tr2 = ref(cc,ac) + ref(cc,ac + ido);             cr2 = ref(cc,ac - ido) + taur*tr2;             ah = i + (k-1)*ido;             ch[ah] = ref(cc,ac - ido) + tr2;             ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);             ci2 = ref(cc,ac - ido + 1) + taur*ti2;             ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;             cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));             ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));             dr2 = cr2 - ci3;             dr3 = cr2 + ci3;             di2 = ci2 + cr3;             di3 = ci2 - cr3;             ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;             ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;             ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;             ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;             }         }         }     } /* passf3 */       static void passf4(int ido, int l1, const Treal cc[], Treal ch[],         const Treal wa1[], const Treal wa2[], const Treal wa3[], int isign)     /* isign == -1 for forward transform and +1 for backward transform */     {         int i, k, ac, ah;         Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;         if (ido == 2) {         for (k=0; k<l1; k++) {             ac = 4*k*ido + 1;             ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);             ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);             tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);             ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);             tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);             tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);             ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);             tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);             ah = k*ido;             ch[ah] = tr2 + tr3;             ch[ah + 2*l1*ido] = tr2 - tr3;             ch[ah + 1] = ti2 + ti3;             ch[ah + 2*l1*ido + 1] = ti2 - ti3;             ch[ah + l1*ido] = tr1 + isign*tr4;             ch[ah + 3*l1*ido] = tr1 - isign*tr4;             ch[ah + l1*ido + 1] = ti1 + isign*ti4;             ch[ah + 3*l1*ido + 1] = ti1 - isign*ti4;         }         } else {         for (k=0; k<l1; k++) {             for (i=0; i<ido-1; i+=2) {             ac = i + 1 + 4*k*ido;             ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);             ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);             ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);             tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);             tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);             tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);             ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);             tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);             ah = i + k*ido;             ch[ah] = tr2 + tr3;             cr3 = tr2 - tr3;             ch[ah + 1] = ti2 + ti3;             ci3 = ti2 - ti3;             cr2 = tr1 + isign*tr4;             cr4 = tr1 - isign*tr4;             ci2 = ti1 + isign*ti4;             ci4 = ti1 - isign*ti4;             ch[ah + l1*ido] = wa1[i]*cr2 - isign*wa1[i + 1]*ci2;             ch[ah + l1*ido + 1] = wa1[i]*ci2 + isign*wa1[i + 1]*cr2;             ch[ah + 2*l1*ido] = wa2[i]*cr3 - isign*wa2[i + 1]*ci3;             ch[ah + 2*l1*ido + 1] = wa2[i]*ci3 + isign*wa2[i + 1]*cr3;             ch[ah + 3*l1*ido] = wa3[i]*cr4 -isign*wa3[i + 1]*ci4;             ch[ah + 3*l1*ido + 1] = wa3[i]*ci4 + isign*wa3[i + 1]*cr4;             }         }         }     } /* passf4 */       static void passf5(int ido, int l1, const Treal cc[], Treal ch[],         const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[], int isign)     /* isign == -1 for forward transform and +1 for backward transform */     {         static const Treal tr11 = 0.309016994374947;         static const Treal ti11 = 0.951056516295154;         static const Treal tr12 = -0.809016994374947;         static const Treal ti12 = 0.587785252292473;         int i, k, ac, ah;         Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,             ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;         if (ido == 2) {         for (k = 1; k <= l1; ++k) {             ac = (5*k - 4)*ido + 1;             ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);             ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);             ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);             ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);             tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);             tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);             tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);             tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);             ah = (k - 1)*ido;             ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;             ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;             cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;             ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;             cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;             ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;             cr5 = isign*(ti11*tr5 + ti12*tr4);             ci5 = isign*(ti11*ti5 + ti12*ti4);             cr4 = isign*(ti12*tr5 - ti11*tr4);             ci4 = isign*(ti12*ti5 - ti11*ti4);             ch[ah + l1*ido] = cr2 - ci5;             ch[ah + 4*l1*ido] = cr2 + ci5;             ch[ah + l1*ido + 1] = ci2 + cr5;             ch[ah + 2*l1*ido + 1] = ci3 + cr4;             ch[ah + 2*l1*ido] = cr3 - ci4;             ch[ah + 3*l1*ido] = cr3 + ci4;             ch[ah + 3*l1*ido + 1] = ci3 - cr4;             ch[ah + 4*l1*ido + 1] = ci2 - cr5;         }         } else {         for (k=1; k<=l1; k++) {             for (i=0; i<ido-1; i+=2) {             ac = i + 1 + (k*5 - 4)*ido;             ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);             ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);             ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);             ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);             tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);             tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);             tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);             tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);             ah = i + (k - 1)*ido;             ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;             ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;             cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;              ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;             cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;              ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;             cr5 = isign*(ti11*tr5 + ti12*tr4);             ci5 = isign*(ti11*ti5 + ti12*ti4);             cr4 = isign*(ti12*tr5 - ti11*tr4);             ci4 = isign*(ti12*ti5 - ti11*ti4);             dr3 = cr3 - ci4;             dr4 = cr3 + ci4;             di3 = ci3 + cr4;             di4 = ci3 - cr4;             dr5 = cr2 + ci5;             dr2 = cr2 - ci5;             di5 = ci2 - cr5;             di2 = ci2 + cr5;             ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;             ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;             ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;             ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;             ch[ah + 3*l1*ido] = wa3[i]*dr4 - isign*wa3[i+1]*di4;             ch[ah + 3*l1*ido + 1] = wa3[i]*di4 + isign*wa3[i+1]*dr4;             ch[ah + 4*l1*ido] = wa4[i]*dr5 - isign*wa4[i+1]*di5;             ch[ah + 4*l1*ido + 1] = wa4[i]*di5 + isign*wa4[i+1]*dr5;             }         }         }     } /* passf5 */       static void passf(int *nac, int ido, int ip, int l1, int idl1,         Treal cc[], Treal ch[],         const Treal wa[], int isign)     /* isign is -1 for forward transform and +1 for backward transform */     {         int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, idj, idl, inc,idp;         Treal wai, war;          idot = ido / 2;         /* nt = ip*idl1;*/         ipph = (ip + 1) / 2;         idp = ip*ido;         if (ido >= l1) {         for (j=1; j<ipph; j++) {             jc = ip - j;             for (k=0; k<l1; k++) {             for (i=0; i<ido; i++) {                 ch[i + (k + j*l1)*ido] =                     ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*ip)*ido);                 ch[i + (k + jc*l1)*ido] =                     ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*ip)*ido);             }             }         }         for (k=0; k<l1; k++)             for (i=0; i<ido; i++)             ch[i + k*ido] = ref(cc,i + k*ip*ido);         } else {         for (j=1; j<ipph; j++) {             jc = ip - j;             for (i=0; i<ido; i++) {             for (k=0; k<l1; k++) {                 ch[i + (k + j*l1)*ido] = ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*                     ip)*ido);                 ch[i + (k + jc*l1)*ido] = ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*                     ip)*ido);             }             }         }         for (i=0; i<ido; i++)             for (k=0; k<l1; k++)             ch[i + k*ido] = ref(cc,i + k*ip*ido);         }          idl = 2 - ido;         inc = 0;         for (l=1; l<ipph; l++) {         lc = ip - l;         idl += ido;         for (ik=0; ik<idl1; ik++) {             cc[ik + l*idl1] = ch[ik] + wa[idl - 2]*ch[ik + idl1];             cc[ik + lc*idl1] = isign*wa[idl-1]*ch[ik + (ip-1)*idl1];         }         idlj = idl;         inc += ido;         for (j=2; j<ipph; j++) {             jc = ip - j;             idlj += inc;             if (idlj > idp) idlj -= idp;             war = wa[idlj - 2];             wai = wa[idlj-1];             for (ik=0; ik<idl1; ik++) {             cc[ik + l*idl1] += war*ch[ik + j*idl1];             cc[ik + lc*idl1] += isign*wai*ch[ik + jc*idl1];             }         }         }         for (j=1; j<ipph; j++)         for (ik=0; ik<idl1; ik++)             ch[ik] += ch[ik + j*idl1];         for (j=1; j<ipph; j++) {         jc = ip - j;         for (ik=1; ik<idl1; ik+=2) {             ch[ik - 1 + j*idl1] = cc[ik - 1 + j*idl1] - cc[ik + jc*idl1];             ch[ik - 1 + jc*idl1] = cc[ik - 1 + j*idl1] + cc[ik + jc*idl1];             ch[ik + j*idl1] = cc[ik + j*idl1] + cc[ik - 1 + jc*idl1];             ch[ik + jc*idl1] = cc[ik + j*idl1] - cc[ik - 1 + jc*idl1];         }         }         *nac = 1;         if (ido == 2) return;         *nac = 0;         for (ik=0; ik<idl1; ik++)         cc[ik] = ch[ik];         for (j=1; j<ip; j++) {         for (k=0; k<l1; k++) {             cc[(k + j*l1)*ido + 0] = ch[(k + j*l1)*ido + 0];             cc[(k + j*l1)*ido + 1] = ch[(k + j*l1)*ido + 1];         }         }         if (idot <= l1) {         idij = 0;         for (j=1; j<ip; j++) {             idij += 2;             for (i=3; i<ido; i+=2) {             idij += 2;             for (k=0; k<l1; k++) {                 cc[i - 1 + (k + j*l1)*ido] =                     wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -                     isign*wa[idij-1]*ch[i + (k + j*l1)*ido];                 cc[i + (k + j*l1)*ido] =                     wa[idij - 2]*ch[i + (k + j*l1)*ido] +                     isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];             }             }         }         } else {         idj = 2 - ido;         for (j=1; j<ip; j++) {             idj += ido;             for (k = 0; k < l1; k++) {             idij = idj;             for (i=3; i<ido; i+=2) {                 idij += 2;                 cc[i - 1 + (k + j*l1)*ido] =                     wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -                     isign*wa[idij-1]*ch[i + (k + j*l1)*ido];                 cc[i + (k + j*l1)*ido] =                     wa[idij - 2]*ch[i + (k + j*l1)*ido] +                     isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];             }             }         }         }     } /* passf */       /* ----------------------------------------------------------------------     radf2,radb2, radf3,radb3, radf4,radb4, radf5,radb5, radfg,radbg.     Treal FFT passes fwd and bwd.     ---------------------------------------------------------------------- */      static void radf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])     {         int i, k, ic;         Treal ti2, tr2;         for (k=0; k<l1; k++) {         ch[2*k*ido] =             ref(cc,k*ido) + ref(cc,(k + l1)*ido);         ch[(2*k+1)*ido + ido-1] =             ref(cc,k*ido) - ref(cc,(k + l1)*ido);         }         if (ido < 2) return;         if (ido != 2) {         for (k=0; k<l1; k++) {             for (i=2; i<ido; i+=2) {             ic = ido - i;             tr2 = wa1[i - 2]*ref(cc, i-1 + (k + l1)*ido) + wa1[i - 1]*ref(cc, i + (k + l1)*ido);             ti2 = wa1[i - 2]*ref(cc, i + (k + l1)*ido) - wa1[i - 1]*ref(cc, i-1 + (k + l1)*ido);             ch[i + 2*k*ido] = ref(cc,i + k*ido) + ti2;             ch[ic + (2*k+1)*ido] = ti2 - ref(cc,i + k*ido);             ch[i - 1 + 2*k*ido] = ref(cc,i - 1 + k*ido) + tr2;             ch[ic - 1 + (2*k+1)*ido] = ref(cc,i - 1 + k*ido) - tr2;             }         }         if (ido % 2 == 1) return;         }         for (k=0; k<l1; k++) {         ch[(2*k+1)*ido] = -ref(cc,ido-1 + (k + l1)*ido);         ch[ido-1 + 2*k*ido] = ref(cc,ido-1 + k*ido);         }     } /* radf2 */       static void radb2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])     {         int i, k, ic;         Treal ti2, tr2;         for (k=0; k<l1; k++) {         ch[k*ido] =             ref(cc,2*k*ido) + ref(cc,ido-1 + (2*k+1)*ido);         ch[(k + l1)*ido] =             ref(cc,2*k*ido) - ref(cc,ido-1 + (2*k+1)*ido);         }         if (ido < 2) return;         if (ido != 2) {         for (k = 0; k < l1; ++k) {             for (i = 2; i < ido; i += 2) {             ic = ido - i;             ch[i-1 + k*ido] =                 ref(cc,i-1 + 2*k*ido) + ref(cc,ic-1 + (2*k+1)*ido);             tr2 = ref(cc,i-1 + 2*k*ido) - ref(cc,ic-1 + (2*k+1)*ido);             ch[i + k*ido] =                 ref(cc,i + 2*k*ido) - ref(cc,ic + (2*k+1)*ido);             ti2 = ref(cc,i + (2*k)*ido) + ref(cc,ic + (2*k+1)*ido);             ch[i-1 + (k + l1)*ido] =                 wa1[i - 2]*tr2 - wa1[i - 1]*ti2;             ch[i + (k + l1)*ido] =                 wa1[i - 2]*ti2 + wa1[i - 1]*tr2;             }         }         if (ido % 2 == 1) return;         }         for (k = 0; k < l1; k++) {         ch[ido-1 + k*ido] = 2*ref(cc,ido-1 + 2*k*ido);         ch[ido-1 + (k + l1)*ido] = -2*ref(cc,(2*k+1)*ido);         }     } /* radb2 */       static void radf3(int ido, int l1, const Treal cc[], Treal ch[],         const Treal wa1[], const Treal wa2[])     {         static const Treal taur = -0.5;         static const Treal taui = 0.866025403784439;         int i, k, ic;         Treal ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;         for (k=0; k<l1; k++) {         cr2 = ref(cc,(k + l1)*ido) + ref(cc,(k + 2*l1)*ido);         ch[3*k*ido] = ref(cc,k*ido) + cr2;         ch[(3*k+2)*ido] = taui*(ref(cc,(k + l1*2)*ido) - ref(cc,(k + l1)*ido));         ch[ido-1 + (3*k + 1)*ido] = ref(cc,k*ido) + taur*cr2;         }         if (ido == 1) return;         for (k=0; k<l1; k++) {         for (i=2; i<ido; i+=2) {             ic = ido - i;             dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) +                 wa1[i - 1]*ref(cc,i + (k + l1)*ido);             di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);             dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + l1*2)*ido) + wa2[i - 1]*ref(cc,i + (k + l1*2)*ido);             di3 = wa2[i - 2]*ref(cc,i + (k + l1*2)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + l1*2)*ido);             cr2 = dr2 + dr3;             ci2 = di2 + di3;             ch[i - 1 + 3*k*ido] = ref(cc,i - 1 + k*ido) + cr2;             ch[i + 3*k*ido] = ref(cc,i + k*ido) + ci2;             tr2 = ref(cc,i - 1 + k*ido) + taur*cr2;             ti2 = ref(cc,i + k*ido) + taur*ci2;             tr3 = taui*(di2 - di3);             ti3 = taui*(dr3 - dr2);             ch[i - 1 + (3*k + 2)*ido] = tr2 + tr3;             ch[ic - 1 + (3*k + 1)*ido] = tr2 - tr3;             ch[i + (3*k + 2)*ido] = ti2 + ti3;             ch[ic + (3*k + 1)*ido] = ti3 - ti2;         }         }     } /* radf3 */       static void radb3(int ido, int l1, const Treal cc[], Treal ch[],         const Treal wa1[], const Treal wa2[])     {         static const Treal taur = -0.5;         static const Treal taui = 0.866025403784439;         int i, k, ic;         Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;         for (k=0; k<l1; k++) {         tr2 = 2*ref(cc,ido-1 + (3*k + 1)*ido);         cr2 = ref(cc,3*k*ido) + taur*tr2;         ch[k*ido] = ref(cc,3*k*ido) + tr2;         ci3 = 2*taui*ref(cc,(3*k + 2)*ido);         ch[(k + l1)*ido] = cr2 - ci3;         ch[(k + 2*l1)*ido] = cr2 + ci3;         }         if (ido == 1) return;         for (k=0; k<l1; k++) {         for (i=2; i<ido; i+=2) {             ic = ido - i;             tr2 = ref(cc,i - 1 + (3*k + 2)*ido) + ref(cc,ic - 1 + (3*k + 1)*ido);             cr2 = ref(cc,i - 1 + 3*k*ido) + taur*tr2;             ch[i - 1 + k*ido] = ref(cc,i - 1 + 3*k*ido) + tr2;             ti2 = ref(cc,i + (3*k + 2)*ido) - ref(cc,ic + (3*k + 1)*ido);             ci2 = ref(cc,i + 3*k*ido) + taur*ti2;             ch[i + k*ido] = ref(cc,i + 3*k*ido) + ti2;             cr3 = taui*(ref(cc,i - 1 + (3*k + 2)*ido) - ref(cc,ic - 1 + (3*k + 1)*ido));             ci3 = taui*(ref(cc,i + (3*k + 2)*ido) + ref(cc,ic + (3*k + 1)*ido));             dr2 = cr2 - ci3;             dr3 = cr2 + ci3;             di2 = ci2 + cr3;             di3 = ci2 - cr3;             ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;             ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;             ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;             ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;         }         }     } /* radb3 */       static void radf4(int ido, int l1, const Treal cc[], Treal ch[],         const Treal wa1[], const Treal wa2[], const Treal wa3[])     {         static const Treal hsqt2 = 0.7071067811865475;         int i, k, ic;         Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;         for (k=0; k<l1; k++) {         tr1 = ref(cc,(k + l1)*ido) + ref(cc,(k + 3*l1)*ido);         tr2 = ref(cc,k*ido) + ref(cc,(k + 2*l1)*ido);         ch[4*k*ido] = tr1 + tr2;         ch[ido-1 + (4*k + 3)*ido] = tr2 - tr1;         ch[ido-1 + (4*k + 1)*ido] = ref(cc,k*ido) - ref(cc,(k + 2*l1)*ido);         ch[(4*k + 2)*ido] = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + l1)*ido);         }         if (ido < 2) return;         if (ido != 2) {         for (k=0; k<l1; k++) {             for (i=2; i<ido; i += 2) {             ic = ido - i;             cr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);             ci2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);             cr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*                 ido);             ci3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*                 ido);             cr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*                 ido);             ci4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*                 ido);             tr1 = cr2 + cr4;             tr4 = cr4 - cr2;             ti1 = ci2 + ci4;             ti4 = ci2 - ci4;             ti2 = ref(cc,i + k*ido) + ci3;             ti3 = ref(cc,i + k*ido) - ci3;             tr2 = ref(cc,i - 1 + k*ido) + cr3;             tr3 = ref(cc,i - 1 + k*ido) - cr3;             ch[i - 1 + 4*k*ido] = tr1 + tr2;             ch[ic - 1 + (4*k + 3)*ido] = tr2 - tr1;             ch[i + 4*k*ido] = ti1 + ti2;             ch[ic + (4*k + 3)*ido] = ti1 - ti2;             ch[i - 1 + (4*k + 2)*ido] = ti4 + tr3;             ch[ic - 1 + (4*k + 1)*ido] = tr3 - ti4;             ch[i + (4*k + 2)*ido] = tr4 + ti3;             ch[ic + (4*k + 1)*ido] = tr4 - ti3;             }         }         if (ido % 2 == 1) return;         }         for (k=0; k<l1; k++) {         ti1 = -hsqt2*(ref(cc,ido-1 + (k + l1)*ido) + ref(cc,ido-1 + (k + 3*l1)*ido));         tr1 = hsqt2*(ref(cc,ido-1 + (k + l1)*ido) - ref(cc,ido-1 + (k + 3*l1)*ido));         ch[ido-1 + 4*k*ido] = tr1 + ref(cc,ido-1 + k*ido);         ch[ido-1 + (4*k + 2)*ido] = ref(cc,ido-1 + k*ido) - tr1;         ch[(4*k + 1)*ido] = ti1 - ref(cc,ido-1 + (k + 2*l1)*ido);         ch[(4*k + 3)*ido] = ti1 + ref(cc,ido-1 + (k + 2*l1)*ido);         }     } /* radf4 */       static void radb4(int ido, int l1, const Treal cc[], Treal ch[],         const Treal wa1[], const Treal wa2[], const Treal wa3[])     {         static const Treal sqrt2 = 1.414213562373095;         int i, k, ic;         Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;         for (k = 0; k < l1; k++) {         tr1 = ref(cc,4*k*ido) - ref(cc,ido-1 + (4*k + 3)*ido);         tr2 = ref(cc,4*k*ido) + ref(cc,ido-1 + (4*k + 3)*ido);         tr3 = ref(cc,ido-1 + (4*k + 1)*ido) + ref(cc,ido-1 + (4*k + 1)*ido);         tr4 = ref(cc,(4*k + 2)*ido) + ref(cc,(4*k + 2)*ido);         ch[k*ido] = tr2 + tr3;         ch[(k + l1)*ido] = tr1 - tr4;         ch[(k + 2*l1)*ido] = tr2 - tr3;         ch[(k + 3*l1)*ido] = tr1 + tr4;         }         if (ido < 2) return;         if (ido != 2) {         for (k = 0; k < l1; ++k) {             for (i = 2; i < ido; i += 2) {             ic = ido - i;             ti1 = ref(cc,i + 4*k*ido) + ref(cc,ic + (4*k + 3)*ido);             ti2 = ref(cc,i + 4*k*ido) - ref(cc,ic + (4*k + 3)*ido);             ti3 = ref(cc,i + (4*k + 2)*ido) - ref(cc,ic + (4*k + 1)*ido);             tr4 = ref(cc,i + (4*k + 2)*ido) + ref(cc,ic + (4*k + 1)*ido);             tr1 = ref(cc,i - 1 + 4*k*ido) - ref(cc,ic - 1 + (4*k + 3)*ido);             tr2 = ref(cc,i - 1 + 4*k*ido) + ref(cc,ic - 1 + (4*k + 3)*ido);             ti4 = ref(cc,i - 1 + (4*k + 2)*ido) - ref(cc,ic - 1 + (4*k + 1)*ido);             tr3 = ref(cc,i - 1 + (4*k + 2)*ido) + ref(cc,ic - 1 + (4*k + 1)*ido);             ch[i - 1 + k*ido] = tr2 + tr3;             cr3 = tr2 - tr3;             ch[i + k*ido] = ti2 + ti3;             ci3 = ti2 - ti3;             cr2 = tr1 - tr4;             cr4 = tr1 + tr4;             ci2 = ti1 + ti4;             ci4 = ti1 - ti4;             ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*cr2 - wa1[i - 1]*ci2;             ch[i + (k + l1)*ido] = wa1[i - 2]*ci2 + wa1[i - 1]*cr2;             ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*cr3 - wa2[i - 1]*ci3;             ch[i + (k + 2*l1)*ido] = wa2[i - 2]*ci3 + wa2[i - 1]*cr3;             ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*cr4 - wa3[i - 1]*ci4;             ch[i + (k + 3*l1)*ido] = wa3[i - 2]*ci4 + wa3[i - 1]*cr4;             }         }         if (ido % 2 == 1) return;         }         for (k = 0; k < l1; k++) {         ti1 = ref(cc,(4*k + 1)*ido) + ref(cc,(4*k + 3)*ido);         ti2 = ref(cc,(4*k + 3)*ido) - ref(cc,(4*k + 1)*ido);         tr1 = ref(cc,ido-1 + 4*k*ido) - ref(cc,ido-1 + (4*k + 2)*ido);         tr2 = ref(cc,ido-1 + 4*k*ido) + ref(cc,ido-1 + (4*k + 2)*ido);         ch[ido-1 + k*ido] = tr2 + tr2;         ch[ido-1 + (k + l1)*ido] = sqrt2*(tr1 - ti1);         ch[ido-1 + (k + 2*l1)*ido] = ti2 + ti2;         ch[ido-1 + (k + 3*l1)*ido] = -sqrt2*(tr1 + ti1);         }     } /* radb4 */       static void radf5(int ido, int l1, const Treal cc[], Treal ch[],         const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])     {         static const Treal tr11 = 0.309016994374947;         static const Treal ti11 = 0.951056516295154;         static const Treal tr12 = -0.809016994374947;         static const Treal ti12 = 0.587785252292473;         int i, k, ic;         Treal ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5,             cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;         for (k = 0; k < l1; k++) {         cr2 = ref(cc,(k + 4*l1)*ido) + ref(cc,(k + l1)*ido);         ci5 = ref(cc,(k + 4*l1)*ido) - ref(cc,(k + l1)*ido);         cr3 = ref(cc,(k + 3*l1)*ido) + ref(cc,(k + 2*l1)*ido);         ci4 = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + 2*l1)*ido);         ch[5*k*ido] = ref(cc,k*ido) + cr2 + cr3;         ch[ido-1 + (5*k + 1)*ido] = ref(cc,k*ido) + tr11*cr2 + tr12*cr3;         ch[(5*k + 2)*ido] = ti11*ci5 + ti12*ci4;         ch[ido-1 + (5*k + 3)*ido] = ref(cc,k*ido) + tr12*cr2 + tr11*cr3;         ch[(5*k + 4)*ido] = ti12*ci5 - ti11*ci4;         }         if (ido == 1) return;         for (k = 0; k < l1; ++k) {         for (i = 2; i < ido; i += 2) {             ic = ido - i;             dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);             di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);             dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*ido);             di3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*ido);             dr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*ido);             di4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*ido);             dr5 = wa4[i - 2]*ref(cc,i - 1 + (k + 4*l1)*ido) + wa4[i - 1]*ref(cc,i + (k + 4*l1)*ido);             di5 = wa4[i - 2]*ref(cc,i + (k + 4*l1)*ido) - wa4[i - 1]*ref(cc,i - 1 + (k + 4*l1)*ido);             cr2 = dr2 + dr5;             ci5 = dr5 - dr2;             cr5 = di2 - di5;             ci2 = di2 + di5;             cr3 = dr3 + dr4;             ci4 = dr4 - dr3;             cr4 = di3 - di4;             ci3 = di3 + di4;             ch[i - 1 + 5*k*ido] = ref(cc,i - 1 + k*ido) + cr2 + cr3;             ch[i + 5*k*ido] = ref(cc,i + k*ido) + ci2 + ci3;             tr2 = ref(cc,i - 1 + k*ido) + tr11*cr2 + tr12*cr3;             ti2 = ref(cc,i + k*ido) + tr11*ci2 + tr12*ci3;             tr3 = ref(cc,i - 1 + k*ido) + tr12*cr2 + tr11*cr3;             ti3 = ref(cc,i + k*ido) + tr12*ci2 + tr11*ci3;             tr5 = ti11*cr5 + ti12*cr4;             ti5 = ti11*ci5 + ti12*ci4;             tr4 = ti12*cr5 - ti11*cr4;             ti4 = ti12*ci5 - ti11*ci4;             ch[i - 1 + (5*k + 2)*ido] = tr2 + tr5;             ch[ic - 1 + (5*k + 1)*ido] = tr2 - tr5;             ch[i + (5*k + 2)*ido] = ti2 + ti5;             ch[ic + (5*k + 1)*ido] = ti5 - ti2;             ch[i - 1 + (5*k + 4)*ido] = tr3 + tr4;             ch[ic - 1 + (5*k + 3)*ido] = tr3 - tr4;             ch[i + (5*k + 4)*ido] = ti3 + ti4;             ch[ic + (5*k + 3)*ido] = ti4 - ti3;         }         }     } /* radf5 */       static void radb5(int ido, int l1, const Treal cc[], Treal ch[],         const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])     {         static const Treal tr11 = 0.309016994374947;         static const Treal ti11 = 0.951056516295154;         static const Treal tr12 = -0.809016994374947;         static const Treal ti12 = 0.587785252292473;         int i, k, ic;         Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,             ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;         for (k = 0; k < l1; k++) {         ti5 = 2*ref(cc,(5*k + 2)*ido);         ti4 = 2*ref(cc,(5*k + 4)*ido);         tr2 = 2*ref(cc,ido-1 + (5*k + 1)*ido);         tr3 = 2*ref(cc,ido-1 + (5*k + 3)*ido);         ch[k*ido] = ref(cc,5*k*ido) + tr2 + tr3;         cr2 = ref(cc,5*k*ido) + tr11*tr2 + tr12*tr3;         cr3 = ref(cc,5*k*ido) + tr12*tr2 + tr11*tr3;         ci5 = ti11*ti5 + ti12*ti4;         ci4 = ti12*ti5 - ti11*ti4;         ch[(k + l1)*ido] = cr2 - ci5;         ch[(k + 2*l1)*ido] = cr3 - ci4;         ch[(k + 3*l1)*ido] = cr3 + ci4;         ch[(k + 4*l1)*ido] = cr2 + ci5;         }         if (ido == 1) return;         for (k = 0; k < l1; ++k) {         for (i = 2; i < ido; i += 2) {             ic = ido - i;             ti5 = ref(cc,i + (5*k + 2)*ido) + ref(cc,ic + (5*k + 1)*ido);             ti2 = ref(cc,i + (5*k + 2)*ido) - ref(cc,ic + (5*k + 1)*ido);             ti4 = ref(cc,i + (5*k + 4)*ido) + ref(cc,ic + (5*k + 3)*ido);             ti3 = ref(cc,i + (5*k + 4)*ido) - ref(cc,ic + (5*k + 3)*ido);             tr5 = ref(cc,i - 1 + (5*k + 2)*ido) - ref(cc,ic - 1 + (5*k + 1)*ido);             tr2 = ref(cc,i - 1 + (5*k + 2)*ido) + ref(cc,ic - 1 + (5*k + 1)*ido);             tr4 = ref(cc,i - 1 + (5*k + 4)*ido) - ref(cc,ic - 1 + (5*k + 3)*ido);             tr3 = ref(cc,i - 1 + (5*k + 4)*ido) + ref(cc,ic - 1 + (5*k + 3)*ido);             ch[i - 1 + k*ido] = ref(cc,i - 1 + 5*k*ido) + tr2 + tr3;             ch[i + k*ido] = ref(cc,i + 5*k*ido) + ti2 + ti3;             cr2 = ref(cc,i - 1 + 5*k*ido) + tr11*tr2 + tr12*tr3;              ci2 = ref(cc,i + 5*k*ido) + tr11*ti2 + tr12*ti3;             cr3 = ref(cc,i - 1 + 5*k*ido) + tr12*tr2 + tr11*tr3;              ci3 = ref(cc,i + 5*k*ido) + tr12*ti2 + tr11*ti3;             cr5 = ti11*tr5 + ti12*tr4;             ci5 = ti11*ti5 + ti12*ti4;             cr4 = ti12*tr5 - ti11*tr4;             ci4 = ti12*ti5 - ti11*ti4;             dr3 = cr3 - ci4;             dr4 = cr3 + ci4;             di3 = ci3 + cr4;             di4 = ci3 - cr4;             dr5 = cr2 + ci5;             dr2 = cr2 - ci5;             di5 = ci2 - cr5;             di2 = ci2 + cr5;             ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;             ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;             ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;             ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;             ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*dr4 - wa3[i - 1]*di4;             ch[i + (k + 3*l1)*ido] = wa3[i - 2]*di4 + wa3[i - 1]*dr4;             ch[i - 1 + (k + 4*l1)*ido] = wa4[i - 2]*dr5 - wa4[i - 1]*di5;             ch[i + (k + 4*l1)*ido] = wa4[i - 2]*di5 + wa4[i - 1]*dr5;         }         }     } /* radb5 */       static void radfg(int ido, int ip, int l1, int idl1,         Treal cc[], Treal ch[], const Treal wa[])     {         static const Treal twopi = 6.28318530717959;         int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;         Treal dc2, ai1, ai2, ar1, ar2, ds2, dcp, arg, dsp, ar1h, ar2h;         arg = twopi / ip;         dcp = cos(arg);         dsp = sin(arg);         ipph = (ip + 1) / 2;         nbd = (ido - 1) / 2;         if (ido != 1) {         for (ik=0; ik<idl1; ik++) ch[ik] = cc[ik];         for (j=1; j<ip; j++)             for (k=0; k<l1; k++)             ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido];         if (nbd <= l1) {             is = -ido;             for (j=1; j<ip; j++) {             is += ido;             idij = is-1;             for (i=2; i<ido; i+=2) {                 idij += 2;                 for (k=0; k<l1; k++) {                 ch[i - 1 + (k + j*l1)*ido] =                     wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];                 ch[i + (k + j*l1)*ido] =                     wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];                 }             }             }         } else {             is = -ido;             for (j=1; j<ip; j++) {             is += ido;             for (k=0; k<l1; k++) {                 idij = is-1;                 for (i=2; i<ido; i+=2) {                 idij += 2;                 ch[i - 1 + (k + j*l1)*ido] =                     wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];                 ch[i + (k + j*l1)*ido] =                     wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];                 }             }             }         }         if (nbd >= l1) {             for (j=1; j<ipph; j++) {             jc = ip - j;             for (k=0; k<l1; k++) {                 for (i=2; i<ido; i+=2) {                 cc[i - 1 + (k + j*l1)*ido] = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];                 cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];                 cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];                 cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];                 }             }             }         } else {             for (j=1; j<ipph; j++) {             jc = ip - j;             for (i=2; i<ido; i+=2) {                 for (k=0; k<l1; k++) {                 cc[i - 1 + (k + j*l1)*ido] =                     ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];                 cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];                 cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];                 cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];                 }             }             }         }         } else {  /* now ido == 1 */         for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];         }         for (j=1; j<ipph; j++) {         jc = ip - j;         for (k=0; k<l1; k++) {             cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido] + ch[(k + jc*l1)*ido];             cc[(k + jc*l1)*ido] = ch[(k + jc*l1)*ido] - ch[(k + j*l1)*ido];         }         }          ar1 = 1;         ai1 = 0;         for (l=1; l<ipph; l++) {         lc = ip - l;         ar1h = dcp*ar1 - dsp*ai1;         ai1 = dcp*ai1 + dsp*ar1;         ar1 = ar1h;         for (ik=0; ik<idl1; ik++) {             ch[ik + l*idl1] = cc[ik] + ar1*cc[ik + idl1];             ch[ik + lc*idl1] = ai1*cc[ik + (ip-1)*idl1];         }         dc2 = ar1;         ds2 = ai1;         ar2 = ar1;         ai2 = ai1;         for (j=2; j<ipph; j++) {             jc = ip - j;             ar2h = dc2*ar2 - ds2*ai2;             ai2 = dc2*ai2 + ds2*ar2;             ar2 = ar2h;             for (ik=0; ik<idl1; ik++) {             ch[ik + l*idl1] += ar2*cc[ik + j*idl1];             ch[ik + lc*idl1] += ai2*cc[ik + jc*idl1];             }         }         }         for (j=1; j<ipph; j++)         for (ik=0; ik<idl1; ik++)             ch[ik] += cc[ik + j*idl1];          if (ido >= l1) {         for (k=0; k<l1; k++) {             for (i=0; i<ido; i++) {             ref(cc,i + k*ip*ido) = ch[i + k*ido];             }         }         } else {         for (i=0; i<ido; i++) {             for (k=0; k<l1; k++) {             ref(cc,i + k*ip*ido) = ch[i + k*ido];             }         }         }         for (j=1; j<ipph; j++) {         jc = ip - j;         j2 = 2*j;         for (k=0; k<l1; k++) {             ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) =                 ch[(k + j*l1)*ido];             ref(cc,(j2 + k*ip)*ido) =                 ch[(k + jc*l1)*ido];         }         }         if (ido == 1) return;         if (nbd >= l1) {         for (j=1; j<ipph; j++) {             jc = ip - j;             j2 = 2*j;             for (k=0; k<l1; k++) {             for (i=2; i<ido; i+=2) {                 ic = ido - i;                 ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];                 ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];                 ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];                 ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];             }             }         }         } else {         for (j=1; j<ipph; j++) {             jc = ip - j;             j2 = 2*j;             for (i=2; i<ido; i+=2) {             ic = ido - i;             for (k=0; k<l1; k++) {                 ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];                 ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];                 ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];                 ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];             }             }         }         }     } /* radfg */       static void radbg(int ido, int ip, int l1, int idl1,         Treal cc[], Treal ch[], const Treal wa[])     {         static const Treal twopi = 6.28318530717959;         int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;         Treal dc2, ai1, ai2, ar1, ar2, ds2;         int nbd;         Treal dcp, arg, dsp, ar1h, ar2h;         arg = twopi / ip;         dcp = cos(arg);         dsp = sin(arg);         nbd = (ido - 1) / 2;         ipph = (ip + 1) / 2;         if (ido >= l1) {         for (k=0; k<l1; k++) {             for (i=0; i<ido; i++) {             ch[i + k*ido] = ref(cc,i + k*ip*ido);             }         }         } else {         for (i=0; i<ido; i++) {             for (k=0; k<l1; k++) {             ch[i + k*ido] = ref(cc,i + k*ip*ido);             }         }         }         for (j=1; j<ipph; j++) {         jc = ip - j;         j2 = 2*j;         for (k=0; k<l1; k++) {             ch[(k + j*l1)*ido] = ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) + ref(cc,ido-1 + (j2 - 1 + k*ip)*                 ido);             ch[(k + jc*l1)*ido] = ref(cc,(j2 + k*ip)*ido) + ref(cc,(j2 + k*ip)*ido);         }         }          if (ido != 1) {         if (nbd >= l1) {             for (j=1; j<ipph; j++) {             jc = ip - j;             for (k=0; k<l1; k++) {                 for (i=2; i<ido; i+=2) {                 ic = ido - i;                 ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,                     ic - 1 + (2*j - 1 + k*ip)*ido);                 ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -                     ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);                 ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic                     + (2*j - 1 + k*ip)*ido);                 ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic                     + (2*j - 1 + k*ip)*ido);                 }             }             }         } else {             for (j=1; j<ipph; j++) {             jc = ip - j;             for (i=2; i<ido; i+=2) {                 ic = ido - i;                 for (k=0; k<l1; k++) {                 ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,                     ic - 1 + (2*j - 1 + k*ip)*ido);                 ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -                     ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);                 ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic                     + (2*j - 1 + k*ip)*ido);                 ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic                     + (2*j - 1 + k*ip)*ido);                 }             }             }         }         }          ar1 = 1;         ai1 = 0;         for (l=1; l<ipph; l++) {         lc = ip - l;         ar1h = dcp*ar1 - dsp*ai1;         ai1 = dcp*ai1 + dsp*ar1;         ar1 = ar1h;         for (ik=0; ik<idl1; ik++) {             cc[ik + l*idl1] = ch[ik] + ar1*ch[ik + idl1];             cc[ik + lc*idl1] = ai1*ch[ik + (ip-1)*idl1];         }         dc2 = ar1;         ds2 = ai1;         ar2 = ar1;         ai2 = ai1;         for (j=2; j<ipph; j++) {             jc = ip - j;             ar2h = dc2*ar2 - ds2*ai2;             ai2 = dc2*ai2 + ds2*ar2;             ar2 = ar2h;             for (ik=0; ik<idl1; ik++) {             cc[ik + l*idl1] += ar2*ch[ik + j*idl1];             cc[ik + lc*idl1] += ai2*ch[ik + jc*idl1];             }         }         }         for (j=1; j<ipph; j++) {         for (ik=0; ik<idl1; ik++) {             ch[ik] += ch[ik + j*idl1];         }         }         for (j=1; j<ipph; j++) {         jc = ip - j;         for (k=0; k<l1; k++) {             ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido] - cc[(k + jc*l1)*ido];             ch[(k + jc*l1)*ido] = cc[(k + j*l1)*ido] + cc[(k + jc*l1)*ido];         }         }          if (ido == 1) return;         if (nbd >= l1) {         for (j=1; j<ipph; j++) {             jc = ip - j;             for (k=0; k<l1; k++) {             for (i=2; i<ido; i+=2) {                 ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];                 ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] + cc[i + (k + jc*l1)*ido];                 ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];                 ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];             }             }         }         } else {         for (j=1; j<ipph; j++) {             jc = ip - j;             for (i=2; i<ido; i+=2) {             for (k=0; k<l1; k++) {                 ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];                 ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j *l1)*ido] + cc[i + (k + jc*l1)*ido];                 ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];                 ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];             }             }         }         }         for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];         for (j=1; j<ip; j++)         for (k=0; k<l1; k++)             cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido];         if (nbd <= l1) {         is = -ido;         for (j=1; j<ip; j++) {             is += ido;             idij = is-1;             for (i=2; i<ido; i+=2) {             idij += 2;             for (k=0; k<l1; k++) {                 cc[i - 1 + (k + j*l1)*ido] = wa[idij - 1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*                     ch[i + (k + j*l1)*ido];                 cc[i + (k + j*l1)*ido] = wa[idij - 1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];             }             }         }         } else {         is = -ido;         for (j=1; j<ip; j++) {             is += ido;             for (k=0; k<l1; k++) {             idij = is - 1;             for (i=2; i<ido; i+=2) {                 idij += 2;                 cc[i - 1 + (k + j*l1)*ido] = wa[idij-1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*                     ch[i + (k + j*l1)*ido];                 cc[i + (k + j*l1)*ido] = wa[idij-1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];             }             }         }         }     } /* radbg */      /* ----------------------------------------------------------------------     cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.     ---------------------------------------------------------------------- */      static void cfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2], int isign)     {         int idot, i;         int k1, l1, l2;         int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;         Treal *cinput, *coutput;         nf = ifac[1];         na = 0;         l1 = 1;         iw = 0;         for (k1=2; k1<=nf+1; k1++) {         ip = ifac[k1];         l2 = ip*l1;         ido = n / l2;         idot = ido + ido;         idl1 = idot*l1;         if (na) {             cinput = ch;             coutput = c;         } else {             cinput = c;             coutput = ch;         }         switch (ip) {         case 4:             ix2 = iw + idot;             ix3 = ix2 + idot;             passf4(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], isign);             na = !na;             break;         case 2:             passf2(idot, l1, cinput, coutput, &wa[iw], isign);             na = !na;             break;         case 3:             ix2 = iw + idot;             passf3(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], isign);             na = !na;             break;         case 5:             ix2 = iw + idot;             ix3 = ix2 + idot;             ix4 = ix3 + idot;             passf5(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);             na = !na;             break;         default:             passf(&nac, idot, ip, l1, idl1, cinput, coutput, &wa[iw], isign);             if (nac != 0) na = !na;         }         l1 = l2;         iw += (ip - 1)*idot;         }         if (na == 0) return;         for (i=0; i<2*n; i++) c[i] = ch[i];     } /* cfftf1 */       void cfftf(int n, Treal c[], Treal wsave[])     {         int iw1, iw2;         if (n == 1) return;         iw1 = 2*n;         iw2 = iw1 + 2*n;         cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), -1);     } /* cfftf */       void cfftb(int n, Treal c[], Treal wsave[])     {         int iw1, iw2;         if (n == 1) return;         iw1 = 2*n;         iw2 = iw1 + 2*n;         cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), +1);     } /* cfftb */       static void factorize(int n, int ifac[MAXFAC+2], const int ntryh[NSPECIAL])     /* Factorize n in factors in ntryh and rest. On exit,     ifac[0] contains n and ifac[1] contains number of factors,     the factors start from ifac[2]. */     {         int ntry=3, i, j=0, ib, nf=0, nl=n, nq, nr;     startloop:         if (j < NSPECIAL)         ntry = ntryh[j];         else         ntry+= 2;         j++;         do {         nq = nl / ntry;         nr = nl - ntry*nq;         if (nr != 0) goto startloop;         nf++;         ifac[nf + 1] = ntry;         nl = nq;         if (ntry == 2 && nf != 1) {             for (i=2; i<=nf; i++) {             ib = nf - i + 2;             ifac[ib + 1] = ifac[ib];             }             ifac[2] = 2;         }         } while (nl != 1);         ifac[0] = n;         ifac[1] = nf;     }       static void cffti1(int n, Treal wa[], int ifac[MAXFAC+2])     {         static const Treal twopi = 6.28318530717959;         Treal arg, argh, argld, fi;         int idot, i, j;         int i1, k1, l1, l2;         int ld, ii, nf, ip;         int ido, ipm;          static const int ntryh[NSPECIAL] = {         3,4,2,5    }; /* Do not change the order of these. */          factorize(n,ifac,ntryh);         nf = ifac[1];         argh = twopi/(Treal)n;         i = 1;         l1 = 1;         for (k1=1; k1<=nf; k1++) {         ip = ifac[k1+1];         ld = 0;         l2 = l1*ip;         ido = n / l2;         idot = ido + ido + 2;         ipm = ip - 1;         for (j=1; j<=ipm; j++) {             i1 = i;             wa[i-1] = 1;             wa[i] = 0;             ld += l1;             fi = 0;             argld = ld*argh;             for (ii=4; ii<=idot; ii+=2) {             i+= 2;             fi+= 1;             arg = fi*argld;             wa[i-1] = cos(arg);             wa[i] = sin(arg);             }             if (ip > 5) {             wa[i1-1] = wa[i-1];             wa[i1] = wa[i];             }         }         l1 = l2;         }     } /* cffti1 */       void cffti(int n, Treal wsave[])     {         int iw1, iw2;         if (n == 1) return;         iw1 = 2*n;         iw2 = iw1 + 2*n;         cffti1(n, wsave+iw1, (int*)(wsave+iw2));     } /* cffti */      /* ----------------------------------------------------------------------     rfftf1, rfftb1, rfftf, rfftb, rffti1, rffti. Treal FFTs.     ---------------------------------------------------------------------- */      static void rfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])     {         int i;         int k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1;         Treal *cinput, *coutput;         nf = ifac[1];         na = 1;         l2 = n;         iw = n-1;         for (k1 = 1; k1 <= nf; ++k1) {         kh = nf - k1;         ip = ifac[kh + 2];         l1 = l2 / ip;         ido = n / l2;         idl1 = ido*l1;         iw -= (ip - 1)*ido;         na = !na;         if (na) {             cinput = ch;             coutput = c;         } else {             cinput = c;             coutput = ch;         }         switch (ip) {         case 4:             ix2 = iw + ido;             ix3 = ix2 + ido;             radf4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);             break;         case 2:             radf2(ido, l1, cinput, coutput, &wa[iw]);             break;         case 3:             ix2 = iw + ido;             radf3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);             break;         case 5:             ix2 = iw + ido;             ix3 = ix2 + ido;             ix4 = ix3 + ido;             radf5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);             break;         default:             if (ido == 1)             na = !na;             if (na == 0) {             radfg(ido, ip, l1, idl1, c, ch, &wa[iw]);             na = 1;             } else {             radfg(ido, ip, l1, idl1, ch, c, &wa[iw]);             na = 0;             }         }         l2 = l1;         }         if (na == 1) return;         for (i = 0; i < n; i++) c[i] = ch[i];     } /* rfftf1 */       void rfftb1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])     {         int i;         int k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;         Treal *cinput, *coutput;         nf = ifac[1];         na = 0;         l1 = 1;         iw = 0;         for (k1=1; k1<=nf; k1++) {         ip = ifac[k1 + 1];         l2 = ip*l1;         ido = n / l2;         idl1 = ido*l1;         if (na) {             cinput = ch;             coutput = c;         } else {             cinput = c;             coutput = ch;         }         switch (ip) {         case 4:             ix2 = iw + ido;             ix3 = ix2 + ido;             radb4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);             na = !na;             break;         case 2:             radb2(ido, l1, cinput, coutput, &wa[iw]);             na = !na;             break;         case 3:             ix2 = iw + ido;             radb3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);             na = !na;             break;         case 5:             ix2 = iw + ido;             ix3 = ix2 + ido;             ix4 = ix3 + ido;             radb5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);             na = !na;             break;         default:             radbg(ido, ip, l1, idl1, cinput, coutput, &wa[iw]);             if (ido == 1) na = !na;         }         l1 = l2;         iw += (ip - 1)*ido;         }         if (na == 0) return;         for (i=0; i<n; i++) c[i] = ch[i];     } /* rfftb1 */       void rfftf(int n, Treal r[], Treal wsave[])     {         if (n == 1) return;         rfftf1(n, r, wsave, wsave+n, (int*)(wsave+2*n));     } /* rfftf */       void rfftb(int n, Treal r[], Treal wsave[])     {         if (n == 1) return;         rfftb1(n, r, wsave, wsave+n, (int*)(wsave+2*n));     } /* rfftb */       static void rffti1(int n, Treal wa[], int ifac[MAXFAC+2])     {         static const Treal twopi = 6.28318530717959;         Treal arg, argh, argld, fi;         int i, j;         int k1, l1, l2;         int ld, ii, nf, ip, is;         int ido, ipm, nfm1;         static const int ntryh[NSPECIAL] = {         4,2,3,5    }; /* Do not change the order of these. */         factorize(n,ifac,ntryh);         nf = ifac[1];         argh = twopi / n;         is = 0;         nfm1 = nf - 1;         l1 = 1;         if (nfm1 == 0) return;         for (k1 = 1; k1 <= nfm1; k1++) {         ip = ifac[k1 + 1];         ld = 0;         l2 = l1*ip;         ido = n / l2;         ipm = ip - 1;         for (j = 1; j <= ipm; ++j) {             ld += l1;             i = is;             argld = (Treal) ld*argh;             fi = 0;             for (ii = 3; ii <= ido; ii += 2) {             i += 2;             fi += 1;             arg = fi*argld;             wa[i - 2] = cos(arg);             wa[i - 1] = sin(arg);             }             is += ido;         }         l1 = l2;         }     } /* rffti1 */       void rffti(int n, Treal wsave[])     {         if (n == 1) return;         rffti1(n, wsave+n, (int*)(wsave+2*n));     } /* rffti */      #ifdef __cplusplus     }     #endif void cffti(int n, Treal wsave[]) # don't even ask me why 4 and why 15. Just took them from C-ext wrapper "as is". void cfftf(int n, Troeal c[], Treal wsave[]);