多项式学习笔记

发表于 更新于 分类于 阅读次数： 本文字数： 68k
苦练计数。

多项式乘法

这里没有很细节的了解原理，几乎就是当黑盒来用的。

FFT

例题是高精度乘法，需要处理进位。

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
#include <bits/stdc++.h>

#define int long long
using namespace std;

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;

using u128 = unsigned __int128;
using i128 = __int128;
using ld = long double;

const double PI = acos(-1.0);

struct Complex {
    double x, y;
    Complex(double x_ = 0, double y_ = 0) {
        x = x_;
        y = y_;
    }
    Complex operator+(const Complex &b) const {
        return Complex(x + b.x, y + b.y);
    }
    Complex operator-(const Complex &b) const {
        return Complex(x - b.x, y - b.y);
    }
    Complex operator*(const Complex &b) const {
        return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
    }
};

void change(Complex y[], int len) {
    for (int i = 1, j = len / 2, k; i < len - 1; i++) {
        if (i < j) {
            std::swap(y[i], y[j]);
        }
        k = len / 2;
        while (j >= k) {
            j -= k;
            k /= 2;
        }
        if (j < k) {
            j += k;
        }
    }
}

void fft(Complex y[], int len, int on) {
    change(y, len);
    for (int h = 2; h <= len; h *= 2) {
        Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
        for (int j = 0; j < len; j += h) {
            Complex w(1, 0);
            for (int k = j; k < j + h / 2; k++) {
                Complex u = y[k];
                Complex t = w * y[k + h / 2];
                y[k] = u + t;
                y[k + h / 2] = u - t;
                w = w * wn;
            }
        }
    }
    if (on == -1) {
        for (int i = 0; i < len; i++) {
            y[i].x /= len;
        }
    }
}

// 注意数组要开到 2^x > n + m
constexpr int N = 2100010;
Complex a[N], b[N], ans[N];

signed main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int n, m;
    cin >> n >> m;
    for (int i = 0; i <= n; i++) {
        cin >> a[i].x;
    }
    for (int i = 0; i <= m; i++) {
        cin >> b[i].x;
    }
    int len = 1;
    while (len <= n + m) {
        len *= 2;
    }
    fft(a, len, 1);
    fft(b, len, 1);
    for (int i = 0; i < len; i++) {
        ans[i] = a[i] * b[i];
    }
    fft(ans, len, -1);
    for (int i = 0; i <= n + m; i++) {
        cout << (int)(ans[i].x + 0.5) << " \n"[i == n + m];
    }

    return 0;
}

固定模数NTT

常见的 NTT 模数有 998244353，1004535809， 469762049。

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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
#include <bits/stdc++.h>

#define int long long
using namespace std;

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;

using u128 = unsigned __int128;
using i128 = __int128;
using ld = long double;

constexpr int P = 998244353;
constexpr int g = 3;

int power(int a, int b) {
    int res = 1;
    for (; b; b /= 2, a = 1LL * a * a % P) {
        if (b % 2) {
            res = 1LL * res * a % P;
        }
    }
    return res;
}

std::vector<int> rev, roots{0, 1};

void dft(std::vector<int> &a) {
    int n = a.size();
    if ((int)(rev.size()) != n) {
        int k = __builtin_ctz(n) - 1;
        rev.resize(n);
        for (int i = 0; i < n; i++) {
            rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
        }
    }
    for (int i = 0; i < n; i++) {
        if (rev[i] < i) {
            std::swap(a[i], a[rev[i]]);
        }
    }
    if (roots.size() < n) {
        int k = __builtin_ctz(roots.size());
        roots.resize(n);
        while ((1 << k) < n) {
            int e = power(g, (P - 1) >> (k + 1));
            for (int i = 1 << (k - 1); i < (1 << k); i++) {
                roots[2 * i] = roots[i];
                roots[2 * i + 1] = 1LL * roots[i] * e % P;
            }
            k++;
        }
    }

    for (int k = 1; k < n; k *= 2) {
        for (int i = 0; i < n; i += 2 * k) {
            for (int j = 0; j < k; j++) {
                int u = a[i + j];
                int v = 1LL * a[i + j + k] * roots[k + j] % P;
                a[i + j] = (u + v) % P;
                a[i + j + k] = (u - v + P) % P;
            }
        }
    }
}

void idft(std::vector<int> &a) {
    int n = a.size();
    std::reverse(a.begin() + 1, a.end());
    dft(a);
    int inv = power(n, P - 2);
    for (int i = 0; i < n; i++) {
        a[i] = 1LL * a[i] * inv % P;
    }
}

std::vector<int> mul(std::vector<int> a, std::vector<int> b) {
    int n = 1, tot = a.size() + b.size() - 1;
    while (n < tot) {
        n *= 2;
    }
    if (tot < 128) {
        std::vector<int> c(a.size() + b.size() - 1);
        for (int i = 0; i < a.size(); i++) {
            for (int j = 0; j < b.size(); j++) {
                c[i + j] = (c[i + j] + 1LL * a[i] * b[j]) % P;
            }
        }
        return c;
    }
    a.resize(n);
    b.resize(n);
    dft(a);
    dft(b);
    for (int i = 0; i < n; i++) {
        a[i] = 1LL * a[i] * b[i] % P;
    }
    idft(a);
    a.resize(tot);
    return a;
}

signed main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int n, m;
    cin >> n >> m;
    vector<int> a(n + 1), b(m + 1);
    for (int i = 0; i <= n; i++) {
        cin >> a[i];
    }
    for (int i = 0; i <= m; i++) {
        cin >> b[i];
    }
    auto c = mul(a, b);
    for (int i = 0; i <= n + m; i++) {
        cout << c[i] << " \n"[i == n + m];
    }

    return 0;
}

任意模数 NTT（MTT）

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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
using i128 = __int128_t;

constexpr i64 mod1 = 998244353;
constexpr i64 mod2 = 1004535809;
constexpr i64 mod3 = 469762049;
constexpr i64 G = 3;

i64 power(i64 a, i64 e, i64 mod) {
    i64 r = 1 % mod;
    a %= mod;
    while (e) {
        if (e & 1) r = (i128)r * a % mod;
        a = (i128)a * a % mod;
        e >>= 1;
    }
    return r;
}

void bit_reverse(vector<i64> &a) {
    int n = (int)a.size();
    int j = 0;
    for (int i = 1; i < n; ++i) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1) j ^= bit;
        j ^= bit;
        if (i < j) swap(a[i], a[j]);
    }
}

void ntt(vector<i64> &a, int sign, i64 mod) {
    int n = (int)a.size();
    bit_reverse(a);
    for (int len = 2; len <= n; len <<= 1) {
        i64 wn = power(G, (mod - 1) / len, mod);
        if (sign == -1) wn = power(wn, mod - 2, mod);
        for (int i = 0; i < n; i += len) {
            i64 w = 1;
            int half = len >> 1;
            for (int j = 0; j < half; ++j) {
                i64 u = a[i + j];
                i64 v = (i128)w * a[i + j + half] % mod;
                a[i + j] = u + v;
                if (a[i + j] >= mod) a[i + j] -= mod;
                a[i + j + half] = u - v;
                if (a[i + j + half] < 0) a[i + j + half] += mod;
                w = (i128)w * wn % mod;
            }
        }
    }
    if (sign == -1) {
        i64 inv = power(n, mod - 2, mod);
        for (int i = 0; i < n; ++i) a[i] = (i128)a[i] * inv % mod;
    }
}

vector<i64> multiply_mod(vector<i64> a, vector<i64> b, i64 mod) {
    if (a.empty() || b.empty()) return {};
    int tot = (int)a.size() + (int)b.size() - 1;
    int sz = 1;
    while (sz < tot) sz <<= 1;
    a.resize(sz);
    b.resize(sz);
    ntt(a, 1, mod);
    ntt(b, 1, mod);
    for (int i = 0; i < sz; ++i) a[i] = (i128)a[i] * b[i] % mod;
    ntt(a, -1, mod);
    a.resize(tot);
    return a;
}

vector<i64> multiply_crt(vector<i64> a, vector<i64> b, i64 p) {
    if (a.empty() || b.empty()) return {};
    auto c1 = multiply_mod(a, b, mod1);
    auto c2 = multiply_mod(a, b, mod2);
    auto c3 = multiply_mod(a, b, mod3);
    int n = (int)c1.size();

    i128 M1 = (i128)mod1;
    i128 M2 = (i128)mod2;
    i128 M3 = (i128)mod3;
    i128 M12 = M1 * M2;
    i64 mod12_ll = (i64)(M12);

    i64 inv1_mod2 = power(mod1 % mod2, mod2 - 2, mod2);
    i64 inv12_mod3 = power((i64)((M12) % mod3), mod3 - 2, mod3);

    vector<i64> res(n);
    for (int i = 0; i < n; ++i) {
        i64 r1 = c1[i];
        i64 r2 = c2[i];
        i64 r3 = c3[i];

        i64 diff = (r2 - r1) % mod2;
        if (diff < 0) diff += mod2;
        i64 k1 = (i128)diff * inv1_mod2 % mod2;
        i128 t12 = (i128)r1 + (i128)k1 * (i128)mod1;

        i64 t12_mod3 = (i64)(t12 % mod3);
        i64 diff3 = (r3 - t12_mod3) % mod3;
        if (diff3 < 0) diff3 += mod3;
        i64 k2 = (i128)diff3 * inv12_mod3 % mod3;

        i128 result128 = t12 + (i128)k2 * M12;

        i64 finalv = (i64)(result128 % p);
        if (finalv < 0) finalv += p;
        res[i] = finalv;
    }
    return res;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int n, m;
    long long p;
    if (!(cin >> n >> m >> p)) return 0;
    vector<i64> A(n + 1), B(m + 1);
    for (int i = 0; i <= n; ++i) cin >> A[i];
    for (int i = 0; i <= m; ++i) cin >> B[i];

    auto C = multiply_crt(A, B, p);
    for (int i = 0; i <= n + m; ++i) {
        cout << C[i] << (i == n + m ? '\n' : ' ');
    }
    return 0;
}

P3338 [ZJOI2014] 力

题意

给定实数 \(q_i\)，定义 \(F_j = \sum_{i=1}^{j - 1}\dfrac{q_i\times q_j}{(i-j)^2} - \sum_{i=j+1}^{n}\dfrac{q_i\times q_j}{(i-j)^2}\)，对于 \(1\le i\le n\)， 求 \(E_i = \dfrac{F_i}{q_i}\)。

题解

学多项式首先要知道多项式拿来干嘛的...上面的多项式乘法学了还是不知道，通过这题可以了解多项式第一个应用——优化卷积运算，通过多项式乘法，我们可以快速求出结果。 \[ E_j =\sum_{i=1}^{j}\dfrac{ q_i}{(i-j)^2} - \sum_{i=j}^{n}\dfrac{q_i}{(i-j)^2}\\\\ \] 我们定义 \(a_i = q_i, b_i = \dfrac{1}{i^2}\)，那么可以写成如下形式： \[ E_j = \sum_{i=1}^ja_ib_{j-i} - \sum_{i=j}^na_ib_{i-j} \] 前者已经是卷积的形式了，继续化简后者，一个常见的做法是定义 \(c_i = a_{n-i}\) \[ \sum_{i=j}^na_ib_{i-j} &= \sum_{i=j}^nc_{n-i}b_{i-j}\\\\ &=\sum_{i=0}^{n-j}c_{n-i-j}b_i\\\\ &=\sum_{i=0}^tc_{t-i}b_i \] 我们设 \(F(x) = \sum_{i=1}^ja_ib_{j-i},G(x) = \sum_{i=0}^tc_{t-i}b_i\)，我们设 \(l_i,r_i\) 分别是 \(F(x),G(x)\) 中 \(x^i\) 项的系数，如果我们直接将多项式 \(a(x)\) 和 \(b(x)\)相乘，第 \(i\) 项的系数就是序列卷积的结果。那么 \(E_i = l_i - r_{n - i}\)​。

