51nod扒下来的蜜汁大数乘法
来源：互联网 发布：淘宝哪家男鞋店好最潮 编辑：程序博客网 时间：2024/05/16 19:39
#include<iostream>#include<stdio.h>#include<string.h>#include<math.h>#include<algorithm>#include<iomanip>#include<stdlib.h>#include <time.h>using namespace std;typedef long long ll;typedef unsigned long long ull;typedef short lint;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 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); } }}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;}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;}const int STRSIZE = 310000;char in1[STRSIZE];char *in = in1 + 4, *tempend;void Input(){ tempend = in + STRSIZE; fread(in, 1, STRSIZE, stdin);}inline void scan(lint a[], int &len){ char c = *(in++); while (!(c >= '0' && c <= '9')) c = *(in++); in--; char *strend = in; c = *(strend++); while (c >= '0' && c <= '9') c = *(strend++); strend--; *(in - 1) = '0'; *(in - 2) = '0'; *(in - 3) = '0'; *(in - 4) = '0'; int i = 0; for (char *pos = strend - 1; pos >= in; pos -= 4) { int x = *pos - '0'; x += (*(pos - 1) - '0') * 10; x += (*(pos - 2) - '0') * 100; x += (*(pos - 3) - '0') * 1000; a[i++] = x; } in = strend; len = i;}char printnum[10000][5];void IniteOutput(){ for (int i = 0; i < 10000; i++) { int number = i; char c[5] = { '0', '0', '0', '0', '0' }; int pos = 4; while (number) { c[pos--] = number % 10 + '0'; number /= 10; } char *pn = printnum[i]; for (int j = 4; j >= 1; j--) pn[j - 1] = c[j]; }}void output(int ans[], int n){ char *out = in1; int number = ans[n - 1]; char c[5] = { '0', '0', '0', '0', '0' }; int pos = 4; while (number) { c[pos--] = number % 10 + '0'; number /= 10; } for (int j = pos + 1; j < 5; j++) *(out++) = c[j]; for (int i = n - 2; i >= 0; i--) { char *posnum = printnum[ans[i]]; while (*posnum) *(out++) = *(posnum++); } fwrite(in1, 1, out - in1, stdout);}void output(lint ans[], int n){ char *out = in1; int number = ans[n - 1]; char c[5] = { '0', '0', '0', '0', '0' }; int pos = 4; while (number) { c[pos--] = number % 10 + '0'; number /= 10; } for (int j = pos + 1; j < 5; j++) *(out++) = c[j]; for (int i = n - 2; i >= 0; i--) { char *posnum = printnum[ans[i]]; while (*posnum) *(out++) = *(posnum++); } fwrite(in1, 1, out - in1, stdout);}void initefft(const int N, double *&a, double *&w, double *&e, int *&ip, lint *&c){ const int n = N >> 1; a = (double*)malloc(sizeof(double)* (4 * n + 2)); e = (double*)malloc(sizeof(double)* 2 * n); w = (double*)malloc(sizeof(double)* n); c = (lint*)malloc(sizeof(lint)* 2 * n); ip = (int*)malloc(sizeof(int)* 8192); ip[0] = 0; makewt(n, ip, w); double pin2 = 2.0*acos(-1.0) / N; int N1 = N >> 1; for (int i = 0; i <= N1; i += 4) { e[i] = cos(pin2*(i >> 1)); e[i + 1] = sin(pin2*(i >> 1)); e[i + 2] = cos(pin2*((i >> 1) + 1)); e[i + 3] = sin(pin2*((i >> 1) + 1)); } for (int i = N1 + 2; i < N; i += 4) { e[i] = -e[N - i]; e[i + 1] = e[N - i + 1]; e[i + 2] = -e[N - i - 2]; e[i + 3] = e[N - i - 1]; }}void mul(lint a[], lint b[], int c[], int lengtha, int lengthb, int &lengthc){ //c[]清零 lengthc = lengtha + lengthb - 1; c[lengthc] = 0; for (int i = 0; i < lengtha; i++) for (int j = 0; j < lengthb; j++) { c[i + j] += a[i] * b[j]; if (c[i + j] >= 10000) { c[i + j + 1] += c[i + j] / 10000; c[i + j] %= 10000; } } if (c[lengthc] != 0) lengthc++;}void mulfft(int N, double *a, double *w, double *e, int *ip, lint *num1, int &lengtha, lint *num2, int &lengthb, lint *c, int &lengthc){ int n = 1; while (n < (lengtha + lengthb)) n <<= 1; n >>= 1; int n1 = n / 2, n2 = 2 * n, n4 = 4 * n; if (lengtha <= lengthb) { for (int i1 = 0, i = 0; i1 < lengtha; i1 += 2, i += 4) { a[i] = num1[i1]; a[i + 1] = num2[i1]; a[i + 2] = num1[i1 + 1]; a[i + 3] = num2[i1 + 1]; } for (int i1 = lengtha, i = (lengtha << 1); i1 < lengthb; i1 += 2, i += 4) { a[i] = 0; a[i + 1] = num2[i1]; a[i + 2] = 0; a[i + 3] = num2[i1 + 1]; } for (int i = (lengthb << 1); i < n4; i += 4) { a[i] = 0; a[i + 1] = 0; a[i + 2] = 0; a[i + 3] = 0; } } else { for (int i1 = 0, i = 0; i1 < lengthb; i1 += 2, i += 4) { a[i] = num1[i1]; a[i + 1] = num2[i1]; a[i + 2] = num1[i1 + 1]; a[i + 3] = num2[i1 + 1]; } for (int i1 = lengthb, i = (lengthb << 1); i1 < lengtha; i1 += 2, i += 4) { a[i] = num1[i1]; a[i + 1] = 0; a[i + 2] = num1[i1 + 1]; a[i + 3] = 0; } for (int i = (lengtha << 1); i < n4; i += 4) { a[i] = 0; a[i + 1] = 0; a[i + 2] = 0; a[i + 3] = 0; } } cdft(4 * n, 1, a, ip, w); double rwi, iwi, rwni, iwni; a[n4] = a[0]; a[n4 + 1] = a[1]; for (int i = 0; i <= n2; i += 2) { rwi = a[i] * a[i] - a[i + 1] * a[i + 1]; iwi = a[i] * a[i + 1] * 2.0; rwni = a[n4 - i] * a[n4 - i] - a[n4 - i + 1] * a[n4 - i + 1]; iwni = a[n4 - i] * a[n4 - i + 1] * 2.0; a[i] = (iwi + iwni)*0.25; a[i + 1] = (rwni - rwi)*0.25; a[n4 - i] = (iwi + iwni)*0.25; a[n4 - i + 1] = (rwi - rwni)*0.25; } int q = N / n; double rw1, iw1, rw2, iw2; for (int i = 0; i < n2; i += 2) { rw1 = a[i] - a[i + n2]; iw1 = a[i + 1] - a[i + n2 + 1]; rw2 = rw1* *e + iw1* *(e + 1); iw2 = iw1 * *e - rw1* *(e + 1); rw1 = a[i] + a[i + n2]; iw1 = a[i + 1] + a[i + n2 + 1]; a[i] = (rw1 - iw2)*0.5; a[i + 1] = (iw1 + rw2)*0.5; e += q; } cdft(2 * n, -1, a, ip, w); lengthc = lengtha + lengthb - 1; double dn1 = 1.0 / n, bit2 = 1.0 / 100000000, bit1 = 1.0 / 10000, temp1; int endn = min(lengthc + 2, n2), temp2; c[lengthc] = 0; a[0] *= dn1; a[0] += 0.5; a[1] *= dn1; a[1] += 0.5; for (int i = 0; i < endn; i++) { temp1 = (int)(a[i] * bit2); a[i + 2] *= dn1; a[i + 2] += 0.5; a[i + 2] += temp1; a[i] -= temp1 * 100000000.0; temp2 = a[i] * bit1; a[i + 1] += temp2; a[i] -= temp2 * 10000; c[i] = a[i]; } if (c[lengthc] != 0) lengthc++;}const static int SIZEN = 131072;lint a[SIZEN], b[SIZEN];int c[SIZEN];int main(){ //freopen("test.in", "r", stdin); //freopen("test.out", "w", stdout); Input(); int lena = 0, lenb = 0; scan(a, lena); scan(b, lenb); int len; if (lena <= 0 || lenb <= 0) { if (lenb < lena) mul(b, a, c, lenb, lena, len); else mul(a, b, c, lena, lenb, len); IniteOutput(); output(c, len); } else { double *a1, *w, *e; int *ip; lint *c; int N = 1; while (N < (lena + lenb)) N <<= 1; initefft(N, a1, w, e, ip, c); mulfft(N, a1, w, e, ip, a, lena, b, lenb, a, len); IniteOutput(); output(a, len); } return 0;}阅读全文
1 0
- 51nod扒下来的蜜汁大数乘法
- 51Nod 大数乘法
- 51nod 大数乘法
- 51nod -1027 . 大数乘法
- 51nod 1027 大数乘法
- 51nod 1027 大数乘法
- 51nod 1027 大数乘法
- 51nod 1027 大数乘法
- 【51nod】1027 大数乘法
- 51Nod 1027 大数乘法
- 51Nod-1027-大数乘法
- 51Nod 1027 大数乘法
- 51nod 1027 大数乘法
- 51Nod 1027 大数乘法
- 51Nod－1027－大数乘法
- 51Nod-1028-大数乘法 V2
- 【51nod 1028 大数乘法 V2】
- 51nod 1028 大数乘法 V2
- CSS选择器
- 分治策略之最大子数组问题
- JS浮点计算问题
- fanxing
- IDEA使用:intellIJ idea 配置 svn
- 51nod扒下来的蜜汁大数乘法
- HDU2222-Keywords Search
- jQuery手风琴图切换特效插件
- Android7.0上PopupWindow的showAsDropDown位置问题 参数无效
- 廖雪峰《python3 基础教程》读书笔记——第四章
- Scala解析Json格式
- do-while语句
- 页面中嵌入帆软报表，帆软报表服务器和项目不是一个，出现跨域的问题，解决方案
- vs2015无法解析外部符号__imp__fprintf 以及imp_iob_func错误
