#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
//#include <stdlib.h>
//#include <math.h>
//#include <stdio.h>
//#include "fft_czt.h"
typedef double REAL;
	/* FFT計算用 数表保持構造体 */
typedef struct {
	int  samples;		/* 標本点の数( 4以上の 2の整数乗に限る) */
	int  *bitrev;		/* ビット反転表 - 要素数は(samples) */
	REAL *sintbl;		/* 三角関数表 - 要素数は(samples*3/4) */
} fft_struct;
	/* CZT計算用 数表保持構造体 */
typedef struct {
	fft_struct fft;		/* 下請け FFT計算用 */
	int  samples;		/* 標本点の数 */
	int  samples_out;	/* 出力する標本点の数 */
	int  samples_ex;	/* 2の整数乗で、(samples + samples_out) <= samples_ex
						 * である最小の数 */
	REAL *wr, *wi;		/* 重みデータ - 要素数は(samples) */
	REAL *vr, *vi;		/* インパルス応答データ - 要素数は(samples_ex) */
	REAL *tr, *ti;		/* 作業用領域 - 要素数は(samples_ex) */
} czt_struct;
	/* fft.c */
int fft_init(fft_struct *fftp, int n);
void fft_end(fft_struct *fftp);
//void fft(fft_struct *fftp, int inv, REAL re[], REAL im[]);
void fft(double*,double*,int,int);
void half_fft(double*,double*,int);
void half_ifft(double*,double*,int);
	/* czt.c */
int czt_init(czt_struct *cztp, int n, int no,int *ip,double *ww);
void czt_end(czt_struct *cztp);
void czt(czt_struct *cztp, int inv, REAL re[], REAL im[],int *ip,double *ww);
    /* czttest.c*/
int charp_FFT(int,double*,double*,double*,double*,int *ip, double *ww);
int charp_IFFT(int,double*,double*,double*,double*,int *ip, double *ww);
/*
Fast Fourier/Cosine/Sine Transform
    dimension   :one
    data length :power of 2
    decimation  :frequency
    radix       :split-radix
    data        :inplace
    table       :use
functions
    cdft: Complex Discrete Fourier Transform
    rdft: Real Discrete Fourier Transform
    ddct: Discrete Cosine Transform
    ddst: Discrete Sine Transform
    dfct: Cosine Transform of RDFT (Real Symmetric DFT)
    dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
function prototypes
    void cdft(int, int, double *, int *, double *);
    void rdft(int, int, double *, int *, double *);
    void ddct(int, int, double *, int *, double *);
    void ddst(int, int, double *, int *, double *);
    void dfct(int, double *, double *, int *, double *);
    void dfst(int, double *, double *, int *, double *);
macro definitions
    USE_CDFT_PTHREADS : default=not defined
        CDFT_THREADS_BEGIN_N  : must be >= 512, default=8192
        CDFT_4THREADS_BEGIN_N : must be >= 512, default=65536
    USE_CDFT_WINTHREADS : default=not defined
        CDFT_THREADS_BEGIN_N  : must be >= 512, default=32768
        CDFT_4THREADS_BEGIN_N : must be >= 512, default=524288


-------- Complex DFT (Discrete Fourier Transform) --------
    [definition]
        <case1>
            X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
        <case2>
            X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
        (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
    [usage]
        <case1>
            ip[0] = 0; // first time only
            cdft(2*n, 1, a, ip, w);
        <case2>
            ip[0] = 0; // first time only
            cdft(2*n, -1, a, ip, w);
    [parameters]
        2*n            :data length (int)
                        n >= 1, n = power of 2
        a[0...2*n-1]   :input/output data (double *)
                        input data
                            a[2*j] = Re(x[j]), 
                            a[2*j+1] = Im(x[j]), 0<=j<n
                        output data
                            a[2*k] = Re(X[k]), 
                            a[2*k+1] = Im(X[k]), 0<=k<n
        ip[0...*]      :work area for bit reversal (int *)
                        length of ip >= 2+sqrt(n)
                        strictly, 
                        length of ip >= 
                            2+(1<<(int)(log(n+0.5)/log(2))/2).
                        ip[0],ip[1] are pointers of the cos/sin table.
        w[0...n/2-1]   :cos/sin table (double *)
                        w[],ip[] are initialized if ip[0] == 0.
    [remark]
        Inverse of 
            cdft(2*n, -1, a, ip, w);
        is 
            cdft(2*n, 1, a, ip, w);
            for (j = 0; j <= 2 * n - 1; j++) {
                a[j] *= 1.0 / n;
            }
        .


-------- Real DFT / Inverse of Real DFT --------
    [definition]
        <case1> RDFT
            R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
            I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
        <case2> IRDFT (excluding scale)
            a[k] = (R[0] + R[n/2]*cos(pi*k))/2 + 
                   sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + 
                   sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
    [usage]
        <case1>
            ip[0] = 0; // first time only
            rdft(n, 1, a, ip, w);
        <case2>
            ip[0] = 0; // first time only
            rdft(n, -1, a, ip, w);
    [parameters]
        n              :data length (int)
                        n >= 2, n = power of 2
        a[0...n-1]     :input/output data (double *)
                        <case1>
                            output data
                                a[2*k] = R[k], 0<=k<n/2
                                a[2*k+1] = I[k], 0<k<n/2
                                a[1] = R[n/2]
                        <case2>
                            input data
                                a[2*j] = R[j], 0<=j<n/2
                                a[2*j+1] = I[j], 0<j<n/2
                                a[1] = R[n/2]
        ip[0...*]      :work area for bit reversal (int *)
                        length of ip >= 2+sqrt(n/2)
                        strictly, 
                        length of ip >= 
                            2+(1<<(int)(log(n/2+0.5)/log(2))/2).
                        ip[0],ip[1] are pointers of the cos/sin table.
        w[0...n/2-1]   :cos/sin table (double *)
                        w[],ip[] are initialized if ip[0] == 0.
    [remark]
        Inverse of 
            rdft(n, 1, a, ip, w);
        is 
            rdft(n, -1, a, ip, w);
            for (j = 0; j <= n - 1; j++) {
                a[j] *= 2.0 / n;
            }
        .


-------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
    [definition]
        <case1> IDCT (excluding scale)
            C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
        <case2> DCT
            C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
    [usage]
        <case1>
            ip[0] = 0; // first time only
            ddct(n, 1, a, ip, w);
        <case2>
            ip[0] = 0; // first time only
            ddct(n, -1, a, ip, w);
    [parameters]
        n              :data length (int)
                        n >= 2, n = power of 2
        a[0...n-1]     :input/output data (double *)
                        output data
                            a[k] = C[k], 0<=k<n
        ip[0...*]      :work area for bit reversal (int *)
                        length of ip >= 2+sqrt(n/2)
                        strictly, 
                        length of ip >= 
                            2+(1<<(int)(log(n/2+0.5)/log(2))/2).
                        ip[0],ip[1] are pointers of the cos/sin table.
        w[0...n*5/4-1] :cos/sin table (double *)
                        w[],ip[] are initialized if ip[0] == 0.
    [remark]
        Inverse of 
            ddct(n, -1, a, ip, w);
        is 
            a[0] *= 0.5;
            ddct(n, 1, a, ip, w);
            for (j = 0; j <= n - 1; j++) {
                a[j] *= 2.0 / n;
            }
        .


-------- DST (Discrete Sine Transform) / Inverse of DST --------
    [definition]
        <case1> IDST (excluding scale)
            S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
        <case2> DST
            S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
    [usage]
        <case1>
            ip[0] = 0; // first time only
            ddst(n, 1, a, ip, w);
        <case2>
            ip[0] = 0; // first time only
            ddst(n, -1, a, ip, w);
    [parameters]
        n              :data length (int)
                        n >= 2, n = power of 2
        a[0...n-1]     :input/output data (double *)
                        <case1>
                            input data
                                a[j] = A[j], 0<j<n
                                a[0] = A[n]
                            output data
                                a[k] = S[k], 0<=k<n
                        <case2>
                            output data
                                a[k] = S[k], 0<k<n
                                a[0] = S[n]
        ip[0...*]      :work area for bit reversal (int *)
                        length of ip >= 2+sqrt(n/2)
                        strictly, 
                        length of ip >= 
                            2+(1<<(int)(log(n/2+0.5)/log(2))/2).
                        ip[0],ip[1] are pointers of the cos/sin table.
        w[0...n*5/4-1] :cos/sin table (double *)
                        w[],ip[] are initialized if ip[0] == 0.
    [remark]
        Inverse of 
            ddst(n, -1, a, ip, w);
        is 
            a[0] *= 0.5;
            ddst(n, 1, a, ip, w);
            for (j = 0; j <= n - 1; j++) {
                a[j] *= 2.0 / n;
            }
        .


-------- Cosine Transform of RDFT (Real Symmetric DFT) --------
    [definition]
        C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
    [usage]
        ip[0] = 0; // first time only
        dfct(n, a, t, ip, w);
    [parameters]
        n              :data length - 1 (int)
                        n >= 2, n = power of 2
        a[0...n]       :input/output data (double *)
                        output data
                            a[k] = C[k], 0<=k<=n
        t[0...n/2]     :work area (double *)
        ip[0...*]      :work area for bit reversal (int *)
                        length of ip >= 2+sqrt(n/4)
                        strictly, 
                        length of ip >= 
                            2+(1<<(int)(log(n/4+0.5)/log(2))/2).
                        ip[0],ip[1] are pointers of the cos/sin table.
        w[0...n*5/8-1] :cos/sin table (double *)
                        w[],ip[] are initialized if ip[0] == 0.
    [remark]
        Inverse of 
            a[0] *= 0.5;
            a[n] *= 0.5;
            dfct(n, a, t, ip, w);
        is 
            a[0] *= 0.5;
            a[n] *= 0.5;
            dfct(n, a, t, ip, w);
            for (j = 0; j <= n; j++) {
                a[j] *= 2.0 / n;
            }
        .


-------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
    [definition]
        S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
    [usage]
        ip[0] = 0; // first time only
        dfst(n, a, t, ip, w);
    [parameters]
        n              :data length + 1 (int)
                        n >= 2, n = power of 2
        a[0...n-1]     :input/output data (double *)
                        output data
                            a[k] = S[k], 0<k<n
                        (a[0] is used for work area)
        t[0...n/2-1]   :work area (double *)
        ip[0...*]      :work area for bit reversal (int *)
                        length of ip >= 2+sqrt(n/4)
                        strictly, 
                        length of ip >= 
                            2+(1<<(int)(log(n/4+0.5)/log(2))/2).
                        ip[0],ip[1] are pointers of the cos/sin table.
        w[0...n*5/8-1] :cos/sin table (double *)
                        w[],ip[] are initialized if ip[0] == 0.
    [remark]
        Inverse of 
            dfst(n, a, t, ip, w);
        is 
            dfst(n, a, t, ip, w);
            for (j = 1; j <= n - 1; j++) {
                a[j] *= 2.0 / n;
            }
        .


Appendix :
    The cos/sin table is recalculated when the larger table required.
    w[] and ip[] are compatible with all routines.
*/
//#include "stdafx.h"

void cdft(int n, int isgn, double *a, int *ip, double *w)
{
    void makewt(int nw, int *ip, double *w);
    void cftfsub(int n, double *a, int *ip, int nw, double *w);
    void cftbsub(int n, double *a, int *ip, int nw, double *w);
    int nw;
    
    nw = ip[0];
    if (n > (nw << 2)) {
        nw = n >> 2;
        makewt(nw, ip, w);
    }
    if (isgn >= 0) {
        cftfsub(n, a, ip, nw, w);
    } else {
        cftbsub(n, a, ip, nw, w);
    }
}


void rdft(int n, int isgn, double *a, int *ip, double *w)
{
    void makewt(int nw, int *ip, double *w);
    void makect(int nc, int *ip, double *c);
    void cftfsub(int n, double *a, int *ip, int nw, double *w);
    void cftbsub(int n, double *a, int *ip, int nw, double *w);
    void rftfsub(int n, double *a, int nc, double *c);
    void rftbsub(int n, double *a, int nc, double *c);
    int nw, nc;
    double xi;
    
    nw = ip[0];
    if (n > (nw << 2)) {
        nw = n >> 2;
        makewt(nw, ip, w);
    }
    nc = ip[1];
    if (n > (nc << 2)) {
        nc = n >> 2;
        makect(nc, ip, w + nw);
    }
    if (isgn >= 0) {
        if (n > 4) {
            cftfsub(n, a, ip, nw, w);
            rftfsub(n, a, nc, w + nw);
        } else if (n == 4) {
            cftfsub(n, a, ip, nw, w);
        }
        xi = a[0] - a[1];
        a[0] += a[1];
        a[1] = xi;
    } else {
        a[1] = 0.5 * (a[0] - a[1]);
        a[0] -= a[1];
        if (n > 4) {
            rftbsub(n, a, nc, w + nw);
            cftbsub(n, a, ip, nw, w);
        } else if (n == 4) {
            cftbsub(n, a, ip, nw, w);
        }
    }
}


void ddct(int n, int isgn, double *a, int *ip, double *w)
{
    void makewt(int nw, int *ip, double *w);
    void makect(int nc, int *ip, double *c);
    void cftfsub(int n, double *a, int *ip, int nw, double *w);
    void cftbsub(int n, double *a, int *ip, int nw, double *w);
    void rftfsub(int n, double *a, int nc, double *c);
    void rftbsub(int n, double *a, int nc, double *c);
    void dctsub(int n, double *a, int nc, double *c);
    int j, nw, nc;
    double xr;
    
    nw = ip[0];
    if (n > (nw << 2)) {
        nw = n >> 2;
        makewt(nw, ip, w);
    }
    nc = ip[1];
    if (n > nc) {
        nc = n;
        makect(nc, ip, w + nw);
    }
    if (isgn < 0) {
        xr = a[n - 1];
        for (j = n - 2; j >= 2; j -= 2) {
            a[j + 1] = a[j] - a[j - 1];
            a[j] += a[j - 1];
        }
        a[1] = a[0] - xr;
        a[0] += xr;
        if (n > 4) {
            rftbsub(n, a, nc, w + nw);
            cftbsub(n, a, ip, nw, w);
        } else if (n == 4) {
            cftbsub(n, a, ip, nw, w);
        }
    }
    dctsub(n, a, nc, w + nw);
    if (isgn >= 0) {
        if (n > 4) {
            cftfsub(n, a, ip, nw, w);
            rftfsub(n, a, nc, w + nw);
        } else if (n == 4) {
            cftfsub(n, a, ip, nw, w);
        }
        xr = a[0] - a[1];
        a[0] += a[1];
        for (j = 2; j < n; j += 2) {
            a[j - 1] = a[j] - a[j + 1];
            a[j] += a[j + 1];
        }
        a[n - 1] = xr;
    }
}