时间复杂度：\(O(n\log 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
#include <bits/stdc++.h>

#define int long long
using namespace std;

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;

using u128 = unsigned __int128;
using i128 = __int128;
using ld = long double;

const double PI = acos(-1.0);

struct Complex {
    double x, y;
    Complex(double x_ = 0, double y_ = 0) {
        x = x_;
        y = y_;
    }
    Complex operator+(const Complex &b) const {
        return Complex(x + b.x, y + b.y);
    }
    Complex operator-(const Complex &b) const {
        return Complex(x - b.x, y - b.y);
    }
    Complex operator*(const Complex &b) const {
        return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
    }
};

void change(Complex y[], int len) {
    for (int i = 1, j = len / 2, k; i < len - 1; i++) {
        if (i < j) {
            std::swap(y[i], y[j]);
        }
        k = len / 2;
        while (j >= k) {
            j -= k;
            k /= 2;
        }
        if (j < k) {
            j += k;
        }
    }
}

void fft(Complex y[], int len, int on) {
    change(y, len);
    for (int h = 2; h <= len; h *= 2) {
        Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
        for (int j = 0; j < len; j += h) {
            Complex w(1, 0);
            for (int k = j; k < j + h / 2; k++) {
                Complex u = y[k];
                Complex t = w * y[k + h / 2];
                y[k] = u + t;
                y[k + h / 2] = u - t;
                w = w * wn;
            }
        }
    }
    if (on == -1) {
        for (int i = 0; i < len; i++) {
            y[i].x /= len;
        }
    }
}

// 注意数组要开到 2^x > n + m
constexpr int N = 270000;
Complex a[N], b[N], c[N], F[N], G[N];

signed main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int n;
    cin >> n;
    for (int i = 1; i <= n; i++) {
        cin >> a[i].x;
        c[n - i].x = a[i].x;
        b[i] = 1.0l / i / i;
    }
    int len = 1;
    while (len <= n + n) {
        len *= 2;
    }
    fft(a, len, 1);
    fft(b, len, 1);
    fft(c, len, 1);
    for (int i = 0; i < len; i++) {
        F[i] = a[i] * b[i];
        G[i] = c[i] * b[i];
    }
    fft(F, len, -1);
    fft(G, len, -1);
    for (int i = 1; i <= n; i++) {
        cout << fixed << setprecision(10) << F[i].x - G[n - i].x << "\n";
    }

    return 0;
}

P3723 [AHOI2017/HNOI2017] 礼物

题意

两个环各有 \(n\) 个珠子，每个珠子写有一个数字 \(a_i,b_i\)，选择 \([-m,m]\) 中的一个数 \(x\)，给一个环的所有数都添加，求 \(\sum_{i=1}^n(a_i - b_i)^2\) 的最小值，环可以旋转。

\(n\le 50000,m\le 100\)。

题解

\[ \sum_{i=1}^n(a_i-b_i + x)^2 = \sum_{i=1}^na_i^2 + \sum_{i=1}^nb_i^2 +nx^2 + 2x(\sum_{i=1}^na_i-\sum_{i=1}^nb_i) - 2\sum_{i=1}^na_ib_i \]

注意我们枚举 \(x\) 时，只有 \(\sum_{i=1}^na_ib_i\) 的值不确定，我们需要求它的最小值。

这是一个经典的卷积方法，我们将 \(a\) 逆序，即 \(a_i \leftarrow a_{n-i+1}\)，就变成了 \(\sum_{i=1}^na_{n-i+1}b_i\)，这个式子可以是卷积形式。对环的处理也非常经典，将 \(b\) 复制一次后，\(\sum_{i=1}^na_ib_i\) 的最小值就是 \(x^{n+1},x^{n+2},...,x^{2n}\) 对应的系数的最小值。

时间复杂度：\(O(n\log n + m)\)

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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
#include <bits/stdc++.h>

#define int long long
using namespace std;

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;

using u128 = unsigned __int128;
using i128 = __int128;
using ld = long double;

constexpr int P = 998244353;

int power(int a, int b) {
    int res = 1;
    for (; b; b /= 2, a = 1LL * a * a % P) {
        if (b % 2) {
            res = 1LL * res * a % P;
        }
    }
    return res;
}

std::vector<int> rev, roots{0, 1};

void dft(std::vector<int> &a) {
    int n = a.size();
    if ((int)(rev.size()) != n) {
        int k = __builtin_ctz(n) - 1;
        rev.resize(n);
        for (int i = 0; i < n; i++) {
            rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
        }
    }
    for (int i = 0; i < n; i++) {
        if (rev[i] < i) {
            std::swap(a[i], a[rev[i]]);
        }
    }
    if (roots.size() < n) {
        int k = __builtin_ctz(roots.size());
        roots.resize(n);
        while ((1 << k) < n) {
            int e = power(31, 1 << (__builtin_ctz(P - 1) - k - 1));
            for (int i = 1 << (k - 1); i < (1 << k); i++) {
                roots[2 * i] = roots[i];
                roots[2 * i + 1] = 1LL * roots[i] * e % P;
            }
            k++;
        }
    }

    for (int k = 1; k < n; k *= 2) {
        for (int i = 0; i < n; i += 2 * k) {
            for (int j = 0; j < k; j++) {
                int u = a[i + j];
                int v = 1LL * a[i + j + k] * roots[k + j] % P;
                a[i + j] = (u + v) % P;
                a[i + j + k] = (u - v + P) % P;
            }
        }
    }
}

void idft(std::vector<int> &a) {
    int n = a.size();
    std::reverse(a.begin() + 1, a.end());
    dft(a);
    int inv = power(n, P - 2);
    for (int i = 0; i < n; i++) {
        a[i] = 1LL * a[i] * inv % P;
    }
}

std::vector<int> mul(std::vector<int> a, std::vector<int> b) {
    int n = 1, tot = a.size() + b.size() - 1;
    while (n < tot) {
        n *= 2;
    }
    if (tot < 128) {
        std::vector<int> c(a.size() + b.size() - 1);
        for (int i = 0; i < a.size(); i++) {
            for (int j = 0; j < b.size(); j++) {
                c[i + j] = (c[i + j] + 1LL * a[i] * b[j]) % P;
            }
        }
        return c;
    }
    a.resize(n);
    b.resize(n);
    dft(a);
    dft(b);
    for (int i = 0; i < n; i++) {
        a[i] = 1LL * a[i] * b[i] % P;
    }
    idft(a);
    a.resize(tot);
    return a;
}

signed main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int n, m;
    cin >> n >> m;
    int suma1 = 0, suma2 = 0, sumb1 = 0, sumb2 = 0;
    vector<int> A(n + 1), a(n + 1), b(2 * n + 1);
    for (int i = 1; i <= n; i++) {
        cin >> A[i];
        suma1 += A[i];
        suma2 += A[i] * A[i];
    }
    for (int i = 1; i <= n; i++) {
        cin >> b[i];
        b[i + n] = b[i];
        sumb1 += b[i];
        sumb2 += b[i] * b[i];
        a[i] = A[n - i + 1];
    }
    auto c = mul(a, b);
    int mx = -2e18;
    for (int i = n + 1; i <= 2 * n; i++) {
        mx = max(mx, c[i]);
    }
    int ans = 2e18;
    for (int x = -m; x <= m; x++) {
        ans = min(ans,
                  suma2 + sumb2 + n * x * x + 2 * x * (suma1 - sumb1) - 2 * mx);
    }
    cout << ans << "\n";

    return 0;
}

P5641 【CSGRound2】开拓者的卓识

题意

记一个子段 \([l,r]\) 的 \(k\) 阶子段和为 \(sum_{k,l,r}\)，有 \(sum_{k,l,r} = \begin{cases}\sum_{i=l}^ra_i,&k=1\\\sum_{i=l}^r\sum_{j=i}^rsum_{k-1,i,j},&k\ge 2\end{cases}\) 对所有 \(i\in [1,n]\)，求 \(sum_{k,1,i}\) 模 998244353。

\(n\le 10^5,k\le 998244353\)

题解

我们不妨考虑单个 \(a_i\) 对答案的贡献：

当 \(k = 1\) 时，只要 \(l\le i\le r\) 就能贡献 1，否则是 0。

当 \(k = 2\)，需要在 \([l,r]\) 选一个子区间 \([l^\prime,r^\prime]\)，满足 \(l\le l^\prime\le i\le r^\prime\le r\)。

当 \(k = 3\)，需要在 \([l,r]\) 里选两层，即 \(l\le l^\prime\le l^{\prime\prime}\le i\le {r^{\prime\prime}}\le r^\prime\le r\)。

考虑枚举 \(r\) 求答案。对于固定的 \(r\)，等价于在 \(i\) 左侧选择 \(k-1\) 个位置，在 \(i\) 右侧选择 \(k - 1\) 个位置，二者相乘就是答案，那么可以推导出 \(sum_{k,1,r} = \sum_{i=1}^na_i\dbinom{i + k - 2}{k-1}\dbinom{r + k - i - 1}{k - 1}\)。因为 \((i + k - 2) + (r + k - i - 1) = r + 2k - 3\)，为常量，可以考虑卷积。

设 \(f_i = a_i\dbinom{i + k - 2}{k-1},g_i = \dbinom{i + k - 2}{k-1}\)，容易发现 \(sum_{k,1,r}\) 就是 \(f_ig_{r-i+1}\) 对应的 \(x^{r+1}\) 项的系数，那么 NTT 即可。

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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#include <bits/stdc++.h>

#define int long long
using namespace std;

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;

using u128 = unsigned __int128;
using i128 = __int128;
using ld = long double;

constexpr int P = 998244353;

int power(int a, int b) {
    int res = 1;
    for (; b; b /= 2, a = 1LL * a * a % P) {
        if (b % 2) {
            res = 1LL * res * a % P;
        }
    }
    return res;
}

std::vector<int> rev, roots{0, 1};

void dft(std::vector<int> &a) {
    int n = a.size();
    if ((int)(rev.size()) != n) {
        int k = __builtin_ctz(n) - 1;
        rev.resize(n);
        for (int i = 0; i < n; i++) {
            rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
        }
    }
    for (int i = 0; i < n; i++) {
        if (rev[i] < i) {
            std::swap(a[i], a[rev[i]]);
        }
    }
    if (roots.size() < n) {
        int k = __builtin_ctz(roots.size());
        roots.resize(n);
        while ((1 << k) < n) {
            int e = power(31, 1 << (__builtin_ctz(P - 1) - k - 1));
            for (int i = 1 << (k - 1); i < (1 << k); i++) {
                roots[2 * i] = roots[i];
                roots[2 * i + 1] = 1LL * roots[i] * e % P;
            }
            k++;
        }
    }

    for (int k = 1; k < n; k *= 2) {
        for (int i = 0; i < n; i += 2 * k) {
            for (int j = 0; j < k; j++) {
                int u = a[i + j];
                int v = 1LL * a[i + j + k] * roots[k + j] % P;
                a[i + j] = (u + v) % P;
                a[i + j + k] = (u - v + P) % P;
            }
        }
    }
}

void idft(std::vector<int> &a) {
    int n = a.size();
    std::reverse(a.begin() + 1, a.end());
    dft(a);
    int inv = power(n, P - 2);
    for (int i = 0; i < n; i++) {
        a[i] = 1LL * a[i] * inv % P;
    }
}

std::vector<int> mul(std::vector<int> a, std::vector<int> b) {
    int n = 1, tot = a.size() + b.size() - 1;
    while (n < tot) {
        n *= 2;
    }
    if (tot < 128) {
        std::vector<int> c(a.size() + b.size() - 1);
        for (int i = 0; i < a.size(); i++) {
            for (int j = 0; j < b.size(); j++) {
                c[i + j] = (c[i + j] + 1LL * a[i] * b[j]) % P;
            }
        }
        return c;
    }
    a.resize(n);
    b.resize(n);
    dft(a);
    dft(b);
    for (int i = 0; i < n; i++) {
        a[i] = 1LL * a[i] * b[i] % P;
    }
    idft(a);
    a.resize(tot);
    return a;
}

