算法模版-数学

发表于 更新于 分类于 阅读次数： 本文字数： 41k
数学模版，慢慢更新。

数学

取模类（模int）

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
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 = 1000000007;
using Z = MInt<P>;  // 静态模数，题目给定模数
using D = MInt<0>;  // 动态模数，可以 D::setMod(x) 改变

线性筛（欧拉筛）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// minp(x) 可以找到 x 大于 1 的最小质因子
std::vector<int> minp, primes;
void sieve(int n) {
    minp.assign(n + 1, 0);
    primes.clear();
    for (int i = 2; i <= n; i++) {
        if (minp[i] == 0) {
            minp[i] = i;
            primes.push_back(i);
        }
        for (auto p : primes) {
            if (i * p > n) {
                break;
            }
            minp[i * p] = p;
            if (p == minp[i]) {
                break;
            }
        }
    }
}

组合数类（搭配取模类）

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
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 类型需要局部初始化

组合数（小范围预处理，逆元+杨辉三角）

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
constexpr int P = 1000000007;
constexpr int L = 10000;

int fac[L + 1], invfac[L + 1];
int sumbinom[L + 1][7];

int binom(int n, int m) {
    if (n < m || m < 0) {
        return 0;
    }
    return 1LL * fac[n] * invfac[m] % P * invfac[n - m] % P;
}

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;
}

int main() {
    fac[0] = 1;
    for (int i = 1; i <= L; i++) {
        fac[i] = 1LL * fac[i - 1] * i % P;
    }
    invfac[L] = power(fac[L], P - 2);
    for (int i = L; i; i--) {
        invfac[i - 1] = 1LL * invfac[i] * i % P;
    }

    sumbinom[0][0] = 1;
    for (int i = 1; i <= L; i++) {
        for (int j = 0; j < 7; j++) {
            sumbinom[i][j] = (sumbinom[i - 1][j] + sumbinom[i - 1][(j + 6) % 7]) % P;
        }
    }
}

矩阵（以斐波那契为例）

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
constexpr int mod = 1e9 + 7;
struct Matrix {
    vector<vector<i64>> a;
    Matrix() { a.assign(3, vector<i64>(3)); }
    Matrix operator*(const Matrix b) const {
        Matrix res;
        for (int i = 1; i <= 2; i++) {
            for (int j = 1; j <= 2; j++) {
                for (int k = 1; k <= 2; k++) {
                    res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j]) % mod;
                }
            }
        }
        return res;
    }
} ans, base;

void qpow(i64 b) {
    while (b) {
        if (b % 2 == 1) {
            ans = ans * base;
        }
        base = base * base;
        b /= 2;
    }
}

void init() {
    base.a[1][1] = base.a[1][2] = base.a[2][1] = 1;
    ans.a[1][1] = ans.a[1][2] = 1;
}

void solve() {
    i64 n;
    cin >> n;
    if (n <= 2) {
        cout << 1 << "\n";
        return;
    }
    init();
    qpow(n - 2);
    cout << ans.a[1][1] % mod << "\n";
}

扩展欧几里得 exgcd，线性同余方程

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
// 返回 gcd(a, b) ax + by = gcd(a, b)
i64 exgcd(i64 a, i64 b, i64 &x, i64 &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    i64 g = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return g;
}

// ax + b = 0 (mod m) 返回最小正整数解，通解的乘数
std::pair<i64, i64> sol(i64 a, i64 b, i64 m) {
    assert(m > 0);
    b *= -1;
  	b = (b % m + m) % m;
    i64 x, y;
    i64 g = exgcd(a, m, x, y);
    if (g < 0) {
        g *= -1;
        x *= -1;
        y *= -1;
    }
    if (b % g != 0) {
        return {-1, -1};
    }
    x = x * (b / g) % (m / g);
    if (x < 0) {
        x += m / g;
    }
    return {x, m / g};
}

std::array<i64, 3> exgcd(i64 a, i64 b) {
    if (!b) {
        return {a, 1, 0};
    }
    auto [g, x, y] = exgcd(b, a % b);
    return {g, y, x - a / b * y};
}

解 xa + yb = n

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// 找 ax + by = n 的正整数解 x, y
bool find_xy(i64 a, i64 b, i64 n, i64 &x, i64 &y) {
    i64 g = exgcd(a, b, x, y);  // ax + by = gcd
    if (n % g != 0) {
        return false;
    }
    i64 k = n / g;
    x *= k;
    y *= k;
    i64 dx = b / g;
    i64 dy = a / g;
    // 通解形式为 x = x + dx * t, y = y - dy * t
    i64 tmin = (x < 0) ? (-x + dx - 1) / dx : 0;
    i64 tmax = (y < 0) ? (-1) : (y / dy);
    if (tmin > tmax) {
        return false;
    }
    x += dx * tmin;
    y -= dy * tmin;
    return true;
}

线性求逆元

1
2
3
4
5
6
7
8
std::vector<int> inv(n + 1);
    inv[1] = 1;
    for (int i = 2; i <= n; i++) {
        inv[i] = 1ll * (p - p / i) * inv[p % i] % p;
    }
    for (int i = 1; i <= n; i++) {
        std::cout << inv[i] << "\n";
    }

中国剩余定理 CRT

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
// 返回 gcd(a, b) ax + by = gcd(a, b)
__int128 exgcd(__int128 a, __int128 b, __int128 &x, __int128 &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    auto g = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return g;
}