void ddst(int n, int isgn, double *a, int *ip, double *w)
{
    void makewt(int nw, int *ip, double *w);
    void makect(int nc, int *ip, double *c);
    void cftfsub(int n, double *a, int *ip, int nw, double *w);
    void cftbsub(int n, double *a, int *ip, int nw, double *w);
    void rftfsub(int n, double *a, int nc, double *c);
    void rftbsub(int n, double *a, int nc, double *c);
    void dstsub(int n, double *a, int nc, double *c);
    int j, nw, nc;
    double xr;
    
    nw = ip[0];
    if (n > (nw << 2)) {
        nw = n >> 2;
        makewt(nw, ip, w);
    }
    nc = ip[1];
    if (n > nc) {
        nc = n;
        makect(nc, ip, w + nw);
    }
    if (isgn < 0) {
        xr = a[n - 1];
        for (j = n - 2; j >= 2; j -= 2) {
            a[j + 1] = -a[j] - a[j - 1];
            a[j] -= a[j - 1];
        }
        a[1] = a[0] + xr;
        a[0] -= xr;
        if (n > 4) {
            rftbsub(n, a, nc, w + nw);
            cftbsub(n, a, ip, nw, w);
        } else if (n == 4) {
            cftbsub(n, a, ip, nw, w);
        }
    }
    dstsub(n, a, nc, w + nw);
    if (isgn >= 0) {
        if (n > 4) {
            cftfsub(n, a, ip, nw, w);
            rftfsub(n, a, nc, w + nw);
        } else if (n == 4) {
            cftfsub(n, a, ip, nw, w);
        }
        xr = a[0] - a[1];
        a[0] += a[1];
        for (j = 2; j < n; j += 2) {
            a[j - 1] = -a[j] - a[j + 1];
            a[j] -= a[j + 1];
        }
        a[n - 1] = -xr;
    }
}


void dfct(int n, double *a, double *t, int *ip, double *w)
{
    void makewt(int nw, int *ip, double *w);
    void makect(int nc, int *ip, double *c);
    void cftfsub(int n, double *a, int *ip, int nw, double *w);
    void rftfsub(int n, double *a, int nc, double *c);
    void dctsub(int n, double *a, int nc, double *c);
    int j, k, l, m, mh, nw, nc;
    double xr, xi, yr, yi;
    
    nw = ip[0];
    if (n > (nw << 3)) {
        nw = n >> 3;
        makewt(nw, ip, w);
    }
    nc = ip[1];
    if (n > (nc << 1)) {
        nc = n >> 1;
        makect(nc, ip, w + nw);
    }
    m = n >> 1;
    yi = a[m];
    xi = a[0] + a[n];
    a[0] -= a[n];
    t[0] = xi - yi;
    t[m] = xi + yi;
    if (n > 2) {
        mh = m >> 1;
        for (j = 1; j < mh; j++) {
            k = m - j;
            xr = a[j] - a[n - j];
            xi = a[j] + a[n - j];
            yr = a[k] - a[n - k];
            yi = a[k] + a[n - k];
            a[j] = xr;
            a[k] = yr;
            t[j] = xi - yi;
            t[k] = xi + yi;
        }
        t[mh] = a[mh] + a[n - mh];
        a[mh] -= a[n - mh];
        dctsub(m, a, nc, w + nw);
        if (m > 4) {
            cftfsub(m, a, ip, nw, w);
            rftfsub(m, a, nc, w + nw);
        } else if (m == 4) {
            cftfsub(m, a, ip, nw, w);
        }
        a[n - 1] = a[0] - a[1];
        a[1] = a[0] + a[1];
        for (j = m - 2; j >= 2; j -= 2) {
            a[2 * j + 1] = a[j] + a[j + 1];
            a[2 * j - 1] = a[j] - a[j + 1];
        }
        l = 2;
        m = mh;
        while (m >= 2) {
            dctsub(m, t, nc, w + nw);
            if (m > 4) {
                cftfsub(m, t, ip, nw, w);
                rftfsub(m, t, nc, w + nw);
            } else if (m == 4) {
                cftfsub(m, t, ip, nw, w);
            }
            a[n - l] = t[0] - t[1];
            a[l] = t[0] + t[1];
            k = 0;
            for (j = 2; j < m; j += 2) {
                k += l << 2;
                a[k - l] = t[j] - t[j + 1];
                a[k + l] = t[j] + t[j + 1];
            }
            l <<= 1;
            mh = m >> 1;
            for (j = 0; j < mh; j++) {
                k = m - j;
                t[j] = t[m + k] - t[m + j];
                t[k] = t[m + k] + t[m + j];
            }
            t[mh] = t[m + mh];
            m = mh;
        }
        a[l] = t[0];
        a[n] = t[2] - t[1];
        a[0] = t[2] + t[1];
    } else {
        a[1] = a[0];
        a[2] = t[0];
        a[0] = t[1];
    }
}


void dfst(int n, double *a, double *t, int *ip, double *w)
{
    void makewt(int nw, int *ip, double *w);
    void makect(int nc, int *ip, double *c);
    void cftfsub(int n, double *a, int *ip, int nw, double *w);
    void rftfsub(int n, double *a, int nc, double *c);
    void dstsub(int n, double *a, int nc, double *c);
    int j, k, l, m, mh, nw, nc;
    double xr, xi, yr, yi;
    
    nw = ip[0];
    if (n > (nw << 3)) {
        nw = n >> 3;
        makewt(nw, ip, w);
    }
    nc = ip[1];
    if (n > (nc << 1)) {
        nc = n >> 1;
        makect(nc, ip, w + nw);
    }
    if (n > 2) {
        m = n >> 1;
        mh = m >> 1;
        for (j = 1; j < mh; j++) {
            k = m - j;
            xr = a[j] + a[n - j];
            xi = a[j] - a[n - j];
            yr = a[k] + a[n - k];
            yi = a[k] - a[n - k];
            a[j] = xr;
            a[k] = yr;
            t[j] = xi + yi;
            t[k] = xi - yi;
        }
        t[0] = a[mh] - a[n - mh];
        a[mh] += a[n - mh];
        a[0] = a[m];
        dstsub(m, a, nc, w + nw);
        if (m > 4) {
            cftfsub(m, a, ip, nw, w);
            rftfsub(m, a, nc, w + nw);
        } else if (m == 4) {
            cftfsub(m, a, ip, nw, w);
        }
        a[n - 1] = a[1] - a[0];
        a[1] = a[0] + a[1];
        for (j = m - 2; j >= 2; j -= 2) {
            a[2 * j + 1] = a[j] - a[j + 1];
            a[2 * j - 1] = -a[j] - a[j + 1];
        }
        l = 2;
        m = mh;
        while (m >= 2) {
            dstsub(m, t, nc, w + nw);
            if (m > 4) {
                cftfsub(m, t, ip, nw, w);
                rftfsub(m, t, nc, w + nw);
            } else if (m == 4) {
                cftfsub(m, t, ip, nw, w);
            }
            a[n - l] = t[1] - t[0];
            a[l] = t[0] + t[1];
            k = 0;
            for (j = 2; j < m; j += 2) {
                k += l << 2;
                a[k - l] = -t[j] - t[j + 1];
                a[k + l] = t[j] - t[j + 1];
            }
            l <<= 1;
            mh = m >> 1;
            for (j = 1; j < mh; j++) {
                k = m - j;
                t[j] = t[m + k] + t[m + j];
                t[k] = t[m + k] - t[m + j];
            }
            t[0] = t[m + mh];
            m = mh;
        }
        a[l] = t[0];
    }
    a[0] = 0;
}


/* -------- initializing routines -------- */


#include <math.h>

void makewt(int nw, int *ip, double *w)
{
    void makeipt(int nw, int *ip);
    int j, nwh, nw0, nw1;
    double delta, wn4r, wk1r, wk1i, wk3r, wk3i;
    
    ip[0] = nw;
    ip[1] = 1;
    if (nw > 2) {
        nwh = nw >> 1;
        delta = atan(1.0) / nwh;
        wn4r = cos(delta * nwh);
        w[0] = 1;
        w[1] = wn4r;
        if (nwh == 4) {
            w[2] = cos(delta * 2);
            w[3] = sin(delta * 2);
        } else if (nwh > 4) {
            makeipt(nw, ip);
            w[2] = 0.5 / cos(delta * 2);
            w[3] = 0.5 / cos(delta * 6);
            for (j = 4; j < nwh; j += 4) {
                w[j] = cos(delta * j);
                w[j + 1] = sin(delta * j);
                w[j + 2] = cos(3 * delta * j);
                w[j + 3] = -sin(3 * delta * j);
            }
        }
        nw0 = 0;
        while (nwh > 2) {
            nw1 = nw0 + nwh;
            nwh >>= 1;
            w[nw1] = 1;
            w[nw1 + 1] = wn4r;
            if (nwh == 4) {
                wk1r = w[nw0 + 4];
                wk1i = w[nw0 + 5];
                w[nw1 + 2] = wk1r;
                w[nw1 + 3] = wk1i;
            } else if (nwh > 4) {
                wk1r = w[nw0 + 4];
                wk3r = w[nw0 + 6];
                w[nw1 + 2] = 0.5 / wk1r;
                w[nw1 + 3] = 0.5 / wk3r;
                for (j = 4; j < nwh; j += 4) {
                    wk1r = w[nw0 + 2 * j];
                    wk1i = w[nw0 + 2 * j + 1];
                    wk3r = w[nw0 + 2 * j + 2];
                    wk3i = w[nw0 + 2 * j + 3];
                    w[nw1 + j] = wk1r;
                    w[nw1 + j + 1] = wk1i;
                    w[nw1 + j + 2] = wk3r;
                    w[nw1 + j + 3] = wk3i;
                }
            }
            nw0 = nw1;
        }
    }
}


void makeipt(int nw, int *ip)
{
    int j, l, m, m2, p, q;
    
    ip[2] = 0;
    ip[3] = 16;
    m = 2;
    for (l = nw; l > 32; l >>= 2) {
        m2 = m << 1;
        q = m2 << 3;
        for (j = m; j < m2; j++) {
            p = ip[j] << 2;
            ip[m + j] = p;
            ip[m2 + j] = p + q;
        }
        m = m2;
    }
}


void makect(int nc, int *ip, double *c)
{
    int j, nch;
    double delta;
    
    ip[1] = nc;
    if (nc > 1) {
        nch = nc >> 1;
        delta = atan(1.0) / nch;
        c[0] = cos(delta * nch);
        c[nch] = 0.5 * c[0];
        for (j = 1; j < nch; j++) {
            c[j] = 0.5 * cos(delta * j);
            c[nc - j] = 0.5 * sin(delta * j);
        }
    }
}


/* -------- child routines -------- */


#ifdef USE_CDFT_PTHREADS
#define USE_CDFT_THREADS
#ifndef CDFT_THREADS_BEGIN_N
#define CDFT_THREADS_BEGIN_N 8192
#endif
#ifndef CDFT_4THREADS_BEGIN_N
#define CDFT_4THREADS_BEGIN_N 65536
#endif
#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#define cdft_thread_t pthread_t
#define cdft_thread_create(thp,func,argp) { \
    if (pthread_create(thp, NULL, func, (void *) argp) != 0) { \
        fprintf(stderr, "cdft thread error\n"); \
        exit(1); \
    } \
}
#define cdft_thread_wait(th) { \
    if (pthread_join(th, NULL) != 0) { \
        fprintf(stderr, "cdft thread error\n"); \
        exit(1); \
    } \
}
#endif /* USE_CDFT_PTHREADS */


#ifdef USE_CDFT_WINTHREADS
#define USE_CDFT_THREADS
#ifndef CDFT_THREADS_BEGIN_N
#define CDFT_THREADS_BEGIN_N 32768
#endif
#ifndef CDFT_4THREADS_BEGIN_N
#define CDFT_4THREADS_BEGIN_N 524288
#endif
#include <windows.h>
#include <stdio.h>
#include <stdlib.h>
#define cdft_thread_t HANDLE
#define cdft_thread_create(thp,func,argp) { \
    DWORD thid; \
    *(thp) = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE) func, (LPVOID) argp, 0, &thid); \
    if (*(thp) == 0) { \
        fprintf(stderr, "cdft thread error\n"); \
        exit(1); \
    } \
}
#define cdft_thread_wait(th) { \
    WaitForSingleObject(th, INFINITE); \
    CloseHandle(th); \
}
#endif /* USE_CDFT_WINTHREADS */


void cftfsub(int n, double *a, int *ip, int nw, double *w)
{
    void bitrv2(int n, int *ip, double *a);
    void bitrv216(double *a);
    void bitrv208(double *a);
    void cftf1st(int n, double *a, double *w);
    void cftrec4(int n, double *a, int nw, double *w);
    void cftleaf(int n, int isplt, double *a, int nw, double *w);
    void cftfx41(int n, double *a, int nw, double *w);
    void cftf161(double *a, double *w);
    void cftf081(double *a, double *w);
    void cftf040(double *a);
    void cftx020(double *a);
#ifdef USE_CDFT_THREADS
    void cftrec4_th(int n, double *a, int nw, double *w);
#endif /* USE_CDFT_THREADS */
    
    if (n > 8) {
        if (n > 32) {
            cftf1st(n, a, &w[nw - (n >> 2)]);
#ifdef USE_CDFT_THREADS
            if (n > CDFT_THREADS_BEGIN_N) {
                cftrec4_th(n, a, nw, w);
            } else 
#endif /* USE_CDFT_THREADS */
            if (n > 512) {
                cftrec4(n, a, nw, w);
            } else if (n > 128) {
                cftleaf(n, 1, a, nw, w);
            } else {
                cftfx41(n, a, nw, w);
            }
            bitrv2(n, ip, a);
        } else if (n == 32) {
            cftf161(a, &w[nw - 8]);
            bitrv216(a);
        } else {
            cftf081(a, w);
            bitrv208(a);
        }
    } else if (n == 8) {
        cftf040(a);
    } else if (n == 4) {
        cftx020(a);
    }
}


void cftbsub(int n, double *a, int *ip, int nw, double *w)
{
    void bitrv2conj(int n, int *ip, double *a);
    void bitrv216neg(double *a);
    void bitrv208neg(double *a);
    void cftb1st(int n, double *a, double *w);
    void cftrec4(int n, double *a, int nw, double *w);
    void cftleaf(int n, int isplt, double *a, int nw, double *w);
    void cftfx41(int n, double *a, int nw, double *w);
    void cftf161(double *a, double *w);
    void cftf081(double *a, double *w);
    void cftb040(double *a);
    void cftx020(double *a);
#ifdef USE_CDFT_THREADS
    void cftrec4_th(int n, double *a, int nw, double *w);
#endif /* USE_CDFT_THREADS */
    
    if (n > 8) {
        if (n > 32) {
            cftb1st(n, a, &w[nw - (n >> 2)]);
#ifdef USE_CDFT_THREADS
            if (n > CDFT_THREADS_BEGIN_N) {
                cftrec4_th(n, a, nw, w);
            } else 
#endif /* USE_CDFT_THREADS */
            if (n > 512) {
                cftrec4(n, a, nw, w);
            } else if (n > 128) {
                cftleaf(n, 1, a, nw, w);
            } else {
                cftfx41(n, a, nw, w);
            }
            bitrv2conj(n, ip, a);
        } else if (n == 32) {
            cftf161(a, &w[nw - 8]);
            bitrv216neg(a);
        } else {
            cftf081(a, w);
            bitrv208neg(a);
        }
    } else if (n == 8) {
        cftb040(a);
    } else if (n == 4) {
        cftx020(a);
    }
}