signed main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int n, k;
    cin >> n >> k;
    vector<int> a(n + 1), f(n + 1), g(n + 1);
    for (int i = 1; i <= n; i++) {
        cin >> a[i];
    }
    vector<int> inv(n + 1);
    inv[1] = 1;
    for (int i = 2; i <= n; i++) {
        inv[i] = (P - P / i) * inv[P % i] % P;
    }
    f[1] = g[1] = 1;
    for (int i = 2; i <= n; i++) {
        f[i] = f[i - 1] * (i + k - 2) % P * inv[i - 1] % P;
        g[i] = f[i];
    }
    for (int i = 1; i <= n; i++) {
        f[i] = f[i] * a[i] % P;
    }
    auto ans = mul(f, g);
    for (int i = 2; i <= n + 1; i++) {
        cout << ans[i] << " \n"[i == n + 1];
    }

    return 0;
}

普通生成函数

之前看着生成函数这四个字就有一种害怕感，现在认真来学总算是慢慢会了一些。

普通生成函数（OGF）跟级数有很大相似之处，定义为 \(F(x) = \sum_{n}a_nx^n\)，可以是有限项，也可以是无穷项。我们其实并不关心 \(x\) 是什么，而只关注 \(a_n\)，这里一般是以 0 为起点的。如果 \(a_n\)​ 有通项公式，那么它的普通生成函数的系数就是通项公式。

一些例子：

  • 序列 \(\{1,2,3\}\) 的普通生成函数是 \(1+2x^2+3x^3\)
  • 序列 \(\{1,1,1,...\}\) 的普通生成函数是 \(\sum_{n\ge 0}x^n\)
  • 序列 \(\{1,2,4,8,16,...\}\) 的普通生成函数是 \(\sum_{n\ge 0}2^nx^n\)
  • 序列 \(\{1,3,5,7,9,...\}\) 的普通生成函数是 \(\sum_{n\ge 0}(2n+1)x^n\)

基本运算

对于两个序列 \(a,b\) 的普通生成函数 \(F(x), G(x)\)，那么 \(F(x) \pm G(x) = \sum_n(a_n\pm b_n)x^n\)。换句话说，\(F(x) \pm G(x)\) 是序列 \(a_n\pm b_n\) 的生成函数。

对于乘法运算，也就是卷积，\(F(x)G(x) = \sum_nx^n\sum_{i=0}^na_ib_{n-i}\)，因此 \(F(x)G(x)\) 是序列 \(\sum_{i=0}^na_ib_{n-i}\)​ 的普通生成函数。因此，我们想求卷积，如果能够求出它们的系数，可以转化为多项式乘法。

封闭形式

我们并不是总是使用求和形式的生成函数，有时又会使用封闭形式，这样更方便化简。这里的化简需要我们利用一些定理或公式，比如二项式定理、等比数列求和公式、求导和积分等。更常见的是，对于无穷项的生成函数，我们得到一个递推公式来解方程。

注意这里牛顿二项式定理的 \(r\) 可以是在复数域下，\(\begin{pmatrix}r\\k\end{pmatrix} = \dfrac{r^{\underline{k}}}{k!}\)，\((1+x)^\alpha = \sum_{n\ge 0}\begin{pmatrix}\alpha\\n\end{pmatrix}x^n\)。

例如，求 \(\{1,1,1,...\}\) 的普通生成函数 \(F(x) = \sum_{n\ge 0}x^n\)。注意到 \(xF(x)\) 就是把每一项乘上 \(x\)，相当于向右移动一位，这样就相对于 \(F(x)\) 差了 \(x^0\) 所在的一项，那么 \(xF(x) + 1 = F(x)\)，我们可以解出 \(F(x) = \dfrac{1}{1-x}\)。

同理，对于等比数列 \(\{1,p,p^2,p^3,...\}\)，生成函数 \(F(x) = \sum_{n\ge 0}p^nx^n\)，\(xpF(x) + 1 = F(x)\)，那么 \(F(x) = \dfrac{1}{1-px}\)。

P5488 差分与前缀和

题意

给定一个长为 \(n\) 的序列 \(a\)，求出其 \(k\) 阶前缀和或差分，对 1004535809 取模。

\(n\le 10^5,a_i\le 10^9,k\le 10^{2333},k\not\equiv 0\pmod {1004535809}\)

题解

考虑序列 \(\{a_0,a_1,...\}\)，它对应的 OGF 是 \(F(x)\)。我们设它的前缀和序列为 \(\{b_0,b_1,...\}\)，我们知道 \(b_i = \sum_{j=0}^ia_j = \sum_{j+k=i}a_jc_{i-j}\)，其中 \(c_i = 1\)，那么我们可以发现，求一次前缀和就是与一个全 1 序列做卷积，而 \(\sum_{i\ge0}x^i = \dfrac{1}{1-x}\)，也就是说，求一次前缀和，就是对 \(F(x)\) 乘一个 \(\dfrac{1}{1-x}\)。

而差分满足 \(b_i = a_i - a_{i-1}\)，写成生成函数就是 \(\sum_{i\ge 0}b_ix^i = \sum_{i\ge 0}(a_i - a_{i-1})x^i = F(x) - xF(x) = (1-x)F(x)\)，那么求一次差分就是对 \(F(x)\) 乘一个 \(1-x\)。

因此前缀和和差分的答案分别是 \(\dfrac{F(x)}{(1-x)^k}\) 和 \(F(x)(1-x)^k\)，可以通过多项式求逆，\(k\) 次幂，多项式乘法计算，时间复杂度 \(O(n(\log n + \log k))\)，不过常数巨大，会 TLE。

考虑换种方式求系数，先考虑差分。我们设 \(G(x) = (1-x)^k = \sum_{i\ge 0}(-1)^i\begin{pmatrix}k\\i\end{pmatrix}x^i\)，虽然 \(k\) 非常大，但是有递推关系 \(\begin{pmatrix}k\\0\end{pmatrix} = 1,\begin{pmatrix}k\\i\end{pmatrix} = \begin{pmatrix}k\\i-1\end{pmatrix}\times \dfrac{k-i+1}{i}\)，可以求出 \(G(x)\) 每一项的系数，再与 \(F(x)\) 进行多项式乘法即可知道卷积的结果。

再考虑前缀和，广义二项式定理有 $(1-x)^{-k} = _{i}(-1)^i \[\begin{pmatrix}-k\\i\end{pmatrix}\] x^i \(，注意组合数学有个重要的性质（上指标反转）：\) \[\begin{pmatrix}r\\k\end{pmatrix}\] = (-1)^k \[\begin{pmatrix}k-r-1\\k\end{pmatrix}\]

$，那么 \((1-x)^{-k} = \sum_{i\ge 0}\begin{pmatrix}i+k-1\\i\end{pmatrix}x^i\)，里面的组合数仍然可以递推，多项式乘法即可。

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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
#include <bits/stdc++.h>

#define int long long
using namespace std;

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;

using u128 = unsigned __int128;
using i128 = __int128;
using ld = long double;

constexpr int P = 1004535809;
constexpr int g = 3;

int power(int a, int b) {
    int res = 1;
    for (; b; b /= 2, a = 1LL * a * a % P) {
        if (b % 2) {
            res = 1LL * res * a % P;
        }
    }
    return res;
}

std::vector<int> rev, roots{0, 1};

void dft(std::vector<int> &a) {
    int n = a.size();
    if ((int)(rev.size()) != n) {
        int k = __builtin_ctz(n) - 1;
        rev.resize(n);
        for (int i = 0; i < n; i++) {
            rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
        }
    }
    for (int i = 0; i < n; i++) {
        if (rev[i] < i) {
            std::swap(a[i], a[rev[i]]);
        }
    }
    if (roots.size() < n) {
        int k = __builtin_ctz(roots.size());
        roots.resize(n);
        while ((1 << k) < n) {
            int e = power(g, (P - 1) >> (k + 1));
            for (int i = 1 << (k - 1); i < (1 << k); i++) {
                roots[2 * i] = roots[i];
                roots[2 * i + 1] = 1LL * roots[i] * e % P;
            }
            k++;
        }
    }

    for (int k = 1; k < n; k *= 2) {
        for (int i = 0; i < n; i += 2 * k) {
            for (int j = 0; j < k; j++) {
                int u = a[i + j];
                int v = 1LL * a[i + j + k] * roots[k + j] % P;
                a[i + j] = (u + v) % P;
                a[i + j + k] = (u - v + P) % P;
            }
        }
    }
}

void idft(std::vector<int> &a) {
    int n = a.size();
    std::reverse(a.begin() + 1, a.end());
    dft(a);
    int inv = power(n, P - 2);
    for (int i = 0; i < n; i++) {
        a[i] = 1LL * a[i] * inv % P;
    }
}

std::vector<int> mul(std::vector<int> a, std::vector<int> b) {
    int n = 1, tot = a.size() + b.size() - 1;
    while (n < tot) {
        n *= 2;
    }
    if (tot < 128) {
        std::vector<int> c(a.size() + b.size() - 1);
        for (int i = 0; i < a.size(); i++) {
            for (int j = 0; j < b.size(); j++) {
                c[i + j] = (c[i + j] + 1LL * a[i] * b[j]) % P;
            }
        }
        return c;
    }
    a.resize(n);
    b.resize(n);
    dft(a);
    dft(b);
    for (int i = 0; i < n; i++) {
        a[i] = 1LL * a[i] * b[i] % P;
    }
    idft(a);
    a.resize(tot);
    return a;
}

int read() {
    string s;
    cin >> s;
    int x = 0;
    for (auto c : s) {
        x *= 10;
        x += (c - '0');
        x %= P;
    }
    return x;
}

signed main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int n, k, t;
    cin >> n;
    k = read();
    cin >> t;
    vector<int> a(n + 1);
    for (int i = 1; i <= n; i++) {
        cin >> a[i];
    }

    vector<int> inv(n + 1);
    inv[1] = 1;
    for (int i = 2; i <= n; i++) {
        inv[i] = (P - P / i) * inv[P % i] % P;
    }

    if (t == 0) {
        vector<int> g(n + 1);
        g[0] = 1;
        for (int i = 1; i <= n; i++) {
            g[i] = g[i - 1] * (i + k - 1) % P * inv[i] % P;
        }
        auto c = mul(a, g);
        for (int i = 1; i <= n; i++) {
            cout << c[i] << " \n"[i == n];
        }
    } else {
        vector<int> g(n + 1);
        g[0] = 1;
        for (int i = 1; i <= n; i++) {
            g[i] = (-g[i - 1] * inv[i] % P * (k - i + 1) % P + P) % P;
        }
        auto c = mul(a, g);
        for (int i = 1; i <= n; i++) {
            cout << c[i] << " \n"[i == n];
        }
    }

    return 0;
}

CF 755 G

题意

\(n\) 个不同的球排成一排，定义一个组可以只包含一个球，或者包含两个相邻的球，一个球只能分在一组里（可以不分配），求将这 \(n\) 个球分成 \(k\) 组的方案数。

\(n\le 10^9,k < 2^{15}\)

题解

先列出朴素的递推关系，设 \(f_{i,j}\) 表示将前 \(i\) 个球分成 \(j\) 组的方案数，那么对于第 \(i\) 个球，它可以自成一组，可以不放入任何一组，也可以和第 \(i - 1\) 个球组成一组，那么 \(f_{i,j} = f_{i-1,j-1} + f_{i-1,j} + f_{i-2,j-1}\)。

考虑用生成函数处理递推关系，我们设 \(F_i(x) = \sum_{j=0}^kf_jx^j\)，那么转移可以写成 \(F_i(x) = xF_{i-1}(x)+F_{i-1}(x) + xF_{i-2}(x)\)。即 \(F_i(x) = (1+x)F_{i-1}(x) + xF_{i-2}(x)\)，这个递推式可以通过矩阵快速幂加速优化，即 \(\begin{bmatrix}F_n\,F_{n-1}\end{bmatrix}\begin{bmatrix}1+x&1\\x&0\end{bmatrix} = \begin{bmatrix}F_{n+1}\,F_{n}\end{bmatrix}\)。我们将矩阵的元素设为多项式就可以 \(O(k\log k\log n)\) 进行计算，要注意多项式乘积时要时刻取模 \(x^{k+1}\)。

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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
#include <bits/stdc++.h>

#define int long long
using namespace std;

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;