__int128 crt(int n, std::vector<i64> &a, std::vector<i64> &b) {
    __int128 M = 1, ans = 0;
    for (int i = 1; i <= n; i++) {
        M = M * b[i];
    }
    for (int i = 1; i <= n; i++) {
        __int128 m = M / b[i];
        __int128 x, y;
        exgcd(m, b[i], x, y);
        ans = (ans + a[i] * m * x) % M;
    }
    ans = (ans % M + M) % M;
    return ans;
};

分解质因数

记住循环最后还有一个 \(n\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int x;
  std::cin >> x;

  for (int i = 2; i * i <= x; i++) {
      if (x % i == 0) {
          while (x % i == 0) {
              std::cout << i << "*";
              x /= i;
          }
      }
  }
  if (x != 1) {
      std::cout << x << "\n";
  }

求所有因数

1
2
3
4
5
6
7
8
9
10
11
12
13
std::vector<int> get_divisors(int n) {
    std::vector<int> res;
    for (int i = 1; i * i <= n; i++) {
        if (n % i == 0) {
            res.push_back(i);
            if (i != n / i) {
                res.push_back(n / i);
            }
        }
    }
    std::sort(res.begin(), res.end());
    return res;
}

欧拉函数

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
// 求解单个数的欧拉函数
int phi(int n) {
  	int res = n;
  	for (int i = 2; i * i <= n; i++) {
    		if(n % i == 0) {
     			 while(n % i == 0) {
        			n /= i;
     			 }
     				res = res / i * (i - 1);
      	}
    }
    if(n > 1) {
      	res = res / n * (n - 1);
    } 
    return res;
}

// 求解全部数的欧拉函数
constexpr int N = 1e7;
constexpr int P = 1000003;

bool isprime[N + 1];
int phi[N + 1];
std::vector<int> primes;

    std::fill(isprime + 2, isprime + N + 1, true);
    phi[1] = 1;
    for (int i = 2; i <= N; i++) {
        if (isprime[i]) {
            primes.push_back(i);
            phi[i] = i - 1;
        }
        for (auto p : primes) {
            if (i * p > N) {
                break;
            }
            isprime[i * p] = false;
            if (i % p == 0) {
                phi[i * p] = phi[i] * p;
                break;
            }
            phi[i * p] = phi[i] * (p - 1);
        }
    }

欧拉定理，扩展欧拉定理

欧拉定理：若 \(\gcd(a,m) = 1, a^{\varphi(m)} = 1\pmod m,a^b = a^{b\pmod{\varphi(m)}}\pmod m\)

扩展欧拉定理：\(a^b = a^{b\pmod {\varphi(m)} + \varphi(m)}\pmod m, b\ge\varphi(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
#include <bits/stdc++.h>

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

// 求解单个数的欧拉函数
int phi(int n) {
    int res = n;
    for (int i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            while (n % i == 0) {
                n /= i;
            }
            res = res / i * (i - 1);
        }
    }
    if (n > 1) {
        res = res / n * (n - 1);
    }
    return res;
}

i64 qpow(i64 a, i64 b, i64 mod) {
    i64 res = 1;
    while (b) {
        if (b & 1) {
            res = res * a % mod;
        }
        a = a * a % mod;
        b /= 2;
    }
    return res;
}

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

    int a, m, b;
    std::cin >> a >> m;
    a %= m;
    int phim = phi(m);
    char ch;
    while ((ch = getchar()) < '0' || ch > '9');
    i64 bm = 0;
    bool flag = false;
    while (bm = bm * 10 + (ch ^ '0'), (ch = getchar()) >= '0' && ch <= '9') {
        if (bm >= phim) {
            flag = true;
            bm %= phim;
        }
    }
    if (bm >= phim) {
        flag = true;
        bm %= phim;
    }
    if (flag) {
        bm += phim;
    }
    std::cout << qpow(a, bm, m) << "\n";

    return 0;
}

数论分块

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// 例：求 \sum_{i=1}^n k \pmod i = nk - \sum_{i=1}^n i(k/i)
i64 ans = 1ll * n * k;
  for (int l = 1, r; l <= n; l = r + 1) {
      int q = k / l;
      if (q == 0) { 
          r = n;
      } else {
          r = std::min(k / q, n);
      }
      ans -= 1ll * (l + r) * (r - l + 1) / 2 * q;
  }
	

// 向上取整的数论分块
for (int l = 1, r; l <= n; l = r + 1) {
      int k = (m + l - 1) / l;
      r = (k == 1) ? n : (m - 1) / (k - 1);
      // 区间 [l, r] \in [1, n] 内，所有 ceil(m / i) 的值都等于 k
      ans = std::min(ans, (k - 1) * l + n - m);
  }

莫比乌斯函数

莫比乌斯函数的定义为： \[ \mu(n) = \left\{ \begin{array}{ll} 1 & \text{if } n = 1 \\ (-1)^k & \text{if } n \text{ 是 } k \text{ 个不同质数的乘积} \\ 0 & \text{if } n \text{ 有平方因子} \end{array} \right. \] 这是一个积性函数（对于 \(\gcd(x,y) = 1,f(x)f(y) = f(xy)\)），具有性质： \[ \sum_{d\mid n}\mu(d) = [n =1] \\\varphi(n) = \sum_{d\mid n}\mu(d)\frac{n}{d} \] 即 \(\mu* I = \epsilon,\mu * id = \varphi\)，其中 \(I\) 代表单位常量函数（\(I(x) = 1\)），\(id\) 代表恒等函数 （\(id(x) = x\)），“*” 代表狄利克雷卷积 \(h(n) = \sum_{d\mid n}f(d)g(\frac{n}{d})\)，记为 \((f*g)(n) = \sum_{d\mid n}f(d)g(\frac{n}{d})\)。

接下来，我们就可以定义莫比乌斯反演了，形象的说就是：莫比乌斯反演 = 用 Möbius 函数做“狄利克雷除法”。 \[ g(n) = \sum_{d\mid n}f(d) \iff f(n) = \sum_{d\mid n}\mu(d)g(\frac{n}{d}) \] 使用条件为：求和必须是整除形式。这样，我们已知 \(f(x)\) 之和 \(g(x)\)，可以通过莫比乌斯反演求出 \(f(x)\)。

例： \[ \begin{aligned} \sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=1] &=\sum_{i=1}^n\sum_{j=1}^m\sum_{d\mid \gcd(i,j)}\mu(d) \\ &=\sum_{i=1}^n\sum_{j=1}^m\sum_{d\mid i,\,d\mid j}\mu(d) \\ &=\sum_{d=1}^n \mu(d)\,\Big\lfloor \tfrac{n}{d}\Big\rfloor \Big\lfloor \tfrac{m}{d}\Big\rfloor \end{aligned} \]

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
// 注意莫比乌斯函数是前缀和的形式
constexpr int N = 1e7;
std::unordered_map<int, Z> fMu;
std::vector<int> minp, primes, phi, mu;
std::vector<i64> sphi;

void sieve(int n) {
    minp.assign(n + 1, 0);
    phi.assign(n + 1, 0);
    sphi.assign(n + 1, 0);
    mu.assign(n + 1, 0);
    primes.clear();
    phi[1] = 1;
    mu[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (minp[i] == 0) {
            minp[i] = i;
            phi[i] = i - 1;
            mu[i] = -1;
            primes.push_back(i);
        }
        for (auto p : primes) {
            if (i * p > n) {
                break;
            }
            minp[i * p] = p;
            if (p == minp[i]) {
                phi[i * p] = phi[i] * p;
                break;
            }
            phi[i * p] = phi[i] * (p - 1);
            mu[i * p] = -mu[i];
        }
    }
    for (int i = 1; i <= n; i++) {
        sphi[i] = sphi[i - 1] + phi[i];
        mu[i] += mu[i - 1];
    }
}

// 快速求 \sum_{i=1}^n \mu(i) 的杜教筛
Z sumMu(int n) {
    if (n <= N) {
        return mu[n];
    }
    if (fMu.count(n)) {
        return fMu[n];
    }
    if (n == 0) {
        return 0;
    }
    Z ans = 1;
    for (int l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        ans -= (r - l + 1) * sumMu(n / l);
    }
    fMu[n] = ans;
    return ans;
}

素数测试，因式分解（Miller-Rabin && Pollard-Rho）

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
i64 mul(i64 a, i64 b, i64 m) {
    return static_cast<__int128>(a) * b % m;
}
i64 power(i64 a, i64 b, i64 m) {
    i64 res = 1 % m;
    for (; b; b >>= 1, a = mul(a, a, m))
        if (b & 1)
            res = mul(res, a, m);
    return res;
}
bool isprime(i64 n) {
    if (n < 2)
        return false;
    static constexpr int A[] = {2, 3, 5, 7, 11, 13, 17, 19, 23};
    int s = __builtin_ctzll(n - 1);
    i64 d = (n - 1) >> s;
    for (auto a : A) {
        if (a == n)
            return true;
        i64 x = power(a, d, n);
        if (x == 1 || x == n - 1)
            continue;
        bool ok = false;
        for (int i = 0; i < s - 1; ++i) {
            x = mul(x, x, n);
            if (x == n - 1) {
                ok = true;
                break;
            }
        }
        if (!ok)
            return false;
    }
    return true;
}
std::vector<i64> factorize(i64 n) {
    std::vector<i64> p;
    std::function<void(i64)> f = [&](i64 n) {
        if (n <= 10000) {
            for (int i = 2; i * i <= n; ++i)
                for (; n % i == 0; n /= i)
                    p.push_back(i);
            if (n > 1)
                p.push_back(n);
            return;
        }
        if (isprime(n)) {
            p.push_back(n);
            return;
        }
        auto g = [&](i64 x) {
            return (mul(x, x, n) + 1) % n;
        };
        i64 x0 = 2;
        while (true) {
            i64 x = x0;
            i64 y = x0;
            i64 d = 1;
            i64 power = 1, lam = 0;
            i64 v = 1;
            while (d == 1) {
                y = g(y);
                ++lam;
                v = mul(v, std::abs(x - y), n);
                if (lam % 127 == 0) {
                    d = std::gcd(v, n);
                    v = 1;
                }
                if (power == lam) {
                    x = y;
                    power *= 2;
                    lam = 0;
                    d = std::gcd(v, n);
                    v = 1;
                }
            }
            if (d != n) {
                f(d);
                f(n / d);
                return;
            }
            ++x0;
        }
    };
    f(n);
    std::sort(p.begin(), p.end());
    return p;
}

O(值域) 预处理的 O(1) gcd

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
struct O1GCD {
    int n, V; // n：值域, V：小范围 gcd 表大小（\sqrt n）
    vector<int> minp, primes;
    vector<array<int, 3>> factors; // 每个数前三个最小质因子
    vector<vector<int>> table; // 小范围 gcd 预处理

    O1GCD(int n) : n(n) {
    	V = sqrt(n) + 1;
        sieve(n);
        pre();
    }

    void sieve(int n) {
        minp.assign(n + 1, 0);
        primes.clear();
        factors.assign(n + 1, {1, 1, 1});

        for (int i = 2; i <= n; i++) {
            if (minp[i] == 0) {
                minp[i] = i;
                primes.push_back(i);
                factors[i] = {1, 1, i};
            }
            for (auto p : primes) {
                if (1LL * i * p > n) {
                	break;
                }
                minp[i * p] = p;
                auto tmp = factors[i];
                tmp[0] *= p;
                sort(tmp.begin(), tmp.end());
                factors[i * p] = tmp;
                if (p == minp[i]) {
                	break;
                }
            }
        }
    }

    void pre() {
        table.assign(V + 1, vector<int>(V + 1));
        for (int i = 0; i <= V; i++) {
            table[i][0] = table[0][i] = i;
        }
        for (int i = 1; i <= V; i++) {
            for (int j = 1; j <= i; j++) {
                table[i][j] = table[j][i] = table[i % j][j];
            }
        }
    }

    int operator()(int a, int b) const {
        int g = 1;
        for (int i = 0; i < 3; i++) {
            int p = factors[a][i];
            if (p == 1) {
            	continue;
            }
            int tmp;
            if (p > V) { 
            	tmp = (b % p == 0 ? p : 1);
            }
            else {
                tmp = table[p][b % p];
            }
            b /= tmp;
            g *= tmp;
        }
        return g;
    }
};

BSGS / exBSGS

注意判不存在是 < 0

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
constexpr int inf = 3e18;
int bsgs(int a, int b, int m, int k = 1) {
	static map<int, int> mp;
	mp.clear();
	int cur = 1, t = sqrt(m) + 1;
	for (int B = 1; B <= t; B++) {
		cur = cur * a % m;
		mp[b * cur % m] = B;
	}
	int now = cur * k % m;
	for (int A = 1; A <= t; A++) {
		auto it = mp.find(now);
		if (it != mp.end()) {
			return A * t - it->second;
		}
		now = now * cur % m;
	}
	return -inf;
}

int exbsgs(int a, int b, int m, int k = 1) {
	int A = a % m, B = b % m, M = m;
	if (b == 1) {
		return 0;
	}
	int cur = 1 % m;
	for (int i = 0;; i++) {
		if (cur == B) {
			return i;
		}
		cur = cur * A % M;
		int d = gcd(a, m);
		if (b % d) {
			return -inf;
		}
		if (d == 1) {
			return bsgs(a, b, m, k * a % m) + i + 1;
		}
		k = k * a / d % m;
		b /= d;
		m /= d;
	}
}

矩阵

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
struct Matrix {
    int n, m;
    vector<vector<int>> a;
    Matrix() : n(0), m(0) {}
    Matrix(int n) : n(n), m(n) { a.assign(n, vector<int>(n)); }
    Matrix(int n, int m) : n(n), m(m) { a.assign(n, vector<int>(m)); }
    vector<int> &operator[](int row) { return a[row]; }
    const vector<int> &operator[](int row) const { return a[row]; }
    Matrix operator*(const Matrix &b) const {
        assert(m == b.n);
        Matrix res(n, b.m);
        for (int k = 0; k < m; k++) {
            for (int i = 0; i < n; i++) {
                for (int j = 0; j < b.m; j++) {
                    res[i][j] = res[i][j] + a[i][k] * b[k][j];
                }
            }
        }
        return res;
    }
    Matrix &operator*=(const Matrix &b) {
        *this = *this * b;
        return *this;
    }
    Matrix power(int b) {
        assert(n == m);
        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;
    }
    int det(int P) {
        assert(n == m);
        int ans = 1, tmp = 1;
        auto b = a;
        for (int i = 0; i < n; i++) {
            for (int j = i + 1; j < n; j++) {
                while (b[i][i]) {
                    int d = b[j][i] / b[i][i];
                    for (int k = i; k < n; k++) {
                        b[j][k] = ((b[j][k] - d * b[i][k]) % P + P) % P;
                    }
                    swap(b[i], b[j]);
                    tmp = -tmp;
                }
                swap(b[i], b[j]);
                tmp = -tmp;
            }
        }
        for (int i = 0; i < n; i++) {
            ans = ans * b[i][i] % P;
        }
        ans = ((ans * tmp) % P + P) % P;
        return ans;
    }
};

求逆矩阵（取模）

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
#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 = 1e9 + 7;
template <class T>
constexpr T power(T a, int b) {
    T res = 1;
    for (; b; b /= 2, a *= a) {
        if (b % 2) {
            res *= a;
        }
    }
    return res;
}

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

bool invMatrix(vector<vector<int>> &a, int P) {
    int n = a.size();
    for (int i = 0, r; i < n; i++) {
        r = i;
        for (int j = i + 1; j < n; j++)
            if (a[j][i] > a[r][i]) {
                r = j;
            }
        if (r != i) {
            swap(a[i], a[r]);
        }
        if (!a[i][i]) {
            return false;
        }

        int inv = power(a[i][i], P - 2, P);
        for (int j = 0; j < 2 * n; j++) {
            a[i][j] = (a[i][j] * inv % P);
        }
        for (int k = 0; k < n; k++) {
            if (k == i || a[k][i] == 0) {
                continue;
            }
            int f = a[k][i];
            for (int j = i; j < n * 2; j++) {
                a[k][j] = (a[k][j] - f * a[i][j]) % P;
                if (a[k][j] < 0) {
                    a[k][j] += P;
                }
            }
        }
    }
    return true;
}

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

    int n;
    cin >> n;
    vector<vector<int>> a(n, vector<int>(2 * n));
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            cin >> a[i][j];
        }
        a[i][i + n] = 1;
    }
    if (!invMatrix(a, P)) {
        cout << "No Solution\n";
        return 0;
    }
    for (int i = 0; i < n; i++) {
        for (int j = n; j < 2 * n; j++) {
            cout << a[i][j] << " \n"[j == 2 * n - 1];
        }
    }

    return 0;
}

