[BZOJ 3527] [ZJOI 2014] 力

卷积 + FFT

Reference:

https://www.cnblogs.com/iwtwiioi/p/4126284.html

FFT 教学:

从多项式乘法到快速傅里叶变换

留一份 BZOJ 谜之 RTE 的代码:

Updated:AC 代码


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
96
97
98
99
100
101
102
103
104
105
106
#include<iostream>
#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' 居然就过了......也是谜......

代码:


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
96
97
98
99
100
101
102
103
104
105
106
#include<iostream>
#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);
}

发表评论

电子邮件地址不会被公开。 必填项已用*标注