using u128 = unsigned __int128;
using i128 = __int128;
using ld = long double;

template <class T>
constexpr T power(T a, i64 b) {
    T res = 1;
    for (; b; b /= 2, a *= a) {
        if (b % 2) {
            res *= a;
        }
    }
    return res;
}

constexpr i64 mul(i64 a, i64 b, i64 p) {
    i64 res = a * b - i64(1.L * a * b / p) * p;
    res %= p;
    if (res < 0) {
        res += p;
    }
    return res;
}
template <int P>
struct MInt {
    int x;
    constexpr MInt() : x{} {}
    constexpr MInt(i64 x) : x{norm(x % getMod())} {}

    static int Mod;
    constexpr static int getMod() {
        if (P > 0) {
            return P;
        } else {
            return Mod;
        }
    }
    constexpr static void setMod(int Mod_) { Mod = Mod_; }
    constexpr int norm(int x) const {
        if (x < 0) {
            x += getMod();
        }
        if (x >= getMod()) {
            x -= getMod();
        }
        return x;
    }
    constexpr int val() const { return x; }
    explicit constexpr operator int() const { return x; }
    constexpr MInt operator-() const {
        MInt res;
        res.x = norm(getMod() - x);
        return res;
    }
    constexpr MInt inv() const {
        assert(x != 0);
        return power(*this, getMod() - 2);
    }
    constexpr MInt &operator*=(MInt rhs) & {
        x = 1LL * x * rhs.x % getMod();
        return *this;
    }
    constexpr MInt &operator+=(MInt rhs) & {
        x = norm(x + rhs.x);
        return *this;
    }
    constexpr MInt &operator-=(MInt rhs) & {
        x = norm(x - rhs.x);
        return *this;
    }
    constexpr MInt &operator/=(MInt rhs) & { return *this *= rhs.inv(); }
    friend constexpr MInt operator*(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res *= rhs;
        return res;
    }
    friend constexpr MInt operator+(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res += rhs;
        return res;
    }
    friend constexpr MInt operator-(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res -= rhs;
        return res;
    }
    friend constexpr MInt operator/(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res /= rhs;
        return res;
    }
    friend constexpr std::istream &operator>>(std::istream &is, MInt &a) {
        i64 v;
        is >> v;
        a = MInt(v);
        return is;
    }
    friend constexpr std::ostream &operator<<(std::ostream &os, const MInt &a) {
        return os << a.val();
    }
    friend constexpr bool operator==(MInt lhs, MInt rhs) {
        return lhs.val() == rhs.val();
    }
    friend constexpr bool operator!=(MInt lhs, MInt rhs) {
        return lhs.val() != rhs.val();
    }
};
template <>
int MInt<0>::Mod = 998244353;

template <int V, int P>
constexpr MInt<P> CInv = MInt<P>(V).inv();

constexpr int P = 998244353;
using Z = MInt<P>;  // 静态模数，题目给定模数
using D = MInt<0>;  // 动态模数，可以 D::setMod(x) 改变

std::vector<int> rev;
std::vector<Z> roots{0, 1};
void dft(std::vector<Z> &a) {
    int n = a.size();

    if ((int)(rev.size()) != n) {
        int k = __builtin_ctz(n) - 1;
        rev.resize(n);
        for (int i = 0; i < n; i++) {
            rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
        }
    }

    for (int i = 0; i < n; i++) {
        if (rev[i] < i) {
            std::swap(a[i], a[rev[i]]);
        }
    }
    if ((int)(roots.size()) < n) {
        int k = __builtin_ctz(roots.size());
        roots.resize(n);
        while ((1 << k) < n) {
            Z e = power(Z(3), (P - 1) >> (k + 1));
            for (int i = 1 << (k - 1); i < (1 << k); i++) {
                roots[2 * i] = roots[i];
                roots[2 * i + 1] = roots[i] * e;
            }
            k++;
        }
    }
    for (int k = 1; k < n; k *= 2) {
        for (int i = 0; i < n; i += 2 * k) {
            for (int j = 0; j < k; j++) {
                Z u = a[i + j];
                Z v = a[i + j + k] * roots[k + j];
                a[i + j] = u + v;
                a[i + j + k] = u - v;
            }
        }
    }
}
void idft(std::vector<Z> &a) {
    int n = a.size();
    std::reverse(a.begin() + 1, a.end());
    dft(a);
    Z inv = (1 - P) / n;
    for (int i = 0; i < n; i++) {
        a[i] *= inv;
    }
}
struct Poly {
    std::vector<Z> a;
    Poly() {}
    explicit Poly(
        int size, std::function<Z(int)> f = [](int) { return 0; })
        : a(size) {
        for (int i = 0; i < size; i++) {
            a[i] = f(i);
        }
    }
    Poly(const std::vector<Z> &a) : a(a) {}
    Poly(const std::initializer_list<Z> &a) : a(a) {}
    int size() const { return a.size(); }
    void resize(int n) { a.resize(n); }
    Z operator[](int idx) const {
        if (idx < size()) {
            return a[idx];
        } else {
            return 0;
        }
    }
    Z &operator[](int idx) { return a[idx]; }
    Poly mulxk(int k) const {
        auto b = a;
        b.insert(b.begin(), k, 0);
        return Poly(b);
    }
    Poly modxk(int k) const {
        k = std::min(k, size());
        return Poly(std::vector<Z>(a.begin(), a.begin() + k));
    }
    Poly divxk(int k) const {
        if (size() <= k) {
            return Poly();
        }
        return Poly(std::vector<Z>(a.begin() + k, a.end()));
    }
    friend Poly operator+(const Poly &a, const Poly &b) {
        std::vector<Z> res(std::max(a.size(), b.size()));
        for (int i = 0; i < (int)(res.size()); i++) {
            res[i] = a[i] + b[i];
        }
        return Poly(res);
    }
    friend Poly operator-(const Poly &a, const Poly &b) {
        std::vector<Z> res(std::max(a.size(), b.size()));
        for (int i = 0; i < (int)(res.size()); i++) {
            res[i] = a[i] - b[i];
        }
        return Poly(res);
    }
    friend Poly operator-(const Poly &a) {
        std::vector<Z> res(a.size());
        for (int i = 0; i < (int)(res.size()); i++) {
            res[i] = -a[i];
        }
        return Poly(res);
    }
    friend Poly operator*(Poly a, Poly b) {
        if (a.size() == 0 || b.size() == 0) {
            return Poly();
        }
        if (a.size() < b.size()) {
            std::swap(a, b);
        }
        if (b.size() < 128) {
            Poly c(a.size() + b.size() - 1);
            for (int i = 0; i < a.size(); i++) {
                for (int j = 0; j < b.size(); j++) {
                    c[i + j] += a[i] * b[j];
                }
            }
            return c;
        }
        int sz = 1, tot = a.size() + b.size() - 1;
        while (sz < tot) {
            sz *= 2;
        }
        a.a.resize(sz);
        b.a.resize(sz);
        dft(a.a);
        dft(b.a);
        for (int i = 0; i < sz; ++i) {
            a.a[i] = a[i] * b[i];
        }
        idft(a.a);
        a.resize(tot);
        return a;
    }
    friend Poly operator*(Z a, Poly b) {
        for (int i = 0; i < (int)(b.size()); i++) {
            b[i] *= a;
        }
        return b;
    }
    friend Poly operator*(Poly a, Z b) {
        for (int i = 0; i < (int)(a.size()); i++) {
            a[i] *= b;
        }
        return a;
    }
    Poly &operator+=(Poly b) { return (*this) = (*this) + b; }
    Poly &operator-=(Poly b) { return (*this) = (*this) - b; }
    Poly &operator*=(Poly b) { return (*this) = (*this) * b; }
    Poly &operator*=(Z b) { return (*this) = (*this) * b; }
    Poly deriv() const {
        if (a.empty()) {
            return Poly();
        }
        std::vector<Z> res(size() - 1);
        for (int i = 0; i < size() - 1; ++i) {
            res[i] = (i + 1) * a[i + 1];
        }
        return Poly(res);
    }
    Poly integr() const {
        std::vector<Z> res(size() + 1);
        for (int i = 0; i < size(); ++i) {
            res[i + 1] = a[i] / (i + 1);
        }
        return Poly(res);
    }
    Poly inv(int m) const {
        Poly x{a[0].inv()};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (Poly{2} - modxk(k) * x)).modxk(k);
        }
        return x.modxk(m);
    }
    Poly log(int m) const { return (deriv() * inv(m)).integr().modxk(m); }
    Poly exp(int m) const {
        Poly x{1};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (Poly{1} - x.log(k) + modxk(k))).modxk(k);
        }
        return x.modxk(m);
    }
    Poly pow(int k, int m) const {
        int i = 0;
        while (i < size() && a[i].val() == 0) {
            i++;
        }
        if (i == size() || 1LL * i * k >= m) {
            return Poly(std::vector<Z>(m));
        }
        Z v = a[i];
        auto f = divxk(i) * v.inv();
        return (f.log(m - i * k) * k).exp(m - i * k).mulxk(i * k) * power(v, k);
    }
    Poly sqrt(int m) const {
        Poly x{1};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x + (modxk(k) * x.inv(k)).modxk(k)) * ((P + 1) / 2);
        }
        return x.modxk(m);
    }
    Poly mulT(Poly b) const {
        if (b.size() == 0) {
            return Poly();
        }
        int n = b.size();
        std::reverse(b.a.begin(), b.a.end());
        return ((*this) * b).divxk(n - 1);
    }
    std::vector<Z> eval(std::vector<Z> x) const {
        if (size() == 0) {
            return std::vector<Z>(x.size(), 0);
        }
        const int n = std::max((int)(x.size()), size());
        std::vector<Poly> q(4 * n);
        std::vector<Z> ans(x.size());
        x.resize(n);
        std::function<void(int, int, int)> build = [&](int p, int l, int r) {
            if (r - l == 1) {
                q[p] = Poly{1, -x[l]};
            } else {
                int m = (l + r) / 2;
                build(2 * p, l, m);
                build(2 * p + 1, m, r);
                q[p] = q[2 * p] * q[2 * p + 1];
            }
        };
        build(1, 0, n);
        std::function<void(int, int, int, const Poly &)> work =
            [&](int p, int l, int r, const Poly &num) {
                if (r - l == 1) {
                    if (l < (int)(ans.size())) {
                        ans[l] = num[0];
                    }
                } else {
                    int m = (l + r) / 2;
                    work(2 * p, l, m, num.mulT(q[2 * p + 1]).modxk(m - l));
                    work(2 * p + 1, m, r, num.mulT(q[2 * p]).modxk(r - m));
                }
            };
        work(1, 0, n, mulT(q[1].inv(n)));
        return ans;
    }
};

int len;

struct Matrix {
    int n;
    vector<vector<Poly>> a;
    Matrix() {}
    Matrix(int n) : n(n) { a.assign(n, vector<Poly>(n)); }
    vector<Poly> &operator[](int row) { return a[row]; }
    const vector<Poly> &operator[](int row) const { return a[row]; }
    Matrix operator*(const Matrix &b) const {
        Matrix res(n);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                for (int k = 0; k < n; k++) {
                    res[i][j] = (res[i][j] + a[i][k] * b[k][j]).modxk(len + 1);
                }
            }
        }
        return res;
    }
    Matrix &operator*=(const Matrix &b) {
        *this = *this * b;
        return *this;
    }
    Matrix power(int b) {
        Matrix res(n);
        auto A = *this;
        for (int i = 0; i < n; i++) {
            res[i][i] = {1};
        }
        while (b) {
            if (b & 1) {
                res = res * A;
            }
            A *= A;
            b /= 2;
        }
        return res;
    }
};

signed main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int n, k;
    cin >> n >> k;
    len = k;
    Matrix A(2), B(2);
    A[0][0] = {1, 1}, A[0][1] = {1};
    B[0][0] = {1, 1};
    B[0][1] = {1};
    B[1][0] = {0, 1};
    B[1][1] = {0};
    auto res = A * B.power(n - 1);
    for (int i = 1; i <= k; i++) {
        cout << (i < res[0][0].size() ? res[0][0][i] : 0) << " \n"[i == k];
    }

    return 0;
}