高斯消元

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
constexpr ld eps = 1e-8;
constexpr ld inf = 2e30;
// ori 为增广矩阵，返回解向量
// -1：无解 0：无穷多解 1:唯一解
pair<vector<ld>, int> gauss(vector<vector<ld>> &ori) {
    int n = ori.size();
    vector<vector<ld>> a(n, vector<ld>(n + 1));
    vector<ld> ans(n, inf);
    for (int i = 0; i < n; i++) {
        for (int j = 0; j <= n; j++) {
            a[i][j] = 1.0 * ori[i][j];
        }
    }
    int r = 0;
    for (int c = 0; c < n; c++) {
        for (int i = r + 1; i < n; i++) {
            if (fabs(a[i][c]) > fabs(a[r][c])) {
                swap(a[r], a[i]);
            }
        }
        if (fabs(a[r][c]) < eps) {
            continue;
        }
        for (int i = n; i >= c; i--) {
            a[r][i] /= a[r][c];
        }
        for (int i = r + 1; i < n; i++) {
            for (int j = n; j >= c; j--) {
                a[i][j] -= a[i][c] * a[r][j];
            }
        }
        r++;
    }
    if (r < n) {
        for (int i = r; i < n; i++) {
            if (fabs(a[i][n]) > eps) {
                return {ans, -1};
            }
        }
        return {ans, 0};
    }
    for (int j = n - 1; j > 0; j--) {
        for (int i = 0; i < j; i++) {
            a[i][n] -= a[j][n] * a[i][j];
            a[i][j] -= a[j][j] * a[i][j];
        }
    }
    for (int i = 0; i < n; i++) {
        ans[i] = a[i][n];
    }
    return {ans, 1};
}

