y3i = a[j1 - 1] - a[j3 - 1];        a[j0] = x0r + x2r;        a[j0 + 1] = x0i + x2i;        a[j0 - 2] = y0r + y2r;        a[j0 - 1] = y0i + y2i;        a[j1] = x0r - x2r;        a[j1 + 1] = x0i - x2i;        a[j1 - 2] = y0r - y2r;        a[j1 - 1] = y0i - y2i;        x0r = x1r - x3i;        x0i = x1i + x3r;        a[j2] = wk1i * x0r - wk1r * x0i;        a[j2 + 1] = wk1i * x0i + wk1r * x0r;        x0r = y1r - y3i;        x0i = y1i + y3r;        a[j2 - 2] = wd1i * x0r - wd1r * x0i;        a[j2 - 1] = wd1i * x0i + wd1r * x0r;        x0r = x1r + x3i;        x0i = x1i - x3r;        a[j3] = wk3i * x0r + wk3r * x0i;        a[j3 + 1] = wk3i * x0i - wk3r * x0r;        x0r = y1r + y3i;        x0i = y1i - y3r;        a[j3 - 2] = wd3i * x0r + wd3r * x0i;        a[j3 - 1] = wd3i * x0i - wd3r * x0r;    }    wk1r = csc1 * (wd1r + wn4r);    wk1i = csc1 * (wd1i + wn4r);    wk3r = csc3 * (wd3r - wn4r);    wk3i = csc3 * (wd3i - wn4r);    j0 = mh;    j1 = j0 + m;    j2 = j1 + m;    j3 = j2 + m;    x0r = a[j0 - 2] + a[j2 - 2];    x0i = a[j0 - 1] + a[j2 - 1];    x1r = a[j0 - 2] - a[j2 - 2];    x1i = a[j0 - 1] - a[j2 - 1];    x2r = a[j1 - 2] + a[j3 - 2];    x2i = a[j1 - 1] + a[j3 - 1];    x3r = a[j1 - 2] - a[j3 - 2];    x3i = a[j1 - 1] - a[j3 - 1];    a[j0 - 2] = x0r + x2r;    a[j0 - 1] = x0i + x2i;    a[j1 - 2] = x0r - x2r;    a[j1 - 1] = x0i - x2i;    x0r = x1r - x3i;    x0i = x1i + x3r;    a[j2 - 2] = wk1r * x0r - wk1i * x0i;    a[j2 - 1] = wk1r * x0i + wk1i * x0r;    x0r = x1r + x3i;    x0i = x1i - x3r;    a[j3 - 2] = wk3r * x0r + wk3i * x0i;    a[j3 - 1] = wk3r * x0i - wk3i * x0r;    x0r = a[j0] + a[j2];    x0i = a[j0 + 1] + a[j2 + 1];    x1r = a[j0] - a[j2];    x1i = a[j0 + 1] - a[j2 + 1];    x2r = a[j1] + a[j3];    x2i = a[j1 + 1] + a[j3 + 1];    x3r = a[j1] - a[j3];    x3i = a[j1 + 1] - a[j3 + 1];    a[j0] = x0r + x2r;    a[j0 + 1] = x0i + x2i;    a[j1] = x0r - x2r;    a[j1 + 1] = x0i - x2i;    x0r = x1r - x3i;    x0i = x1i + x3r;    a[j2] = wn4r * (x0r - x0i);    a[j2 + 1] = wn4r * (x0i + x0r);    x0r = x1r + x3i;    x0i = x1i - x3r;    a[j3] = -wn4r * (x0r + x0i);    a[j3 + 1] = -wn4r * (x0i - x0r);    x0r = a[j0 + 2] + a[j2 + 2];    x0i = a[j0 + 3] + a[j2 + 3];    x1r = a[j0 + 2] - a[j2 + 2];    x1i = a[j0 + 3] - a[j2 + 3];    x2r = a[j1 + 2] + a[j3 + 2];    x2i = a[j1 + 3] + a[j3 + 3];    x3r = a[j1 + 2] - a[j3 + 2];    x3i = a[j1 + 3] - a[j3 + 3];    a[j0 + 2] = x0r + x2r;    a[j0 + 3] = x0i + x2i;    a[j1 + 2] = x0r - x2r;    a[j1 + 3] = x0i - x2i;    x0r = x1r - x3i;    x0i = x1i + x3r;    a[j2 + 2] = wk1i * x0r - wk1r * x0i;    a[j2 + 3] = wk1i * x0i + wk1r * x0r;    x0r = x1r + x3i;    x0i = x1i - x3r;    a[j3 + 2] = wk3i * x0r + wk3r * x0i;    a[j3 + 3] = wk3i * x0i - wk3r * x0r;}void cftb1st(int n, double *a, double *w){    int j, j0, j1, j2, j3, k, m, mh;    double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i,         wd1r, wd1i, wd3r, wd3i;    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;        mh = n >> 3;    m = 2 * mh;    j1 = m;    j2 = j1 + m;    j3 = j2 + m;    x0r = a[0] + a[j2];    x0i = -a[1] - a[j2 + 1];    x1r = a[0] - a[j2];    x1i = -a[1] + a[j2 + 1];    x2r = a[j1] + a[j3];    x2i = a[j1 + 1] + a[j3 + 1];    x3r = a[j1] - a[j3];    x3i = a[j1 + 1] - a[j3 + 1];    a[0] = x0r + x2r;    a[1] = x0i - x2i;    a[j1] = x0r - x2r;    a[j1 + 1] = x0i + x2i;    a[j2] = x1r + x3i;    a[j2 + 1] = x1i + x3r;    a[j3] = x1r - x3i;    a[j3 + 1] = x1i - x3r;    wn4r = w[1];    csc1 = w[2];    csc3 = w[3];    wd1r = 1;    wd1i = 0;    wd3r = 1;    wd3i = 0;    k = 0;    for (j = 2; j < mh - 2; j += 4) {        k += 4;        wk1r = csc1 * (wd1r + w[k]);        wk1i = csc1 * (wd1i + w[k + 1]);        wk3r = csc3 * (wd3r + w[k + 2]);        wk3i = csc3 * (wd3i + w[k + 3]);        wd1r = w[k];        wd1i = w[k + 1];        wd3r = w[k + 2];        wd3i = w[k + 3];        j1 = j + m;        j2 = j1 + m;        j3 = j2 + m;        x0r = a[j] + a[j2];        x0i = -a[j + 1] - a[j2 + 1];        x1r = a[j] - a[j2];        x1i = -a[j + 1] + a[j2 + 1];        y0r = a[j + 2] + a[j2 + 2];        y0i = -a[j + 3] - a[j2 + 3];        y1r = a[j + 2] - a[j2 + 2];        y1i = -a[j + 3] + a[j2 + 3];        x2r = a[j1] + a[j3];        x2i = a[j1 + 1] + a[j3 + 1];        x3r = a[j1] - a[j3];        x3i = a[j1 + 1] - a[j3 + 1];        y2r = a[j1 + 2] + a[j3 + 2];        y2i = a[j1 + 3] + a[j3 + 3];        y3r = a[j1 + 2] - a[j3 + 2];        y3i = a[j1 + 3] - a[j3 + 3];        a[j] = x0r + x2r;        a[j + 1] = x0i - x2i;        a[j + 2] = y0r + y2r;        a[j + 3] = y0i - y2i;        a[j1] = x0r - x2r;        a[j1 + 1] = x0i + x2i;        a[j1 + 2] = y0r - y2r;        a[j1 + 3] = y0i + y2i;        x0r = x1r + x3i;        x0i = x1i + x3r;        a[j2] = wk1r * x0r - wk1i * x0i;        a[j2 + 1] = wk1r * x0i + wk1i * x0r;        x0r = y1r + y3i;        x0i = y1i + y3r;        a[j2 + 2] = wd1r * x0r - wd1i * x0i;        a[j2 + 3] = wd1r * x0i + wd1i * x0r;        x0r = x1r - x3i;        x0i = x1i - x3r;        a[j3] = wk3r * x0r + wk3i * x0i;        a[j3 + 1] = wk3r * x0i - wk3i * x0r;        x0r = y1r - y3i;        x0i = y1i - y3r;        a[j3 + 2] = wd3r * x0r + wd3i * x0i;        a[j3 + 3] = wd3r * x0i - wd3i * x0r;        j0 = m - j;        j1 = j0 + m;        j2 = j1 + m;        j3 = j2 + m;        x0r = a[j0] + a[j2];        x0i = -a[j0 + 1] - a[j2 + 1];        x1r = a[j0] - a[j2];        x1i = -a[j0 + 1] + a[j2 + 1];        y0r = a[j0 - 2] + a[j2 - 2];        y0i = -a[j0 - 1] - a[j2 - 1];        y1r = a[j0 - 2] - a[j2 - 2];        y1i = -a[j0 - 1] + a[j2 - 1];        x2r = a[j1] + a[j3];        x2i = a[j1 + 1] + a[j3 + 1];        x3r = a[j1] - a[j3];        x3i = a[j1 + 1] - a[j3 + 1];        y2r = a[j1 - 2] + a[j3 - 2];        y2i = a[j1 - 1] + a[j3 - 1];        y3r = a[j1 - 2] - a[j3 - 2];        y3i = a[j1 - 1] - a[j3 - 1];        a[j0] = x0r + x2r;        a[j0 + 1] = x0i - x2i;        a[j0 - 2] = y0r + y2r;        a[j0 - 1] = y0i - y2i;        a[j1] = x0r - x2r;        a[j1 + 1] = x0i + x2i;        a[j1 - 2] = y0r - y2r;        a[j1 - 1] = y0i + y2i;        x0r = x1r + x3i;        x0i = x1i + x3r;        a[j2] = wk1i * x0r - wk1r * x0i;        a[j2 + 1] = wk1i * x0i + wk1r * x0r;        x0r = y1r + y3i;        x0i = y1i + y3r;        a[j2 - 2] = wd1i * x0r - wd1r * x0i;        a[j2 - 1] = wd1i * x0i + wd1r * x0r;        x0r = x1r - x3i;        x0i = x1i - x3r;        a[j3] = wk3i * x0r + wk3r * x0i;        a[j3 + 1] = wk3i * x0i - wk3r * x0r;        x0r = y1r - y3i;        x0i = y1i - y3r;        a[j3 - 2] = wd3i * x0r + wd3r * x0i;        a[j3 - 1] = wd3i * x0i - wd3r * x0r;    }    wk1r = csc1 * (wd1r + wn4r);    wk1i = csc1 * (wd1i + wn4r);    wk3r = csc3 * (wd3r - wn4r);    wk3i = csc3 * (wd3i - wn4r);    j0 = mh;    j1 = j0 + m;    j2 = j1 + m;    j3 = j2 + m;    x0r = a[j0 - 2] + a[j2 - 2];    x0i = -a[j0 - 1] - a[j2 - 1];    x1r = a[j0 - 2] - a[j2 - 2];    x1i = -a[j0 - 1] + a[j2 - 1];    x2r = a[j1 - 2] + a[j3 - 2];    x2i = a[j1 - 1] + a[j3 - 1];    x3r = a[j1 - 2] - a[j3 - 2];    x3i = a[j1 - 1] - a[j3 - 1];    a[j0 - 2] = x0r + x2r;    a[j0 - 1] = x0i - x2i;    a[j1 - 2] = x0r - x2r;    a[j1 - 1] = x0i + x2i;    x0r = x1r + x3i;    x0i = x1i + x3r;    a[j2 - 2] = wk1r * x0r - wk1i * x0i;    a[j2 - 1] = wk1r * x0i + wk1i * x0r;    x0r = x1r - x3i;    x0i = x1i - x3r;    a[j3 - 2] = wk3r * x0r + wk3i * x0i;    a[j3 - 1] = wk3r * x0i - wk3i * x0r;    x0r = a[j0] + a[j2];    x0i = -a[j0 + 1] - a[j2 + 1];    x1r = a[j0] - a[j2];    x1i = -a[j0 + 1] + a[j2 + 1];    x2r = a[j1] + a[j3];    x2i = a[j1 + 1] + a[j3 + 1];    x3r = a[j1] - a[j3];    x3i = a[j1 + 1] - a[j3 + 1];    a[j0] = x0r + x2r;    a[j0 + 1] = x0i - x2i;    a[j1] = x0r - x2r;    a[j1 + 1] = x0i + x2i;    x0r = x1r + x3i;    x0i = x1i + x3r;    a[j2] = wn4r * (x0r - x0i);    a[j2 + 1] = wn4r * (x0i + x0r);    x0r = x1r - x3i;    x0i = x1i - x3r;    a[j3] = -wn4r * (x0r + x0i);    a[j3 + 1] = -wn4r * (x0i - x0r);    x0r = a[j0 + 2] + a[j2 + 2];    x0i = -a[j0 + 3] - a[j2 + 3];    x1r = a[j0 + 2] - a[j2 + 2];    x1i = -a[j0 + 3] + a[j2 + 3];    x2r = a[j1 + 2] + a[j3 + 2];    x2i = a[j1 + 3] + a[j3 + 3];    x3r = a[j1 + 2] - a[j3 + 2];    x3i = a[j1 + 3] - a[j3 + 3];    a[j0 + 2] = x0r + x2r;    a[j0 + 3] = x0i - x2i;    a[j1 + 2] = x0r - x2r;    a[j1 + 3] = x0i + x2i;    x0r = x1r + x3i;    x0i = x1i + x3r;    a[j2 + 2] = wk1i * x0r - wk1r * x0i;    a[j2 + 3] = wk1i * x0i + wk1r * x0r;    x0r = x1r - x3i;    x0i = x1i - x3r;    a[j3 + 2] = wk3i * x0r + wk3r * x0i;    a[j3 + 3] = wk3i * x0i - wk3r * x0r;}#ifdef USE_CDFT_THREADSstruct cdft_arg_st {    int n0;    int n;    double *a;    int nw;    double *w;};typedef struct cdft_arg_st cdft_arg_t;void cftrec4_th(int n, double *a, int nw, double *w){    void *cftrec1_th(void *p);    void *cftrec2_th(void *p);    int i, idiv4, m, nthread;    cdft_thread_t th[4];    cdft_arg_t ag[4];        nthread = 2;    idiv4 = 0;    m = n >> 1;    if (n > CDFT_4THREADS_BEGIN_N) {        nthread = 4;        idiv4 = 1;        m >>= 1;    }    for (i = 0; i < nthread; i++) {        ag[i].n0 = n;        ag[i].n = m;        ag[i].a = &a[i * m];        ag[i].nw = nw;        ag[i].w = w;        if (i != idiv4) {            cdft_thread_create(&th[i], cftrec1_th, &ag[i]);        } else {            cdft_thread_create(&th[i], cftrec2_th, &ag[i]);        }    }    for (i = 0; i < nthread; i++) {        cdft_thread_wait(th[i]);    }}void *cftrec1_th(void *p){    int cfttree(int n, int j, int k, double *a, int nw, double *w);    void cftleaf(int n, int isplt, double *a, int nw, double *w);    void cftmdl1(int n, double *a, double *w);    int isplt, j, k, m, n, n0, nw;    double *a, *w;        n0 = ((cdft_arg_t *) p)->n0;    n = ((cdft_arg_t *) p)->n;    a = ((cdft_arg_t *) p)->a;    nw = ((cdft_arg_t *) p)->nw;    w = ((cdft_arg_t *) p)->w;    m = n0;    while (m > 512) {        m >>= 2;        cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);    }    cftleaf(m, 1, &a[n - m], nw, w);    k = 0;    for (j = n - m; j > 0; j -= m) {        k++;        isplt = cfttree(m, j, k, a, nw, w);        cftleaf(m, isplt, &a[j - m], nw, w);    }    return (void *) 0;}void *cftrec2_th(void *p){    int cfttree(int n, int j, int k, double *a, int nw, double *w);    void cftleaf(int n, int isplt, double *a, int nw, double *w);    void cftmdl2(int n, double *a, double *w);    int isplt, j, k, m, n, n0, nw;    double *a, *w;        n0 = ((cdft_arg_t *) p)->n0;    n = ((cdft_arg_t *) p)->n;    a = ((cdft_arg_t *) p)->a;    nw = ((cdft_arg_t *) p)->nw;    w = ((cdft_arg_t *) p)->w;    k = 1;    m = n0;    while (m > 512) {        m >>= 2;        k <<= 2;        cftmdl2(m, &a[n - m], &w[nw - m]);    }    cftleaf(m, 0, &a[n - m], nw, w);    k >>= 1;    for (j = n - m; j > 0; j -= m) {        k++;        isplt = cfttree(m, j, k, a, nw, w);        cftleaf(m, isplt, &a[j - m], nw, w);    }    return (void *) 0;}#endif /* USE_CDFT_THREADS */void cftrec4(int n, double *a, int nw, double *w){    int cfttree(int n, int j, int k, double *a, int nw, double *w);    void cftleaf(int n, int isplt, double *a, int nw, double *w);    void cftmdl1(int n, double *a, double *w);    int isplt, j, k, m;        m = n;    while (m > 512) {        m >>= 2;        cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);    }    cftleaf(m, 1, &a[n - m], nw, w);    k = 0;    for (j = n - m; j > 0; j -= m) {        k++;        isplt = cfttree(m, j, k, a, nw, w);        cftleaf(m, isplt, &a[j - m], nw, w);    }}int cfttree(int n, int j, int k, double *a, int nw, double *w){    void cftmdl1(int n, double *a, double *w);    void cftmdl2(int n, double *a, double *w);    int i, isplt, m;        if ((k & 3) != 0) {        isplt = k & 1;        if (isplt != 0) {            cftmdl1(n, &a[j - n], &w[nw - (n >> 1)]);        } else {            cftmdl2(n, &a[j - n], &w[nw - n]);        }    } else {        m = n;        for (i = k; (i & 3) == 0; i >>= 2) {            m <<= 2;        }        isplt = i & 1;        if (isplt != 0) {            while (m > 128) {                cftmdl1(m, &a[j - m], &w[nw - (m >> 1)]);                m >>= 2;            }        } else {            while (m > 128) {                cftmdl2(m, &a[j - m], &w[nw - m]);                m >>= 2;            }        }    }    return isplt;}void cftleaf(int n, int isplt, double *a, int nw, double *w){    void cftmdl1(int n, double *a, double *w);    void cftmdl2(int n, double *a, double *w);