P4389 付公主的背包

题意

\(n\) 个物品体积是 \(v_i\)。给定 \(m\)，对于 \(s\in[1,m]\) 回答用这些物品恰好装 \(s\) 体积的方案数，模 998244353。

\(n,m\le 10^5\)

题解

数据范围限制了我们使用 dp。对于单个体积为 \(v\) 的物品我们可以构造生成函数 \(A(x) = \sum\limits_{i=0}^\infin[i\equiv0\pmod v]x^i\)，含义是只有 \(i\) 为 \(k\) 的倍数时才有 1 个方案，否则没有。那么 \(n\) 个物品就是 \(F(x) = \prod\limits_{i=1}^nA_i(x)\)，复杂度 \(O(nm\log m)\)。

考虑优化，上面是多项式的乘积，转化为对数是很自然的想法，有一个重要的结论：\(\ln(1-x^v) = -\sum\limits_{i\ge 1}\dfrac{x^{vi}}{i}\)，在末尾给出证明。

那么 \(F(x) = e^{\ln F(x)} = e^{\sum_\limits{i=1}^n\sum\limits_{j\ge 1}\frac{x^{v_ij}}{j}}\)。指数作为多项式求 exp 即可。然而指数直接进行加法也是 \(O(nm)\) 的，所以我们把 \(v_i\) 相同的打包在一起计算，这样的复杂度是调和级数 \(O(m\log m)\)​。

给出 \(\ln(1-x^v) = -\sum\limits_{i\ge 1}\dfrac{x^{vi}}{i}\) 的证明：

设 \(F(x) = 1-x^v\)，\(\ln F(x) = G(x)\)。两边求导有 \(\dfrac{F^\prime(x)}{F(x)} = \dfrac{-vx^{v-1}}{1-x^v} =G^\prime(x)\)，将 \(\dfrac{1}{1-x^v}\) 展开为级数有：

\(G^\prime(x) = \sum\limits_{i=0}^\infin-vx^{v-1}x^{vi}= -\sum\limits_{i=0}^\infin vx^NaN\)，那么积分后得到\(G(x) = -\sum_{i=0}\limits^\infty\dfrac{x^{v(1+i)}}{1+i} = -\sum_\limits{i\ge 1}\dfrac{x^{vi}}{i}\)。

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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
#include <bits/stdc++.h>

#define int long long
using namespace std;

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;

using u128 = unsigned __int128;
using i128 = __int128;
using ld = long double;

template <class T>
constexpr T power(T a, i64 b) {
    T res = 1;
    for (; b; b /= 2, a *= a) {
        if (b % 2) {
            res *= a;
        }
    }
    return res;
}

constexpr i64 mul(i64 a, i64 b, i64 p) {
    i64 res = a * b - i64(1.L * a * b / p) * p;
    res %= p;
    if (res < 0) {
        res += p;
    }
    return res;
}
template <int P>
struct MInt {
    int x;
    constexpr MInt() : x{} {}
    constexpr MInt(i64 x) : x{norm(x % getMod())} {}

    static int Mod;
    constexpr static int getMod() {
        if (P > 0) {
            return P;
        } else {
            return Mod;
        }
    }
    constexpr static void setMod(int Mod_) { Mod = Mod_; }
    constexpr int norm(int x) const {
        if (x < 0) {
            x += getMod();
        }
        if (x >= getMod()) {
            x -= getMod();
        }
        return x;
    }
    constexpr int val() const { return x; }
    explicit constexpr operator int() const { return x; }
    constexpr MInt operator-() const {
        MInt res;
        res.x = norm(getMod() - x);
        return res;
    }
    constexpr MInt inv() const {
        assert(x != 0);
        return power(*this, getMod() - 2);
    }
    constexpr MInt &operator*=(MInt rhs) & {
        x = 1LL * x * rhs.x % getMod();
        return *this;
    }
    constexpr MInt &operator+=(MInt rhs) & {
        x = norm(x + rhs.x);
        return *this;
    }
    constexpr MInt &operator-=(MInt rhs) & {
        x = norm(x - rhs.x);
        return *this;
    }
    constexpr MInt &operator/=(MInt rhs) & { return *this *= rhs.inv(); }
    friend constexpr MInt operator*(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res *= rhs;
        return res;
    }
    friend constexpr MInt operator+(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res += rhs;
        return res;
    }
    friend constexpr MInt operator-(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res -= rhs;
        return res;
    }
    friend constexpr MInt operator/(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res /= rhs;
        return res;
    }
    friend constexpr std::istream &operator>>(std::istream &is, MInt &a) {
        i64 v;
        is >> v;
        a = MInt(v);
        return is;
    }
    friend constexpr std::ostream &operator<<(std::ostream &os, const MInt &a) {
        return os << a.val();
    }
    friend constexpr bool operator==(MInt lhs, MInt rhs) {
        return lhs.val() == rhs.val();
    }
    friend constexpr bool operator!=(MInt lhs, MInt rhs) {
        return lhs.val() != rhs.val();
    }
};
template <>
int MInt<0>::Mod = 998244353;

template <int V, int P>
constexpr MInt<P> CInv = MInt<P>(V).inv();

constexpr int P = 998244353;
using Z = MInt<P>;  // 静态模数，题目给定模数
using D = MInt<0>;  // 动态模数，可以 D::setMod(x) 改变

std::vector<int> rev;
std::vector<Z> roots{0, 1};
void dft(std::vector<Z> &a) {
    int n = a.size();

    if ((int)(rev.size()) != n) {
        int k = __builtin_ctz(n) - 1;
        rev.resize(n);
        for (int i = 0; i < n; i++) {
            rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
        }
    }

    for (int i = 0; i < n; i++) {
        if (rev[i] < i) {
            std::swap(a[i], a[rev[i]]);
        }
    }
    if ((int)(roots.size()) < n) {
        int k = __builtin_ctz(roots.size());
        roots.resize(n);
        while ((1 << k) < n) {
            Z e = power(Z(3), (P - 1) >> (k + 1));
            for (int i = 1 << (k - 1); i < (1 << k); i++) {
                roots[2 * i] = roots[i];
                roots[2 * i + 1] = roots[i] * e;
            }
            k++;
        }
    }
    for (int k = 1; k < n; k *= 2) {
        for (int i = 0; i < n; i += 2 * k) {
            for (int j = 0; j < k; j++) {
                Z u = a[i + j];
                Z v = a[i + j + k] * roots[k + j];
                a[i + j] = u + v;
                a[i + j + k] = u - v;
            }
        }
    }
}
void idft(std::vector<Z> &a) {
    int n = a.size();
    std::reverse(a.begin() + 1, a.end());
    dft(a);
    Z inv = (1 - P) / n;
    for (int i = 0; i < n; i++) {
        a[i] *= inv;
    }
}
struct Poly {
    std::vector<Z> a;
    Poly() {}
    explicit Poly(
        int size, std::function<Z(int)> f = [](int) { return 0; })
        : a(size) {
        for (int i = 0; i < size; i++) {
            a[i] = f(i);
        }
    }
    Poly(const std::vector<Z> &a) : a(a) {}
    Poly(const std::initializer_list<Z> &a) : a(a) {}
    int size() const { return a.size(); }
    void resize(int n) { a.resize(n); }
    Z operator[](int idx) const {
        if (idx < size()) {
            return a[idx];
        } else {
            return 0;
        }
    }
    Z &operator[](int idx) { return a[idx]; }
    Poly mulxk(int k) const {
        auto b = a;
        b.insert(b.begin(), k, 0);
        return Poly(b);
    }
    Poly modxk(int k) const {
        k = std::min(k, size());
        return Poly(std::vector<Z>(a.begin(), a.begin() + k));
    }
    Poly divxk(int k) const {
        if (size() <= k) {
            return Poly();
        }
        return Poly(std::vector<Z>(a.begin() + k, a.end()));
    }
    friend Poly operator+(const Poly &a, const Poly &b) {
        std::vector<Z> res(std::max(a.size(), b.size()));
        for (int i = 0; i < (int)(res.size()); i++) {
            res[i] = a[i] + b[i];
        }
        return Poly(res);
    }
    friend Poly operator-(const Poly &a, const Poly &b) {
        std::vector<Z> res(std::max(a.size(), b.size()));
        for (int i = 0; i < (int)(res.size()); i++) {
            res[i] = a[i] - b[i];
        }
        return Poly(res);
    }
    friend Poly operator-(const Poly &a) {
        std::vector<Z> res(a.size());
        for (int i = 0; i < (int)(res.size()); i++) {
            res[i] = -a[i];
        }
        return Poly(res);
    }
    friend Poly operator*(Poly a, Poly b) {
        if (a.size() == 0 || b.size() == 0) {
            return Poly();
        }
        if (a.size() < b.size()) {
            std::swap(a, b);
        }
        if (b.size() < 128) {
            Poly c(a.size() + b.size() - 1);
            for (int i = 0; i < a.size(); i++) {
                for (int j = 0; j < b.size(); j++) {
                    c[i + j] += a[i] * b[j];
                }
            }
            return c;
        }
        int sz = 1, tot = a.size() + b.size() - 1;
        while (sz < tot) {
            sz *= 2;
        }
        a.a.resize(sz);
        b.a.resize(sz);
        dft(a.a);
        dft(b.a);
        for (int i = 0; i < sz; ++i) {
            a.a[i] = a[i] * b[i];
        }
        idft(a.a);
        a.resize(tot);
        return a;
    }
    friend Poly operator*(Z a, Poly b) {
        for (int i = 0; i < (int)(b.size()); i++) {
            b[i] *= a;
        }
        return b;
    }
    friend Poly operator*(Poly a, Z b) {
        for (int i = 0; i < (int)(a.size()); i++) {
            a[i] *= b;
        }
        return a;
    }
    Poly &operator+=(Poly b) { return (*this) = (*this) + b; }
    Poly &operator-=(Poly b) { return (*this) = (*this) - b; }
    Poly &operator*=(Poly b) { return (*this) = (*this) * b; }
    Poly &operator*=(Z b) { return (*this) = (*this) * b; }
    Poly deriv() const {
        if (a.empty()) {
            return Poly();
        }
        std::vector<Z> res(size() - 1);
        for (int i = 0; i < size() - 1; ++i) {
            res[i] = (i + 1) * a[i + 1];
        }
        return Poly(res);
    }
    Poly integr() const {
        std::vector<Z> res(size() + 1);
        for (int i = 0; i < size(); ++i) {
            res[i + 1] = a[i] / (i + 1);
        }
        return Poly(res);
    }
    Poly inv(int m) const {
        Poly x{a[0].inv()};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (Poly{2} - modxk(k) * x)).modxk(k);
        }
        return x.modxk(m);
    }
    Poly log(int m) const { return (deriv() * inv(m)).integr().modxk(m); }
    Poly exp(int m) const {
        Poly x{1};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (Poly{1} - x.log(k) + modxk(k))).modxk(k);
        }
        return x.modxk(m);
    }
    Poly pow(int k, int m) const {
        int i = 0;
        while (i < size() && a[i].val() == 0) {
            i++;
        }
        if (i == size() || 1LL * i * k >= m) {
            return Poly(std::vector<Z>(m));
        }
        Z v = a[i];
        auto f = divxk(i) * v.inv();
        return (f.log(m - i * k) * k).exp(m - i * k).mulxk(i * k) * power(v, k);
    }
    Poly sqrt(int m) const {
        Poly x{1};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x + (modxk(k) * x.inv(k)).modxk(k)) * ((P + 1) / 2);
        }
        return x.modxk(m);
    }
    Poly mulT(Poly b) const {
        if (b.size() == 0) {
            return Poly();
        }
        int n = b.size();
        std::reverse(b.a.begin(), b.a.end());
        return ((*this) * b).divxk(n - 1);
    }
    std::vector<Z> eval(std::vector<Z> x) const {
        if (size() == 0) {
            return std::vector<Z>(x.size(), 0);
        }
        const int n = std::max((int)(x.size()), size());
        std::vector<Poly> q(4 * n);
        std::vector<Z> ans(x.size());
        x.resize(n);
        std::function<void(int, int, int)> build = [&](int p, int l, int r) {
            if (r - l == 1) {
                q[p] = Poly{1, -x[l]};
            } else {
                int m = (l + r) / 2;
                build(2 * p, l, m);
                build(2 * p + 1, m, r);
                q[p] = q[2 * p] * q[2 * p + 1];
            }
        };
        build(1, 0, n);
        std::function<void(int, int, int, const Poly &)> work =
            [&](int p, int l, int r, const Poly &num) {
                if (r - l == 1) {
                    if (l < (int)(ans.size())) {
                        ans[l] = num[0];
                    }
                } else {
                    int m = (l + r) / 2;
                    work(2 * p, l, m, num.mulT(q[2 * p + 1]).modxk(m - l));
                    work(2 * p + 1, m, r, num.mulT(q[2 * p]).modxk(r - m));
                }
            };
        work(1, 0, n, mulT(q[1].inv(n)));
        return ans;
    }
};