void bitrv2(int n, int *ip, double *a)
{
    int j, j1, k, k1, l, m, nh, nm;
    double xr, xi, yr, yi;
    
    m = 1;
    for (l = n >> 2; l > 8; l >>= 2) {
        m <<= 1;
    }
    nh = n >> 1;
    nm = 4 * m;
    if (l == 8) {
        for (k = 0; k < m; k++) {
            for (j = 0; j < k; j++) {
                j1 = 4 * j + 2 * ip[m + k];
                k1 = 4 * k + 2 * ip[m + j];
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 += 2 * nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 -= nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 += 2 * nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nh;
                k1 += 2;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 -= 2 * nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 += nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 -= 2 * nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += 2;
                k1 += nh;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 += 2 * nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 -= nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 += 2 * nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nh;
                k1 -= 2;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 -= 2 * nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 += nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 -= 2 * nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
            }
            k1 = 4 * k + 2 * ip[m + k];
            j1 = k1 + 2;
            k1 += nh;
            xr = a[j1];
            xi = a[j1 + 1];
            yr = a[k1];
            yi = a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
            j1 += nm;
            k1 += 2 * nm;
            xr = a[j1];
            xi = a[j1 + 1];
            yr = a[k1];
            yi = a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
            j1 += nm;
            k1 -= nm;
            xr = a[j1];
            xi = a[j1 + 1];
            yr = a[k1];
            yi = a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
            j1 -= 2;
            k1 -= nh;
            xr = a[j1];
            xi = a[j1 + 1];
            yr = a[k1];
            yi = a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
            j1 += nh + 2;
            k1 += nh + 2;
            xr = a[j1];
            xi = a[j1 + 1];
            yr = a[k1];
            yi = a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
            j1 -= nh - nm;
            k1 += 2 * nm - 2;
            xr = a[j1];
            xi = a[j1 + 1];
            yr = a[k1];
            yi = a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
        }
    } else {
        for (k = 0; k < m; k++) {
            for (j = 0; j < k; j++) {
                j1 = 4 * j + ip[m + k];
                k1 = 4 * k + ip[m + j];
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 += nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nh;
                k1 += 2;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 -= nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += 2;
                k1 += nh;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 += nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nh;
                k1 -= 2;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 -= nm;
                xr = a[j1];
                xi = a[j1 + 1];
                yr = a[k1];
                yi = a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
            }
            k1 = 4 * k + ip[m + k];
            j1 = k1 + 2;
            k1 += nh;
            xr = a[j1];
            xi = a[j1 + 1];
            yr = a[k1];
            yi = a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
            j1 += nm;
            k1 += nm;
            xr = a[j1];
            xi = a[j1 + 1];
            yr = a[k1];
            yi = a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
        }
    }
}


void bitrv2conj(int n, int *ip, double *a)
{
    int j, j1, k, k1, l, m, nh, nm;
    double xr, xi, yr, yi;
    
    m = 1;
    for (l = n >> 2; l > 8; l >>= 2) {
        m <<= 1;
    }
    nh = n >> 1;
    nm = 4 * m;
    if (l == 8) {
        for (k = 0; k < m; k++) {
            for (j = 0; j < k; j++) {
                j1 = 4 * j + 2 * ip[m + k];
                k1 = 4 * k + 2 * ip[m + j];
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 += 2 * nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 -= nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 += 2 * nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nh;
                k1 += 2;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 -= 2 * nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 += nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 -= 2 * nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += 2;
                k1 += nh;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 += 2 * nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 -= nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 += 2 * nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nh;
                k1 -= 2;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 -= 2 * nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 += nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 -= 2 * nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
            }
            k1 = 4 * k + 2 * ip[m + k];
            j1 = k1 + 2;
            k1 += nh;
            a[j1 - 1] = -a[j1 - 1];
            xr = a[j1];
            xi = -a[j1 + 1];
            yr = a[k1];
            yi = -a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
            a[k1 + 3] = -a[k1 + 3];
            j1 += nm;
            k1 += 2 * nm;
            xr = a[j1];
            xi = -a[j1 + 1];
            yr = a[k1];
            yi = -a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
            j1 += nm;
            k1 -= nm;
            xr = a[j1];
            xi = -a[j1 + 1];
            yr = a[k1];
            yi = -a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
            j1 -= 2;
            k1 -= nh;
            xr = a[j1];
            xi = -a[j1 + 1];
            yr = a[k1];
            yi = -a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
            j1 += nh + 2;
            k1 += nh + 2;
            xr = a[j1];
            xi = -a[j1 + 1];
            yr = a[k1];
            yi = -a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
            j1 -= nh - nm;
            k1 += 2 * nm - 2;
            a[j1 - 1] = -a[j1 - 1];
            xr = a[j1];
            xi = -a[j1 + 1];
            yr = a[k1];
            yi = -a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
            a[k1 + 3] = -a[k1 + 3];
        }
    } else {
        for (k = 0; k < m; k++) {
            for (j = 0; j < k; j++) {
                j1 = 4 * j + ip[m + k];
                k1 = 4 * k + ip[m + j];
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 += nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nh;
                k1 += 2;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 -= nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += 2;
                k1 += nh;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += nm;
                k1 += nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nh;
                k1 -= 2;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 -= nm;
                k1 -= nm;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
            }
            k1 = 4 * k + ip[m + k];
            j1 = k1 + 2;
            k1 += nh;
            a[j1 - 1] = -a[j1 - 1];
            xr = a[j1];
            xi = -a[j1 + 1];
            yr = a[k1];
            yi = -a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
            a[k1 + 3] = -a[k1 + 3];
            j1 += nm;
            k1 += nm;
            a[j1 - 1] = -a[j1 - 1];
            xr = a[j1];
            xi = -a[j1 + 1];
            yr = a[k1];
            yi = -a[k1 + 1];
            a[j1] = yr;
            a[j1 + 1] = yi;
            a[k1] = xr;
            a[k1 + 1] = xi;
            a[k1 + 3] = -a[k1 + 3];
        }
    }
}


void bitrv216(double *a)
{
    double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
        x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i, 
        x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i;
    
    x1r = a[2];
    x1i = a[3];
    x2r = a[4];
    x2i = a[5];
    x3r = a[6];
    x3i = a[7];
    x4r = a[8];
    x4i = a[9];
    x5r = a[10];
    x5i = a[11];
    x7r = a[14];
    x7i = a[15];
    x8r = a[16];
    x8i = a[17];
    x10r = a[20];
    x10i = a[21];
    x11r = a[22];
    x11i = a[23];
    x12r = a[24];
    x12i = a[25];
    x13r = a[26];
    x13i = a[27];
    x14r = a[28];
    x14i = a[29];
    a[2] = x8r;
    a[3] = x8i;
    a[4] = x4r;
    a[5] = x4i;
    a[6] = x12r;
    a[7] = x12i;
    a[8] = x2r;
    a[9] = x2i;
    a[10] = x10r;
    a[11] = x10i;
    a[14] = x14r;
    a[15] = x14i;
    a[16] = x1r;
    a[17] = x1i;
    a[20] = x5r;
    a[21] = x5i;
    a[22] = x13r;
    a[23] = x13i;
    a[24] = x3r;
    a[25] = x3i;
    a[26] = x11r;
    a[27] = x11i;
    a[28] = x7r;
    a[29] = x7i;
}


void bitrv216neg(double *a)
{
    double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
        x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i, 
        x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i, 
        x13r, x13i, x14r, x14i, x15r, x15i;
    
    x1r = a[2];
    x1i = a[3];
    x2r = a[4];
    x2i = a[5];
    x3r = a[6];
    x3i = a[7];
    x4r = a[8];
    x4i = a[9];
    x5r = a[10];
    x5i = a[11];
    x6r = a[12];
    x6i = a[13];
    x7r = a[14];
    x7i = a[15];
    x8r = a[16];
    x8i = a[17];
    x9r = a[18];
    x9i = a[19];
    x10r = a[20];
    x10i = a[21];
    x11r = a[22];
    x11i = a[23];
    x12r = a[24];
    x12i = a[25];
    x13r = a[26];
    x13i = a[27];
    x14r = a[28];
    x14i = a[29];
    x15r = a[30];
    x15i = a[31];
    a[2] = x15r;
    a[3] = x15i;
    a[4] = x7r;
    a[5] = x7i;
    a[6] = x11r;
    a[7] = x11i;
    a[8] = x3r;
    a[9] = x3i;
    a[10] = x13r;
    a[11] = x13i;
    a[12] = x5r;
    a[13] = x5i;
    a[14] = x9r;
    a[15] = x9i;
    a[16] = x1r;
    a[17] = x1i;
    a[18] = x14r;
    a[19] = x14i;
    a[20] = x6r;
    a[21] = x6i;
    a[22] = x10r;
    a[23] = x10i;
    a[24] = x2r;
    a[25] = x2i;
    a[26] = x12r;
    a[27] = x12i;
    a[28] = x4r;
    a[29] = x4i;
    a[30] = x8r;
    a[31] = x8i;
}


void bitrv208(double *a)
{
    double x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
    
    x1r = a[2];
    x1i = a[3];
    x3r = a[6];
    x3i = a[7];
    x4r = a[8];
    x4i = a[9];
    x6r = a[12];
    x6i = a[13];
    a[2] = x4r;
    a[3] = x4i;
    a[6] = x6r;
    a[7] = x6i;
    a[8] = x1r;
    a[9] = x1i;
    a[12] = x3r;
    a[13] = x3i;
}


void bitrv208neg(double *a)
{
    double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
        x5r, x5i, x6r, x6i, x7r, x7i;
    
    x1r = a[2];
    x1i = a[3];
    x2r = a[4];
    x2i = a[5];
    x3r = a[6];
    x3i = a[7];
    x4r = a[8];
    x4i = a[9];
    x5r = a[10];
    x5i = a[11];
    x6r = a[12];
    x6i = a[13];
    x7r = a[14];
    x7i = a[15];
    a[2] = x7r;
    a[3] = x7i;
    a[4] = x3r;
    a[5] = x3i;
    a[6] = x5r;
    a[7] = x5i;
    a[8] = x1r;
    a[9] = x1i;
    a[10] = x6r;
    a[11] = x6i;
    a[12] = x2r;
    a[13] = x2i;
    a[14] = x4r;
    a[15] = x4i;
}


void cftf1st(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;
}


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_THREADS
struct 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);
    void cftf161(double *a, double *w);
    void cftf162(double *a, double *w);
    void cftf081(double *a, double *w);
    void cftf082(double *a, double *w);
    
    if (n == 512) {
        cftmdl1(128, a, &w[nw - 64]);
        cftf161(a, &w[nw - 8]);
        cftf162(&a[32], &w[nw - 32]);
        cftf161(&a[64], &w[nw - 8]);
        cftf161(&a[96], &w[nw - 8]);
        cftmdl2(128, &a[128], &w[nw - 128]);
        cftf161(&a[128], &w[nw - 8]);
        cftf162(&a[160], &w[nw - 32]);
        cftf161(&a[192], &w[nw - 8]);
        cftf162(&a[224], &w[nw - 32]);
        cftmdl1(128, &a[256], &w[nw - 64]);
        cftf161(&a[256], &w[nw - 8]);
        cftf162(&a[288], &w[nw - 32]);
        cftf161(&a[320], &w[nw - 8]);
        cftf161(&a[352], &w[nw - 8]);
        if (isplt != 0) {
            cftmdl1(128, &a[384], &w[nw - 64]);
            cftf161(&a[480], &w[nw - 8]);
        } else {
            cftmdl2(128, &a[384], &w[nw - 128]);
            cftf162(&a[480], &w[nw - 32]);
        }
        cftf161(&a[384], &w[nw - 8]);
        cftf162(&a[416], &w[nw - 32]);
        cftf161(&a[448], &w[nw - 8]);
    } else {
        cftmdl1(64, a, &w[nw - 32]);
        cftf081(a, &w[nw - 8]);
        cftf082(&a[16], &w[nw - 8]);
        cftf081(&a[32], &w[nw - 8]);
        cftf081(&a[48], &w[nw - 8]);
        cftmdl2(64, &a[64], &w[nw - 64]);
        cftf081(&a[64], &w[nw - 8]);
        cftf082(&a[80], &w[nw - 8]);
        cftf081(&a[96], &w[nw - 8]);
        cftf082(&a[112], &w[nw - 8]);
        cftmdl1(64, &a[128], &w[nw - 32]);
        cftf081(&a[128], &w[nw - 8]);
        cftf082(&a[144], &w[nw - 8]);
        cftf081(&a[160], &w[nw - 8]);
        cftf081(&a[176], &w[nw - 8]);
        if (isplt != 0) {
            cftmdl1(64, &a[192], &w[nw - 32]);
            cftf081(&a[240], &w[nw - 8]);
        } else {
            cftmdl2(64, &a[192], &w[nw - 64]);
            cftf082(&a[240], &w[nw - 8]);
        }
        cftf081(&a[192], &w[nw - 8]);
        cftf082(&a[208], &w[nw - 8]);
        cftf081(&a[224], &w[nw - 8]);
    }
}


void cftmdl1(int n, double *a, double *w)
{
    int j, j0, j1, j2, j3, k, m, mh;
    double wn4r, wk1r, wk1i, wk3r, wk3i;
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
    
    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];
    k = 0;
    for (j = 2; j < mh; j += 2) {
        k += 4;
        wk1r = w[k];
        wk1i = w[k + 1];
        wk3r = w[k + 2];
        wk3i = 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];
        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[j] = x0r + x2r;
        a[j + 1] = x0i + x2i;
        a[j1] = x0r - x2r;
        a[j1 + 1] = x0i - x2i;
        x0r = x1r - x3i;
        x0i = x1i + x3r;
        a[j2] = wk1r * x0r - wk1i * x0i;
        a[j2 + 1] = wk1r * x0i + wk1i * x0r;
        x0r = x1r + x3i;
        x0i = x1i - x3r;
        a[j3] = wk3r * x0r + wk3i * x0i;
        a[j3 + 1] = wk3r * x0i - wk3i * 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];
        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] = wk1i * x0r - wk1r * x0i;
        a[j2 + 1] = wk1i * x0i + wk1r * x0r;
        x0r = x1r + x3i;
        x0i = x1i - x3r;
        a[j3] = wk3i * x0r + wk3r * x0i;
        a[j3 + 1] = wk3i * x0i - wk3r * x0r;
    }
    j0 = mh;
    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];
    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);
}