constexpr int N = 2010;

高斯消元（取模）

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
template <class T>
constexpr T power(T a, int b) {
    T res = 1;
    for (; b; b /= 2, a *= a) {
        if (b % 2) {
            res *= a;
        }
    }
    return res;
}

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

// -1：无解 0：无穷多解 1:唯一解
// 注意读入矩阵时要给系数取模
pair<vector<int>, int> gauss(vector<vector<int>> &a, int P) {
    int n = a.size();
    vector<int> ans(n, -1);
    int r = 0;
    for (int c = 0; c < n; c++) {
        for (int i = r + 1; i < n; i++) {
            if (a[i][c] > a[r][c]) {
                swap(a[r], a[i]);
            }
        }
        if (a[r][c] == 0) {
            continue;
        }
        int inv = power(a[r][c], P - 2, P);
        inv = (inv % P + P) % P;
        for (int i = n; i >= c; i--) {
            a[r][i] *= inv;
            a[r][i] %= P;
        }
        for (int i = r + 1; i < n; i++) {
            for (int j = n; j >= c; j--) {
                a[i][j] -= a[i][c] * a[r][j] % P;
                a[i][j] = (a[i][j] % P + P) % P;
            }
        }
        r++;
    }
    if (r < n) {
        for (int i = r; i < n; i++) {
            if (a[i][n]) {
                return {ans, -1};
            }
        }
        return {ans, 0};
    }
    for (int j = n - 1; j > 0; j--) {
        for (int i = 0; i < j; i++) {
            a[i][n] -= a[j][n] * a[i][j] % P;
            a[i][n] = (a[i][n] % P + P) % P;
            a[i][j] -= a[j][j] * a[i][j] % P;
            a[i][j] = (a[i][j] % P + P) % P;
        }
    }
    for (int i = 0; i < n; i++) {
        ans[i] = a[i][n];
    }
    return {ans, 1};
}