signed main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int n, m;
    cin >> n >> m;
    vector<int> buk(m + 1);
    for (int i = 1; i <= n; i++) {
        int v;
        cin >> v;
        buk[v]++;
    }
    vector<int> inv(m + 1);
    inv[1] = 1;
    for (int i = 2; i <= m; i++) {
        inv[i] = (P - P / i) * inv[P % i] % P;
    }
    vector<Z> a(m + 1);
    for (int i = 1; i <= m; i++) {
        if (buk[i]) {
            for (int j = i; j <= m; j += i) {
                a[j] += Z(inv[j / i]) * buk[i];
            }
        }
    }
    Poly p(a);
    p = p.exp(m + 1);
    for (int i = 1; i <= m; i++) {
        cout << p[i] << "\n";
    }

    return 0;
}

P5824 十二重计数法

题意

\(n\) 个球，\(m\) 个盒子，在以下条件下计数，模 998244353:

I：球之间互不相同，盒子之间互不相同。 II：球之间互不相同，盒子之间互不相同，每个盒子至多装一个球。 III：球之间互不相同，盒子之间互不相同，每个盒子至少装一个球。

IV：球之间互不相同，盒子全部相同。 V：球之间互不相同，盒子全部相同，每个盒子至多装一个球。 VI：球之间互不相同，盒子全部相同，每个盒子至少装一个球。

VII：球全部相同，盒子之间互不相同。 VIII：球全部相同，盒子之间互不相同，每个盒子至多装一个球。 IX：球全部相同，盒子之间互不相同，每个盒子至少装一个球。

X：球全部相同，盒子全部相同。 XI：球全部相同，盒子全部相同，每个盒子至多装一个球。 XII：球全部相同，盒子全部相同，每个盒子至少装一个球。

\(n,m\le 2\times 10^5\)

题解

  • 球不同，盒子不同，无限制：每个球都能放到任意盒子中，\(m^n\)。
  • 球不同，盒子不同，每个盒子至多 1 个球：每个球找一个没放过的盒子放进去，\(m^{\underline{n}}\)。
  • 球不同，盒子不同，每个盒子至少 1 个球：考虑容斥，设 \(f(S)\) 代表集合 \(S\) 中的盒子不放球的方案数，那么有 \(f(S) = (m - |S|)^n\)，且 \(f(S) =\sum_{S\subseteq T}\limits g(T)\)，根据二项式反演，\(g(S) = \sum_{S\subseteq T}\limits(-1)^{|T|-|S|}f(T)\)，答案是 \(g(\varnothing) = \sum_{i=0}^m(-1)^i\dbinom{m}{i}(m-i)^n\)。
  • 球不同，盒子同，无限制：盒子相同要想到斯特林数，答案就是 \(\sum_{i=1}^m\left\{\begin{matrix}n\\i\end{matrix}\right\}\)。求一行斯特林数通过斯特林数的卷积公式来求。
  • 球不同，盒子同，每个盒子至多 1 个球：每个球放哪个盒子都是一样的，答案就是 \([n\le m]\)。
  • 球不同，盒子同，每个盒子至少 1 个球：根据第二类斯特林数的定义，答案就是 \(\left\{\begin{matrix}n\\m\end{matrix}\right\}\)​。
  • 球同，盒子不同，无限制：球相同要想到插板法，这里分成的 \(m\) 个集合可以为空，那么答案就是 \(\dbinom{n + m - 1}{m - 1}\)。
  • 球同，盒子不同，每个盒子至多 1 个球：取 \(n\) 个盒子放球，\(\dbinom{m}{n}\)
  • 球同，盒子不同，每个盒子至少 1 个球：插板法，\(\dbinom{n - 1}{m - 1}\)。
  • 球同，盒子同，无限制：这里需要引入分拆数——将自然数 \(n\) 写成递降正整数的和的方案数 \(p_n\)，如果将 \(n\) 恰好分成 \(k\) 部分，就叫做 \(k\) 部分拆数 \(p_{n,k}\)。显然 \(k\) 部分拆数也是方程 \(n - k = y_1 + y_2+...+y_k,y_1\ge y_2\ge ...\ge y_k\ge 0\) 的解的数量。这里我们定义 \(f(n,m)\) 为将整数 \(n\) 划分成不超过 \(m\) 个非负整数的方案数，那么 \(f\) 的转移可以向后面加一个 0，也可以全体加 1，有 \(f(n,m) = f(n - m,m) + f(n, m - 1)\)。我们令 \(F_k(x) = \sum_{i=0}^\infty\limits x^if(i,k)\)，那么 \(F_k(x) = \sum_{i=0}^\infty\limits (x^if(i - k, k) + x^if(i, k - 1)) = F_{k-1} + x^kF_{k}\)。因此 \(F_k(x) = \dfrac{F_{k-1}(x)}{1-x^k}\)，且 \(F_1(x) = 1 + x + x^2+... = \dfrac{1}{1 - x}\)，因此可以推导出 \(F_k(x) = \prod_{i=1}^k\limits\dfrac{1}{1 - x^i}\)，类似付公主的背包求多项式 exp 即可。答案是 \([x^n]F_m(x)\)。
  • 球同，盒子同，每个盒子至多 1 个球：\([n\le m]\)。
  • 球同，盒子同，每个盒子至少 1 个球：先往每个盒子放一个球，那么答案就是 \([x^{n-m}]F_m(x)\)。
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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
#include <bits/stdc++.h>

#define int long long
using namespace std;

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;

using u128 = unsigned __int128;
using i128 = __int128;
using ld = long double;

template <class T>
constexpr T power(T a, i64 b) {
    T res = 1;
    for (; b; b /= 2, a *= a) {
        if (b % 2) {
            res *= a;
        }
    }
    return res;
}

constexpr i64 mul(i64 a, i64 b, i64 p) {
    i64 res = a * b - i64(1.L * a * b / p) * p;
    res %= p;
    if (res < 0) {
        res += p;
    }
    return res;
}
template <int P>
struct MInt {
    int x;
    constexpr MInt() : x{} {}
    constexpr MInt(i64 x) : x{norm(x % getMod())} {}

    static int Mod;
    constexpr static int getMod() {
        if (P > 0) {
            return P;
        } else {
            return Mod;
        }
    }
    constexpr static void setMod(int Mod_) { Mod = Mod_; }
    constexpr int norm(int x) const {
        if (x < 0) {
            x += getMod();
        }
        if (x >= getMod()) {
            x -= getMod();
        }
        return x;
    }
    constexpr int val() const { return x; }
    explicit constexpr operator int() const { return x; }
    constexpr MInt operator-() const {
        MInt res;
        res.x = norm(getMod() - x);
        return res;
    }
    constexpr MInt inv() const {
        assert(x != 0);
        return power(*this, getMod() - 2);
    }
    constexpr MInt &operator*=(MInt rhs) & {
        x = 1LL * x * rhs.x % getMod();
        return *this;
    }
    constexpr MInt &operator+=(MInt rhs) & {
        x = norm(x + rhs.x);
        return *this;
    }
    constexpr MInt &operator-=(MInt rhs) & {
        x = norm(x - rhs.x);
        return *this;
    }
    constexpr MInt &operator/=(MInt rhs) & { return *this *= rhs.inv(); }
    friend constexpr MInt operator*(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res *= rhs;
        return res;
    }
    friend constexpr MInt operator+(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res += rhs;
        return res;
    }
    friend constexpr MInt operator-(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res -= rhs;
        return res;
    }
    friend constexpr MInt operator/(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res /= rhs;
        return res;
    }
    friend constexpr std::istream &operator>>(std::istream &is, MInt &a) {
        i64 v;
        is >> v;
        a = MInt(v);
        return is;
    }
    friend constexpr std::ostream &operator<<(std::ostream &os, const MInt &a) {
        return os << a.val();
    }
    friend constexpr bool operator==(MInt lhs, MInt rhs) {
        return lhs.val() == rhs.val();
    }
    friend constexpr bool operator!=(MInt lhs, MInt rhs) {
        return lhs.val() != rhs.val();
    }
};
template <>
int MInt<0>::Mod = 998244353;

template <int V, int P>
constexpr MInt<P> CInv = MInt<P>(V).inv();

constexpr int P = 998244353;
using Z = MInt<P>;  // 静态模数，题目给定模数
using D = MInt<0>;  // 动态模数，可以 D::setMod(x) 改变

struct Comb {
    int n;
    std::vector<Z> _fac;
    std::vector<Z> _invfac;
    std::vector<Z> _inv;

    Comb() : n{0}, _fac{1}, _invfac{1}, _inv{0} {}
    Comb(int n) : Comb() { init(n); }

    void init(int m) {
        if (m <= n) return;
        _fac.resize(m + 1);
        _invfac.resize(m + 1);
        _inv.resize(m + 1);

        for (int i = n + 1; i <= m; i++) {
            _fac[i] = _fac[i - 1] * i;
        }
        _invfac[m] = _fac[m].inv();
        for (int i = m; i > n; i--) {
            _invfac[i - 1] = _invfac[i] * i;
            _inv[i] = _invfac[i] * _fac[i - 1];
        }
        n = m;
    }

    Z fac(int m) {
        if (m > n) init(2 * m);
        return _fac[m];
    }
    Z invfac(int m) {
        if (m > n) init(2 * m);
        return _invfac[m];
    }
    Z inv(int m) {
        if (m > n) init(2 * m);
        return _inv[m];
    }
    Z binom(int n, int m) {
        if (n < m || m < 0) return 0;
        return fac(n) * invfac(m) * invfac(n - m);
    }
} comb(200000);  // Z 类型全局初始化，D 类型需要局部初始化