void cftmdl2(int n, double *a, double *w)
{
    int j, j0, j1, j2, j3, k, kr, m, mh;
    double wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i;
    
    mh = n >> 3;
    m = 2 * mh;
    wn4r = w[1];
    j1 = m;
    j2 = j1 + m;
    j3 = j2 + m;
    x0r = a[0] - a[j2 + 1];
    x0i = a[1] + a[j2];
    x1r = a[0] + a[j2 + 1];
    x1i = a[1] - a[j2];
    x2r = a[j1] - a[j3 + 1];
    x2i = a[j1 + 1] + a[j3];
    x3r = a[j1] + a[j3 + 1];
    x3i = a[j1 + 1] - a[j3];
    y0r = wn4r * (x2r - x2i);
    y0i = wn4r * (x2i + x2r);
    a[0] = x0r + y0r;
    a[1] = x0i + y0i;
    a[j1] = x0r - y0r;
    a[j1 + 1] = x0i - y0i;
    y0r = wn4r * (x3r - x3i);
    y0i = wn4r * (x3i + x3r);
    a[j2] = x1r - y0i;
    a[j2 + 1] = x1i + y0r;
    a[j3] = x1r + y0i;
    a[j3 + 1] = x1i - y0r;
    k = 0;
    kr = 2 * m;
    for (j = 2; j < mh; j += 2) {
        k += 4;
        wk1r = w[k];
        wk1i = w[k + 1];
        wk3r = w[k + 2];
        wk3i = w[k + 3];
        kr -= 4;
        wd1i = w[kr];
        wd1r = w[kr + 1];
        wd3i = w[kr + 2];
        wd3r = w[kr + 3];
        j1 = j + m;
        j2 = j1 + m;
        j3 = j2 + m;
        x0r = a[j] - a[j2 + 1];
        x0i = a[j + 1] + a[j2];
        x1r = a[j] + a[j2 + 1];
        x1i = a[j + 1] - a[j2];
        x2r = a[j1] - a[j3 + 1];
        x2i = a[j1 + 1] + a[j3];
        x3r = a[j1] + a[j3 + 1];
        x3i = a[j1 + 1] - a[j3];
        y0r = wk1r * x0r - wk1i * x0i;
        y0i = wk1r * x0i + wk1i * x0r;
        y2r = wd1r * x2r - wd1i * x2i;
        y2i = wd1r * x2i + wd1i * x2r;
        a[j] = y0r + y2r;
        a[j + 1] = y0i + y2i;
        a[j1] = y0r - y2r;
        a[j1 + 1] = y0i - y2i;
        y0r = wk3r * x1r + wk3i * x1i;
        y0i = wk3r * x1i - wk3i * x1r;
        y2r = wd3r * x3r + wd3i * x3i;
        y2i = wd3r * x3i - wd3i * x3r;
        a[j2] = y0r + y2r;
        a[j2 + 1] = y0i + y2i;
        a[j3] = y0r - y2r;
        a[j3 + 1] = y0i - y2i;
        j0 = m - j;
        j1 = j0 + m;
        j2 = j1 + m;
        j3 = j2 + m;
        x0r = a[j0] - a[j2 + 1];
        x0i = a[j0 + 1] + a[j2];
        x1r = a[j0] + a[j2 + 1];
        x1i = a[j0 + 1] - a[j2];
        x2r = a[j1] - a[j3 + 1];
        x2i = a[j1 + 1] + a[j3];
        x3r = a[j1] + a[j3 + 1];
        x3i = a[j1 + 1] - a[j3];
        y0r = wd1i * x0r - wd1r * x0i;
        y0i = wd1i * x0i + wd1r * x0r;
        y2r = wk1i * x2r - wk1r * x2i;
        y2i = wk1i * x2i + wk1r * x2r;
        a[j0] = y0r + y2r;
        a[j0 + 1] = y0i + y2i;
        a[j1] = y0r - y2r;
        a[j1 + 1] = y0i - y2i;
        y0r = wd3i * x1r + wd3r * x1i;
        y0i = wd3i * x1i - wd3r * x1r;
        y2r = wk3i * x3r + wk3r * x3i;
        y2i = wk3i * x3i - wk3r * x3r;
        a[j2] = y0r + y2r;
        a[j2 + 1] = y0i + y2i;
        a[j3] = y0r - y2r;
        a[j3 + 1] = y0i - y2i;
    }
    wk1r = w[m];
    wk1i = w[m + 1];
    j0 = mh;
    j1 = j0 + m;
    j2 = j1 + m;
    j3 = j2 + m;
    x0r = a[j0] - a[j2 + 1];
    x0i = a[j0 + 1] + a[j2];
    x1r = a[j0] + a[j2 + 1];
    x1i = a[j0 + 1] - a[j2];
    x2r = a[j1] - a[j3 + 1];
    x2i = a[j1 + 1] + a[j3];
    x3r = a[j1] + a[j3 + 1];
    x3i = a[j1 + 1] - a[j3];
    y0r = wk1r * x0r - wk1i * x0i;
    y0i = wk1r * x0i + wk1i * x0r;
    y2r = wk1i * x2r - wk1r * x2i;
    y2i = wk1i * x2i + wk1r * x2r;
    a[j0] = y0r + y2r;
    a[j0 + 1] = y0i + y2i;
    a[j1] = y0r - y2r;
    a[j1 + 1] = y0i - y2i;
    y0r = wk1i * x1r - wk1r * x1i;
    y0i = wk1i * x1i + wk1r * x1r;
    y2r = wk1r * x3r - wk1i * x3i;
    y2i = wk1r * x3i + wk1i * x3r;
    a[j2] = y0r - y2r;
    a[j2 + 1] = y0i - y2i;
    a[j3] = y0r + y2r;
    a[j3 + 1] = y0i + y2i;
}


void cftfx41(int n, double *a, int nw, double *w)
{
    void cftf161(double *a, double *w);
    void cftf162(double *a, double *w);
    void cftf081(double *a, double *w);
    void cftf082(double *a, double *w);
    
    if (n == 128) {
        cftf161(a, &w[nw - 8]);
        cftf162(&a[32], &w[nw - 32]);
        cftf161(&a[64], &w[nw - 8]);
        cftf161(&a[96], &w[nw - 8]);
    } else {
        cftf081(a, &w[nw - 8]);
        cftf082(&a[16], &w[nw - 8]);
        cftf081(&a[32], &w[nw - 8]);
        cftf081(&a[48], &w[nw - 8]);
    }
}


void cftf161(double *a, double *w)
{
    double wn4r, wk1r, wk1i, 
        x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
        y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
        y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, 
        y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, 
        y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
    
    wn4r = w[1];
    wk1r = w[2];
    wk1i = w[3];
    x0r = a[0] + a[16];
    x0i = a[1] + a[17];
    x1r = a[0] - a[16];
    x1i = a[1] - a[17];
    x2r = a[8] + a[24];
    x2i = a[9] + a[25];
    x3r = a[8] - a[24];
    x3i = a[9] - a[25];
    y0r = x0r + x2r;
    y0i = x0i + x2i;
    y4r = x0r - x2r;
    y4i = x0i - x2i;
    y8r = x1r - x3i;
    y8i = x1i + x3r;
    y12r = x1r + x3i;
    y12i = x1i - x3r;
    x0r = a[2] + a[18];
    x0i = a[3] + a[19];
    x1r = a[2] - a[18];
    x1i = a[3] - a[19];
    x2r = a[10] + a[26];
    x2i = a[11] + a[27];
    x3r = a[10] - a[26];
    x3i = a[11] - a[27];
    y1r = x0r + x2r;
    y1i = x0i + x2i;
    y5r = x0r - x2r;
    y5i = x0i - x2i;
    x0r = x1r - x3i;
    x0i = x1i + x3r;
    y9r = wk1r * x0r - wk1i * x0i;
    y9i = wk1r * x0i + wk1i * x0r;
    x0r = x1r + x3i;
    x0i = x1i - x3r;
    y13r = wk1i * x0r - wk1r * x0i;
    y13i = wk1i * x0i + wk1r * x0r;
    x0r = a[4] + a[20];
    x0i = a[5] + a[21];
    x1r = a[4] - a[20];
    x1i = a[5] - a[21];
    x2r = a[12] + a[28];
    x2i = a[13] + a[29];
    x3r = a[12] - a[28];
    x3i = a[13] - a[29];
    y2r = x0r + x2r;
    y2i = x0i + x2i;
    y6r = x0r - x2r;
    y6i = x0i - x2i;
    x0r = x1r - x3i;
    x0i = x1i + x3r;
    y10r = wn4r * (x0r - x0i);
    y10i = wn4r * (x0i + x0r);
    x0r = x1r + x3i;
    x0i = x1i - x3r;
    y14r = wn4r * (x0r + x0i);
    y14i = wn4r * (x0i - x0r);
    x0r = a[6] + a[22];
    x0i = a[7] + a[23];
    x1r = a[6] - a[22];
    x1i = a[7] - a[23];
    x2r = a[14] + a[30];
    x2i = a[15] + a[31];
    x3r = a[14] - a[30];
    x3i = a[15] - a[31];
    y3r = x0r + x2r;
    y3i = x0i + x2i;
    y7r = x0r - x2r;
    y7i = x0i - x2i;
    x0r = x1r - x3i;
    x0i = x1i + x3r;
    y11r = wk1i * x0r - wk1r * x0i;
    y11i = wk1i * x0i + wk1r * x0r;
    x0r = x1r + x3i;
    x0i = x1i - x3r;
    y15r = wk1r * x0r - wk1i * x0i;
    y15i = wk1r * x0i + wk1i * x0r;
    x0r = y12r - y14r;
    x0i = y12i - y14i;
    x1r = y12r + y14r;
    x1i = y12i + y14i;
    x2r = y13r - y15r;
    x2i = y13i - y15i;
    x3r = y13r + y15r;
    x3i = y13i + y15i;
    a[24] = x0r + x2r;
    a[25] = x0i + x2i;
    a[26] = x0r - x2r;
    a[27] = x0i - x2i;
    a[28] = x1r - x3i;
    a[29] = x1i + x3r;
    a[30] = x1r + x3i;
    a[31] = x1i - x3r;
    x0r = y8r + y10r;
    x0i = y8i + y10i;
    x1r = y8r - y10r;
    x1i = y8i - y10i;
    x2r = y9r + y11r;
    x2i = y9i + y11i;
    x3r = y9r - y11r;
    x3i = y9i - y11i;
    a[16] = x0r + x2r;
    a[17] = x0i + x2i;
    a[18] = x0r - x2r;
    a[19] = x0i - x2i;
    a[20] = x1r - x3i;
    a[21] = x1i + x3r;
    a[22] = x1r + x3i;
    a[23] = x1i - x3r;
    x0r = y5r - y7i;
    x0i = y5i + y7r;
    x2r = wn4r * (x0r - x0i);
    x2i = wn4r * (x0i + x0r);
    x0r = y5r + y7i;
    x0i = y5i - y7r;
    x3r = wn4r * (x0r - x0i);
    x3i = wn4r * (x0i + x0r);
    x0r = y4r - y6i;
    x0i = y4i + y6r;
    x1r = y4r + y6i;
    x1i = y4i - y6r;
    a[8] = x0r + x2r;
    a[9] = x0i + x2i;
    a[10] = x0r - x2r;
    a[11] = x0i - x2i;
    a[12] = x1r - x3i;
    a[13] = x1i + x3r;
    a[14] = x1r + x3i;
    a[15] = x1i - x3r;
    x0r = y0r + y2r;
    x0i = y0i + y2i;
    x1r = y0r - y2r;
    x1i = y0i - y2i;
    x2r = y1r + y3r;
    x2i = y1i + y3i;
    x3r = y1r - y3r;
    x3i = y1i - y3i;
    a[0] = x0r + x2r;
    a[1] = x0i + x2i;
    a[2] = x0r - x2r;
    a[3] = x0i - x2i;
    a[4] = x1r - x3i;
    a[5] = x1i + x3r;
    a[6] = x1r + x3i;
    a[7] = x1i - x3r;
}


void cftf162(double *a, double *w)
{
    double wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, 
        x0r, x0i, x1r, x1i, x2r, x2i, 
        y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
        y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, 
        y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, 
        y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
    
    wn4r = w[1];
    wk1r = w[4];
    wk1i = w[5];
    wk3r = w[6];
    wk3i = -w[7];
    wk2r = w[8];
    wk2i = w[9];
    x1r = a[0] - a[17];
    x1i = a[1] + a[16];
    x0r = a[8] - a[25];
    x0i = a[9] + a[24];
    x2r = wn4r * (x0r - x0i);
    x2i = wn4r * (x0i + x0r);
    y0r = x1r + x2r;
    y0i = x1i + x2i;
    y4r = x1r - x2r;
    y4i = x1i - x2i;
    x1r = a[0] + a[17];
    x1i = a[1] - a[16];
    x0r = a[8] + a[25];
    x0i = a[9] - a[24];
    x2r = wn4r * (x0r - x0i);
    x2i = wn4r * (x0i + x0r);
    y8r = x1r - x2i;
    y8i = x1i + x2r;
    y12r = x1r + x2i;
    y12i = x1i - x2r;
    x0r = a[2] - a[19];
    x0i = a[3] + a[18];
    x1r = wk1r * x0r - wk1i * x0i;
    x1i = wk1r * x0i + wk1i * x0r;
    x0r = a[10] - a[27];
    x0i = a[11] + a[26];
    x2r = wk3i * x0r - wk3r * x0i;
    x2i = wk3i * x0i + wk3r * x0r;
    y1r = x1r + x2r;
    y1i = x1i + x2i;
    y5r = x1r - x2r;
    y5i = x1i - x2i;
    x0r = a[2] + a[19];
    x0i = a[3] - a[18];
    x1r = wk3r * x0r - wk3i * x0i;
    x1i = wk3r * x0i + wk3i * x0r;
    x0r = a[10] + a[27];
    x0i = a[11] - a[26];
    x2r = wk1r * x0r + wk1i * x0i;
    x2i = wk1r * x0i - wk1i * x0r;
    y9r = x1r - x2r;
    y9i = x1i - x2i;
    y13r = x1r + x2r;
    y13i = x1i + x2i;
    x0r = a[4] - a[21];
    x0i = a[5] + a[20];
    x1r = wk2r * x0r - wk2i * x0i;
    x1i = wk2r * x0i + wk2i * x0r;
    x0r = a[12] - a[29];
    x0i = a[13] + a[28];
    x2r = wk2i * x0r - wk2r * x0i;
    x2i = wk2i * x0i + wk2r * x0r;
    y2r = x1r + x2r;
    y2i = x1i + x2i;
    y6r = x1r - x2r;
    y6i = x1i - x2i;
    x0r = a[4] + a[21];
    x0i = a[5] - a[20];
    x1r = wk2i * x0r - wk2r * x0i;
    x1i = wk2i * x0i + wk2r * x0r;
    x0r = a[12] + a[29];
    x0i = a[13] - a[28];
    x2r = wk2r * x0r - wk2i * x0i;
    x2i = wk2r * x0i + wk2i * x0r;
    y10r = x1r - x2r;
    y10i = x1i - x2i;
    y14r = x1r + x2r;
    y14i = x1i + x2i;
    x0r = a[6] - a[23];
    x0i = a[7] + a[22];
    x1r = wk3r * x0r - wk3i * x0i;
    x1i = wk3r * x0i + wk3i * x0r;
    x0r = a[14] - a[31];
    x0i = a[15] + a[30];
    x2r = wk1i * x0r - wk1r * x0i;
    x2i = wk1i * x0i + wk1r * x0r;
    y3r = x1r + x2r;
    y3i = x1i + x2i;
    y7r = x1r - x2r;
    y7i = x1i - x2i;
    x0r = a[6] + a[23];
    x0i = a[7] - a[22];
    x1r = wk1i * x0r + wk1r * x0i;
    x1i = wk1i * x0i - wk1r * x0r;
    x0r = a[14] + a[31];
    x0i = a[15] - a[30];
    x2r = wk3i * x0r - wk3r * x0i;
    x2i = wk3i * x0i + wk3r * x0r;
    y11r = x1r + x2r;
    y11i = x1i + x2i;
    y15r = x1r - x2r;
    y15i = x1i - x2i;
    x1r = y0r + y2r;
    x1i = y0i + y2i;
    x2r = y1r + y3r;
    x2i = y1i + y3i;
    a[0] = x1r + x2r;
    a[1] = x1i + x2i;
    a[2] = x1r - x2r;
    a[3] = x1i - x2i;
    x1r = y0r - y2r;
    x1i = y0i - y2i;
    x2r = y1r - y3r;
    x2i = y1i - y3i;
    a[4] = x1r - x2i;
    a[5] = x1i + x2r;
    a[6] = x1r + x2i;
    a[7] = x1i - x2r;
    x1r = y4r - y6i;
    x1i = y4i + y6r;
    x0r = y5r - y7i;
    x0i = y5i + y7r;
    x2r = wn4r * (x0r - x0i);
    x2i = wn4r * (x0i + x0r);
    a[8] = x1r + x2r;
    a[9] = x1i + x2i;
    a[10] = x1r - x2r;
    a[11] = x1i - x2i;
    x1r = y4r + y6i;
    x1i = y4i - y6r;
    x0r = y5r + y7i;
    x0i = y5i - y7r;
    x2r = wn4r * (x0r - x0i);
    x2i = wn4r * (x0i + x0r);
    a[12] = x1r - x2i;
    a[13] = x1i + x2r;
    a[14] = x1r + x2i;
    a[15] = x1i - x2r;
    x1r = y8r + y10r;
    x1i = y8i + y10i;
    x2r = y9r - y11r;
    x2i = y9i - y11i;
    a[16] = x1r + x2r;
    a[17] = x1i + x2i;
    a[18] = x1r - x2r;
    a[19] = x1i - x2i;
    x1r = y8r - y10r;
    x1i = y8i - y10i;
    x2r = y9r + y11r;
    x2i = y9i + y11i;
    a[20] = x1r - x2i;
    a[21] = x1i + x2r;
    a[22] = x1r + x2i;
    a[23] = x1i - x2r;
    x1r = y12r - y14i;
    x1i = y12i + y14r;
    x0r = y13r + y15i;
    x0i = y13i - y15r;
    x2r = wn4r * (x0r - x0i);
    x2i = wn4r * (x0i + x0r);
    a[24] = x1r + x2r;
    a[25] = x1i + x2i;
    a[26] = x1r - x2r;
    a[27] = x1i - x2i;
    x1r = y12r + y14i;
    x1i = y12i - y14r;
    x0r = y13r - y15i;
    x0i = y13i + y15r;
    x2r = wn4r * (x0r - x0i);
    x2i = wn4r * (x0i + x0r);
    a[28] = x1r - x2i;
    a[29] = x1i + x2r;
    a[30] = x1r + x2i;
    a[31] = x1i - x2r;
}


