x0i = y1i + y5i;            a[j + 2] = wk1r * x0r - wk1i * x0i;            a[j + 3] = wk1r * x0i + wk1i * x0r;            x0r = y1r - y5r;            x0i = y1i - y5i;            a[j + 10] = wk5r * x0r - wk5i * x0i;            a[j + 11] = wk5r * x0i + wk5i * x0r;            x0r = y3r - y7i;            x0i = y3i + y7r;            a[j + 6] = wk3r * x0r - wk3i * x0i;            a[j + 7] = wk3r * x0i + wk3i * x0r;            x0r = y3r + y7i;            x0i = y3i - y7r;            a[j + 14] = wk7r * x0r - wk7i * x0i;            a[j + 15] = wk7r * x0i + wk7i * x0r;            a[j] = y0r + y4r;            a[j + 1] = y0i + y4i;            x0r = y0r - y4r;            x0i = y0i - y4i;            a[j + 8] = wk4r * x0r - wk4i * x0i;            a[j + 9] = wk4r * x0i + wk4i * x0r;            x0r = y2r - y6i;            x0i = y2i + y6r;            a[j + 4] = wk2r * x0r - wk2i * x0i;            a[j + 5] = wk2r * x0i + wk2i * x0r;            x0r = y2r + y6i;            x0i = y2i - y6r;            a[j + 12] = wk6r * x0r - wk6i * x0i;            a[j + 13] = wk6r * x0i + wk6i * x0r;        }    }}void cftmdl(int n, int l, double *a){    int j, j1, j2, j3, j4, j5, j6, j7, k, kj, kr, m;    double ew, wn4r, wtmp, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i,         wk4r, wk4i, wk5r, wk5i, wk6r, wk6i, wk7r, wk7i;    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;        m = l << 3;    wn4r = WR5000;    for (j = 0; j < l; j += 2) {        j1 = j + l;        j2 = j1 + l;        j3 = j2 + l;        j4 = j3 + l;        j5 = j4 + l;        j6 = j5 + l;        j7 = j6 + l;        x0r = a[j] + a[j1];        x0i = a[j + 1] + a[j1 + 1];        x1r = a[j] - a[j1];        x1i = a[j + 1] - a[j1 + 1];        x2r = a[j2] + a[j3];        x2i = a[j2 + 1] + a[j3 + 1];        x3r = a[j2] - a[j3];        x3i = a[j2 + 1] - a[j3 + 1];        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[j4] + a[j5];        x0i = a[j4 + 1] + a[j5 + 1];        x1r = a[j4] - a[j5];        x1i = a[j4 + 1] - a[j5 + 1];        x2r = a[j6] + a[j7];        x2i = a[j6 + 1] + a[j7 + 1];        x3r = a[j6] - a[j7];        x3i = a[j6 + 1] - a[j7 + 1];        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[j1] = y1r + y5r;        a[j1 + 1] = y1i + y5i;        a[j5] = y1r - y5r;        a[j5 + 1] = y1i - y5i;        a[j3] = y3r - y7i;        a[j3 + 1] = y3i + y7r;        a[j7] = y3r + y7i;        a[j7 + 1] = y3i - y7r;        a[j] = y0r + y4r;        a[j + 1] = y0i + y4i;        a[j4] = y0r - y4r;        a[j4 + 1] = y0i - y4i;        a[j2] = y2r - y6i;        a[j2 + 1] = y2i + y6r;        a[j6] = y2r + y6i;        a[j6 + 1] = y2i - y6r;    }    if (m < n) {        wk1r = WR2500;        wk1i = WI2500;        for (j = m; j < l + m; j += 2) {            j1 = j + l;            j2 = j1 + l;            j3 = j2 + l;            j4 = j3 + l;            j5 = j4 + l;            j6 = j5 + l;            j7 = j6 + l;            x0r = a[j] + a[j1];            x0i = a[j + 1] + a[j1 + 1];            x1r = a[j] - a[j1];            x1i = a[j + 1] - a[j1 + 1];            x2r = a[j2] + a[j3];            x2i = a[j2 + 1] + a[j3 + 1];            x3r = a[j2] - a[j3];            x3i = a[j2 + 1] - a[j3 + 1];            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[j4] + a[j5];            x0i = a[j4 + 1] + a[j5 + 1];            x1r = a[j4] - a[j5];            x1i = a[j4 + 1] - a[j5 + 1];            x2r = a[j6] + a[j7];            x2i = a[j6 + 1] + a[j7 + 1];            x3r = a[j6] - a[j7];            x3i = a[j6 + 1] - a[j7 + 1];            y4r = x0r + x2r;            y4i = x0i + x2i;            y6r = x0r - x2r;            y6i = x0i - x2i;            x0r = x1r - x3i;            x0i = x1i + x3r;            x2r = x1r + x3i;            x2i = x3r - x1i;            y5r = wk1i * x0r - wk1r * x0i;            y5i = wk1i * x0i + wk1r * x0r;            y7r = wk1r * x2r + wk1i * x2i;            y7i = wk1r * x2i - wk1i * x2r;            x0r = wk1r * y1r - wk1i * y1i;            x0i = wk1r * y1i + wk1i * y1r;            a[j1] = x0r + y5r;            a[j1 + 1] = x0i + y5i;            a[j5] = y5i - x0i;            a[j5 + 1] = x0r - y5r;            x0r = wk1i * y3r - wk1r * y3i;            x0i = wk1i * y3i + wk1r * y3r;            a[j3] = x0r - y7r;            a[j3 + 1] = x0i + y7i;            a[j7] = y7i - x0i;            a[j7 + 1] = x0r + y7r;            a[j] = y0r + y4r;            a[j + 1] = y0i + y4i;            a[j4] = y4i - y0i;            a[j4 + 1] = y0r - y4r;            x0r = y2r - y6i;            x0i = y2i + y6r;            a[j2] = wn4r * (x0r - x0i);            a[j2 + 1] = wn4r * (x0i + x0r);            x0r = y6r - y2i;            x0i = y2r + y6i;            a[j6] = wn4r * (x0r - x0i);            a[j6 + 1] = wn4r * (x0i + x0r);        }        ew = M_PI_2 / n;        kr = n >> 2;        for (k = 2 * m; k < n; k += m) {            for (kj = n >> 2; kj > (kr ^= kj); kj >>= 1);            wk1r = cos(ew * kr);            wk1i = sin(ew * kr);            wk2r = 1 - 2 * wk1i * wk1i;            wk2i = 2 * wk1i * wk1r;            wtmp = 2 * wk2i;            wk3r = wk1r - wtmp * wk1i;            wk3i = wtmp * wk1r - wk1i;            wk4r = 1 - wtmp * wk2i;            wk4i = wtmp * wk2r;            wtmp = 2 * wk4i;            wk5r = wk3r - wtmp * wk1i;            wk5i = wtmp * wk1r - wk3i;            wk6r = wk2r - wtmp * wk2i;            wk6i = wtmp * wk2r - wk2i;            wk7r = wk1r - wtmp * wk3i;            wk7i = wtmp * wk3r - wk1i;            for (j = k; j < l + k; j += 2) {                j1 = j + l;                j2 = j1 + l;                j3 = j2 + l;                j4 = j3 + l;                j5 = j4 + l;                j6 = j5 + l;                j7 = j6 + l;                x0r = a[j] + a[j1];                x0i = a[j + 1] + a[j1 + 1];                x1r = a[j] - a[j1];                x1i = a[j + 1] - a[j1 + 1];                x2r = a[j2] + a[j3];                x2i = a[j2 + 1] + a[j3 + 1];                x3r = a[j2] - a[j3];                x3i = a[j2 + 1] - a[j3 + 1];                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[j4] + a[j5];                x0i = a[j4 + 1] + a[j5 + 1];                x1r = a[j4] - a[j5];                x1i = a[j4 + 1] - a[j5 + 1];                x2r = a[j6] + a[j7];                x2i = a[j6 + 1] + a[j7 + 1];                x3r = a[j6] - a[j7];                x3i = a[j6 + 1] - a[j7 + 1];                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);                x0r = y1r + y5r;                x0i = y1i + y5i;                a[j1] = wk1r * x0r - wk1i * x0i;                a[j1 + 1] = wk1r * x0i + wk1i * x0r;                x0r = y1r - y5r;                x0i = y1i - y5i;                a[j5] = wk5r * x0r - wk5i * x0i;                a[j5 + 1] = wk5r * x0i + wk5i * x0r;                x0r = y3r - y7i;                x0i = y3i + y7r;                a[j3] = wk3r * x0r - wk3i * x0i;                a[j3 + 1] = wk3r * x0i + wk3i * x0r;                x0r = y3r + y7i;                x0i = y3i - y7r;                a[j7] = wk7r * x0r - wk7i * x0i;                a[j7 + 1] = wk7r * x0i + wk7i * x0r;                a[j] = y0r + y4r;                a[j + 1] = y0i + y4i;                x0r = y0r - y4r;                x0i = y0i - y4i;                a[j4] = wk4r * x0r - wk4i * x0i;                a[j4 + 1] = wk4r * x0i + wk4i * x0r;                x0r = y2r - y6i;                x0i = y2i + y6r;                a[j2] = wk2r * x0r - wk2i * x0i;                a[j2 + 1] = wk2r * x0i + wk2i * x0r;                x0r = y2r + y6i;                x0i = y2i - y6r;                a[j6] = wk6r * x0r - wk6i * x0i;                a[j6 + 1] = wk6r * x0i + wk6i * x0r;            }        }    }}void rftfsub(int n, double *a){    int i, i0, j, k;    double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;        ec = 2 * M_PI_2 / n;    wkr = 0;    wki = 0;    wdi = cos(ec);    wdr = sin(ec);    wdi *= wdr;    wdr *= wdr;    w1r = 1 - 2 * wdr;    w1i = 2 * wdi;    ss = 2 * w1i;    i = n >> 1;    for (;;) {        i0 = i - 4 * RDFT_LOOP_DIV;        if (i0 < 4) {            i0 = 4;        }        for (j = i - 4; j >= i0; j -= 4) {            k = n - j;            xr = a[j + 2] - a[k - 2];            xi = a[j + 3] + a[k - 1];            yr = wdr * xr - wdi * xi;            yi = wdr * xi + wdi * xr;            a[j + 2] -= yr;            a[j + 3] -= yi;            a[k - 2] += yr;            a[k - 1] -= yi;            wkr += ss * wdi;            wki += ss * (0.5 - wdr);            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;            wdr += ss * wki;            wdi += ss * (0.5 - wkr);        }        if (i0 == 4) {            break;        }        wkr = 0.5 * sin(ec * i0);        wki = 0.5 * cos(ec * i0);        wdr = 0.5 - (wkr * w1r - wki * w1i);        wdi = wkr * w1i + wki * w1r;        wkr = 0.5 - wkr;        i = i0;    }    xr = a[2] - a[n - 2];    xi = a[3] + a[n - 1];    yr = wdr * xr - wdi * xi;    yi = wdr * xi + wdi * xr;    a[2] -= yr;    a[3] -= yi;    a[n - 2] += yr;    a[n - 1] -= yi;}void rftbsub(int n, double *a){    int i, i0, j, k;    double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;        ec = 2 * M_PI_2 / n;    wkr = 0;    wki = 0;    wdi = cos(ec);    wdr = sin(ec);    wdi *= wdr;    wdr *= wdr;    w1r = 1 - 2 * wdr;    w1i = 2 * wdi;    ss = 2 * w1i;    i = n >> 1;    a[i + 1] = -a[i + 1];    for (;;) {        i0 = i - 4 * RDFT_LOOP_DIV;        if (i0 < 4) {            i0 = 4;        }        for (j = i - 4; j >= i0; j -= 4) {            k = n - j;            xr = a[j + 2] - a[k - 2];            xi = a[j + 3] + a[k - 1];            yr = wdr * xr + wdi * xi;            yi = wdr * xi - wdi * xr;            a[j + 2] -= yr;            a[j + 3] = yi - a[j + 3];            a[k - 2] += yr;            a[k - 1] = yi - a[k - 1];            wkr += ss * wdi;            wki += ss * (0.5 - wdr);            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[j + 1];            a[k] += yr;            a[k + 1] = yi - a[k + 1];            wdr += ss * wki;            wdi += ss * (0.5 - wkr);        }        if (i0 == 4) {            break;        }        wkr = 0.5 * sin(ec * i0);        wki = 0.5 * cos(ec * i0);        wdr = 0.5 - (wkr * w1r - wki * w1i);        wdi = wkr * w1i + wki * w1r;        wkr = 0.5 - wkr;        i = i0;    }    xr = a[2] - a[n - 2];    xi = a[3] + a[n - 1];    yr = wdr * xr + wdi * xi;    yi = wdr * xi - wdi * xr;    a[2] -= yr;    a[3] = yi - a[3];    a[n - 2] += yr;    a[n - 1] = yi - a[n - 1];    a[1] = -a[1];}void dctsub(int n, double *a){    int i, i0, j, k, m;    double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;        ec = M_PI_2 / n;    wkr = 0.5;    wki = 0.5;    w1r = cos(ec);    w1i = sin(ec);    wdr = 0.5 * (w1r - w1i);    wdi = 0.5 * (w1r + w1i);    ss = 2 * w1i;    m = n >> 1;    i = 0;    for (;;) {        i0 = i + 2 * DCST_LOOP_DIV;        if (i0 > m - 2) {            i0 = m - 2;        }        for (j = i + 2; j <= i0; j += 2) {            k = n - j;            xr = wdi * a[j - 1] - wdr * a[k + 1];            xi = wdr * a[j - 1] + wdi * a[k + 1];            wkr -= ss * wdi;            wki += ss * wdr;            yr = wki * a[j] - wkr * a[k];            yi = wkr * a[j] + wki * a[k];            wdr -= ss * wki;            wdi += ss * wkr;            a[k + 1] = xr;            a[k] = yr;            a[j - 1] = xi;            a[j] = yi;        }        if (i0 == m - 2) {            break;        }        wdr = cos(ec * i0);        wdi = sin(ec * i0);        wkr = 0.5 * (wdr - wdi);        wki = 0.5 * (wdr + wdi);        wdr = wkr * w1r - wki * w1i;        wdi = wkr * w1i + wki * w1r;        i = i0;    }    xr = wdi * a[m - 1] - wdr * a[m + 1];    a[m - 1] = wdr * a[m - 1] + wdi * a[m + 1];    a[m + 1] = xr;    a[m] *= wki + ss * wdr;}void dstsub(int n, double *a){    int i, i0, j, k, m;    double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;        ec = M_PI_2 / n;    wkr = 0.5;    wki = 0.5;    w1r = cos(ec);    w1i = sin(ec);    wdr = 0.5 * (w1r - w1i);    wdi = 0.5 * (w1r + w1i);    ss = 2 * w1i;    m = n >> 1;    i = 0;    for (;;) {        i0 = i + 2 * DCST_LOOP_DIV;        if (i0 > m - 2) {            i0 = m - 2;        }        for (j = i + 2; j <= i0; j += 2) {            k = n - j;            xr = wdi * a[k + 1] - wdr * a[j - 1];            xi = wdr * a[k + 1] + wdi * a[j - 1];            wkr -= ss * wdi;            wki += ss * wdr;            yr = wki * a[k] - wkr * a[j];            yi = wkr * a[k] + wki * a[j];            wdr -= ss * wki;            wdi += ss * wkr;            a[j - 1] = xr;            a[j] = yr;            a[k + 1] = xi;            a[k] = yi;        }        if (i0 == m - 2) {            break;        }        wdr = cos(ec * i0);        wdi = sin(ec * i0);        wkr = 0.5 * (wdr - wdi);        wki = 0.5 * (wdr + wdi);        wdr = wkr * w1r - wki * w1i;        wdi = wkr * w1i + wki * w1r;        i = i0;    }    xr = wdi * a[m + 1] - wdr * a[m - 1];    a[m + 1] = wdr * a[m + 1] + wdi * a[m - 1];    a[m - 1] = xr;    a[m] *= wki + ss * wdr;}void dctsub4(int n, double *a){    int m;    double wki, wdr, wdi, xr;        wki = WR5000;    m = n >> 1;    if (m == 2) {        wdr = wki * WI2500;        wdi = wki * WR2500;        xr = wdi * a[1] - wdr * a[3];        a[1] = wdr * a[1] + wdi * a[3];        a[3] = xr;    }    a[m] *= wki;}void dstsub4(int n, double *a){    int m;    double wki, wdr, wdi, xr;        wki = WR5000;    m = n >> 1;    if (m == 2) {        wdr = wki * WI2500;        wdi = wki * WR2500;        xr = wdi * a[3] - wdr * a[1];        a[3] = wdr * a[3] + wdi * a[1];        a[1] = xr;    }    a[m] *= wki;}