std::vector<int> rev;
std::vector<Z> roots{0, 1};
void dft(std::vector<Z> &a) {
    int n = a.size();

    if ((int)(rev.size()) != n) {
        int k = __builtin_ctz(n) - 1;
        rev.resize(n);
        for (int i = 0; i < n; i++) {
            rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
        }
    }

    for (int i = 0; i < n; i++) {
        if (rev[i] < i) {
            std::swap(a[i], a[rev[i]]);
        }
    }
    if ((int)(roots.size()) < n) {
        int k = __builtin_ctz(roots.size());
        roots.resize(n);
        while ((1 << k) < n) {
            Z e = power(Z(3), (P - 1) >> (k + 1));
            for (int i = 1 << (k - 1); i < (1 << k); i++) {
                roots[2 * i] = roots[i];
                roots[2 * i + 1] = roots[i] * e;
            }
            k++;
        }
    }
    for (int k = 1; k < n; k *= 2) {
        for (int i = 0; i < n; i += 2 * k) {
            for (int j = 0; j < k; j++) {
                Z u = a[i + j];
                Z v = a[i + j + k] * roots[k + j];
                a[i + j] = u + v;
                a[i + j + k] = u - v;
            }
        }
    }
}
void idft(std::vector<Z> &a) {
    int n = a.size();
    std::reverse(a.begin() + 1, a.end());
    dft(a);
    Z inv = (1 - P) / n;
    for (int i = 0; i < n; i++) {
        a[i] *= inv;
    }
}
struct Poly {
    std::vector<Z> a;
    Poly() {}
    explicit Poly(
        int size, std::function<Z(int)> f = [](int) { return 0; })
        : a(size) {
        for (int i = 0; i < size; i++) {
            a[i] = f(i);
        }
    }
    Poly(const std::vector<Z> &a) : a(a) {}
    Poly(const std::initializer_list<Z> &a) : a(a) {}
    int size() const { return a.size(); }
    void resize(int n) { a.resize(n); }
    Z operator[](int idx) const {
        if (idx < size()) {
            return a[idx];
        } else {
            return 0;
        }
    }
    Z &operator[](int idx) { return a[idx]; }
    Poly mulxk(int k) const {
        auto b = a;
        b.insert(b.begin(), k, 0);
        return Poly(b);
    }
    Poly modxk(int k) const {
        k = std::min(k, size());
        return Poly(std::vector<Z>(a.begin(), a.begin() + k));
    }
    Poly divxk(int k) const {
        if (size() <= k) {
            return Poly();
        }
        return Poly(std::vector<Z>(a.begin() + k, a.end()));
    }
    friend Poly operator+(const Poly &a, const Poly &b) {
        std::vector<Z> res(std::max(a.size(), b.size()));
        for (int i = 0; i < (int)(res.size()); i++) {
            res[i] = a[i] + b[i];
        }
        return Poly(res);
    }
    friend Poly operator-(const Poly &a, const Poly &b) {
        std::vector<Z> res(std::max(a.size(), b.size()));
        for (int i = 0; i < (int)(res.size()); i++) {
            res[i] = a[i] - b[i];
        }
        return Poly(res);
    }
    friend Poly operator-(const Poly &a) {
        std::vector<Z> res(a.size());
        for (int i = 0; i < (int)(res.size()); i++) {
            res[i] = -a[i];
        }
        return Poly(res);
    }
    friend Poly operator*(Poly a, Poly b) {
        if (a.size() == 0 || b.size() == 0) {
            return Poly();
        }
        if (a.size() < b.size()) {
            std::swap(a, b);
        }
        if (b.size() < 128) {
            Poly c(a.size() + b.size() - 1);
            for (int i = 0; i < a.size(); i++) {
                for (int j = 0; j < b.size(); j++) {
                    c[i + j] += a[i] * b[j];
                }
            }
            return c;
        }
        int sz = 1, tot = a.size() + b.size() - 1;
        while (sz < tot) {
            sz *= 2;
        }
        a.a.resize(sz);
        b.a.resize(sz);
        dft(a.a);
        dft(b.a);
        for (int i = 0; i < sz; ++i) {
            a.a[i] = a[i] * b[i];
        }
        idft(a.a);
        a.resize(tot);
        return a;
    }
    friend Poly operator*(Z a, Poly b) {
        for (int i = 0; i < (int)(b.size()); i++) {
            b[i] *= a;
        }
        return b;
    }
    friend Poly operator*(Poly a, Z b) {
        for (int i = 0; i < (int)(a.size()); i++) {
            a[i] *= b;
        }
        return a;
    }
    Poly &operator+=(Poly b) { return (*this) = (*this) + b; }
    Poly &operator-=(Poly b) { return (*this) = (*this) - b; }
    Poly &operator*=(Poly b) { return (*this) = (*this) * b; }
    Poly &operator*=(Z b) { return (*this) = (*this) * b; }
    Poly deriv() const {
        if (a.empty()) {
            return Poly();
        }
        std::vector<Z> res(size() - 1);
        for (int i = 0; i < size() - 1; ++i) {
            res[i] = (i + 1) * a[i + 1];
        }
        return Poly(res);
    }
    Poly integr() const {
        std::vector<Z> res(size() + 1);
        for (int i = 0; i < size(); ++i) {
            res[i + 1] = a[i] / (i + 1);
        }
        return Poly(res);
    }
    Poly inv(int m) const {
        Poly x{a[0].inv()};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (Poly{2} - modxk(k) * x)).modxk(k);
        }
        return x.modxk(m);
    }
    Poly log(int m) const { return (deriv() * inv(m)).integr().modxk(m); }
    Poly exp(int m) const {
        Poly x{1};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (Poly{1} - x.log(k) + modxk(k))).modxk(k);
        }
        return x.modxk(m);
    }
    Poly pow(int k, int m) const {
        int i = 0;
        while (i < size() && a[i].val() == 0) {
            i++;
        }
        if (i == size() || 1LL * i * k >= m) {
            return Poly(std::vector<Z>(m));
        }
        Z v = a[i];
        auto f = divxk(i) * v.inv();
        return (f.log(m - i * k) * k).exp(m - i * k).mulxk(i * k) * power(v, k);
    }
    Poly sqrt(int m) const {
        Poly x{1};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x + (modxk(k) * x.inv(k)).modxk(k)) * ((P + 1) / 2);
        }
        return x.modxk(m);
    }
    Poly mulT(Poly b) const {
        if (b.size() == 0) {
            return Poly();
        }
        int n = b.size();
        std::reverse(b.a.begin(), b.a.end());
        return ((*this) * b).divxk(n - 1);
    }
    std::vector<Z> eval(std::vector<Z> x) const {
        if (size() == 0) {
            return std::vector<Z>(x.size(), 0);
        }
        const int n = std::max((int)(x.size()), size());
        std::vector<Poly> q(4 * n);
        std::vector<Z> ans(x.size());
        x.resize(n);
        std::function<void(int, int, int)> build = [&](int p, int l, int r) {
            if (r - l == 1) {
                q[p] = Poly{1, -x[l]};
            } else {
                int m = (l + r) / 2;
                build(2 * p, l, m);
                build(2 * p + 1, m, r);
                q[p] = q[2 * p] * q[2 * p + 1];
            }
        };
        build(1, 0, n);
        std::function<void(int, int, int, const Poly &)> work =
            [&](int p, int l, int r, const Poly &num) {
                if (r - l == 1) {
                    if (l < (int)(ans.size())) {
                        ans[l] = num[0];
                    }
                } else {
                    int m = (l + r) / 2;
                    work(2 * p, l, m, num.mulT(q[2 * p + 1]).modxk(m - l));
                    work(2 * p + 1, m, r, num.mulT(q[2 * p]).modxk(r - m));
                }
            };
        work(1, 0, n, mulT(q[1].inv(n)));
        return ans;
    }
};

signed main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    vector<Z> inv(2e5 + 1);
    inv[1] = 1;
    for (int i = 2; i <= 2e5; i++) {
        inv[i] = (P - P / i) * inv[P % i];
    }

    int n, m;
    cin >> n >> m;
    // 球不同，盒子不同，无限制
    cout << power(Z(m), n) << "\n";

    // 球不同，盒子不同，每个盒子至多 1 个球
    cout << (n <= m ? comb.fac(m) / comb.fac(m - n) : 0) << "\n";
    // 球不同，盒子不同，每个盒子至少 1 个球

    Z ans = 0;
    for (int i = 0; i <= m; i++) {
        ans += comb.binom(m, i) * power(Z(m - i), n) * (i & 1 ? -1 : 1);
    }
    cout << ans << "\n";
    // 球不同，盒子同，无限制

    ans = 0;
    vector<Z> f(m + 1), g(m + 1);
    for (int i = 0; i <= m; i++) {
        f[i] = power(Z(i), n) * comb.invfac(i);
        g[i] = comb.invfac(i) * (i & 1 ? -1 : 1);
    }
    Poly F(f), G(g);
    auto strling = F * G;
    for (int i = 1; i <= m; i++) {
        ans += strling[i];
    }
    cout << ans << "\n";
    // 球不同，盒子同，每个盒子至多 1 个球
    cout << (n <= m) << "\n";
    // 球不同，盒子同，每个盒子至少 1 个球
    cout << strling[m] << "\n";

    // 球同，盒子不同，无限制
    cout << comb.binom(n + m - 1, m - 1) << "\n";
    // 球同，盒子不同，每个盒子至多 1 个球
    cout << comb.binom(m, n) << "\n";
    // 球同，盒子不同，每个盒子至少 1 个球
    cout << comb.binom(n - 1, m - 1) << "\n";
    // 球同，盒子同，无限制
    vector<Z> a(n + 1);
    for (int i = 1; i <= m; i++) {
        for (int j = i; j <= n; j += i) {
			a[j] += inv[j / i];
		}
  }
	Poly A(a);
	A = A.exp(n + 1);
	cout << A[n] << "\n";
    // 球同，盒子同，每个盒子至多 1 个球
	cout << (n <= m) << "\n";
    // 球同，盒子同，每个盒子至少 1 个球
	cout << (n >= m ? A[n - m] : 0) << "\n";

    return 0;
}

指数生成函数

指数生成函数的形式是：\(\hat F(x) = \sum_{n}\limits a_n\dfrac{x^n}{n!}\)，这么定义的原因是用来表示不可区分的对象，需要结合运算来理解。我们定义序列 \(a\) 的指数生成函数是 \(F(x) = \sum_{n\ge 0}\limits a_n\dfrac{x^n}{n!}\)，那么 \(F(x)\) 实际也是序列 \(\dfrac{a_n}{n!}\)​ 的普通生成函数，二者没有区别，只是理解问题的方式。

我们通过组合数学来理解一下这样的定理，我们知道多重集排列数，假设共有 \(n\) 个球，其中有 \(k\) 类不同种类的球，第 \(i\) 类球有 \(i\) 个，则这 \(n\) 个球的全排列为 \(\dfrac{n!}{n_1!n_2!...n_k!}\) 。我们尝试在集合中选择去拿取多少个物品的时候，由于每类物品是无序的，我们尝试去掉顺序，因此定义系数为 \(\dfrac{a_i}{i!}\)，而最终生成函数的乘积中，\(x^n\) 前面的系数并不是我们的答案，还需要乘上 \(n!\)。

基本运算

指数生成函数相加和普通生成函数是一致的。

重点是指数生成函数的乘法运算，设两个序列 \(a,b\) 分别对应的指数生成函数是 \(\hat F(x),\hat G(x)\)，那么有： \[ \begin{align*} \hat F(x) \hat G(x) &= \sum_{i \ge 0} a_i \frac{x^i}{i!} \sum_{j \ge 0} b_j \frac{x^j}{j!} \\ &= \sum_{n \ge 0} x^n \sum_{i=0}^{n} a_i b_{n-i} \frac{1}{i!(n-i)!} \\ &= \sum_{n \ge 0} \frac{x^n}{n!} \sum_{i=0}^{n} \binom{n}{i} a_i b_{n-i} \end{align*} \] 那么 \(\hat F(x)\hat G(x)\) 是序列 \(\sum_{i=0}^n\limits\dbinom{n}{i}a_ib_{n-i}\) 的生成函数。

封闭形式

序列 \(\{1,1,1,...\}\) 的指数生成函数是 \(\hat F(x) = \sum_{n\ge 0}\dfrac{x^n}{n!} = e^x\)。

同理，\(\{1,p,p^2,...\}\) 的指数生成函数是 \(\hat F(x) = \sum_{n\ge 0}\dfrac{p^nx^n}{n!} = e^{px}\)。

ABC 422 G

题意

有三个盒子，将 \(n\) 个物品放入这三个盒子中，要求这三个盒子中放的球的数量必须分别是 \(A,B,C\) 的倍数。解决如下两种问题，模 998244353：

  1. 假设物品是不可区分的，求放 \(n\) 个物品的方案数。
  2. 假设物品是是可区分的，求放 \(n\) 个物品的方案数。

题解

第一问是常规生成函数，和付公主的背包一题类似，不过按照那题的方法求多项式 exp 会被卡常。由于盒子数量不多，可以尝试直接系数相加。

第二问考虑到有序性，在单独考虑每个盒子时除以 \(n!\) 去掉顺序在求生成函数的乘法，答案就是 \([x^n]G_1(x)G_2(x)G_3(x)\times n!\)，其中 \([x^i]\) 符号代表从生成函数中取 \(x^i\) 项的系数。

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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
#include <bits/stdc++.h>

#define int long long
using namespace std;

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;

using u128 = unsigned __int128;
using i128 = __int128;
using ld = long double;