void cftf081(double *a, double *w)
{
    double wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
        y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
        y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
    
    wn4r = w[1];
    x0r = a[0] + a[8];
    x0i = a[1] + a[9];
    x1r = a[0] - a[8];
    x1i = a[1] - a[9];
    x2r = a[4] + a[12];
    x2i = a[5] + a[13];
    x3r = a[4] - a[12];
    x3i = a[5] - a[13];
    y0r = x0r + x2r;
    y0i = x0i + x2i;
    y2r = x0r - x2r;
    y2i = x0i - x2i;
    y1r = x1r - x3i;
    y1i = x1i + x3r;
    y3r = x1r + x3i;
    y3i = x1i - x3r;
    x0r = a[2] + a[10];
    x0i = a[3] + a[11];
    x1r = a[2] - a[10];
    x1i = a[3] - a[11];
    x2r = a[6] + a[14];
    x2i = a[7] + a[15];
    x3r = a[6] - a[14];
    x3i = a[7] - a[15];
    y4r = x0r + x2r;
    y4i = x0i + x2i;
    y6r = x0r - x2r;
    y6i = x0i - x2i;
    x0r = x1r - x3i;
    x0i = x1i + x3r;
    x2r = x1r + x3i;
    x2i = x1i - x3r;
    y5r = wn4r * (x0r - x0i);
    y5i = wn4r * (x0r + x0i);
    y7r = wn4r * (x2r - x2i);
    y7i = wn4r * (x2r + x2i);
    a[8] = y1r + y5r;
    a[9] = y1i + y5i;
    a[10] = y1r - y5r;
    a[11] = y1i - y5i;
    a[12] = y3r - y7i;
    a[13] = y3i + y7r;
    a[14] = y3r + y7i;
    a[15] = y3i - y7r;
    a[0] = y0r + y4r;
    a[1] = y0i + y4i;
    a[2] = y0r - y4r;
    a[3] = y0i - y4i;
    a[4] = y2r - y6i;
    a[5] = y2i + y6r;
    a[6] = y2r + y6i;
    a[7] = y2i - y6r;
}


void cftf082(double *a, double *w)
{
    double wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i, 
        y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
        y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
    
    wn4r = w[1];
    wk1r = w[2];
    wk1i = w[3];
    y0r = a[0] - a[9];
    y0i = a[1] + a[8];
    y1r = a[0] + a[9];
    y1i = a[1] - a[8];
    x0r = a[4] - a[13];
    x0i = a[5] + a[12];
    y2r = wn4r * (x0r - x0i);
    y2i = wn4r * (x0i + x0r);
    x0r = a[4] + a[13];
    x0i = a[5] - a[12];
    y3r = wn4r * (x0r - x0i);
    y3i = wn4r * (x0i + x0r);
    x0r = a[2] - a[11];
    x0i = a[3] + a[10];
    y4r = wk1r * x0r - wk1i * x0i;
    y4i = wk1r * x0i + wk1i * x0r;
    x0r = a[2] + a[11];
    x0i = a[3] - a[10];
    y5r = wk1i * x0r - wk1r * x0i;
    y5i = wk1i * x0i + wk1r * x0r;
    x0r = a[6] - a[15];
    x0i = a[7] + a[14];
    y6r = wk1i * x0r - wk1r * x0i;
    y6i = wk1i * x0i + wk1r * x0r;
    x0r = a[6] + a[15];
    x0i = a[7] - a[14];
    y7r = wk1r * x0r - wk1i * x0i;
    y7i = wk1r * x0i + wk1i * x0r;
    x0r = y0r + y2r;
    x0i = y0i + y2i;
    x1r = y4r + y6r;
    x1i = y4i + y6i;
    a[0] = x0r + x1r;
    a[1] = x0i + x1i;
    a[2] = x0r - x1r;
    a[3] = x0i - x1i;
    x0r = y0r - y2r;
    x0i = y0i - y2i;
    x1r = y4r - y6r;
    x1i = y4i - y6i;
    a[4] = x0r - x1i;
    a[5] = x0i + x1r;
    a[6] = x0r + x1i;
    a[7] = x0i - x1r;
    x0r = y1r - y3i;
    x0i = y1i + y3r;
    x1r = y5r - y7r;
    x1i = y5i - y7i;
    a[8] = x0r + x1r;
    a[9] = x0i + x1i;
    a[10] = x0r - x1r;
    a[11] = x0i - x1i;
    x0r = y1r + y3i;
    x0i = y1i - y3r;
    x1r = y5r + y7r;
    x1i = y5i + y7i;
    a[12] = x0r - x1i;
    a[13] = x0i + x1r;
    a[14] = x0r + x1i;
    a[15] = x0i - x1r;
}


void cftf040(double *a)
{
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
    
    x0r = a[0] + a[4];
    x0i = a[1] + a[5];
    x1r = a[0] - a[4];
    x1i = a[1] - a[5];
    x2r = a[2] + a[6];
    x2i = a[3] + a[7];
    x3r = a[2] - a[6];
    x3i = a[3] - a[7];
    a[0] = x0r + x2r;
    a[1] = x0i + x2i;
    a[2] = x1r - x3i;
    a[3] = x1i + x3r;
    a[4] = x0r - x2r;
    a[5] = x0i - x2i;
    a[6] = x1r + x3i;
    a[7] = x1i - x3r;
}


void cftb040(double *a)
{
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
    
    x0r = a[0] + a[4];
    x0i = a[1] + a[5];
    x1r = a[0] - a[4];
    x1i = a[1] - a[5];
    x2r = a[2] + a[6];
    x2i = a[3] + a[7];
    x3r = a[2] - a[6];
    x3i = a[3] - a[7];
    a[0] = x0r + x2r;
    a[1] = x0i + x2i;
    a[2] = x1r + x3i;
    a[3] = x1i - x3r;
    a[4] = x0r - x2r;
    a[5] = x0i - x2i;
    a[6] = x1r - x3i;
    a[7] = x1i + x3r;
}


void cftx020(double *a)
{
    double x0r, x0i;
    
    x0r = a[0] - a[2];
    x0i = a[1] - a[3];
    a[0] += a[2];
    a[1] += a[3];
    a[2] = x0r;
    a[3] = x0i;
}


void rftfsub(int n, double *a, int nc, double *c)
{
    int j, k, kk, ks, m;
    double wkr, wki, xr, xi, yr, yi;
    
    m = n >> 1;
    ks = 2 * nc / m;
    kk = 0;
    for (j = 2; j < m; j += 2) {
        k = n - j;
        kk += ks;
        wkr = 0.5 - c[nc - kk];
        wki = c[kk];
        xr = a[j] - a[k];
        xi = a[j + 1] + a[k + 1];
        yr = wkr * xr - wki * xi;
        yi = wkr * xi + wki * xr;
        a[j] -= yr;
        a[j + 1] -= yi;
        a[k] += yr;
        a[k + 1] -= yi;
    }
}


void rftbsub(int n, double *a, int nc, double *c)
{
    int j, k, kk, ks, m;
    double wkr, wki, xr, xi, yr, yi;
    
    m = n >> 1;
    ks = 2 * nc / m;
    kk = 0;
    for (j = 2; j < m; j += 2) {
        k = n - j;
        kk += ks;
        wkr = 0.5 - c[nc - kk];
        wki = c[kk];
        xr = a[j] - a[k];
        xi = a[j + 1] + a[k + 1];
        yr = wkr * xr + wki * xi;
        yi = wkr * xi - wki * xr;
        a[j] -= yr;
        a[j + 1] -= yi;
        a[k] += yr;
        a[k + 1] -= yi;
    }
}


void dctsub(int n, double *a, int nc, double *c)
{
    int j, k, kk, ks, m;
    double wkr, wki, xr;
    
    m = n >> 1;
    ks = nc / n;
    kk = 0;
    for (j = 1; j < m; j++) {
        k = n - j;
        kk += ks;
        wkr = c[kk] - c[nc - kk];
        wki = c[kk] + c[nc - kk];
        xr = wki * a[j] - wkr * a[k];
        a[j] = wkr * a[j] + wki * a[k];
        a[k] = xr;
    }
    a[m] *= c[0];
}


void dstsub(int n, double *a, int nc, double *c)
{
    int j, k, kk, ks, m;
    double wkr, wki, xr;
    
    m = n >> 1;
    ks = nc / n;
    kk = 0;
    for (j = 1; j < m; j++) {
        k = n - j;
        kk += ks;
        wkr = c[kk] - c[nc - kk];
        wki = c[kk] + c[nc - kk];
        xr = wki * a[k] - wkr * a[j];
        a[k] = wkr * a[k] + wki * a[j];
        a[j] = xr;
    }
    a[m] *= c[0];
}
/*
**	fft.c -- Fast Fourier Tramsform (FFT)
**	時間軸間引きFFT (Cooley--Tukeyのアルゴリズム)
**
**	『Ｃ言語による最新アルゴリズム事典』（奥村晴彦・著）に掲載の
**	プログラムを改変して利用させて頂きました.
**
**	Revised by MIYASAKA Masaru <alkaid@coral.ocn.ne.jp> (Sep 15, 2003)
**
**	[TAB4]
*/

//#include <stdlib.h>
//#include <math.h>
//#include <stdio.h>

//#include "fft_czt.h"

//#undef PI
//#define PI	3.1415926535897932384626433832795L	/* 円周率(π) */


/*
**	三角関数表を作る.
*/
static void make_sintbl(int n, REAL sintbl[])
{
	/* ↓精度が重要でないなら double でもいいかも */
	double PI;
	PI =(double)atan(1.0)*4.0;
	long double r, d = (2 * PI) / n;
	int i, n2 = n / 2, n4 = n / 4, n8 = n / 8;

	for (i = 0; i <= n8; i++) {
		r = d * i;
		sintbl[i]      = sin(r);			/* sin(2*PI*i/n) */
		sintbl[n4 - i] = cos(r);
	}
	for (i = 0; i < n4; i++) {
		sintbl[n2 - i] =  sintbl[i];
		sintbl[i + n2] = -sintbl[i];
	}
}


/*
**	ビット反転表を作る.
*/
static void make_bitrev(int n, int bitrev[])
{
	int i = 0, j = 0, k;
	int n2 = n / 2;

	for (;;) {
		bitrev[i] = j;
		if (++i >= n) break;
		k = n2;
		while (k <= j) { j -= k; k /= 2; }
		j += k;
	}
}


/*
**	FFT計算用構造体に対し、標本数 n 用の数表データを作成する。
**
**	fftp	= FFT計算用構造体へのポインタ
**	n		= 標本点の数 ( 4以上の 2の整数乗に限る)
**	return	= 0:正常終了 1:nが無効な数 2:メモリ不足
*/
int fft_init(fft_struct *fftp, int n)
{
	int i;

	for (i = 4; i < n; i *= 2) ;
	if (i != n) return 1;

	fftp->samples = n;
	fftp->sintbl  = (REAL*)malloc((n / 4 * 3) * sizeof(REAL));
	fftp->bitrev  = (int*)malloc(n * sizeof(int));
	if (fftp->sintbl == NULL || fftp->bitrev == NULL) {
		fft_end(fftp);
		return 2;
	}
//	make_sintbl(n, fftp->sintbl);
//	make_bitrev(n, fftp->bitrev);

	return 0;
}


/*
**	FFT計算用構造体の数表データを消去してそのメモリ領域を開放する。
**
**	fftp	= FFT計算用構造体へのポインタ
*/
void fft_end(fft_struct *fftp)
{
	if (fftp->sintbl != NULL) { free(fftp->sintbl); fftp->sintbl = NULL; }
	if (fftp->bitrev != NULL) { free(fftp->bitrev); fftp->bitrev = NULL; }

//	fftp->samples = 0;
}


