快速傅里叶变换（FFT）

Posted 2021-02-12 ac-evil

 

 1 #include <bits/stdc++.h>
 2 
 3 using namespace std;
 4 
 5 #define re register
 6 #define rep(i, a, b) for (re int i = a; i <= b; ++i)
 7 #define repd(i, a, b) for (re int i = a; i >= b; --i)
 8 #define For(i, a, b, s) for (re int i = a; i <= b; s)
 9 #define maxx(a, b) a = max(a, b)
10 #define minn(a, b) a = min(a, b)
11 #define LL long long
12 #define INF (1 << 30)
13 
14 inline int read() {
15     int w = 0, f = 1; char c = getchar();
16     while (!isdigit(c)) f = c == ‘-‘ ? -1 : f, c = getchar();
17     while (isdigit(c)) w = (w << 3) + (w << 1) + (c ^ ‘0‘), c = getchar();
18     return w * f;
19 }
20 
21 const double PI = acos(-1);
22 const int maxn = 1e6 + 5;
23 
24 struct comp {
25     double x, y;
26     comp(double x = 0, double y = 0) : x(x), y(y) {}
27 };
28 comp operator + (comp a, comp b) { return comp(a.x+b.x, a.y+b.y); }
29 comp operator - (comp a, comp b) { return comp(a.x-b.x, a.y-b.y); }
30 comp operator * (comp a, comp b) { return comp(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x); }    
31 
32 int n, m, N, lg = 0, rev[maxn << 2];
33 comp a[maxn << 2], b[maxn << 2];
34 
35 void fft(comp *a, int opt) {
36     rep(i, 0, N-1) if (i < rev[i]) swap(a[i], a[rev[i]]);
37     rep(s, 1, lg) {
38         int q = 1<<(s-1); comp Wn(cos(2*PI/(q<<1)), opt * sin(2*PI/(q<<1)));
39         for (register int p = 0; p < N; p += q<<1) {
40             comp W(1, 0);
41             rep(i, 0, q-1) {
42                 comp tmp = W * a[p+i+q];
43                 a[p+i+q] = a[p+i] - tmp;
44                 a[p+i] = a[p+i] + tmp;
45                 W = W * Wn;
46             }
47         }
48     }
49 }
50 
51 int main() {
52     n = read(), m = read();
53     N = 1; lg = 0;
54     while (N <= n+m) N <<= 1, lg++;
55     rep(i, 0, N-1) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(lg-1));
56     rep(i, 0, n) a[i] = comp(read(), 0);
57     rep(i, 0, m) b[i] = comp(read(), 0);
58     fft(a, 1); fft(b, 1);
59     rep(i, 0, N-1) a[i] = a[i] * b[i];
60     fft(a, -1);
61     rep(i, 0, n+m) printf("%d ", (int)(a[i].x / N + 0.5));
62     return 0;
63 }