template <class T>
constexpr T power(T a, i64 b) {
    T res = 1;
    for (; b; b /= 2, a *= a) {
        if (b % 2) {
            res *= a;
        }
    }
    return res;
}

constexpr i64 mul(i64 a, i64 b, i64 p) {
    i64 res = a * b - i64(1.L * a * b / p) * p;
    res %= p;
    if (res < 0) {
        res += p;
    }
    return res;
}
template <int P>
struct MInt {
    int x;
    constexpr MInt() : x{} {}
    constexpr MInt(i64 x) : x{norm(x % getMod())} {}

    static int Mod;
    constexpr static int getMod() {
        if (P > 0) {
            return P;
        } else {
            return Mod;
        }
    }
    constexpr static void setMod(int Mod_) { Mod = Mod_; }
    constexpr int norm(int x) const {
        if (x < 0) {
            x += getMod();
        }
        if (x >= getMod()) {
            x -= getMod();
        }
        return x;
    }
    constexpr int val() const { return x; }
    explicit constexpr operator int() const { return x; }
    constexpr MInt operator-() const {
        MInt res;
        res.x = norm(getMod() - x);
        return res;
    }
    constexpr MInt inv() const {
        assert(x != 0);
        return power(*this, getMod() - 2);
    }
    constexpr MInt &operator*=(MInt rhs) & {
        x = 1LL * x * rhs.x % getMod();
        return *this;
    }
    constexpr MInt &operator+=(MInt rhs) & {
        x = norm(x + rhs.x);
        return *this;
    }
    constexpr MInt &operator-=(MInt rhs) & {
        x = norm(x - rhs.x);
        return *this;
    }
    constexpr MInt &operator/=(MInt rhs) & { return *this *= rhs.inv(); }
    friend constexpr MInt operator*(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res *= rhs;
        return res;
    }
    friend constexpr MInt operator+(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res += rhs;
        return res;
    }
    friend constexpr MInt operator-(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res -= rhs;
        return res;
    }
    friend constexpr MInt operator/(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res /= rhs;
        return res;
    }
    friend constexpr std::istream &operator>>(std::istream &is, MInt &a) {
        i64 v;
        is >> v;
        a = MInt(v);
        return is;
    }
    friend constexpr std::ostream &operator<<(std::ostream &os, const MInt &a) {
        return os << a.val();
    }
    friend constexpr bool operator==(MInt lhs, MInt rhs) {
        return lhs.val() == rhs.val();
    }
    friend constexpr bool operator!=(MInt lhs, MInt rhs) {
        return lhs.val() != rhs.val();
    }
};
template <>
int MInt<0>::Mod = 998244353;

template <int V, int P>
constexpr MInt<P> CInv = MInt<P>(V).inv();

constexpr int P = 998244353;
using Z = MInt<P>;  // 静态模数，题目给定模数
using D = MInt<0>;  // 动态模数，可以 D::setMod(x) 改变

std::vector<int> rev;
std::vector<Z> roots{0, 1};
void dft(std::vector<Z> &a) {
    int n = a.size();

    if ((int)(rev.size()) != n) {
        int k = __builtin_ctz(n) - 1;
        rev.resize(n);
        for (int i = 0; i < n; i++) {
            rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
        }
    }

    for (int i = 0; i < n; i++) {
        if (rev[i] < i) {
            std::swap(a[i], a[rev[i]]);
        }
    }
    if ((int)(roots.size()) < n) {
        int k = __builtin_ctz(roots.size());
        roots.resize(n);
        while ((1 << k) < n) {
            Z e = power(Z(3), (P - 1) >> (k + 1));
            for (int i = 1 << (k - 1); i < (1 << k); i++) {
                roots[2 * i] = roots[i];
                roots[2 * i + 1] = roots[i] * e;
            }
            k++;
        }
    }
    for (int k = 1; k < n; k *= 2) {
        for (int i = 0; i < n; i += 2 * k) {
            for (int j = 0; j < k; j++) {
                Z u = a[i + j];
                Z v = a[i + j + k] * roots[k + j];
                a[i + j] = u + v;
                a[i + j + k] = u - v;
            }
        }
    }
}
void idft(std::vector<Z> &a) {
    int n = a.size();
    std::reverse(a.begin() + 1, a.end());
    dft(a);
    Z inv = (1 - P) / n;
    for (int i = 0; i < n; i++) {
        a[i] *= inv;
    }
}
struct Poly {
    std::vector<Z> a;
    Poly() {}
    explicit Poly(
        int size, std::function<Z(int)> f = [](int) { return 0; })
        : a(size) {
        for (int i = 0; i < size; i++) {
            a[i] = f(i);
        }
    }
    Poly(const std::vector<Z> &a) : a(a) {}
    Poly(const std::initializer_list<Z> &a) : a(a) {}
    int size() const { return a.size(); }
    void resize(int n) { a.resize(n); }
    Z operator[](int idx) const {
        if (idx < size()) {
            return a[idx];
        } else {
            return 0;
        }
    }
    Z &operator[](int idx) { return a[idx]; }
    Poly mulxk(int k) const {
        auto b = a;
        b.insert(b.begin(), k, 0);
        return Poly(b);
    }
    Poly modxk(int k) const {
        k = std::min(k, size());
        return Poly(std::vector<Z>(a.begin(), a.begin() + k));
    }
    Poly divxk(int k) const {
        if (size() <= k) {
            return Poly();
        }
        return Poly(std::vector<Z>(a.begin() + k, a.end()));
    }
    friend Poly operator+(const Poly &a, const Poly &b) {
        std::vector<Z> res(std::max(a.size(), b.size()));
        for (int i = 0; i < (int)(res.size()); i++) {
            res[i] = a[i] + b[i];
        }
        return Poly(res);
    }
    friend Poly operator-(const Poly &a, const Poly &b) {
        std::vector<Z> res(std::max(a.size(), b.size()));
        for (int i = 0; i < (int)(res.size()); i++) {
            res[i] = a[i] - b[i];
        }
        return Poly(res);
    }
    friend Poly operator-(const Poly &a) {
        std::vector<Z> res(a.size());
        for (int i = 0; i < (int)(res.size()); i++) {
            res[i] = -a[i];
        }
        return Poly(res);
    }
    friend Poly operator*(Poly a, Poly b) {
        if (a.size() == 0 || b.size() == 0) {
            return Poly();
        }
        if (a.size() < b.size()) {
            std::swap(a, b);
        }
        if (b.size() < 128) {
            Poly c(a.size() + b.size() - 1);
            for (int i = 0; i < a.size(); i++) {
                for (int j = 0; j < b.size(); j++) {
                    c[i + j] += a[i] * b[j];
                }
            }
            return c;
        }
        int sz = 1, tot = a.size() + b.size() - 1;
        while (sz < tot) {
            sz *= 2;
        }
        a.a.resize(sz);
        b.a.resize(sz);
        dft(a.a);
        dft(b.a);
        for (int i = 0; i < sz; ++i) {
            a.a[i] = a[i] * b[i];
        }
        idft(a.a);
        a.resize(tot);
        return a;
    }
    friend Poly operator*(Z a, Poly b) {
        for (int i = 0; i < (int)(b.size()); i++) {
            b[i] *= a;
        }
        return b;
    }
    friend Poly operator*(Poly a, Z b) {
        for (int i = 0; i < (int)(a.size()); i++) {
            a[i] *= b;
        }
        return a;
    }
    Poly &operator+=(Poly b) { return (*this) = (*this) + b; }
    Poly &operator-=(Poly b) { return (*this) = (*this) - b; }
    Poly &operator*=(Poly b) { return (*this) = (*this) * b; }
    Poly &operator*=(Z b) { return (*this) = (*this) * b; }
    Poly deriv() const {
        if (a.empty()) {
            return Poly();
        }
        std::vector<Z> res(size() - 1);
        for (int i = 0; i < size() - 1; ++i) {
            res[i] = (i + 1) * a[i + 1];
        }
        return Poly(res);
    }
    Poly integr() const {
        std::vector<Z> res(size() + 1);
        for (int i = 0; i < size(); ++i) {
            res[i + 1] = a[i] / (i + 1);
        }
        return Poly(res);
    }
    Poly inv(int m) const {
        Poly x{a[0].inv()};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (Poly{2} - modxk(k) * x)).modxk(k);
        }
        return x.modxk(m);
    }
    Poly log(int m) const { return (deriv() * inv(m)).integr().modxk(m); }
    Poly exp(int m) const {
        Poly x{1};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (Poly{1} - x.log(k) + modxk(k))).modxk(k);
        }
        return x.modxk(m);
    }
    Poly pow(int k, int m) const {
        int i = 0;
        while (i < size() && a[i].val() == 0) {
            i++;
        }
        if (i == size() || 1LL * i * k >= m) {
            return Poly(std::vector<Z>(m));
        }
        Z v = a[i];
        auto f = divxk(i) * v.inv();
        return (f.log(m - i * k) * k).exp(m - i * k).mulxk(i * k) * power(v, k);
    }
    Poly sqrt(int m) const {
        Poly x{1};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x + (modxk(k) * x.inv(k)).modxk(k)) * ((P + 1) / 2);
        }
        return x.modxk(m);
    }
    Poly mulT(Poly b) const {
        if (b.size() == 0) {
            return Poly();
        }
        int n = b.size();
        std::reverse(b.a.begin(), b.a.end());
        return ((*this) * b).divxk(n - 1);
    }
    std::vector<Z> eval(std::vector<Z> x) const {
        if (size() == 0) {
            return std::vector<Z>(x.size(), 0);
        }
        const int n = std::max((int)(x.size()), size());
        std::vector<Poly> q(4 * n);
        std::vector<Z> ans(x.size());
        x.resize(n);
        std::function<void(int, int, int)> build = [&](int p, int l, int r) {
            if (r - l == 1) {
                q[p] = Poly{1, -x[l]};
            } else {
                int m = (l + r) / 2;
                build(2 * p, l, m);
                build(2 * p + 1, m, r);
                q[p] = q[2 * p] * q[2 * p + 1];
            }
        };
        build(1, 0, n);
        std::function<void(int, int, int, const Poly &)> work =
            [&](int p, int l, int r, const Poly &num) {
                if (r - l == 1) {
                    if (l < (int)(ans.size())) {
                        ans[l] = num[0];
                    }
                } else {
                    int m = (l + r) / 2;
                    work(2 * p, l, m, num.mulT(q[2 * p + 1]).modxk(m - l));
                    work(2 * p + 1, m, r, num.mulT(q[2 * p]).modxk(r - m));
                }
            };
        work(1, 0, n, mulT(q[1].inv(n)));
        return ans;
    }
};

signed main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int a, b, c, n;
    cin >> n >> a >> b >> c;
    vector<Z> inv(n + 1);
    inv[1] = 1;
    for (int i = 2; i <= n; i++) {
        inv[i] = (P - P / i) * inv[P % i];
    }
    vector<Z> f(n + 1);
    for (int i = 1; i <= n; i++) {
        if (i % a == 0) {
            f[i] += inv[i / a];
        }
        if (i % b == 0) {
            f[i] += inv[i / b];
        }
        if (i % c == 0) {
            f[i] += inv[i / c];
        }
    }
    Poly ans1(f);
    ans1 = ans1.exp(n + 1);
    cout << ans1[n] << "\n";

    vector<Z> fac(n + 1), invfac(n + 1);
    fac[0] = 1;
    for (int i = 1; i <= n; i++) {
        fac[i] = fac[i - 1] * i;
    }
    invfac[n] = power(fac[n], P - 2);
    for (int i = n - 1; i >= 0; i--) {
        invfac[i] = invfac[i + 1] * (i + 1);
    }
    vector<Z> g1(n + 1), g2(n + 1), g3(n + 1);
    for (int i = 0; i <= n; i++) {
        if (i % a == 0) {
            g1[i] = invfac[i];
        }
        if (i % b == 0) {
            g2[i] = invfac[i];
        }
        if (i % c == 0) {
            g3[i] = invfac[i];
        }
    }
    Poly G1(g1), G2(g2), G3(g3);
    auto mid = G1 * G2;
    auto G = mid * G3;
    cout << G[n] * fac[n] << "\n";

    return 0;
}