/*Fast Fourier/Cosine/Sine Transform    dimension   :one    data length :power of 2    decimation  :frequency    radix       :4, 2    data        :inplace    table       :usefunctions    cdft: Complex Discrete Fourier Transform    rdft: Real Discrete Fourier Transform    function prototypes    void cdft(int, int, double *, int *, double *);    void rdft(int, int, double *, int *, double *);    -------- 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;            }        .Appendix :    The cos/sin table is recalculated when the larger table required.    w[] and ip[] are compatible with all routines.*/void cdft(int n, int isgn, double *a, int *ip, double *w){    void makewt(int nw, int *ip, double *w);    void bitrv2(int n, int *ip, double *a);    void bitrv2conj(int n, int *ip, double *a);    void cftfsub(int n, double *a, double *w);    void cftbsub(int n, double *a, double *w);        if (n > (ip[0] << 2))
	{        makewt(n >> 2, ip, w);    }    if (n > 4) 
	{        if (isgn >= 0) 
		{            bitrv2(n, ip + 2, a);            cftfsub(n, a, w);        } 
		else 
		{            bitrv2conj(n, ip + 2, a);            cftbsub(n, a, w);        }    } 
	else if (n == 4)
	{        cftfsub(n, a, 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 bitrv2(int n, int *ip, double *a);    void cftfsub(int n, double *a, double *w);    void cftbsub(int n, double *a, 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) 
		{            bitrv2(n, ip + 2, a);            cftfsub(n, a, w);            rftfsub(n, a, nc, w + nw);        } 
		else if (n == 4) 
		{            cftfsub(n, a, 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);            bitrv2(n, ip + 2, a);            cftbsub(n, a, w);        } 
		else if (n == 4) 
		{            cftfsub(n, a, w);        }    }}
/* -------- initializing routines -------- */#include <math.h>void makewt(int nw, int *ip, double *w){    void bitrv2(int n, int *ip, double *a);    int j, nwh;    double delta, x, y;        ip[0] = nw;    ip[1] = 1;    if (nw > 2) 
	{        nwh = nw >> 1;        delta = atan(1.0) / nwh;        w[0] = 1;        w[1] = 0;        w[nwh] = cos(delta * nwh);        w[nwh + 1] = w[nwh];        if (nwh > 2) 
		{            for (j = 2; j < nwh; j += 2) 
			{                x = cos(delta * j);                y = sin(delta * j);                w[j] = x;                w[j + 1] = y;                w[nw - j] = y;                w[nw - j + 1] = x;            }            bitrv2(nw, ip + 2, w);        }    }}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 -------- */void bitrv2(int n, int *ip, double *a){    int j, j1, k, k1, l, m, m2;    double xr, xi, yr, yi;        ip[0] = 0;    l = n;    m = 1;    while ((m << 3) < l) 
	{        l >>= 1;        for (j = 0; j < m; j++) 
		{            ip[m + j] = ip[j] + l;        }        m <<= 1;    }    m2 = 2 * m;    if ((m << 3) == l) 
	{        for (k = 0; k < m; k++) 
		{            for (j = 0; j < k; j++) 
			{                j1 = 2 * j + ip[k];                k1 = 2 * k + ip[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 += m2;                k1 += 2 * m2;                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 += m2;                k1 -= m2;                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 += m2;                k1 += 2 * m2;                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 * k + m2 + ip[k];            k1 = j1 + m2;            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 = 1; k < m; k++) 
		{            for (j = 0; j < k; j++) 
			{                j1 = 2 * j + ip[k];                k1 = 2 * k + ip[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 += m2;                k1 += m2;                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, m2;    double xr, xi, yr, yi;        ip[0] = 0;    l = n;    m = 1;    while ((m << 3) < l) 
	{        l >>= 1;        for (j = 0; j < m; j++) 
		{            ip[m + j] = ip[j] + l;        }        m <<= 1;    }    m2 = 2 * m;    if ((m << 3) == l) 
	{        for (k = 0; k < m; k++) 
		{            for (j = 0; j < k; j++) 
			{                j1 = 2 * j + ip[k];                k1 = 2 * k + ip[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 += m2;                k1 += 2 * m2;                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 += m2;                k1 -= m2;                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 += m2;                k1 += 2 * m2;                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 = 2 * k + ip[k];            a[k1 + 1] = -a[k1 + 1];            j1 = k1 + m2;            k1 = j1 + m2;            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;