forked from Dev-XYS/Algorithms
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Fast-Fourier-Transform(Iterative).cpp
95 lines (87 loc) · 1.65 KB
/
Fast-Fourier-Transform(Iterative).cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
#include <cstdio>
#include <cmath>
using namespace std;
struct complex
{
double re, im;
complex(double _re = 0, double _im = 0) : re(_re), im(_im) { }
complex operator + (const complex &x) { return complex(re + x.re, im + x.im); }
complex operator - (const complex &x) { return complex(re - x.re, im - x.im); }
complex operator * (const complex &x) { return complex(re * x.re - im * x.im, re * x.im + im * x.re); }
};
int n, m, N, k;
complex a[1 << 18], b[1 << 18], res_a[1 << 18], res_b[1 << 18];
int log2(int x)
{
int res = -1;
while (x > 0)
{
res++;
x >>= 1;
}
return res;
}
inline int reverse(int x)
{
int res = 0;
for (int i = 0; i <= k; i++)
{
if ((x & (1 << i)) != 0)
{
res |= (1 << (k - i));
}
}
return res;
}
void init(complex *a, complex *res)
{
for (int i = 0; i < N; i++)
{
res[reverse(i)] = a[i];
}
}
void dft(complex *a, complex *res, int inv)
{
init(a, res);
for (int i = 2; i <= N; i <<= 1)
{
complex w0(cos(M_PI * 2 / i), inv * sin(M_PI * 2 / i));
for (int j = 0; j < N; j += i)
{
complex w(1, 0);
for (int k = j; k < j + (i >> 1); k++)
{
complex temp = res[k];
res[k] = temp + res[k + (i >> 1)] * w;
res[k + (i >> 1)] = temp - res[k + (i >> 1)] * w;
w = w * w0;
}
}
}
}
int main()
{
scanf("%d%d", &n, &m);
for (int i = 0; i <= n; i++)
{
scanf("%lf", &a[i].re);
}
for (int i = 0; i <= m; i++)
{
scanf("%lf", &b[i].re);
}
k = log2(m + n);
N = 1 << (k + 1);
dft(a, res_a, 1);
dft(b, res_b, 1);
for (int i = 0; i < N; i++)
{
a[i] = res_a[i] * res_b[i];
}
dft(a, res_a, -1);
for (int i = 0; i <= n + m; i++)
{
printf("%d ", (int)((res_a[i].re / N) + 0.5));
}
return 0;
}