DSP

递归FFT

2019-07-13 15:05发布

站内文章 / DSP
14814 0 1793
#include #include #include #include const int mx_n = 4e5 + 10; const double pi = acos(-1.0); #define rep(i, s, t) for(register int i = s; i <= t; ++i) using namespace std; typedef complex<double> LF; void DFT(LF *p, int n, int f) { if(n == 1) return ; n >>= 1; LF p0[n], p1[n]; rep(i, 0, n-1) p0[i] = p[i<<1], p1[i] = p[i<<1|1]; DFT(p0, n, f), DFT(p1, n, f); LF e(1, 0), w(cos(2*pi/(n<<1)), f*sin(2*pi/(n<<1))); rep(i, 0, n-1) { p[i] = p0[i] + e * p1[i]; p[i+n] = p0[i] - e * p1[i]; e = e * w; } } int n, m; LF a[mx_n], b[mx_n], c[mx_n]; int main() { #ifndef ONLINE_JUDGE freopen("input.in", "r", stdin); freopen("res.out", "w", stdout); #endif scanf("%d%d", &n, &m); rep(i, 0, n) scanf("%lf", &a[i]); rep(i, 0, m) scanf("%lf", &b[i]); int N = 1; for(; N <= n + m; N <<= 1); DFT(a, N, 1), DFT(b, N, 1); rep(i, 0, N-1) c[i] = a[i] * b[i]; DFT(c, N, -1); rep(i, 0, n + m) printf("%d%c", (int)(c[i].real()/N + 0.5), i ^ (n+m)? ' ':' '); return 0; }

热门文章