高斯消元（01矩阵 异或）

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
// get xor rank
int gauss(vector<vector<int>> &ori) {
    int n = ori.size();
    int m = ori[0].size();
    vector<bitset<N>> a(n);
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < m; j++) {
            if (ori[i][j] == 1) {
                a[i][j] = 1;
            }
        }
    }
    int r = 0;
    for (int c = 0; c < m && r < n; c++) {
        int p = r;
        while (p < n && !a[p][c]) {
            p++;
        }
        if (p == n) {
            continue;
        }
        swap(a[r], a[p]);
        for (int i = r + 1; i < n; i++) {
            if (a[i][c]) {
                a[i] ^= a[r];
            }
        }
        r++;
    }
    return n - r;
}

组合数学相关

\(\begin{pmatrix}n\\m\end{pmatrix}\) 读作 “\(n\) 选 \(m\)”，\(\begin{pmatrix}n\\m\end{pmatrix} =\displaystyle \frac{n(n-1)...(n - m + 1)}{m(m - 1)...1} = \frac{n!}{m!(n - m)!}\)。

可以定义广义二项式系数：\(\begin{pmatrix}r\\k\end{pmatrix}=\begin{cases}\frac{r(r-1)...(r-k+1)}{k(k-1)...1} = \frac{r^{\underline{k}}}{k!},&整数k\ge 0\\0, &整数 k < 0\end{cases}\)，这里 \(r\) 为非负整数时，就退化为普通二项式系数，\(r\) 为负数或非整数时，仍然可以按照公式计算。

  • 预处理阶乘及逆元：预处理 \(O(n)\)，查询 \(O(1)\)
  • 当 \(n\) 很小，或者题目没有给定模数，用加法公式预处理 （可能用到高精度），预处理 \(O(n^2)\)，查询 \(O(1)\)。
  • 当 \(m\) 很小，通过 \(\begin{pmatrix}n\\m\end{pmatrix} = \frac{n^{\underline{m}}}{m!}\) 暴力计算（可以预处理分母），预处理 \(O(m)\)，查询 \(O(m)\)
  • 当模数 \(P\) 很小，根据 Lucas 定理，\(\begin{pmatrix}n\\m\end{pmatrix}\equiv \begin{pmatrix}n\mod p\\m\mod p\end{pmatrix}\begin{pmatrix}\lfloor\frac{n}{p}\rfloor\\\lfloor\frac{m}{p}\rfloor\end{pmatrix}\pmod P\)，预处理 \(O(P)\)，查询 \(O(1)\)。

