卷积 + FFT
Reference:
https://www.cnblogs.com/iwtwiioi/p/4126284.html
FFT 教学：
留一份 BZOJ 谜之 RTE 的代码：
Updated：AC 代码
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
96
97
98
99
100
101
102
103
104
105
106
#include<algorithm>
#include<cstdio>
#include<iomanip>
#include<complex>
#include<cmath>
using namespace std;
inline void set_file_IO(string);
inline void close_IO(void);
inline void work(void);
int main(void) {
#ifndef ONLINE_JUDGE
set_file_IO("force");
#endif
ios::sync_with_stdio(false);
work();
#ifndef ONLINE_JUDGE
close_IO();
#endif
return 0;
}
typedef complex<double> cp;
const int nSize = 1 << 18;
cp f[nSize], g[nSize], eps[nSize], inv_eps[nSize];
inline void init_eps(const int p) {
const double Pi = M_PI;
// r = 1 -> e^(i*2*Pi*k/n) = cos(2*Pi*k/n) + i * sin(2*pi*k/n)
for (int i=0; i<p; ++i) {
eps[i] = cp(cos(2.*Pi*i/p), sin(2.*Pi*i/p));
}
for (int i=0; i<p; ++i) {
inv_eps[i] = conj(eps[i]);
}
}
inline void transform(const int n, cp *x, const cp *w) {
for (int i=0, j=0; i<n; ++i) {
if (i > j) {
swap(x[i], x[j]);
}
for (int k=n>>1; (j^=k)<k; k>>=1) ;
}
for (int i=2; i<=n; i<<=1) {
for (int j=0, m=i>>1; j<n; j+=i) {
for (int k=0; k<m; ++k) {
const cp z = x[j+m+k] * w[n/i*k];
x[j+m+k] = x[j+k] - z;
x[j+k] += z;
}
}
}
}
inline void work(void) {
int n, p = 1;
cin >> n;
while (p < n) {
p <<= 1;
}
p <<= 1;
for (int i=0; i<p; ++i) {
f[i] = g[i] = 0.;
}
for (int i=0; i<n; ++i) {
cin >> f[i];
}
for (int i=0; i<n-1; ++i) {
g[i] = 1. / (n - i - 1.) / (n - i - 1.);
g[2 * n - i - 2] = -g[i];
}
init_eps(p);
reverse(g, g+p);
transform(p, f, eps);
transform(p, g, eps);
for (int i=0; i<p; ++i) {
f[i] *= g[i];
}
transform(p, f, inv_eps);
cout << setiosflags(ios::fixed) << setprecision(3);
for (int i=p-n; i<p; ++i) {
cout << f[i].real() / p << endl;
}
}
inline void set_file_IO(string name) {
freopen((name + ".in" ).c_str(), "r", stdin );
freopen((name + ".out").c_str(), "w", stdout);
}
inline void close_IO(void) {
fclose(stdin);
fclose(stdout);
}
在不关闭流同步的情况下把 endl 换成 '\n' 居然就过了......也是谜......
代码：
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
96
97
98
99
100
101
102
103
104
105
106
#include<algorithm>
#include<cstdio>
#include<iomanip>
#include<complex>
#include<cmath>
using namespace std;
inline void set_file_IO(string);
inline void close_IO(void);
inline void work(void);
int main(void) {
#ifndef ONLINE_JUDGE
set_file_IO("force");
#endif
//ios::sync_with_stdio(false);
work();
#ifndef ONLINE_JUDGE
close_IO();
#endif
return 0;
}
typedef complex<double> cp;
const int nSize = 1 << 18;
cp f[nSize], g[nSize], eps[nSize], inv_eps[nSize];
inline void init_eps(const int p) {
const double Pi = M_PI;
// r = 1 -> e^(i*2*Pi*k/n) = cos(2*Pi*k/n) + i * sin(2*pi*k/n)
for (int i=0; i<p; ++i) {
eps[i] = cp(cos(2.*Pi*i/p), sin(2.*Pi*i/p));
}
for (int i=0; i<p; ++i) {
inv_eps[i] = conj(eps[i]);
}
}
inline void transform(const int n, cp *x, const cp *w) {
for (int i=0, j=0; i<n; ++i) {
if (i > j) {
swap(x[i], x[j]);
}
for (int k=n>>1; (j^=k)<k; k>>=1) ;
}
for (int i=2; i<=n; i<<=1) {
for (int j=0, m=i>>1; j<n; j+=i) {
for (int k=0; k<m; ++k) {
const cp z = x[j+m+k] * w[n/i*k];
x[j+m+k] = x[j+k] - z;
x[j+k] += z;
}
}
}
}
inline void work(void) {
int n, p = 1;
cin >> n;
while (p < n) {
p <<= 1;
}
p <<= 1;
for (int i=0; i<p; ++i) {
f[i] = g[i] = 0.;
}
for (int i=0; i<n; ++i) {
cin >> f[i];
}
for (int i=0; i<n-1; ++i) {
g[i] = 1. / (n - i - 1.) / (n - i - 1.);
g[2 * n - i - 2] = -g[i];
}
init_eps(p);
reverse(g, g+p);
transform(p, f, eps);
transform(p, g, eps);
for (int i=0; i<p; ++i) {
f[i] *= g[i];
}
transform(p, f, inv_eps);
cout << setiosflags(ios::fixed) << setprecision(3);
for (int i=p-n; i<p; ++i) {
cout << f[i].real() / p << '\n';
}
}
inline void set_file_IO(string name) {
freopen((name + ".in" ).c_str(), "r", stdin );
freopen((name + ".out").c_str(), "w", stdout);
}
inline void close_IO(void) {
fclose(stdin);
fclose(stdout);
}