/*
**	高速Fourier変換 (FFT) (Cooley--Tukeyのアルゴリズム).
**	re[] が実部, im[] が虚部. 結果は re[],im[] に上書きされる.
**	inv!=0 (=TRUE) なら逆変換を行う. fftp には, 計算用データが
**	入っている構造体を指定する.
*/
void fft_tukawan(fft_struct *fftp, int inv, REAL re[], REAL im[])
{
	REAL si, co, dr, di, t;
	int i, j, k, ik, h, d, k2, n, n4;

	n  = fftp->samples;
	n4 = n / 4;

	for (i = 0; i < n; i++) {		/* ビット反転 */
		j = fftp->bitrev[i];
		if (i < j) {
			t = re[i];  re[i] = re[j];  re[j] = t;
			t = im[i];  im[i] = im[j];  im[j] = t;
		}
	}
	for (k = 1; k < n; k = k2) {	/* 変換 */
		h = 0;  k2 = k + k;  d = n / k2;
		for (j = 0; j < k; j++) {
			co = fftp->sintbl[h + n4];
			si = fftp->sintbl[h];
			if (inv) si = -si;
			for (i = j; i < n; i += k2) {
				ik = i + k;
				dr = si * im[ik] + co * re[ik];
				di = co * im[ik] - si * re[ik];
				re[ik] = re[i] - dr;  re[i] += dr;
				im[ik] = im[i] - di;  im[i] += di;
			}
			h += d;
		}
	}
	double t1;
	if (!inv) {						/* 逆変換でないならnで割る */
		t = 1.0 / n;			/* 逆数をかける(除算は遅いので) */
        t1 = 0.5 / n;
//		t = (4.3775) / (double)n ;			/* 逆数をかける(除算は遅いので) */
		for (i = 0; i < n; i++) {
			re[i] *= t;
			im[i] *= t;
		}
//		for (i = 0; i < n; i++) {
//			if(i==0){
//			re[i] *= t1;
//			im[i] *= t1;
//			} else { 
//			re[i] *= t;
//			im[i] *= t;
//			}
//		}
	}
}
/*
void fft(double *real,double *imag, int power, int flg)
{
double c,s,t1,t2;
int arg,inc,i,j,k,lmax,lix,j1,j2,pnum,pnum1,nv2;
pnum = (int)((pow(2.0,(float)(power)))+0.5);
double PI = atan(1.0) * 4.0;

lmax = pnum;
inc = 1;
for(i=0;i<power;i++){
lix=lmax;
lmax /= 2;
arg = 0;
for(j=0;j<lmax;j++){
c=cos(2.0*PI*((double)(arg))/((double)(pnum)));
s=flg*sin(2.0*PI*((double)(arg))/((double)(pnum)));
arg += inc;
for(k=lix;k<=pnum;k+=lix){
j1=k-lix+j;
j2=j1+lmax;
t1=(*(real+j1))-(*(real+j2));
t2=(*(imag+j1))-(*(imag+j2));
*(real+j1) += *(real+j2);
*(imag+j1) += *(imag+j2);
*(real+j2) = c*t1+s*t2;
*(imag+j2) = c*t2-s*t1;
}
}
inc *= 2;
}
j=0;
nv2=pnum/2;
pnum1=pnum-1;
for(i=0;i<pnum1;i++){
if(i<j){
t1=*(real+j);
t2=*(imag+j);
*(real+j)=*(real+i);
*(imag+j)=*(imag+i);
*(real+i)=t1;
*(imag+i)=t2;
}
k=nv2;
while(k<j+1){
j -= k;
k /= 2;
}
j+=k;
}
	double sc=1.0/(double)(pnum);
	if(flg==-1){
 		for(i=0;i<pnum;i++){
			*(real+i) *= sc;
			*(imag+i) *= sc;
		}
	}
}
*/

/*
**	czt.c -- Chirp z-Tramsform (CZT)
**	任意データ数での計算が可能な高速 DFT アルゴリズム
**
**	『トランジスタ技術』1993年2月号 P.363〜368 の解説を参考にしています.
**
**	Public domain by MIYASAKA Masaru <alkaid@coral.ocn.ne.jp> (Sep 15, 2003)
**
**	[TAB4]
*/


/*
**	fft_czt.h -- header file for fft.c and czt.c
**
**	Public domain by MIYASAKA Masaru <alkaid@coral.ocn.ne.jp> (Sep 15, 2003)
**
**	[TAB4]
*/

//#ifndef FFT_CZT_H
//#define FFT_CZT_H

	/* データ実数型 REAL の定義 (float or double) */
	/* この型を long double にする場合は、czt.c / fft.c で使われている
	 * 三角関数 sin() / cos() を long double 版 (C99 対応コンパイラでは
	 * sinl / cosl) に変更する必要あり */



//#endif	/* FFT_CZT_H */


//#undef PI
//#define PI	3.1415926535897932384626433832795L	/* 円周率(π) */


/*
**	重みデータ・インパルス応答データを作る
*/
static void make_cztdata(int n, int no, int nx, REAL wr[], REAL wi[],
                         REAL vr[], REAL vi[])
{
	/* ↓精度が重要でないなら double でもいいかも */
	double PI;
	PI = (double)atan(1.0)*4.0;
	long double r, d = PI / n;
	int i, j, n2 = n * 2;

	for (i = j = 0; i < n; i++, j = (j + (2 * i - 1)) % n2) {
		r = d * j;
		wr[i] =  cos(r);					/*  cos(i*i*PI/n) */
		wi[i] = -sin(r);					/* -sin(i*i*PI/n) */
	}

	vr[0] = nx;					/* nx: FFT順変換のスケーリングに対応する補正 */
	vi[0] = 0;
	for (i = 1; i < n; i++) {
		vr[i] = vr[nx - i] =  nx * wr[i];	/* nx * cos(i*i*PI/n) */
		vi[i] = vi[nx - i] = -nx * wi[i];	/* nx * sin(i*i*PI/n) */
	}
	for (i = no, j = nx - n; i <= j; i++, j--) {
		vr[j] = vr[i] = 0;					/* 残りの拡張部には 0 を */
		vi[j] = vi[i] = 0;
	}
}


/*
**	CZT計算用構造体に対し、標本数 n, 出力データ数 no 用の数表データを
**	作成する。
**
**	cztp	= CZT計算用構造体へのポインタ
**	n		= 標本点の数
**	no		= 出力するデータの数
**	return	= 0:正常終了 1:nが無効な数 2:メモリ不足
*/
int czt_init(czt_struct *cztp, int n, int no,int *ip,double *ww)
{
	int i, nx,powe;

	if (n <= 1) return 1;
	if (no <= 1 || n < no) no = n;

	nx = n + no;		/* (n + no)を2の整数乗まで拡張する(nx) */
	for (i = 1; i < nx; i *= 2) ;
	nx = i;//64
	//debug start
	//printf("nx=%d\n",nx);
	//debug end

	if (fft_init(&cztp->fft, nx) != 0) return 2;

	cztp->samples     = n;
	cztp->samples_out = no;
	cztp->samples_ex  = nx;
	cztp->wr = (REAL*)malloc(2 * n  * sizeof(REAL));
	cztp->vr = (REAL*)malloc(2 * nx * sizeof(REAL));
	cztp->tr = (REAL*)malloc(2 * nx * sizeof(REAL));
	if (cztp->wr == NULL || cztp->vr == NULL || cztp->tr == NULL) {
		czt_end(cztp);
		return 2;
	}
	cztp->wi = cztp->wr + n;
	cztp->vi = cztp->vr + nx;
	cztp->ti = cztp->tr + nx;

	make_cztdata(n, no, nx, cztp->wr, cztp->wi, cztp->vr, cztp->vi);
	double *temp/*,*ww*/;
	temp = (double *)malloc(sizeof(double)*nx*2);
//	ww = (double *)malloc(sizeof(double)*nx);
	for(i=0;i<nx;i++){
		temp[2*i] = cztp->vr[i];
		temp[2*i+1] = cztp->vi[i];
	}
//	int *ip;
//	ip = (int *)malloc(sizeof(int)*nx);
//	ip[0] = 0;
//debug s
//debug e
	cdft(2*nx,-1,temp,ip,ww);
	for(i=0;i<nx;i++){
		cztp->vr[i] = temp[2*i]/nx ;
		cztp->vi[i] = temp[2*i+1]/nx ;
	}
	free(temp);

//	fft(&cztp->fft, 0, cztp->vr, cztp->vi);
//	powe = (int)((log((double)(cztp->fft.samples))/log(2.0))+0.5);
	//debug start
//	printf("cztp->fft.samples=%d,powe=%d\n",cztp->fft.samples,powe);
	//debug end
//debug start	
//	double *zzzr,*zzzi;
//	zzzr=(double*)malloc(sizeof(double)*cztp->fft.samples);
//	zzzi=(double*)malloc(sizeof(double)*cztp->fft.samples);
//	for(i=0;i<cztp->fft.samples;i++){
//		zzzr[i] = cztp->vr[i];
//		zzzi[i] = cztp->vi[i];
//	}	
//	half_fft(cztp->vr, cztp->vi,powe);
//	half_fft(zzzr, zzzi,powe);
//	fft(cztp->vr, cztp->vi,powe,-1);
//	for(i=0;i<1024;i++){
//		printf("N= %d %d half %f full %f\n",cztp->fft.samples,i,zzzi[i],cztp->vi[i]);
//	}
//	while(1);
//debug end	
	return 0;
}


/*
**	CZT計算用構造体の数表データを消去してそのメモリ領域を開放する。
**
**	cztp	= CZT計算用構造体へのポインタ
*/
void czt_end(czt_struct *cztp)
{
	if (cztp->wr != NULL) { free(cztp->wr); cztp->wr = NULL; }
	if (cztp->vr != NULL) { free(cztp->vr); cztp->vr = NULL; }
	if (cztp->tr != NULL) { free(cztp->tr); cztp->tr = NULL; }

	cztp->samples     = 0;
	cztp->samples_out = 0;
	cztp->samples_ex  = 0;

	fft_end(&cztp->fft);
}


/*
**	CZT (Chirp z-Tramsform) アルゴリズムによる高速 DFT.
**	re[] が実部, im[] が虚部. 結果は re[],im[] に上書きされる.
**	inv!=0 (=TRUE) なら逆変換を行う. cztp には, 計算用データが
**	入っている構造体を指定する.
*/
void czt(czt_struct *cztp, int inv, REAL re[], REAL im[],int *ip,double *ww)
{
	REAL xr, xi, yr, yi, t;
	int i, n, no, nx, powe;

	n  = cztp->samples;
	no = cztp->samples_out;
	nx = cztp->samples_ex;

	for (i = 0; i < n; i++) {		/* 入力の重みデータ乗算 */
		yr = cztp->wr[i];
		yi = cztp->wi[i];
		if (inv) yi = -yi;
		cztp->tr[i] = re[i] * yr - im[i] * yi;
		cztp->ti[i] = im[i] * yr + re[i] * yi;
	}
	for (; i < nx; i++) {			/* 残りの拡張部には 0 を */
		cztp->tr[i] = 0;
		cztp->ti[i] = 0;
	}
	double *temp/*,*ww*/;
	temp = (double *)malloc(sizeof(double)*nx*2);
//	ww = (double *)malloc(sizeof(double)*nx);
	for(i=0;i<nx;i++){
		temp[2*i] = cztp->tr[i];
		temp[2*i+1] = cztp->ti[i];
	}
//	fft(&cztp->fft, 0, cztp->tr, cztp->ti);
	cdft(2*nx,-1,temp,ip,ww);
	for(i=0;i<nx;i++){
		cztp->tr[i] = temp[2*i]/nx ;
		cztp->ti[i] = temp[2*i+1]/nx ;
	}

//	fft(&cztp->fft, 0, cztp->tr, cztp->ti);
//	powe = (int)((log((double)(cztp->fft.samples))/log(2.0))+0.5);
	//debug start
//	printf("cztp->fft.samples=%d,powe=%d\n",cztp->fft.samples,powe);
	//debug end
//	half_fft(cztp->tr, cztp->ti,powe);
//	fft(cztp->tr, cztp->ti,powe,-1);

	for (i = 0; i < nx; i++) {		/* コンボリューション演算 */
		xr = cztp->tr[i];
		xi = cztp->ti[i];
		yr = cztp->vr[i];
		yi = cztp->vi[i];
		if (inv) yi = -yi;
		cztp->tr[i] = xr * yr - xi * yi;
		cztp->ti[i] = xi * yr + xr * yi;
	}

//	fft(&cztp->fft, 1, cztp->tr, cztp->ti);
//	powe = (int)((log((double)(cztp->fft.samples))/log(2.0))+0.5);
	//debug start
//	printf("cztp->fft.samples=%d,powe=%d\n",cztp->fft.samples,powe);
	//debug end
//	half_ifft(cztp->tr, cztp->ti,powe);
//	fft(cztp->tr, cztp->ti,powe,1);
	for(i=0;i<nx;i++){
		temp[2*i] = cztp->tr[i];
		temp[2*i+1] = cztp->ti[i];
	}
	cdft(2*nx,1,temp,ip,ww);
	for(i=0;i<nx;i++){
		cztp->tr[i] = temp[2*i]/*/nx*/ ;
		cztp->ti[i] = temp[2*i+1]/*/nx*/ ;
	}

	for (i = 0; i < no; i++) {		/* 出力の重みデータ乗算 */
		yr = cztp->wr[i];
		yi = cztp->wi[i];
		if (inv) yi = -yi;
		re[i] = cztp->tr[i] * yr - cztp->ti[i] * yi;
		im[i] = cztp->ti[i] * yr + cztp->tr[i] * yi;
	}

	if (!inv) {						/* 逆変換でないならnで割る */
		t = 1.0 / n;			/* 逆数をかける(除算は遅いので) */
		for (i = 0; i < no; i++) {
			re[i] *= t;
			im[i] *= t;
		}
	}
}
/*
**	czttest.c -- czt.c のテストルーチン
**
**	『Ｃ言語による最新アルゴリズム事典』（奥村晴彦・著）に掲載の
**	fft.c のテストルーチンを参考にさせて頂きました.
**
**	Public domain by MIYASAKA Masaru <alkaid@coral.ocn.ne.jp> (Sep 15, 2003)
**
**	[TAB4]
*/

//#include <stdio.h>
//#include <stdlib.h>
//#include <math.h>

//#include "fft_czt.h"

//#undef PI
//#define PI	3.1415926535897932384626433832795	/* 円周率(π) */

//#define N 8000	/* テスト用データのサンプル数 */


int charp_FFT(int N,double *R,double *I,double *OR, double *OI,int *ip, double *ww)
{
	static REAL  r2[8192], i2[8192];
	czt_struct cztd;
	int i;
	for (i=0;i<N;i++){
		r2[i]=*(R+i);
		i2[i]=*(I+i);
	}
		if (czt_init(&cztd, N, N,ip,ww) != 0) {printf("failure!!!!\n"); while(1); return EXIT_FAILURE;}
	czt(&cztd, 0, r2, i2,ip,ww);		/* 順変換 */
	for(i=0;i<N;i++){
		*(OR+i) = r2[i];
		*(OI+i) = i2[i];
	}
//	czt(&cztd, 1, r3, i3);		/* 逆変換 */
	czt_end(&cztd);


	return EXIT_SUCCESS;
}
		
int charp_IFFT(int N,double *R,double *I,double *OR,double *OI,int *ip, double *ww)
{
	static REAL r2[8192], i2[8192];
	czt_struct cztd;
	int i;
	for (i=0;i<N;i++){
		r2[i]=*(R+i);
		i2[i]=*(I+i);
	}

	if (czt_init(&cztd, N, N,ip,ww) != 0) return EXIT_FAILURE;
//	czt(&cztd, 0, r2, i2);		/* 順変換 */
	for(i=0;i<N;i++){
		*(OR+i) = r2[i];
		*(OI+i) = i2[i];
	}

	czt(&cztd, 1, r2, i2,ip,ww);		/* 逆変換 */
	for(i=0;i<N;i++){
		*(OR+i) = r2[i];
		*(OI+i) = i2[i];
	}
	czt_end(&cztd);


	return EXIT_SUCCESS;
}



/***************************************************************************************/
/*								  	  	       */
/*      Program Name : bmp_dct_up_resize_ver**.c	  	  	               */
/*								  	  	       */
/*	Author       : zuruneko		  	  	                               */
/*								  	  	       */
/*	Function     : EVEN FUNCTIONIZED DFT RESIZE				       */	
/*										       */
/*	Limitation   : bitmap 24bit true color only				       */
/*		       						  	  	       */
/*      ver01  06.04.05             New Create			   	  	       */
/*      ver02  06.04.06             Counter add                                        */
/*								  	               */
/***************************************************************************************/
				