组合恒等变换

  • 对称性

\(\begin{pmatrix}n\\m\end{pmatrix} = \begin{pmatrix}n\\n - m\end{pmatrix}\)，整数 \(n\ge 0\)，\(m\) 是整数。

  • 二项式定理

当 \(n\le 0\) 且 \(n\) 为整数时，是最常见的情况：\((x+y)^n = \sum_{k=0}^n\begin{pmatrix}n\\k\end{pmatrix}x^ky^{n-k}\)。

当 \(|\frac{y}{x}|\le 1\) 时，\(n\) 可以是任意实数或复数，这是广义二项式展开的条件，即 Newton 二项式公式：\((x+y)^n=\sum_{k=0}^\infty\begin{pmatrix}n\\k\end{pmatrix}x^{n-k}y^k\)，其中 \(\begin{pmatrix}n\\k\end{pmatrix}\) 也是广义二项式系数。

  • 吸收恒等式

\(k\begin{pmatrix}r\\k\end{pmatrix} = r\begin{pmatrix}r-1\\k-1\end{pmatrix}\)，\(k\) 是整数。将组合数展开为阶乘即可推导。

推论：\((r - k)\begin{pmatrix}r\\k\end{pmatrix} = r\begin{pmatrix}r - 1\\k\end{pmatrix}\)，\(k\) 是整数。因为 \(r\begin{pmatrix}r-1\\k\end{pmatrix}=(k+1)\begin{pmatrix}r\\k+1\end{pmatrix} = (r-k)\begin{pmatrix}r\\k\end{pmatrix}\)。

  • 加法公式

\(\begin{pmatrix}r\\k\end{pmatrix} = \begin{pmatrix}r - 1\\k - 1\end{pmatrix} + \begin{pmatrix}r - 1\\k\end{pmatrix}\)，\(k\) 是整数。通过杨辉三角的意义即可推导。

  • 平行求和法

\(\sum_{k\le n}\begin{pmatrix}r+k\\k\end{pmatrix} = \begin{pmatrix}r\\0\end{pmatrix} + \begin{pmatrix}r+1\\1\end{pmatrix} + ... + \begin{pmatrix}r+n\\n\end{pmatrix} = \begin{pmatrix}r + n + 1\\n\end{pmatrix}\)​，\(n\) 是整数。注意没有 \(r\) 的限制，可以推广到广义二项式系数。

  • 上指标求和法

\(\sum_{0\le k\le n}\begin{pmatrix}k\\m\end{pmatrix} = \begin{pmatrix}0\\m\end{pmatrix} + \begin{pmatrix}1\\m\end{pmatrix} + ... + \begin{pmatrix}n\\m\end{pmatrix} = \begin{pmatrix}n + 1\\m + 1\end{pmatrix}\)，整数 \(m,n\ge 0\)。

  • 上指标反转

\(\begin{pmatrix}r\\k\end{pmatrix} = (-1)^k\begin{pmatrix}k - r - 1\\k\end{pmatrix}\)，\(k\)​ 是整数。

根据广义二项式系数，\(\begin{pmatrix}k - r - 1\\k\end{pmatrix} = \displaystyle\frac{(k-r-1)(k-r-2)...(-r)}{k!} = \frac{(-1)^kr(r-1)...(r-k+1)}{k!}\)。

  • 三项式版恒等式

\(\begin{pmatrix}r\\m\end{pmatrix} \begin{pmatrix}m\\k\end{pmatrix}= \begin{pmatrix}r\\k\end{pmatrix}\begin{pmatrix}r-k\\m-k\end{pmatrix}\)​。

考虑组合意义，即左边的意义是先选 \(m\) 个元素，再从这 \(m\) 个元素里选 \(k\) 个。右边的意义是先选 \(k\) 个元素，再从剩下的 \(r-k\) 个元素中选 \(m-k\) 个。

  • 三项式系数

\(\begin{pmatrix}a+b+c\\a,b,c\end{pmatrix} = \displaystyle\frac{(a+b+c)!}{a!b!c!} = \begin{pmatrix}a+b+c\\b+c\end{pmatrix}\begin{pmatrix}b+c\\c\end{pmatrix}\)​。

表示从 \(a+b+c\) 个元素中选出 \(a\) 个放在第一类，\(b\) 个放在第二类，\(c\) 个放在第三类的方案数。

同理，可以推广到多项式系数。

  • 范德蒙徳卷积

\(\sum_{k}\begin{pmatrix}r\\k\end{pmatrix}\begin{pmatrix}s\\n-k\end{pmatrix} = \begin{pmatrix}r+s\\n\end{pmatrix}\)。考虑组合意义，从两个集合中一共选 \(n\) 个的方案数，就是枚举在第一个集合中选 \(k\) 个的方案数，乘上在第二个集合中选 \(n - k\)​ 个的方案数。

特殊的数

错位排列

定义 \(D_n\) 为有多少个长度为 \(n\) 的排列 \(P\) 满足 \(P_i \ne i\)。

有递推公式：\(D_n = (n - 1)(D_{n - 1} + D_{n - 2})\)。含义是：在 \(D_{n-1}\) 的基础上加第 \(n\) 个数，那么把 \(P_n\) 跟前面任意一个位置交换即可，或者前面 \(n - 1\) 个数有一个位置满足 \(P_i = i\)，那么新增一个数后，让 \(P_i\) 和 \(P_n\) 互换位置即可。

初值 \(D_1 = 0, D_2 = 1\)。

卡特兰数

卡特兰数对应序列为：\(H_0 = 1, H_1 = 1, H_2 = 2, H_3 = 5,H_4 = 14, H_5=42,H_6=132,...\)。

应用于以下问题：

  1. \(2n\) 个人排成一行进入剧场，入场费 5 元，其中 \(n\) 个人只有一张 5 元，另外 \(n\) 人只有一张 10 元，剧院无其他钞票，求有多少种方法使得只要有 10 元钞票，售票处就有 5 元的钞票找零。
  2. 有一个大小为 \(n\times n\) 的方格图，左下角 \((0,0)\)，右下角 \((n,n)\)。从左下角开始每次都只能向右或者向上一个单位，不走到对角线 \(y = x\) 上方（但可以触碰）的情况下到达右上角的方案数。
  3. 在圆上选择 \(2n\) 个点，将这些点成对连接起来使得得到的 \(n\) 条线段不相交的方案数。
  4. 对角线不相交的情况下，将一个凸多边形分成三角形区域的方案数。
  5. 一个栈的进栈序列为 \(1,2,3,...,n\)，求有多少种不同的出栈序列。
  6. \(n\) 个节点可以构造多少个不同的二叉树。
  7. 由 \(n\) 个 +1 和 \(n\) 个 -1 组成的 \(2n\) 个数 \(a_1,a_2,...,a_{2n}\)，其部分和满足 \(a_1+a_2+...+a_k\ge 0,k=1,2,3,...,2n\)​，有多少个满足条件的数列。

以上问题的答案均是卡特兰数 \(H_i\)。 卡特兰数有几个常见公式，可以对应上述情况：

\(H_n = \begin{cases}\sum_{i=1}^nH_{i-1}H_{n-i}&n\ge 2,n\in N_+\\1&n=0,1\end{cases}\)

\(H_n = \displaystyle\frac{H_{n-1}(4n-2)}{n+1}\)

\(H_n = \begin{pmatrix}2n\\n\end{pmatrix} - \begin{pmatrix}2n\\n - 1\end{pmatrix}\)

例如，我们处理问题 3，我们将圆顺时针编号 \(1,2,...,2n\)，选择一个固定点 1，他必须与一个偶数编号的点配对，我们将这条线段画出后，圆被分成上下两部分，也就是拆成两个子问题。若弦一侧含有 \(2k\) 个点，就能组成 \(k\) 对，枚举 \(k = 0,1,2,...,n-1\) 得到递推 \(C_n = \sum_{k=0}^{n-1}C_{k}C_{n-1-k}\)，符合卡特兰数的形式，初值 \(C_0 = 1\)​。

我们还可以给卡特兰数抽象出这样一种数学模型：对于一个只含有两种操作（设为操作1，操作2）的问题，问题的解是由这两种操作构成的排列，并且操作1的总数等于操作2的总数，任取前 \(k\) 项（\(k\ge1\)），若满足操作 1 的个数大于等于操作 2 的个数才是问题的一种解，则问题解的个数就是卡特兰数。