#define uchar           unsigned char
#define ushort          unsigned short
#define uint            unsigned int
#define ulong           unsigned long
#define rchar           register char
#define rshort          register short
#define rint            register int
#define rlong           register long
#define ruchar          register uchar
#define rushort         register ushort
#define ruint           register uint
#define rulong          register ulong
#define schar           static char
#define sshhort         static short
#define sint            static int
#define slong           static long
#define suchar          static uchar
#define sushort         static ushort
#define suint           static uint
#define sulong          static ulong
#define ORG_IMAGE_X     720
#define ORG_IMAGE_Y     480

#define BUNBOMAX 32
#define BUNSHIMAX 31
#define TAPMAXY 50
#define TAPMAXC 50
                                                                                                                                                         
typedef struct control_data {
  uint          sizex;
  uint          sizey;
  uint          startx;
  uint          starty;
  char         InputFileName[128];
  char         InputFileName2[128];
  char         OutputFileName[128];
  uint          hv;
  uint		sm;
  uint          smy;	
  uint          bunshi;
  uint          bunbo;
  uint          norm;
  uint          mado;
  uint          tapnum_y;
  uint          tapnum_c;
  double        Fs;		
} CNTL;
        char inputfile[50];
        char inputfile2[50];
        char outputfile[50];
        int hsize;
        int vsize;
        int org_image_x;
        int org_image_y;
	void coeffhairetsuculc(uint , uint, uint , uint, uint, uint, uint, uint, double, double, double, double, double *,double *, uint,double);
        double INITPHASECULC(uint , uint , uint , uint , uint );
        uint FILPATCULC(uint , uint , uint , double , double *);
        void COEFCULC( uint , uint , uint , uint , uint , uint , uint, double );
  	void normarize(uint ,uint ,uint ,uint , double *, uint  );
	void DFT1DV(int , double *, double *, double *, double *,int *ip,double *ww);
	void DFT1DH(int , int  , double *, double*, double *, double *,int *ip,double *ww);
	void IDFT1DV(int , double *, double *, double *, double *,int *ip,double *ww);
	void IDFT1DH(int , int  , double *, double*, double *, double *,int *ip,double *ww);
void resize( int hsize, int vsize, int sizeminus, int smy,  unsigned char *vin_buffer , unsigned char *vout);
void    print_help( void )
{
  fprintf( stderr, "Usage: EVEN FUNCTIONIZED DFT UP RESIZE [options]\n" );
  fprintf( stderr, "\toptions: \n" );
  fprintf( stderr, "\t  -help          : help message display\n" );
  fprintf( stderr, "\t  -i             : input bitmap file name 24 bit true color only\n" );
  fprintf( stderr, "\t  -o             : output bitmap file name (auto add .bmp extension)\n" );
//  fprintf( stderr, "\t  -x             : \n" );
//  fprintf( stderr, "\t  -orgx          : 水平画素サイズ(任意の自然数)\n" );
//  fprintf( stderr, "\t  -orgy          : 垂直画素サイズ(任意の自然数)\n" );
  fprintf( stderr, "\t  -sizexminus    : difference between input horizontal picture size and output one \n" );
  fprintf( stderr, "\t  -sizexminusy   : difference between input vertical picture size and output one\n" );
//  fprintf( stderr, "\t  -hv       : 1:水平 2:プログレッシブ垂直 3:インタレース垂直\n" );
//  fprintf( stderr, "\t  -bunshi   : 縮小率分子 max31\n" );
//  fprintf( stderr, "\t  -bunbo    : 縮小率分母 max32\n" );
//  fprintf( stderr, "\t  -tapnum_y : 輝度フィルタタップ数  推奨 21 max50\n" );
//  fprintf( stderr, "\t  -tapnum_c : 色差フィルタタップ数  推奨 11 max50\n" );
//  fprintf( stderr, "\t  -norm     : 総係数合計 推奨4096 サブサンプルプログラムで使う場合\n" );
//  fprintf( stderr, "\t  -mado     : 1:矩形 2:Lanczos 3:hanning(両端が0) 4:hamming\n" );
//  fprintf( stderr, "\t  -Fsc      : サブサンプル後のサンプリング周波数を１としたときカットオフ周波数を指定(推奨0.5)\n" );

}
void    option_set(short argc,char *argv[], CNTL *ptr)
{
  short i;
                                                                                                                                                         
  for(i=1;i < argc;i++) {
    if(argv[i][0] != '-')
      continue;

    if ( !strcmp( &argv[i][1], "h" ) ){
      /* \245\330\245\353\245\327\311\275\274\250 */
      print_help();
      exit(0);
    } else if ( !strcmp( &argv[i][1], "o" ) ){
      /* \275\320\316\317\245\325\245\241\245\244\245\353\314\276 */
      strcpy(ptr->OutputFileName, argv[++i]);
    } else if ( !strcmp( &argv[i][1], "i" ) ){
      strcpy(ptr->InputFileName, argv[++i]);
    } else if ( !strcmp( &argv[i][1], "i2" ) ){
      strcpy(ptr->InputFileName2, argv[++i]);
                                                                                                                                                         
    } else if ( !strcmp( &argv[i][1], "x" ) ){
      ptr->sizex = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "y" ) ){
      ptr->sizey = atoi(argv[++i]);
                                                                                                                                                         
    } else if ( !strcmp( &argv[i][1], "startx" ) ){
      ptr->startx = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "starty" ) ){
      ptr->starty = atoi(argv[++i]);
                                                                                                                                                         
    } else if ( !strcmp( &argv[i][1], "hv" ) ){
      ptr->hv = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "bunshi" ) ){
      ptr->bunshi = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "bunbo" ) ){
      ptr->bunbo = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "norm" ) ){
      ptr->norm = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "mado" ) ){
      ptr->mado = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "tapnum_y" ) ){
      ptr->tapnum_y = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "tapnum_c" ) ){
      ptr->tapnum_c = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "Fsc" ) ){
      ptr->Fs = atof(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "sizexminus" ) ){
      ptr->sm = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "sizexminusy" ) ){
      ptr->smy = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "orgx" ) ){
      /* \277\345\312\277\312\375\270\376\306\376\316\317\262\350\301\374\245\265\245\244\245\272 */
      org_image_x = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "orgy" ) ){
      /* \277\342\304\276\312\375\270\376\306\376\316\317\262\350\301\374\245\265\245\244\245\272 */
      org_image_y = atoi(argv[++i]);
    } else {
      print_help();
      exit(0);
    }
  }
                                                                                                                                                         
}

int main(int argc, char *argv[]){
FILE *fpr1;
FILE *fpr2;
FILE *fpw;
unsigned char *vin;
unsigned char *vout;
unsigned char *vin_buffer;
#define HSIZE 360
#define VSIZE 288
#define tmp HSIZE*VSIZE 
int i,j,k,l,m,n,o ;
CNTL  *cnt;
int center0;
int sizeminus;
int filesizebyte;
unsigned char *B;
unsigned char *G;
unsigned char *R;
unsigned char *BIN;
unsigned char *GIN;
unsigned char *RIN;
unsigned char *RUPOUT;
unsigned char *GUPOUT;
unsigned char *BUPOUT;
unsigned char temp;
unsigned char temp2;
unsigned int filesize;


  org_image_x = ORG_IMAGE_X;
  org_image_y = ORG_IMAGE_Y;

        if( argc < 2 ) {
                print_help();
                exit(0);
        }
        cnt = (CNTL *)malloc( sizeof(CNTL)*sizeof(uchar) );
                                                                                                                                                         
        /* \245\307\245\325\245\251\245\353\245\310\303\315\300\337\304\352 */
        cnt->sizex = ORG_IMAGE_X;    
        cnt->sizey = ORG_IMAGE_Y;     
        cnt->startx = 0;
        cnt->starty = 0;
                                                                                                                                                         
        /* \245\252\245\327\245\267\245\347\245\363\245\301\245\247\245\303\245\257 */
        option_set(argc, argv, cnt);
                                                                                                                                                         
                                                                                                                                                         
        hsize = org_image_x;
        vsize = org_image_y;
        strcpy(inputfile,cnt->InputFileName);
        strcpy(inputfile2,cnt->InputFileName2);
        strcpy(outputfile,cnt->OutputFileName);
	sizeminus = cnt->sm;
                                                                                                                                                         
/********************************/ 
/*				*/
/*	  MAIN ROUTINE		*/
/*				*/
/********************************/

//	vin_buffer = (unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);
//	vout = (unsigned char *)malloc(sizeof(unsigned char)*(hsize+sizeminus)*(vsize+smy)*2);

        if((fpr1 = fopen(inputfile,"rb"))==NULL){
                exit(0);
        }
/***************************************************/
/* picture size & file size get from bitmap header */
/***************************************************/
	fseek(fpr1,18,SEEK_SET);
	fread(&hsize,sizeof(int),1,fpr1);
	fread(&vsize,sizeof(int),1,fpr1);
	fseek(fpr1,2,SEEK_SET);
	fread(&filesizebyte,sizeof(int),1,fpr1);

/***************/
/* picture get */
/***************/
	vin = (unsigned char *)malloc(sizeof(unsigned char)*(filesizebyte-54));
	B = (unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);
	G = (unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);
	R = (unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);
	fseek(fpr1,54,SEEK_SET);
	if(fread(vin,sizeof(unsigned char),filesizebyte - 54,fpr1) != filesizebyte -54 ){
		printf("input file size err!!\n");
		exit(1);
	}
	fclose(fpr1);

/***************************************/
/* set to B,G,R array with RGB PICTURE */
/***************************************/
	for(i=0,k=0;i<vsize;i++){
		for(j=0;j<hsize;j++){
			*(B+j*vsize+i) = *(vin+k);
			k++;
			*(G+j*vsize+i) = *(vin+k);
			k++;
			*(R+j*vsize+i) = *(vin+k);
			k++;
		}
		if(((hsize*3)%4) == 1){
			k=k+3;
		}
		if(((hsize*3)%4) == 2){
			k=k+2;
		}
		if(((hsize*3)%4) == 3){
			k=k+1;
		}
		if(((hsize*3)%4) == 0){
			;
		}
	}
//debug start
//for(;;);
//debug end
	BIN = (unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);
	GIN = (unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);
	RIN = (unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);

/*******************/
/* RGB arrangement */	
/*******************/

	for(i=0;i<vsize;i++){
		for(j=0;j<hsize;j++){
			*(RIN+j*vsize+i) = *(R+j*vsize+vsize-1-i);
			*(GIN+j*vsize+i) = *(G+j*vsize+vsize-1-i);
			*(BIN+j*vsize+i) = *(B+j*vsize+vsize-1-i);	
		}
	}

	RUPOUT = (unsigned char *)malloc(sizeof(unsigned char)*(hsize+sizeminus)*(vsize+cnt->smy));
	GUPOUT = (unsigned char *)malloc(sizeof(unsigned char)*(hsize+sizeminus)*(vsize+cnt->smy));
	BUPOUT = (unsigned char *)malloc(sizeof(unsigned char)*(hsize+sizeminus)*(vsize+cnt->smy));
	printf("RED START\n");
	resize(  hsize,  vsize,  sizeminus,  cnt->smy,   RIN ,  RUPOUT);
	printf("GREEN START\n");
	resize(  hsize,  vsize,  sizeminus,  cnt->smy,   GIN ,  GUPOUT);
	printf("BLUE START\n");
	resize(  hsize,  vsize,  sizeminus,  cnt->smy,   BIN ,  BUPOUT);

	sprintf(cnt->OutputFileName,"%s.bmp",outputfile);
	fpw = fopen(cnt->OutputFileName,"wb");
	hsize = hsize + sizeminus;
	vsize = vsize + cnt->smy;
/************************/
/* HEADER 54 byte write */
/************************/
	if(((hsize*3)%4) == 0){
		filesize = 0x36 + hsize*vsize*3;
	} else {
		filesize = 0x36 + hsize*vsize*3+(4-(hsize*3)%4)*vsize;
	}

	temp = 0x42; //'B'
	fwrite(&temp,sizeof(unsigned char),1,fpw);
	temp = 0x4D; //'M'
	fwrite(&temp,sizeof(unsigned char),1,fpw);
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	filesize = 0;
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	filesize = 0x00000036;
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	filesize = 0x28;
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	fwrite(&hsize,sizeof(int),1,fpw);	
	fwrite(&vsize,sizeof(int),1,fpw);
	filesize = 0x00180001;
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	filesize = 0;
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	filesize = 0x00000EC4;
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	filesize = 0;
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	temp2 =0;

	for(i=0;i<vsize;i++){
		for(j=0;j<hsize;j++){
			fwrite((BUPOUT+j*vsize+vsize-1-i),sizeof(unsigned char),1,fpw);
			fwrite((GUPOUT+j*vsize+vsize-1-i),sizeof(unsigned char),1,fpw);
			fwrite((RUPOUT+j*vsize+vsize-1-i),sizeof(unsigned char),1,fpw);
		}
		temp = 4-((hsize*3)%4) ;
		if(temp != 4 ){
			for(k=0;k<temp;k++){
				fwrite( &temp2,sizeof(unsigned char),1,fpw);
			}
		}
	}	

	fclose(fpw);
        return 0;
}