我们定义非降路径是只能向上或向右走的路径。

  1. 从 \((0,0)\) 走到 \((m,n)\) 的非降路径数等于 \(m\) 个 \(x\) 和 \(n\) 个 \(y\) 的排列数，即 \(\begin{pmatrix}n+m\\n\end{pmatrix}\)
  2. 从 \((0,0)\) 走到 \((n,n)\) 的除端点外不接触直线 \(y = x\) 的非降路径数。我们先考虑 \(y = x\) 下方的路径，此时第一步一定是走到 \((1,0)\)，最后一步一定是从 \((n,n-1)\) 走到 \((n,n)\)，那么等价于从 \((1,0)\) 到 \((n,n-1)\) 全程不触碰 \(y=x\) 的路径数，这里可以采用折线法，无限制的话方案就是 \(\begin{pmatrix}2n - 2\\n - 1\end{pmatrix}\)，再考虑触碰到 \(y = x\)，在将一个交点之后的路径对称，等价于走到 \((n-1,n)\) 的方案数，即 \(\begin{pmatrix}2n-2\\n - 2\end{pmatrix}\)。因此 \(y = x\) 下方的方案数就是 \(\begin{pmatrix}2n - 2\\n - 1\end{pmatrix} - \begin{pmatrix}2n-2\\n - 2\end{pmatrix}\)。而 \(y = x\) 上方完全对称，总方案数就是 \(2\begin{pmatrix}2n - 2\\n - 1\end{pmatrix} - 2\begin{pmatrix}2n-2\\n - 2\end{pmatrix}\)。
  3. 从 \((0,0)\) 到 \((n,n)\) 的除端点外不穿过直线 \(y = x\) 的方案数。和上一题类似，先考虑 \(y = x\) 下方的路径，此时“不穿过”等价于不接触 \(y = x + 1\)，同理可以推出总方案数就是 \(2\begin{pmatrix}2n \\n \end{pmatrix} - 2\begin{pmatrix}2n\\n - 1\end{pmatrix} = \displaystyle \frac{2}{n+1}\begin{pmatrix}2n \\n \end{pmatrix}\)​。

第二类斯特林数

第二类斯特林数 \(\left\{\begin{matrix}n\\k\end{matrix} \right\}\)，也可记作 \(S(n,k)\)，表示将 \(n\) 个两两不同的元素，划分为 \(k\) 个互不区分的非空子集的方案数。

递推式：\(\left\{\begin{matrix}n\\k\end{matrix} \right\} = \left\{\begin{matrix}n - 1\\k - 1\end{matrix} \right\} + k\left\{\begin{matrix}n- 1\\k\end{matrix} \right\}\)，边界是 \(\left\{\begin{matrix}n\\0\end{matrix} \right\} = [n = 0]\)。

通过组合意义来证明：插入一个新元素时，有两种方案：将新元素单独放入一个子集，有 \(\left\{\begin{matrix}n-1\\k-1\end{matrix} \right\}\) 种方案；将新元素放入一个现有的非空子集，有 \(k\left\{\begin{matrix}n - 1\\k\end{matrix} \right\}\) 种方案。根据加法原理，二者相加即可得到递推式。

根据这个递推式，可以 \(O(n^2)\) 预处理。

同时，第二类斯特林数具有通项公式：\(\left\{\begin{matrix}n\\m\end{matrix} \right\} = \sum_{i=0}^m\displaystyle\frac{(-1)^{m-i}i^n}{i!(m-i)!}\)。预处理阶乘后，可以 \(O(n\log n)\) 求出单个第二类斯特林数。

注意到通项公式中含有 \(i^n\)，这启示我们当题目所求和式带有 \(i^n\) 项时，可以考虑将其展开为第二类斯特林数，即 \(i^n = \sum_{k=0}^n\left\{\begin{matrix}n\\k\end{matrix} \right\}i^{\underline{k}}\)​。

容斥及反演

容斥原理

就是将一般的“\(x\) 人喜欢语文，\(y\) 人喜欢数学，\(z\) 人喜欢英语，求至少喜欢一门课的人数” 推广一下。

容斥原理常用于集合的计数问题，对于两个集合的函数 \(f(S),g(S)\)，若 \(f(S) = \sum_{T\subseteq S}g(T)\)，则 \(g(S) = \sum_{T\subseteq S}(-1)^{|S| - |T|}f(T)\)。还有一种推论是若 \(f(S) = \sum_{S\subseteq T}g(T)\)，则 \(g(S) = \sum_{S\subseteq T}(-1)^{|T|-|S|}g(T)\)​。

反演

简单理解一下反演，如果能够从 \(f\) 推出 \(g\)，则称之为正演，那么反过来通过 \(g\) 推出 \(f\) 就是反演。

二项式反演

如果有 \(g_n = \sum_{i=0}^n\begin{pmatrix}n\\i\end{pmatrix}f_i\)，那么 \(f_n = \sum_{i=0}^n\begin{pmatrix}n\\i\end{pmatrix}(-1)^{n-i}g_i\)。

上面这对式子就是二项式反演。还有变式：若 \(F(n) = \sum_{i=n}^m\begin{pmatrix}i\\n\end{pmatrix}G(i)\)，则 $ G(n) = _{i=n}m(-1){i-n} \[\begin{pmatrix}i\\n\end{pmatrix}\]

F(i)$

Min-Max 反演

对于满足全序关系并且其中元素满足可加减性的序列 \(\{x_i\}\)，设其长度为 \(n\)，并设 \(S = \{1,2,3,...,n\}\)，则有：

\(\displaystyle \max_{i\in S}x_i = \sum_{T\subseteq S}(-1)^{|T| - 1}\min_{j\in T}x_j\)

\(\displaystyle \min_{i\in S}x_i = \sum_{T\subseteq S}(-1)^{|T| - 1}\max_{j\in T}x_j\)

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
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(3, (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;
}

任意模数 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;
}

多项式类（Z）

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
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;
    }
};