void resize( int hsize, int vsize, int sizeminus, int smy,  unsigned char *vin_buffer , unsigned char *vout)
{ 

double *VINR;
double *VINI;
double *VOUTR;
double *VOUTI;
double *VINRptr;
double *VINIptr;
double *VOUTRptr;
double *VOUTIptr;
double *kariR;
double *kariI;
int i,j,k,l,m,n,o ;
double PI;
double *fourieR;
double *fourieI;
	PI = atan(1)*4.0;
	kariR = (double *)malloc(sizeof(double)*(hsize+sizeminus)*(vsize+smy)*4);
	kariI = (double *)malloc(sizeof(double)*(hsize+sizeminus)*(vsize+smy)*4);
	VINR = (double *)malloc(sizeof(double)*(hsize+sizeminus)*(vsize+smy)*4);
	VINI = (double *)malloc(sizeof(double)*(hsize+sizeminus)*(vsize+smy)*4);
	VOUTR = (double *)malloc(sizeof(double)*(hsize+sizeminus)*(vsize+smy)*4);
	VOUTI = (double *)malloc(sizeof(double)*(hsize+sizeminus)*(vsize+smy)*4);
	VINRptr = VINR;
	VINIptr = VINI;
	VOUTRptr = VOUTR;
	VOUTIptr = VOUTI;
/************************/
/*			*/
/*      BUFFERING	*/
/*			*/
/************************/
//	for(i=0,k=0;i<vsize;i++){
//		for(j=0;j<hsize;j++){
//			*(vin_buffer+j*vsize+i) = *(vin+k);
//			k++;
//		}
//	}
/*************************/
/*			 */
/*     第一象限代入	 */
/*			 */
/*************************/
	for(i=0;i<vsize;i++){
		for(j=0;j<hsize;j++){
			*(VINR+j*2*vsize+i) = (double)(*(vin_buffer+j*vsize+i));
			*(VINI+j*2*vsize+i) = 0.0;
		}
	}
/*************************/
/*			 */
/*     第二象限代入	 */
/*			 */
/*************************/
	for(i=0;i<vsize;i++){
		for(j=0;j<hsize;j++){
			*(VINR+(j+hsize)*2*vsize+i) = (double)(*(vin_buffer+(hsize-1-j)*vsize+i));
			*(VINI+(j+hsize)*2*vsize+i) = 0.0; 
		}
	}
/*************************/
/*			 */
/*     第三象限代入	 */
/*			 */
/*************************/
	for(i=0;i<vsize;i++){
		for(j=0;j<hsize;j++){
			*(VINR+j*2*vsize+i+vsize) = (double)(*(vin_buffer+j*vsize+(vsize-1-i)));
			*(VINI+j*2*vsize+i+vsize) = 0.0;
		}
	}
/*************************/
/*			 */
/*     第四象限代入	 */
/*			 */
/*************************/
	for(i=0;i<vsize;i++){
		for(j=0;j<hsize;j++){
			*(VINR+(j+hsize)*2*vsize+i+vsize) = (double)(*(vin_buffer+(hsize-1-j)*vsize+(vsize-1-i)));
			*(VINI+(j+hsize)*2*vsize+i+vsize) = 0.0;
		}
	}

//debug
//for(;;);
//------
	fprintf(stderr,"horizontal DFT start!!\n");
	double *ww;
	ww = (double *)malloc(sizeof(double)*8192*4);
	if(ww==NULL){printf("ww no memory!!\n");exit(1);}
	int *ip;
	ip = (int *)malloc(sizeof(int)*8192*4);
	if(ip==NULL){printf("ip no memory!!\n");exit(1);}

	for(i=0;i<vsize*2;i++){
		if(i==0){
			ip[0]=0;
			DFT1DH (hsize*2,vsize*2,VINRptr+i,VINIptr+i,VOUTRptr+i,VOUTIptr+i,ip,ww);
		} else {
			DFT1DH (hsize*2,vsize*2,VINRptr+i,VINIptr+i,VOUTRptr+i,VOUTIptr+i,ip,ww);
		}
		//printf("DFTH %d\n",i);
	}

	fprintf(stderr,"horizontal DFT end\n");
	for(i=0;i<vsize*hsize*4;i++){
		*(VINR+i) = *(VOUTR+i);
		*(VINI+i) = *(VOUTI+i);
		//debug start
//		fprintf(stderr,"------%d=%f\n",i,*(VOUTR+i));
		//debug end
	}
//debug
//for(;;);
//------

	fprintf(stderr,"vertical DFT start!!\n");

	for(i=0;i<hsize*2;i++){
		if(i==0){
		ip[0] =0;
			DFT1DV (vsize*2,VINRptr+(vsize*2)*i,VINIptr+(vsize*2)*i,VOUTRptr+(vsize*2)*i,VOUTIptr+(vsize*2)*i,ip,ww);
		} else {
			DFT1DV (vsize*2,VINRptr+(vsize*2)*i,VINIptr+(vsize*2)*i,VOUTRptr+(vsize*2)*i,VOUTIptr+(vsize*2)*i,ip,ww);
		}
		//printf("DFTV %d\n",i);
	}
	fprintf(stderr,"vertical DFT end!!\n");


//	fourieR = (double *)malloc(sizeof(double)*(hsize-sizeminus)*vsize);
//	fourieI = (double *)malloc(sizeof(double)*(hsize-sizeminus)*vsize);
	fourieR = (double *)malloc(sizeof(double)*(hsize+sizeminus)*(vsize+smy)*4);
	fourieI = (double *)malloc(sizeof(double)*(hsize+sizeminus)*(vsize+smy)*4);

/*縮小した画像のフーリエ変換像の配列作成*/

	fprintf(stderr,"high frequency 0 staffing start!!\n");
//	for(i=0;i<vsize*2;i++){
//		*(VOUTRptr+((hsize+sizeminus))*vsize*2+i) = *(VOUTRptr+((hsize+sizeminus))*vsize*2+i) * 2.0;
//		*(VOUTIptr+((hsize+sizeminus))*vsize*2+i) = 0.0;	
//	}

	for(i=0;i<vsize*2;i++){
		for(j=0;j<hsize;j++){
			*(fourieR+j*vsize*2+i)=*(VOUTRptr+j*vsize*2+i);
			*(fourieI+j*vsize*2+i)=*(VOUTIptr+j*vsize*2+i);
		}
		for(j=0;j<sizeminus*2;j++){
			*(fourieR+(j+hsize)*vsize*2+i) = 0.0;
			*(fourieI+(j+hsize)*vsize*2+i) = 0.0;
		}
		for(j=0;j<hsize;j++){
			*(fourieR+((hsize+sizeminus*2)+j)*vsize*2+i) = *(VOUTRptr+(hsize+j)*vsize*2+i);
			*(fourieI+((hsize+sizeminus*2)+j)*vsize*2+i) = *(VOUTIptr+(hsize+j)*vsize*2+i);
		}
	}

//	for(j=0;j<(hsize-sizeminus);j++){
//		*(fourieR+j*vsize*2+(vsize+smy)) = *(fourieR+j*vsize*2+(vsize+smy)) * 2.0;
//		*(fourieI+j*vsize*2+(vsize+smy)) = 0.0;
//	}

	for(j=0;j<2*(hsize+sizeminus);j++){
		for(i=0;i<vsize;i++){
			*(kariR+j*2*(vsize+(smy))+i)=*(fourieR+j*2*vsize+i);
			*(kariI+j*2*(vsize+(smy))+i)=*(fourieI+j*2*vsize+i);
		}
		for(i=0;i<(smy)*2;i++){
			*(kariR+j*2*(vsize+(smy))+i+vsize) = 0.0;
			*(kariI+j*2*(vsize+(smy))+i+vsize) = 0.0;
		}
		for(i=0;i<vsize;i++){
			*(kariR+j*2*(vsize+(smy))+i+(vsize+(smy)*2)) =*(fourieR+j*2*vsize+i+vsize);
			*(kariI+j*2*(vsize+(smy))+i+(vsize+(smy)*2)) =*(fourieI+j*2*vsize+i+vsize);
		}
	}
//debug start
//for(i=0;i<2*(vsize-smy);i++){
//	for(j=0;j<2*(hsize-sizeminus);j++){
//		printf("%f\n",*(kariR+j*2*(vsize-smy)+i));
//	}
//}	
//debug end
	for(i=0;i<2*(vsize+smy);i++){
		for(j=0;j<2*(hsize+sizeminus);j++){
			*(fourieR+j*2*(vsize+smy)+i) = *(kariR+j*2*(vsize+smy)+i);
			*(fourieI+j*2*(vsize+smy)+i) = *(kariI+j*2*(vsize+smy)+i);
		}
	}
	fprintf(stderr,"high frequency 0 staffing end!!\n");
//debug start
//for(i=0;i<2*(vsize-smy);i++){
//	for(j=0;j<2*(hsize-sizeminus);j++){
//		printf("%f\n",*(fourieR+j*2*(vsize-smy)+i));
//	}
//}	
//debug end
	fprintf(stderr,"horizontal IDFT start!!\n");
	
	for(i=0;i<2*(vsize+smy);i++){
		if(i==0){
		ip[0]=0;
		IDFT1DH (2*(hsize+sizeminus),2*(vsize+smy),fourieR+i,fourieI+i,VOUTRptr+i,VOUTIptr+i,ip,ww);
		} else {
		IDFT1DH (2*(hsize+sizeminus),2*(vsize+smy),fourieR+i,fourieI+i,VOUTRptr+i,VOUTIptr+i,ip,ww);
		}
		//printf("IDFTH %d\n",i);
	}
	fprintf(stderr,"horizontal IDFT end!!\n");
//debug start
//for(i=0;i<2*(vsize-smy);i++){
//	for(j=0;j<2*(hsize-sizeminus);j++){
//		printf("%f\n",*(VOUTRptr+j*2*(vsize-smy)+i));
//	}
//}	
//debug end
	fprintf(stderr,"vertical IDFT start!!\n");
	for(i=0;i<2*(hsize+sizeminus);i++){
		if(i==0){
			ip[0]=0;
		IDFT1DV (2*(vsize+smy),VOUTRptr+2*(vsize+smy)*i,VOUTIptr+2*(vsize+smy)*i,fourieR+2*(vsize+smy)*i,fourieI+2*(vsize+smy)*i,ip,ww);
		} else {
		IDFT1DV (2*(vsize+smy),VOUTRptr+2*(vsize+smy)*i,VOUTIptr+2*(vsize+smy)*i,fourieR+2*(vsize+smy)*i,fourieI+2*(vsize+smy)*i,ip,ww);
		}
		//printf("IDFTV %d\n",i);
	}
	fprintf(stderr,"vertical IDFT end!!\n");
	
//	fprintf(stderr,"出力中。\n");
	
	for(i=0,k=0;i<vsize+smy;i++){
		for(j=0;j<(hsize+sizeminus);j++){
				if(*(fourieR+j*2*(vsize+smy)+i) >= 255.0){
				*(vout+j*(vsize+smy)+i) = 255;
				} else if(*(fourieR+j*2*(vsize+smy)+i) <= 0.0){
				*(vout+j*(vsize+smy)+i) = 0;
			} else {
				*(vout+j*(vsize+smy)+i) = (unsigned char)(*(fourieR+j*2*(vsize+smy)+i));///**(double)(300*240/*hsize*vsize*/)/*/(16.0*4096.0)*//**((double)(hsize+sizeminus))*(double)(vsize+smy)*/);

				
				
				
				
				
				
				//				*(vout+k) = (unsigned char)(*(fourieR+j*(vsize-smy)+i)/(double)(hsize*vsize)/((double)(vsize-smy) * (double)(hsize-sizeminus)));
//				*(vout+k) = (unsigned char)(*(fourieR+j*(vsize-smy)+i)/(double)(hsize*vsize)/(double)(hsize*vsize)*((double)(vsize-smy) * (double)(hsize-sizeminus)));
//				*(vout+k) = (unsigned char)((*(fourieR+j*(vsize-smy)+i))/(sqrt((double)(hsize*vsize)*((double)(vsize-smy) * (double)(hsize-sizeminus)))));
//				*(vout+k) = (unsigned char)((*(fourieR+j*(vsize-smy)+i))*((double)(hsize*vsize)));
			}
//debug
//printf("%f\n",((*(fourieR+j*2*(vsize-smy)+i))/(double)(hsize*vsize*4)));
//
			k++;
		}
	}
}// resize kansuu end	

void DFT1DV(int size, double *VINR, double *VINI, double *VOUTR, double *VOUTI,int *ip,double *ww)
{
	int i,j;
	double PI = atan(1)*4.0;
	for(i=0;i<size;i++){
		*(VOUTR+i) = 0.0;
		*(VOUTI+i) = 0.0;
	}
	charp_FFT(size,VINR,VINI,VOUTR,VOUTI,ip,ww);
//	*(VOUTR) /= 1.0;
//	for(i=0;i<size;i++){
//		for(j=0;j<size;j++){
//			*(VOUTR+i) += (*(VINR+j)*cos(2.0*PI*(double)i*(double)j/(double)size) + *(VINI+j)*sin(2.0*PI*(double)i*(double)j/(double)size));
//			*(VOUTI+i) += (*(VINI+j)*cos(2.0*PI*(double)i*(double)j/(double)size) - *(VINR+j)*sin(2.0*PI*(double)i*(double)j/(double)size));
//		}
//	}	
		
}

void DFT1DH(int hsize, int vsize , double *VINR, double*VINI, double *VOUTR, double *VOUTI,int *ip,double *ww)
{
	int i,j;
        double PI = atan(1)*4.0;
        for(i=0;i<hsize;i++){
                *(VOUTR+i*vsize) = 0.0;
                *(VOUTI+i*vsize) = 0.0;
        }
	double *R,*I,*OR,*OI;
	R= (double*)malloc(sizeof(double)*hsize);
	I= (double*)malloc(sizeof(double)*hsize);
	OR= (double*)malloc(sizeof(double)*hsize);
	OI= (double*)malloc(sizeof(double)*hsize);
	for(i=0;i<hsize;i++){
		*(R+i)= *(VINR+i*vsize);
		*(I+i)= *(VINI+i*vsize);
	}
	charp_FFT(hsize,R,I,OR,OI,ip,ww);
	for(i=0;i<hsize;i++){
		*(VOUTR+i*vsize)=*(OR+i);
		*(VOUTI+i*vsize)=*(OI+i);
//debug start
//		fprintf(stderr,"%d======%f\n",i,*(VOUTR+i*vsize));
//debug end
	}
//	*(VOUTR) /= 1.0;
	free(R);free(I);free(OR);free(OI);
//	for(i=0;i<hsize;i++){
  //              for(j=0;j<hsize;j++){
    //                    *(VOUTR+i*vsize) += (*(VINR+j*vsize)*cos(2.0*PI*(double)i*(double)j/(double)hsize) + *(VINI+j*vsize)*sin(2.0*PI*(double)i*(double)j/(double)hsize));
      //                  *(VOUTI+i*vsize) += (*(VINI+j*vsize)*cos(2.0*PI*(double)i*(double)j/(double)hsize) - *(VINR+j*vsize)*sin(2.0*PI*(double)i*(double)j/(double)hsize));
        //        }
//        }
}
void IDFT1DV(int size, double *VINR, double *VINI, double *VOUTR, double *VOUTI,int *ip,double *ww)
{
	int i,j;
	double PI = atan(1)*4.0;
	for(i=0;i<size;i++){
		*(VOUTR+i) = 0.0;
		*(VOUTI+i) = 0.0;
	}
	charp_IFFT(size,VINR,VINI,VOUTR,VOUTI,ip,ww);
//	for(i=0;i<size;i++){
//		for(j=0;j<size;j++){
//			*(VOUTR+i) += (*(VINR+j)*cos(2.0*PI*(double)i*(double)j/(double)size) - *(VINI+j)*sin(2.0*PI*(double)i*(double)j/(double)size));
//			*(VOUTI+i) += (*(VINI+j)*cos(2.0*PI*(double)i*(double)j/(double)size) + *(VINR+j)*sin(2.0*PI*(double)i*(double)j/(double)size));
//		}
//	}	
		
}

void IDFT1DH(int hsize, int vsize , double *VINR, double*VINI, double *VOUTR, double *VOUTI,int *ip,double *ww)
{
	int i,j;
        double PI = atan(1)*4.0;
        for(i=0;i<hsize;i++){
                *(VOUTR+i*vsize) = 0.0;
                *(VOUTI+i*vsize) = 0.0;
        }
	double *R,*I,*OR,*OI;
	R= (double*)malloc(sizeof(double)*hsize);
	I= (double*)malloc(sizeof(double)*hsize);
	OR= (double*)malloc(sizeof(double)*hsize);
	OI= (double*)malloc(sizeof(double)*hsize);
	for(i=0;i<hsize;i++){
		*(R+i)= *(VINR+i*vsize);
		*(I+i)= *(VINI+i*vsize);
	}
	charp_IFFT(hsize,R,I,OR,OI,ip,ww);
	for(i=0;i<hsize;i++){
		*(VOUTR+i*vsize)=*(OR+i);
		*(VOUTI+i*vsize)=*(OI+i);
		
	}
	free(R);free(I);free(OR);free(OI);
//	for(i=0;i<hsize;i++){
  //              for(j=0;j<hsize;j++){
    //                    *(VOUTR+i*vsize) += (*(VINR+j*vsize)*cos(2.0*PI*(double)i*(double)j/(double)hsize) - *(VINI+j*vsize)*sin(2.0*PI*(double)i*(double)j/(double)hsize));
      //                  *(VOUTI+i*vsize) += (*(VINI+j*vsize)*cos(2.0*PI*(double)i*(double)j/(double)hsize) + *(VINR+j*vsize)*sin(2.0*PI*(double)i*(double)j/(double)hsize));
        //        }
  //}
}