拓扑排序#

1
#include <bits/stdc++.h>
2
using namespace std;
3

4
/**
5
 * @brief 拓扑排序封装
6
 * 适用场景：有向无环图（DAG）的先后顺序处理，检测环。
7
 */
8
struct TopoSort {
9
    int n;
10
    vector<int> in;          // 入度数组
11
    vector<vector<int>> g;   // 邻接表
12

13
    TopoSort(int n) : n(n), in(n + 1, 0), g(n + 1) {}
14

15
    void add_edge(int u, int v) {
16
        g[u].push_back(v);
17
        in[v]++;
18
    }
19

20
    vector<int> solve() {
21
        queue<int> q;
22
        vector<int> res;
23
        for (int i = 1; i <= n; i++)
24
            if (in[i] == 0) q.push(i);
25

26
        while (!q.empty()) {
27
            int u = q.front(); q.pop();
28
            res.push_back(u);
29
            for (int v : g[u]) {
30
                if (--in[v] == 0) q.push(v);
31
            }
32
        }
33
        // 如果结果集大小不足 n，说明图中存在环
34
        return (res.size() == n) ? res : vector<int>();
35
    }
36
};
37

38
int main() {
39
    ios::sync_with_stdio(false); cin.tie(nullptr);
40
    int n, m;
41
    cin >> n >> m;
42
    TopoSort ts(n);
43
    for (int i = 0; i < m; i++) {
44
        int u, v; cin >> u >> v;
45
        ts.add_edge(u, v);
46
    }
47
    vector<int> ans = ts.solve();
48
    if (ans.empty()) cout << "IMPOSSIBLE\n";
49
    else {
50
        for (int i = 0; i < ans.size(); i++)
51
            cout << ans[i] << (i == ans.size() - 1 ? "" : " ");
52
        cout << "\n";
53
    }
54
    return 0;
55
}

树#

树的直径#

1
#include <bits/stdc++.h>
2
using namespace std;
3

4
/**
5
 * @brief 树的直径 (两次 DFS 法)
6
 * 注意：仅适用于边权为正的树。
7
 */
8
struct TreeDiameter {
9
    int n, far_node;
10
    vector<vector<int>> g;
11
    vector<int> dist;
12

13
    TreeDiameter(int n) : n(n), g(n + 1), dist(n + 1), far_node(0) {}
14

15
    void add_edge(int u, int v) {
16
        g[u].push_back(v);
17
        g[v].push_back(u);
18
    }
19

20
    void dfs(int u, int fa) {
21
        if (dist[u] > dist[far_node]) far_node = u;
22
        for (int v : g[u]) {
23
            if (v == fa) continue;
24
            dist[v] = dist[u] + 1;
25
            dfs(v, u);
26
        }
27
    }
28

29
    int solve() {
30
        far_node = 1;
31
        fill(dist.begin(), dist.end(), 0);
32
        dfs(1, 0); // 第一次 DFS 找到距离起点最远的点
33

34
        int start_node = far_node;
35
        fill(dist.begin(), dist.end(), 0);
36
        dfs(start_node, 0); // 第二次 DFS 从最远点出发找到直径
37
        return dist[far_node];
38
    }
39
};
40

41
int main() {
42
    int n; cin >> n;
43
    TreeDiameter td(n);
44
    for (int i = 1; i < n; i++) {
45
        int u, v; cin >> u >> v;
46
        td.add_edge(u, v);
47
    }
48
    cout << td.solve() << endl;
49
    return 0;
50
}

树的重心#

1
#include <bits/stdc++.h>
2
using namespace std;
3

4
/**
5
 * @brief 树的重心
6
 * 性质：删除重心后，最大连通分量的节点数最小，且该值 <= n/2。
7
 */
8
struct TreeCentroid {
9
    int n;
10
    vector<vector<int>> g;
11
    vector<int> siz, wei, centroids;
12

13
    TreeCentroid(int n) : n(n), g(n + 1), siz(n + 1, 0), wei(n + 1, 0) {}
14

15
    void add_edge(int u, int v) {
16
        g[u].push_back(v);
17
        g[v].push_back(u);
18
    }
19

20
    void dfs(int u, int fa) {
21
        siz[u] = 1;
22
        wei[u] = 0;
23
        for (int v : g[u]) {
24
            if (v == fa) continue;
25
            dfs(v, u);
26
            siz[u] += siz[v];
27
            wei[u] = max(wei[u], siz[v]); // 最大的子树大小
28
        }
29
        wei[u] = max(wei[u], n - siz[u]); // 向上连接的部分
30
        if (wei[u] <= n / 2) centroids.push_back(u);
31
    }
32

33
    int get_one() {
34
        dfs(1, 0);
35
        return centroids.empty() ? -1 : centroids[0];
36
    }
37
};
38

39
int main() {
40
    int n; cin >> n;
41
    TreeCentroid tc(n);
42
    for (int i = 1; i < n; i++) {
43
        int u, v; cin >> u >> v;
44
        tc.add_edge(u, v);
45
    }
46
    cout << tc.get_one() << endl;
47
    return 0;
48
}

倍增LCA#

1
#include <bits/stdc++.h>
2
using namespace std;
3

4
/**
5
 * @brief 倍增法求最近公共祖先 (LCA)
6
 * 预处理 O(n log n)，单次查询 O(log n)。
7
 */
8
struct BinaryLCA {
9
    int n, LOG;
10
    vector<vector<int>> g, anc;
11
    vector<int> depth;
12

13
    BinaryLCA(int n) : n(n), LOG(20), g(n + 1), anc(n + 1, vector<int>(21, 0)), depth(n + 1, 0) {}
14

15
    void add_edge(int u, int v) {
16
        g[u].push_back(v);
17
        g[v].push_back(u);
18
    }
19

20
    void dfs(int u, int fa, int d) {
21
        depth[u] = d;
22
        anc[u][0] = fa;
23
        for (int i = 1; i <= LOG; i++) {
24
            anc[u][i] = anc[anc[u][i - 1]][i - 1];
25
        }
26
        for (int v : g[u]) {
27
            if (v != fa) dfs(v, u, d + 1);
28
        }
29
    }
30

31
    int get_lca(int a, int b) {
32
        if (depth[a] < depth[b]) swap(a, b);
33
        for (int i = LOG; i >= 0; i--) {
34
            if (depth[a] - (1 << i) >= depth[b]) a = anc[a][i];
35
        }
36
        if (a == b) return a;
37
        for (int i = LOG; i >= 0; i--) {
38
            if (anc[a][i] != anc[b][i]) {
39
                a = anc[a][i];
40
                b = anc[b][i];
41
            }
42
        }
43
        return anc[a][0];
44
    }
45
};
46

47
int main() {
48
    ios::sync_with_stdio(false); cin.tie(nullptr);
49
    int n, q; cin >> n >> q;
50
    BinaryLCA lca(n);
51
    for (int i = 2; i <= n; i++) {
52
        int p; cin >> p; // 假设给出的是父节点
53
        lca.add_edge(p, i);
54
    }
55
    lca.dfs(1, 0, 1);
56
    while (q--) {
57
        int a, b; cin >> a >> b;
58
        cout << lca.get_lca(a, b) << "\n";
59
    }
60
    return 0;
61
}

树上启发式合并 (DSU on Tree)#

1
#include <bits/stdc++.h>
2
using namespace std;
3

4
/**
5
 * @brief 树上启发式合并 (DSU on Tree)
6
 * 复杂度 O(n log n)。常用于统计子树内颜色数量等问题。
7
 */
8
struct DSUonTree {
9
    int n;
10
    vector<vector<int>> g;
11
    vector<int> color, siz, hson, ans, cnt;
12
    int dist; // 当前统计的颜色种类数
13

14
    DSUonTree(int n) : n(n), g(n + 1), color(n + 1), siz(n + 1, 0),
15
                       hson(n + 1, 0), ans(n + 1), cnt(1000005, 0), dist(0) {}
16

17
    void add_edge(int u, int v) {
18
        g[u].push_back(v);
19
        g[v].push_back(u);
20
    }
21

22
    void dfs_siz(int u, int fa) {
23
        siz[u] = 1;
24
        for (int v : g[u]) {
25
            if (v == fa) continue;
26
            dfs_siz(v, u);
27
            siz[u] += siz[v];
28
            if (siz[v] > siz[hson[u]]) hson[u] = v;
29
        }
30
    }
31

32
    void update(int u, int fa, int val) {
33
        if (val == 1) {
34
            if (cnt[color[u]] == 0) dist++;
35
            cnt[color[u]]++;
36
        } else {
37
            cnt[color[u]]--;
38
            if (cnt[color[u]] == 0) dist--;
39
        }
40
        for (int v : g[u]) {
41
            if (v != fa) update(v, u, val);
42
        }
43
    }
44

45
    // 核心逻辑：先处理轻儿子（不保留结果），再处理重儿子（保留结果），最后合并
46
    void dfs_solve(int u, int fa, bool keep, bool is_hson = false) {
47
        for (int v : g[u]) {
48
            if (v != fa && v != hson[u]) dfs_solve(v, u, false);
49
        }
50
        if (hson[u]) dfs_solve(hson[u], u, true, true);
51

52
        for (int v : g[u]) {
53
            if (v != fa && v != hson[u]) update(v, u, 1);
54
        }
55

56
        if (cnt[color[u]] == 0) dist++;
57
        cnt[color[u]]++;
58

59
        ans[u] = dist;
60

61
        if (!keep) update(u, fa, -1);
62
    }
63
};
64

65
int main() {
66
    int n; cin >> n;
67
    DSUonTree dsu(n);
68
    for (int i = 1; i <= n; i++) cin >> dsu.color[i];
69
    for (int i = 1; i < n; i++) {
70
        int u, v; cin >> u >> v;
71
        dsu.add_edge(u, v);
72
    }
73
    dsu.dfs_siz(1, 0);
74
    dsu.dfs_solve(1, 0, false);
75
    for (int i = 1; i <= n; i++) cout << dsu.ans[i] << " ";
76
    return 0;
77
}

Tarjan#

SCC (强连通分量)#

1
#include <bits/stdc++.h>
2
using namespace std;
3

4
/**
5
 * @brief Tarjan 算法求强连通分量 (SCC)
6
 * 适用于有向图。
7
 */
8
struct TarjanSCC {
9
    int n, timer, scc_cnt;
10
    vector<vector<int>> g, sccs;
11
    vector<int> dfn, low, scc_id;
12
    vector<bool> instk;
13
    stack<int> stk;
14

15
    TarjanSCC(int n) : n(n), timer(0), scc_cnt(0), g(n + 1),
16
                       dfn(n + 1, 0), low(n + 1, 0), scc_id(n + 1, 0), instk(n + 1, false) {}
17

18
    void add_edge(int u, int v) { g[u].push_back(v); }
19

20
    void tarjan(int u) {
21
        dfn[u] = low[u] = ++timer;
22
        stk.push(u);
23
        instk[u] = true;
24

25
        for (int v : g[u]) {
26
            if (!dfn[v]) {
27
                tarjan(v);
28
                low[u] = min(low[u], low[v]);
29
            } else if (instk[v]) {
30
                low[u] = min(low[u], dfn[v]);
31
            }
32
        }
33

34
        if (low[u] == dfn[u]) {
35
            scc_cnt++;
36
            vector<int> current_scc;
37
            while (true) {
38
                int v = stk.top(); stk.pop();
39
                instk[v] = false;
40
                scc_id[v] = scc_cnt;
41
                current_scc.push_back(v);
42
                if (u == v) break;
43
            }
44
            sccs.push_back(current_scc);
45
        }
46
    }
47

48
    void solve() {
49
        for (int i = 1; i <= n; i++) if (!dfn[i]) tarjan(i);
50
        reverse(sccs.begin(), sccs.end()); // 拓扑序反转
51
    }
52
};
53

54
int main() {
55
    int n, m; cin >> n >> m;
56
    TarjanSCC ts(n);
57
    for (int i = 0; i < m; i++) {
58
        int u, v; cin >> u >> v;
59
        ts.add_edge(u, v);
60
    }
61
    ts.solve();
62
    cout << ts.scc_cnt << endl;
63
    for (auto &scc : ts.sccs) {
64
        cout << scc.size() << " ";
65
        for (int v : scc) cout << v << " ";
66
        cout << "\n";
67
    }
68
    return 0;
69
}

割边 (Bridges)#

1
#include <bits/stdc++.h>
2
using namespace std;
3

4
/**
5
 * @brief Tarjan 算法求割边（桥）
6
 * 适用于无向图。
7
 */
8
struct TarjanBridge {
9
    int n, m, timer;
10
    struct Edge { int to, id; };
11
    vector<vector<Edge>> g;
12
    vector<int> dfn, low;
13
    vector<bool> is_bridge;
14

15
    TarjanBridge(int n, int m) : n(n), m(m), timer(0), g(n + 1),
16
                                 dfn(n + 1, 0), low(n + 1, 0), is_bridge(m + 1, false) {}
17

18
    void add_edge(int u, int v, int id) {
19
        g[u].push_back({v, id});
20
        g[v].push_back({u, id});
21
    }
22

23
    void tarjan(int u, int from_edge) {
24
        dfn[u] = low[u] = ++timer;
25
        for (auto &e : g[u]) {
26
            if (!dfn[e.to]) {
27
                tarjan(e.to, e.id);
28
                low[u] = min(low[u], low[e.to]);
29
                if (low[e.to] > dfn[u]) is_bridge[e.id] = true;
30
            } else if (e.id != from_edge) {
31
                low[u] = min(low[u], dfn[e.to]);
32
            }
33
        }
34
    }
35
};
36

37
int main() {
38
    int n, m; cin >> n >> m;
39
    TarjanBridge tb(n, m);
40
    for (int i = 1; i <= m; i++) {
41
        int u, v; cin >> u >> v;
42
        tb.add_edge(u, v, i);
43
    }
44
    for (int i = 1; i <= n; i++) if (!tb.dfn[i]) tb.tarjan(i, -1);
45
    // 后续可根据 is_bridge 统计或 DFS 求双连通分量
46
    return 0;
47
}

割点 (Cut Points)#

1
#include <bits/stdc++.h>
2
using namespace std;
3

4
/**
5
 * @brief Tarjan 算法求割点
6
 * 适用于无向图。
7
 */
8
struct TarjanCutPoint {
9
    int n, timer;
10
    vector<vector<int>> g;
11
    vector<int> dfn, low;
12
    vector<bool> is_cut;
13

14
    TarjanCutPoint(int n) : n(n), timer(0), g(n + 1),
15
                            dfn(n + 1, 0), low(n + 1, 0), is_cut(n + 1, false) {}
16

17
    void add_edge(int u, int v) {
18
        g[u].push_back(v);
19
        g[v].push_back(u);
20
    }
21

22
    void tarjan(int u, int fa) {
23
        dfn[u] = low[u] = ++timer;
24
        int child = 0;
25
        for (int v : g[u]) {
26
            if (!dfn[v]) {
27
                child++;
28
                tarjan(v, u);
29
                low[u] = min(low[u], low[v]);
30
                if (fa != -1 && low[v] >= dfn[u]) is_cut[u] = true;
31
            } else if (v != fa) {
32
                low[u] = min(low[u], dfn[v]);
33
            }
34
        }
35
        if (fa == -1 && child >= 2) is_cut[u] = true;
36
    }
37

38
    int count() {
39
        int cnt = 0;
40
        for (int i = 1; i <= n; i++) if (is_cut[i]) cnt++;
41
        return cnt;
42
    }
43
};
44

45
int main() {
46
    int n, m; cin >> n >> m;
47
    TarjanCutPoint tcp(n);
48
    for (int i = 0; i < m; i++) {
49
        int u, v; cin >> u >> v;
50
        tcp.add_edge(u, v);
51
    }
52
    for (int i = 1; i <= n; i++) if (!tcp.dfn[i]) tcp.tarjan(i, -1);
53
    cout << tcp.count() << endl;
54
    return 0;
55
}

网络流#

最大流/最小割 (Dinic)#

1
#include <bits/stdc++.h>
2
using namespace std;
3
typedef long long ll;
4

5
/**
6
 * @brief Dinic 算法求最大流
7
 * 时间复杂度：一般图 O(V^2 E)，二分图 O(E sqrt(V))。
8
 */
9
struct Dinic {
10
    struct Edge {
11
        int to;
12
        ll cap;
13
        int rev;
14
    };
15
    int n;
16
    vector<vector<Edge>> g;
17
    vector<int> level, cur;
18

19
    Dinic(int n) : n(n), g(n + 1), level(n + 1), cur(n + 1) {}
20

21
    void add_edge(int u, int v, ll cap) {
22
        g[u].push_back({v, cap, (int)g[v].size()});
23
        g[v].push_back({u, 0, (int)g[u].size() - 1});
24
    }
25

26
    bool bfs(int s, int t) {
27
        fill(level.begin(), level.end(), -1);
28
        level[s] = 0;
29
        queue<int> q;
30
        q.push(s);
31
        while (!q.empty()) {
32
            int u = q.front(); q.pop();
33
            for (auto &e : g[u]) {
34
                if (e.cap > 0 && level[e.to] == -1) {
35
                    level[e.to] = level[u] + 1;
36
                    q.push(e.to);
37
                }
38
            }
39
        }
40
        return level[t] != -1;
41
    }
42

43
    ll dfs(int u, int t, ll f) {
44
        if (u == t || f == 0) return f;
45
        for (int &i = cur[u]; i < g[u].size(); i++) {
46
            Edge &e = g[u][i];
47
            if (e.cap > 0 && level[e.to] == level[u] + 1) {
48
                ll d = dfs(e.to, t, min(f, e.cap));
49
                if (d > 0) {
50
                    e.cap -= d;
51
                    g[e.to][e.rev].cap += d;
52
                    return d;
53
                }
54
            }
55
        }
56
        return 0;
57
    }
58

59
    ll max_flow(int s, int t) {
60
        ll flow = 0;
61
        while (bfs(s, t)) {
62
            fill(cur.begin(), cur.end(), 0);
63
            while (ll f = dfs(s, t, LLONG_MAX)) flow += f;
64
        }
65
        return flow;
66
    }
67

68
    // 求最小割边：BFS 找到从 S 可达的点，跨越可达与不可达点集的边即为割边
69
    vector<pair<int, int>> get_min_cut_edges(const vector<pair<int, int>>& original_edges) {
70
        vector<int> vis(n + 1, 0);
71
        queue<int> q; q.push(1); vis[1] = 1;
72
        while(!q.empty()){
73
            int u = q.front(); q.pop();
74
            for(auto &e : g[u]) if(e.cap > 0 && !vis[e.to]) { vis[e.to] = 1; q.push(e.to); }
75
        }
76
        vector<pair<int, int>> res;
77
        for(auto &p : original_edges) if(vis[p.first] ^ vis[p.second]) res.push_back(p);
78
        return res;
79
    }
80
};
81

82
int main() {
83
    int n, m; cin >> n >> m;
84
    Dinic dinic(n);
85
    for (int i = 0; i < m; i++) {
86
        int u, v; ll c; cin >> u >> v >> c;
87
        dinic.add_edge(u, v, c);
88
    }
89
    cout << dinic.max_flow(1, n) << endl;
90
    return 0;
91
}

数据结构模板库#

树状数组#

1
template <typename T>
2
class FenWick
3
{
4
private:
5
    vector<T> BIT;
6
    T lowbit(T x);
7
    int n;
8

9
public:
10
    FenWick(T n) : n(n)
11
    {
12
        this->BIT.resize(n, 0);
13
    }
14
    FenWick(vector<T> &ori) : FenWick(ori.size())
15
    {
16
        this->BIT = ori;
17
    }
18
    void set(int poi, T value);
19
    T query(int poi);
20
    T query(int l, int r);
21
};
22
template <typename T>
23
T FenWick<T>::lowbit(T x)
24
{
25
    return x & (-x);
26
}
27
template <typename T>
28
void FenWick<T>::set(int poi, T value)
29
{
30
    for (; poi < n; poi += (lowbit(poi)))
31
        (this->BIT[poi]) += value;
32
    return;
33
}
34
template <typename T>
35
T FenWick<T>::query(int poi)
36
{
37
    T ans = 0;
38
    for (; poi > 0; poi -= lowbit(poi))
39
    {
40
        ans += (this->BIT)[poi];
41
    }
42
    return ans;
43
}
44
template <typename T>
45
T FenWick<T>::query(int l, int r)
46
{
47
    return (this->query(r) - this->query(l - 1));
48
}

并查集#

基础并查集#

1
//这是并查集的一个极简实现
2
struct Disjoint_set_unoin{
3
    vector<int>fa;
4
    Disjoint_set_unoin(int n){fa.resize(n);iota(fa.begin(),fa.end(),0);};
5
    int find(int n){return n==fa[n]?n:fa[n] = find(fa[n]);}
6
    void unite(int x,int y){fa[find(x)] = find(y);}
7
};

可回滚并查集#

1
struct Disjoint_set_union_with_rollback
2
{
3
    //由于要保留路径信息，因此不能再使用路径压缩
4
    vector<int>fa;
5
    vector<int>siz;
6
    stack<pair<int,int>> undost;
7
    Disjoint_set_union_with_rollback(int n){
8
        fa.resize(n);
9
        iota(fa.begin(),fa.end(),0);
10
        siz.resize(n,1);
11
    }
12
    int find(int x){
13
        return x==fa[x]?x:find(fa[x]);
14
    }
15
    void unite(int x,int y){
16
        int fx = find(x) ,fy = find(y);
17
        if(fx==fy){
18
            //虽然没有实际操作，但是为了确保timer一致性，push空操作
19
            undost.push({-1,-1});
20
            return;
21
        }
22
        if(siz[fx]<=siz[fy]){
23
            siz[fy] += siz[fx];
24
            fa[fx] = fy;
25
            undost.push({fx,fy});
26
        }else
27
        {
28
            siz[fx] += siz[fy];
29
            fa[fy] = fx;
30
            undost.push({fy,fx});
31
        }
32
        return;
33
    }
34
    void undo(){
35
        auto [x,y] = undost.top();
36
        undost.pop();
37
        if(x==-1) return;
38
        else{
39
            siz[y] -= siz[x];
40
            fa[x] = x;
41
            return;
42
        }
43
    }
44
};

带删除并查集#

1
struct Disjoint_set_union_with_delete
2
{
3
    vector<int> fa;
4
    vector<int> siz;
5
    vector<int> id;
6
    int tot;
7

8
    // n 为初始逻辑节点数，q 为预计的最大删除次数（用于预留物理节点空间）
9
    Disjoint_set_union_with_delete(int n, int q){
10
        int max_nodes = n + q + 10;
11
        fa.resize(max_nodes);
12
        iota(fa.begin(), fa.end(), 0);
13
        siz.resize(max_nodes, 1);
14

15
        id.resize(n);
16
        iota(id.begin(), id.end(), 0);
17
        tot = n - 1;
18
    }
19

20
    int find(int x){
21
        return x == fa[x] ? x : fa[x] = find(fa[x]);
22
    }
23

24
    void unite(int x, int y){
25
        int fx = find(id[x]), fy = find(id[y]);
26
        if(fx == fy){
27
            return;
28
        }
29
        if(siz[fx] <= siz[fy]){
30
            siz[fy] += siz[fx];
31
            fa[fx] = fy;
32
        }else
33
        {
34
            siz[fx] += siz[fy];
35
            fa[fy] = fx;
36
        }
37
        return;
38
    }
39

40
    void del(int x){
41
        int fx = find(id[x]);
42
        // 从原集合的 size 中扣除自己 (视具体题目需求，通常删除节点意味着原集合大小-1)
43
        siz[fx]--;
44

45
        // 为逻辑节点 x 分配一个新的独立物理节点
46
        tot++;
47
        id[x] = tot;
48
        fa[tot] = tot;
49
        siz[tot] = 1;
50
        return;
51
    }
52

53
    // 获取逻辑节点 x 所在集合的大小
54
    int get_size(int x){
55
        return siz[find(id[x])];
56
    }
57
};

线段树#

基础线段树#

1
struct D{
2
    int val = 0;//按需配置
3
};
4
D operator+(const D& left, const D& right){return {left.val+right.val};} //按需重载
5
struct segment_tree{
6
    vector<D>data;
7
    D e; //单位元
8
    segment_tree(int n){
9
        data.resize(4*n+10);
10
    }
11
    void modify(int nl,int nr ,int tp ,D val ,int p){
12
        if(nl==nr){
13
            data[p] = val;
14
            return;
15
        }
16
        else{
17
            int mid = (nl+nr)>>1;
18
            if(tp<=mid) modify(nl,mid,tp,val,2*p);
19
            if(tp>mid) modify(mid+1,nr,tp,val,2*p+1);
20
            //push up
21
            data[p] = data[2*p] + data[2*p+1];
22
            return;
23
        }
24
    }
25
    D query(int nl,int nr,int tl,int tr,int p) const{
26
        if(nl>=tl&&nr<=tr) return data[p];
27
        else{
28
            int mid = (nl+nr)>>1;
29
            D lt = e ,rt = e;
30
            if(tl<=mid) lt = query(nl,mid,tl,tr,2*p);
31
            if(tr>mid) rt = query(mid+1,nr,tl,tr,2*p+1);
32
            return lt+rt;
33
        }
34
    }
35
};

懒标记线段树#

1
//懒标记线段树的实现较为简单，但本身由于TAG的存在，非常灵活
2
//特别是涉及多种操作的时候，要仔细想清楚TAG的与DATA的合并逻辑与优先级
3
//基础逻辑就一句话，在递归之前，必须下传TAG，在结束之后，必须根据子区间信息push up
4

5
//这里给出一个区间赋值/加和/查询 操作的线段树实现，实际实现时对照修改即可
6
#include<bits/stdc++.h>
7
using namespace std;
8
using ll = long long;
9
struct D{
10
    ll val = 0;
11
};
12
D operator+(const D&l_op,const D& r_op){return {l_op.val+r_op.val};};
13
struct TAG{
14
    ll add = 0;
15
    ll set = -1;
16
};
17

18
TAG operator+(const TAG&l_op,const TAG&r_op){ //r_op是优先级更高的tag
19
    if(r_op.set!=-1){
20
        return {r_op.add,r_op.set};
21
    }
22
    else{
23
        return {l_op.add+ r_op.add,l_op.set};
24
    }
25
}
26
TAG& operator+=(TAG& l_op,const TAG& r_op){
27
    l_op = l_op + r_op;
28
    return l_op;
29
}
30
struct Segment_tree_with_lazy_tag{
31
    vector<D> data;
32
    vector<TAG> lz_tag;
33
    D e; //单位元
34
    Segment_tree_with_lazy_tag(int n){
35
        data.resize(4*n);
36
        lz_tag.resize(4*n);
37
    }
38
    //计算区间信息与tag合并的函数
39
    D merge(const D& data,const TAG& tag,ll length){
40
        const auto[add,set] = tag;
41
        if(set!=-1){
42
            return {(set+add)*length};
43
        }
44
        else{
45
            return {data.val + add*length};
46
        }
47
    }
48
    //区间编辑
49
    void modify(int nl,int nr,int tl,int tr,TAG val,int p){
50
        if(nl>=tl&&nr<=tr){
51
            lz_tag[p] += val;
52
            return;
53
        }
54
        else{
55
            //递归前push down
56
            lz_tag[2*p] += lz_tag[p];
57
            lz_tag[2*p+1] += lz_tag[p];
58
            lz_tag[p] = TAG();
59

60
            int mid = nl +((nr-nl)>>1);
61
            if(tl<=mid) modify (nl,mid,tl,tr,val,2*p);
62
            if(tr>mid) modify (mid+1,nr,tl,tr,val,2*p+1);
63

64
            //push up
65
            data[p] = merge(data[2*p],lz_tag[2*p],(mid - nl+1)) + merge(data[2*p+1],lz_tag[2*p+1],(nr-mid));
66
            return ;
67
        }
68
    }
69
    //区间查询
70
    D query(int nl,int nr ,int tl,int tr,int p){
71
        if(nl>=tl&&nr<=tr){
72
            return merge(data[p],lz_tag[p],nr-nl+1);
73
        }
74
        else{
75
            //递归前push down
76
            lz_tag[2*p] += lz_tag[p];
77
            lz_tag[2*p+1] += lz_tag[p];
78
            lz_tag[p] = TAG();
79

80
            D lt = e;
81
            D rt = e;
82
            int mid = nl +((nr-nl)>>1);
83
            if(tl<=mid) lt = query(nl,mid,tl,tr,2*p);
84
            if(tr>mid) rt = query(mid+1,nr,tl,tr,2*p+1);
85

86
            //push up
87
            data[p] = merge(data[2*p],lz_tag[2*p],(mid - nl+1)) + merge(data[2*p+1],lz_tag[2*p+1],(nr-mid));
88
            return lt+rt;
89
        }
90
    }
91
};

动态开点线段树#

1
struct D{
2
    int val = 0; //按需配置
3
};
4
D operator+(const D& left, const D& right){return {left.val+right.val};} //按需重载
5
struct segment_tree_with_dynamic_points{
6
    vector<D> data;
7
    vector<pair<int,int>> nxt;
8
    D e; //单位元
9
    segment_tree_with_dynamic_points(){
10
        data.resize(2);
11
        nxt.resize(2);
12
    }
13
    int newp(){
14
        data.emplace_back();
15
        nxt.emplace_back();
16
        return data.size() -1;
17
    }
18
    void modify(int nl,int nr ,int tp,D val,int p){
19
        if(nl==nr) data[p] = val;
20
        else{
21
            auto [lp,rp] = nxt[p]; //注意此处不能引用，避免扩容后悬垂
22
            int mid = (nl+nr)>>1;
23
            if(tp<=mid){
24
                if(lp==0) lp = newp();
25
                modify(nl,mid,tp,val,lp);
26
            }
27
            if(tp>mid){
28
                if(rp==0) rp =newp();
29
                modify(mid+1,nr,tp,val,rp);
30
            }
31
            nxt[p].first = lp;
32
            nxt[p].second = rp;
33
            //push up;
34
            D lt =e ,rt = e;
35
            if(lp!=0) lt = data[lp];
36
            if(rp!=0) rt = data[rp];
37
            data[p] = lt + rt;
38

39
            return;
40
        }
41
    }
42
    D query(int nl,int nr ,int tl ,int tr ,int p) const{
43
        if(nl>=tl&&nr<=tr) return data[p];
44
        else{
45
            const auto&[lp,rp] = nxt[p];
46
            int mid = (nl+nr)>>1;
47
            D lt =e ,rt = e;
48
            if(lp!=0&&tl<=mid) lt = query(nl,mid,tl,tr,lp);
49
            if(rp!=0&&tr>mid) rt = query(mid+1,nr,tl,tr,rp);
50
            return lt+rt;
51
        }
52
    }
53
};

可持久化线段树（主席树）#

1
struct D{
2
    int val = 0;
3
};
4
D operator+ (const D& left ,const D& right){return {left.val+right.val};}
5
struct persistent_segment_tree{
6
    int cnt_v;
7
    vector<int> h;
8
    vector<D> data;
9
    D e;    //单位元
10
    vector<pair<int,int>>nxt;
11
    persistent_segment_tree(D e) : e(e){
12
        data.resize(2,e);
13
        nxt.resize(2);
14
        h.resize(1);
15
        h[0] = 1;
16
    }
17
    int newp(){
18
        data.push_back(e);
19
        nxt.push_back({});
20
        return data.size() - 1;
21
    }
22
    //每次修改都创建新节点，不修改的部分复用原来的节点，因此递归新旧同时进行
23
    void modify(int nl,int nr ,int tp ,D val,int p ,int old_p){
24
        if(nl==nr){
25
            data[p] = val;
26
            return;
27
        }else{
28
            const auto [o_lp,o_rp] = nxt[old_p]; //不允许使用引用，以免悬垂指针
29
            ll mid  = (nl+nr)>>1;
30
            if(tp<=mid){
31
                nxt[p].second = o_rp;
32
                nxt[p].first = newp();
33
                modify(nl,mid,tp,val,nxt[p].first,o_lp);
34
            }
35
            else{
36
                nxt[p].first = o_lp;
37
                nxt[p].second = newp();
38
                modify(mid+1,nr,tp,val,nxt[p].second,o_rp);
39
            }
40
            //push up;
41
            const auto[lp,rp] = nxt[p];
42
            data[p] = data[lp] + data[rp];
43
        }
44
    }
45
    //查询和正常线段树完全一致
46
    D query(int nl,int nr ,int tl ,int tr ,int p) const{
47
        if(nl>=tl&&nr<=tr){
48
            return data[p];
49
        }
50
        else{
51
            const auto [lp,rp] = nxt[p];
52
            ll mid = (nl+nr)>>1;
53
            D lt = e ,rt = e;
54
            if(tl<=mid&&lp)
55
                lt = query(nl,mid,tl,tr,lp);
56
            if(tr>mid&&rp)
57
                rt = query(mid+1,nr,tl,tr,rp);
58
            return lt+rt;
59
        }
60
    }
61
    void update(int tk,int nl,int nr ,int tp ,D val){
62
        int old_head = h[tk];
63
        int new_head = newp();
64
        modify(nl ,nr ,tp,val,new_head,old_head);
65
        h[tk] = new_head;
66
    }
67
    int backup(int tk){
68
        h.push_back(h[tk]);
69
        return h.size()-1;
70
    }
71
};

标记永久化线段树#

1
//在线段树区间查询的时候，除了打懒标记并下传，还有一种写法
2
//那就是我们将标记"保留"在区间中，查询的时候累计路径上的所有标记并计算结果
3
//坏处很明显，这要求标记满足结合律和交换律，这显然是阿贝尔群的性质。
4
//好处是什么？由于这样实现不要求下传标记，因此可以很方便的进行 "持久化"
5

6
//下给出一个普通线段树的标记永久化区间和实现
7
//支持区间增加及查询
8
struct D
9
{
10
    int val = 0;
11
};
12
D operator+(const D&l_op,const D&r_op) {return {l_op.val+r_op.val};}
13
struct TAG
14
{
15
    int add = 0;
16
};
17
TAG operator+(const TAG&l_op,const TAG&r_op){ //r_op是优先级更高的tag
18
   return {l_op.add+r_op.add};
19
}
20
TAG& operator+=(TAG& l_op,const TAG& r_op){
21
    l_op = l_op + r_op;
22
    return l_op;
23
}
24
struct Segment_tree_with_persistent_marks
25
{
26
    vector<D>data;
27
    vector<TAG>tags;
28
    Segment_tree_with_persistent_marks(int n){
29
        data.resize(4*n+10);
30
        tags.resize(4*n+10);
31
    }
32
    D e;    //单位元
33
    D merge(const D&data, const TAG&tag,int length){
34
        return {data.val+tag.add*length};
35
    }
36
    void modify(int nl,int nr,int tl,int tr,TAG add,int p){
37
        //每次途经节点，累加实际受影响的交集长度贡献
38
        int overlap = min(nr, tr) - max(nl, tl) + 1;
39
        data[p] = merge(data[p],add,overlap);
40
        if(nl>=tl&&nr<=tr){
41

42
            tags[p] += add; //标记留在此处，不再递归
43
            return;
44
        }
45
        else{
46
            int mid = (nl+((nr-nl)>>1));
47
            if(tl<=mid) modify(nl,mid,tl,tr,add,2*p);
48
            if(tr>mid) modify(mid+1,nr,tl,tr,add,2*p+1);
49
            return;
50
        }
51
    }
52

53
    //查询时，只计算父节点的标记和！
54
    D query(int nl,int nr,int tl,int tr,int p){
55
        if(nl>=tl&&nr<=tr){
56
            return data[p];
57
        }
58
        else{
59
            int mid = (nl + ((nr-nl)>>1));
60
            D lt = e;
61
            D rt = e;
62
            if(tl<=mid) lt = query(nl,mid,tl,tr,2*p);
63
            if(tr>mid) rt = query(mid+1,nr,tl,tr,2*p+1);
64

65
            int overlap = min(nr, tr) - max(nl, tl) + 1;
66
            return merge(lt + rt, tags[p], overlap);
67
        }
68
    }
69
};

字典树#

基础字典树#

1
struct D{
2
    int nxt[26] ={0};
3
    int cnt = 0;
4
};
5
struct Trie{
6
    vector<D>node;
7
    Trie(){node.resize(2);}
8
    int newp(){node.push_back({}); return node.size()-1;}
9
    //只给出insert的基本实现，简单来说，不能使用递归实现，查询时不断挪动指针向下
10
    //其余操作的查询操作均以此为基础，不再赘述
11
    void insert(const string& s){
12
        int n = s.size();
13
        int cur = 1;
14
        for(int i = 0;i<n;i++){
15
            auto val = s[i];
16
            if(!node[cur].nxt[val-'a']) {
17
                node[cur].nxt[val-'a'] = newp();
18
            }
19
            cur = node[cur].nxt[val-'a'] ;
20
        }
21
        node[cur].cnt++;
22
        return;
23
    }
24
};

01字典树#

1
//01字典树一般用来解决位运算相关的问题
2
//如果原题可以按位贪心，可以从01字典树的角度思考
3
struct D{
4
    int nxt[2] = {0};
5
};
6
struct Trie_01{
7
    vector<D> node;
8
    Trie_01(){
9
        node.resize(2);
10
    }
11
    int newp(){
12
        node.push_back({});
13
        return node.size() - 1;
14
    }
15
    //插入一个数
16
    void insert(int num){
17
        int cur = 1;
18
        for(int i = 30;i>=0;i--){
19
            int bit = (num>>i) &1;
20
            if(node[cur].nxt[bit]==0){
21
                node[cur].nxt[bit]  = newp();
22
            }
23
            cur = node[cur].nxt[bit];
24
        }
25
    }
26
    //返回异或值最大的结果
27
    int query(int num) const{
28
        int cur = 1;
29
        int ret = 0;
30
        for(int i = 30;i>=0;i--){
31
            int bit = (num>>i)&1;
32
            if(node[cur].nxt[!bit] ==0){
33
                cur = node[cur].nxt[bit];
34
            }
35
            else{
36
                cur = node[cur].nxt[!bit];
37
                ret += (1<<i);
38
            }
39
        }
40
        return ret;
41
    }
42
};

可持久化01字典树#

1
//为01字典树添加了持久化操作
2
//思路与主席树类似，只要操作，就修改附近所有节点，然后复用之前节点
3
#include<bits/stdc++.h>
4
using namespace std;
5
struct D{
6
    int nxt[2] = {0};
7
    int cnt = 0;
8
};
9
struct Trie_01{
10
    vector<D> node;
11
    vector<int> h;
12
    Trie_01(){
13
        node.resize(2);
14
        h.resize(1);
15
        h[0] = 1;
16
    }
17
    int newp(){
18
        node.push_back({});
19
        return node.size() - 1;
20
    }
21
    //给版本k,插入一个数
22
    void insert(int tk,int num){
23
        int cur = newp();
24
        int o_cur = h[tk];
25
        h[tk] = cur;    //cur 更换新的头节点
26
        for(int i = 30;i>=0;i--){
27
            int bit = (num>>i)&1;
28
            node[cur].nxt[bit] = newp();
29
            node[cur].cnt = node[o_cur].cnt + 1;
30

31
            node[cur].nxt[!bit] = node[o_cur].nxt[!bit];
32

33
            cur = node[cur].nxt[bit];
34
            o_cur = node[o_cur].nxt[bit];
35
        }
36
        node[cur].cnt = node[o_cur].cnt + 1;
37
    }
38
    //复制版本k，并返回新的版本编号
39
    int copy(int tk){
40
        h.push_back(h[tk]);
41
        return h.size()-1;
42
    }
43
    //查询版本k的最大值
44
    int query(int tk,int num) const{
45
        int cur = h[tk];
46
        int ret = 0;
47
        for(int i = 30;i>=0;i--){
48
            int bit = (num>>i)&1;
49
            if(node[cur].nxt[!bit] ==0){
50
                cur = node[cur].nxt[bit];
51
            }
52
            else{
53
                cur = node[cur].nxt[!bit];
54
                ret += (1<<i);
55
            }
56
        }
57
        return ret;
58
    }
59
    //查询版本l 与版本r之间的最大值  [l,r] 双闭区间
60
    //往往这个更有效果，因为可以处理区间
61
    int query(int tl,int tr,int num) const{
62
        if(tl==0) return query(tr,num);
63
        int curl = h[tl - 1];
64
        int curr = h[tr];
65
        int ret = 0;
66
        for(int i = 30;i>=0;i--){
67
            int bit = (num>>i)&1;
68
            int cnt = node[node[curr].nxt[!bit]].cnt - node[node[curl].nxt[!bit]].cnt;
69
            if(cnt>0){
70
                //共同进入 !bit 分支
71
                ret += (1<<i);
72
                curl = node[curl].nxt[!bit];
73
                curr = node[curr].nxt[!bit];
74
            }
75
            else{
76
                //共同进入 bit分支
77
                curl = node[curl].nxt[bit];
78
                curr = node[curr].nxt[bit];
79
            }
80
        }
81
        return ret;
82
    }
83
};

平衡树#

pbds实现#

1
//pbds自带平衡树，用这个非常方便，一下给出示例（Generated by Gemini 3.1 pro)
2
#include <iostream>
3
// 1. 引入 pb_ds 的核心头文件
4
#include <ext/pb_ds/assoc_container.hpp>
5
#include <ext/pb_ds/tree_policy.hpp>
6

7
using namespace std;
8
using namespace __gnu_pbds; // 必须引入此命名空间
9

10
// ==========================================
11
// 1. 普通平衡树（Set：元素不可重复）模板
12
// ==========================================
13
typedef tree<
14
    int,                                  // Key 类型
15
    null_type,                            // Value 类型（null_type 表示作为 set 使用）
16
    less<int>,                            // 比较函数（less 为升序，greater 为降序）
17
    rb_tree_tag,                          // 底层结构（红黑树）
18
    tree_order_statistics_node_update     // 开启名次树功能（必须带上这个）
19
> ordered_set;
20

21
// ==========================================
22
// 2. 多重集平衡树（Multiset：元素可重复）模板
23
// 竞赛避坑：直接用 less_equal 会导致 erase 行为异常
24
// 解决方案：通过插入 pair<值, 唯一ID> 强行变回单集
25
// ==========================================
26
typedef pair<int, int> pii;
27
typedef tree<
28
    pii,
29
    null_type,
30
    less<pii>,
31
    rb_tree_tag,
32
    tree_order_statistics_node_update
33
> ordered_multiset;
34

35
int main() {
36
    // --------------------------------------
37
    // 基础单集演示
38
    // --------------------------------------
39
    ordered_set tr;
40

41
    // 插入元素 (O(log N))
42
    tr.insert(10);
43
    tr.insert(30);
44
    tr.insert(20);
45
    tr.insert(50);
46
    tr.insert(40);
47

48
    // 【核心操作 1】order_of_key(x) -> 查询排名
49
    // 作用：返回严格小于 x 的元素个数（排名从 0 开始）
50
    int rank = tr.order_of_key(30);
51
    cout << "严格小于 30 的元素个数: " << rank << "\n"; // 输出 2
52

53
    // 【核心操作 2】find_by_order(k) -> 按排名查询
54
    // 作用：返回第 k 小元素的迭代器（k 从 0 开始计算）
55
    auto it = tr.find_by_order(2);
56
    cout << "第 2 小（0-indexed）的元素: " << *it << "\n\n"; // 输出 30
57

58
    // --------------------------------------
59
    // 多重集演示 (竞赛实战推荐)
60
    // --------------------------------------
61
    ordered_multiset mtr;
62
    int id_counter = 0; // 全局或局部的唯一时间戳/ID分配器
63

64
    // 插入多个相同元素
65
    mtr.insert({15, ++id_counter});
66
    mtr.insert({15, ++id_counter});
67
    mtr.insert({15, ++id_counter});
68
    mtr.insert({5, ++id_counter});
69

70
    // 1. 查询排名：直接查 {x, 0}
71
    // 因为 second=0 比所有真实插入的同值元素（ID>=1）都小
72
    cout << "严格小于 15 的元素个数: " << mtr.order_of_key({15, 0}) << "\n"; // 输出 1 (只有 5)
73

74
    // 2. 按排名查询：k_th_val = find_by_order(k)->first
75
    cout << "第 1 小的元素的值: " << mtr.find_by_order(1)->first << "\n"; // 输出 15
76

77
    // 3. 安全删除操作：需要先用 lower_bound 找到具体的 pair，然后再 erase
78
    auto del_it = mtr.lower_bound({15, 0});
79
    if (del_it != mtr.end() && del_it->first == 15) {
80
        mtr.erase(del_it); // 只删除了其中一个 15
81
    }
82

83
    // 验证删除后 15 还剩下几个 (原本3个，删了1个，应该剩2个，加上前置的5，总size为3)
84
    cout << "删除一个 15 后，树的大小: " << mtr.size() << "\n"; // 输出 3
85

86
    return 0;
87
}

隐式Treap（无旋Treap / FHQ-Treap）#

1
struct Implicit_treap {
2
    vector<int> val;
3
    vector<int> siz;
4
    vector<unsigned int> pri;
5
    vector<pair<int, int>> nxt;
6
    vector<bool> rev;   // 区间翻转懒标记
7
    int root;
8
    mt19937 rng;
9

10
    Implicit_treap() : root(0), rng(random_device{}()) {
11
        val.push_back(0);
12
        siz.push_back(0);
13
        pri.push_back(0);
14
        nxt.push_back({0, 0});
15
        rev.push_back(false);
16
    }
17

18
    int newp(int v) {
19
        val.push_back(v);
20
        siz.push_back(1);
21
        pri.push_back(rng());
22
        nxt.push_back({0, 0});
23
        rev.push_back(false);
24
        return siz.size() - 1;
25
    }
26

27
    void push_up(int p) {
28
        if (p) {
29
            const auto [lp, rp] = nxt[p];
30
            siz[p] = siz[lp] + siz[rp] + 1;
31
        }
32
    }
33

34
    void push_down(int p) {
35
        if (p && rev[p]) {
36
            // 1. 实际执行翻转：交换左右儿子
37
            swap(nxt[p].first, nxt[p].second);
38
            // 2. 交换后，再取出当前真正的左右儿子去下传标记
39
            const auto [lp, rp] = nxt[p];
40
            if (lp) rev[lp] = !rev[lp];
41
            if (rp) rev[rp] = !rev[rp];
42
            // 3. 清除当前节点的标记
43
            rev[p] = false;
44
        }
45
    }
46

47
    pair<int, int> split(int p, int k) {
48
        if (!p) return {0, 0};
49

50
        push_down(p);
51

52
        const auto [lp, rp] = nxt[p];
53
        if (siz[lp] + 1 <= k) {
54
            auto [rx, ry] = split(rp, k - siz[lp] - 1);
55
            nxt[p].second = rx;
56
            push_up(p);
57
            return {p, ry};
58
        } else {
59
            auto [lx, ly] = split(lp, k);
60
            nxt[p].first = ly;
61
            push_up(p);
62
            return {lx, p};
63
        }
64
    }
65

66
    int merge(int x, int y) {
67
        if (!x || !y) return x | y;
68
        push_down(x);
69
        push_down(y);
70

71
        if (pri[x] > pri[y]) {
72
            nxt[x].second = merge(nxt[x].second, y);
73
            push_up(x);
74
            return x;
75
        } else {
76
            nxt[y].first = merge(x, nxt[y].first);
77
            push_up(y);
78
            return y;
79
        }
80
    }
81

82
    // --- 0-based 对外接口 ---
83

84
    // 在下标 pos 处插入值 v (插入后它将成为新的 array[pos])
85
    // 允许 pos 取 0 到 siz[root] (即支持在首尾插入)
86
    void insert(int pos, int v) {
87
        auto [x, y] = split(root, pos);
88
        root = merge(merge(x, newp(v)), y);
89
    }
90

91
    // 删除下标为 pos 的节点 (0 <= pos < siz[root])
92
    void erase(int pos) {
93
        auto [x, z] = split(root, pos + 1);
94
        auto [lx, del] = split(x, pos);
95
        root = merge(lx, z);
96
    }
97

98
    // 反转区间 [l, r] (0 <= l <= r < siz[root])
99
    void reverse_range(int l, int r) {
100
        if(l >= r) return;
101
        auto [x, z] = split(root, r + 1);
102
        auto [lx, y] = split(x, l);
103
        rev[y] = !rev[y];
104
        root = merge(merge(lx, y), z);
105
    }
106
};

其他数据结构与算法#

稀疏表#

1
#include <bits/stdc++.h>
2
using namespace std;
3
template <typename S>
4
struct SpTable {
5
    int n;
6
    int log;
7
    vector<vector<S>> st;
8

9
    SpTable(const vector<S>& a) {
10
        build(a);
11
    }
12
    void build(const vector<S>& a) {
13
        n = a.size();
14
        log = 32 - __builtin_clz(n);
15
        st.assign(log, vector<S>(n));
16
        for (int i = 0; i < n; i++) {
17
            st[0][i] = a[i];
18
        }
19
        for (int j = 1; j < log; j++) {
20
            for (int i = 0; i + (1 << j) <= n; i++) {
21
                st[j][i] = op(st[j - 1][i], st[j - 1][i + (1 << (j - 1))]);
22
            }
23
        }
24
    }
25
    S query(int l, int r) {
26
        int j = 31 - __builtin_clz(r - l);
27
        return op(st[j][l], st[j][r - (1 << j)]);
28
    }
29
};

离散化#

1
// 离散化
2
// 许多题目在处理数据的时候，并不关心值本身，只关心值之间的大小关系
3
// 这个时候进行离散化是比较正确的决定
4
template <typename T>
5
vector<T> Discretization(vector<T> o_array)
6
{
7
    vector<T> knum = o_array;
8
    sort(knum.begin(), knum.end());
9
    knum.erase(unique(knum.begin(), knum.end()),knum.end());
10
    for(int i = 0;i<o_array.size();i++){
11
        o_array[i] = lower_bound(knum.begin(),knum.end(),o_array[i]) - knum.begin();
12
    }
13
    return o_array;
14
}
15
//传入一个数组，传出离散化之后的数组 如 num = Normalization(num);

哈希表#

1
#include <bits/stdc++.h>
2
using namespace std;
3
using ull = unsigned long long;
4
#include <ext/pb_ds/assoc_container.hpp>
5
#include <ext/pb_ds/tree_policy.hpp>
6
struct chash {
7
    static ull splitmix64(ull x) {
8
        x += 0x9e3779b97f4a7c15;
9
        x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;
10
        x = (x ^ (x >> 27)) * 0x94d049bb133111eb;
11
        return x ^ (x >> 31);
12
    }
13

14
    size_t operator()(ull x) const {
15
        static const ull F =
16
            chrono::steady_clock::now().time_since_epoch().count();
17
        return splitmix64(x + F);
18
    }
19
};
20
template<typename K, typename V>
21
using HashMap = __gnu_pbds::gp_hash_table<K, V, chash>;

欧拉序#

1
//许多处理树上问题，特别是对 子树 进行操作的时候 ，往往利用欧拉序将树拍平后，再利用线段树维护
2
//实现了一个结构体，初始化n为节点数，并传入边（需人为确保是一棵树）
3
//调用成员函数.init后，会生成两个个长度为n的数组
4
//第一个为树的dfn序，注意默认视1为根节点
5
//第二个为子树索引，tp[p] 代表 p节点所在子树在dfn上对应的连续区间 ，请使用auto[lp,rp] = pt[p] 调用
6
struct GenerateEulerTourOrder
7
{
8
    vector<int>dfn;
9
    vector<pair<int,int>>tp;
10
    vector<vector<int>> g;
11
    GenerateEulerTourOrder(int n){
12
        dfn.reserve(n);
13
        tp.resize(n+1);
14
        g.resize(n+1);
15
    }
16
    void add(int u,int v){
17
        g[u].push_back(v);
18
        g[v].push_back(u);
19
    }
20
    void init(int head = 1){
21
        auto dfs = [this](auto&&self,int n,int p) ->void {
22
            auto & [lp,rp] = tp[n];
23
            dfn.emplace_back(n);
24
            lp = dfn.size()-1;
25
            for(auto nextp:g[n]){
26
                if(nextp== p) continue;
27
                else{
28
                    self(self,nextp,n);
29
                }
30
            }
31
            rp = dfn.size() - 1;
32
        };
33
        dfs(dfs,head,-1);
34
        return;
35
    }
36
};

莫队#

1
struct Mo {
2
    // 存储左闭右开区间[L, R)
3
    vector<pair<int, int>> queries;
4

5
    // 增加一个左闭右开的查询
6
    void add_query(int lt, int rt) {
7
        queries.emplace_back(lt, rt);
8
    }
9

10
    // 计算最优执行顺序
11
    vector<int> get_mo_order() {
12
        int qsiz = queries.size();
13
        if (qsiz == 0) return {};
14

15
        int tmp = 0;
16
        for (auto&& [lt, rt] : queries) {
17
            if (lt > tmp) tmp = lt;
18
            if (rt > tmp) tmp = rt;
19
        }
20

21
        // 理论最优块大小为 sqrt(3) * N / sqrt(2 * Q)
22
        int blocksiz = max(1, (int)(1.732 * tmp / sqrt(max(1.0, 2.0 * qsiz))));
23

24
        vector<int> ret(qsiz);
25
        iota(ret.begin(), ret.end(), 0);
26

27
        sort(ret.begin(), ret.end(), [&](int a, int b) {
28
            int ab = queries[a].first / blocksiz, bb = queries[b].first / blocksiz;
29
            if (ab != bb) return ab < bb;
30
            return (ab & 1) ? (queries[a].second > queries[b].second) : (queries[a].second < queries[b].second);
31
        });
32
        // 如果卡常, 则启用下面这段注释
33

34
        // auto get_cost = [&](int a, int b) -> int {
35
        //     return abs(queries[ret[a]].first - queries[ret[b]].first) + abs(queries[ret[a]].second - queries[ret[b]].second);
36
        // };
37

38
        // for (int k = 0; k + 4 < qsiz; ++k) {
39
        //     if (get_cost(k, k + 2) + get_cost(k + 1, k + 3) < get_cost(k, k + 1) + get_cost(k + 2, k + 3)) {
40
        //         swap(ret[k + 1], ret[k + 2]);
41
        //     }
42
        //     if (get_cost(k, k + 3) + get_cost(k + 1, k + 4) < get_cost(k, k + 1) + get_cost(k + 3, k + 4)) {
43
        //         swap(ret[k + 1], ret[k + 3]);
44
        //     }
45
        // }
46
        return ret;
47
    }
48

49
    // 非对称增删
50
    template <class F1, class F2, class F3, class F4, class Q>
51
    void solve(F1 add_l, F2 add_r, F3 rm_l, F4 rm_r, Q query) {
52
        auto I = get_mo_order();
53
        int l = 0, r = 0;
54
        for (int idx : I) {
55
            while (l > queries[idx].first) add_l(--l);
56
            while (r < queries[idx].second) add_r(r++);
57
            while (l < queries[idx].first) rm_l(l++);
58
            while (r > queries[idx].second) rm_r(--r);
59
            query(idx);
60
        }
61
    }
62

63
    // 对称增删
64
    template <class F1, class F2, class Q>
65
    void solve(F1 add, F2 rm, Q query) {
66
        solve(add, add, rm, rm, query);
67
    }
68
};

回滚莫队#

此模板仅包含思路注释，无具体实现代码。

1
//常规莫队在挪动指针时，总是双指针交替挪动，O(1)更新答案
2
//回滚莫队的思想是：
3
//对于一些可离线的查询，我们进行add操作比较简单，del操作比较复杂
4
//我们在查询的时候，先确保右指针单调挪动，每次右指针挪动前，把左指针回溯到最初始的状态（区块右边界）
5
//再将指针左移以匹配区间
6
//注意到这样就只有add操作，我们只需要维护一个操作栈即可正确回溯，注意到左指针最多挪动sqrt(n)次

几何#

库实数类实现（双精度）#

1
using Real = int;
2
using Point = complex<Real>;
3

4
Real cross(const Point &a, const Point &b) {
5
    return (conj(a) * b).imag();
6
}
7
Real dot(const Point &a, const Point &b) {
8
    return (conj(a) * b).real();
9
}

平面几何必要初始化#

字符串读入浮点数#

1
const int Knum = 4;
2
int read(int k = Knum) {
3
    string s;
4
    cin >> s;
5

6
    int num = 0;
7
    int it = s.find('.');
8
    if (it != -1) { // 存在小数点
9
        num = s.size() - it - 1; // 计算小数位数
10
        s.erase(s.begin() + it); // 删除小数点
11
    }
12
    for (int i = 1; i <= k - num; i++) { // 补全小数位数
13
        s += '0';
14
    }
15
    return stoi(s);
16
}

预置函数#

1
using ld = long double;
2
const ld PI = acos(-1);
3
const ld EPS = 1e-7;
4
const ld INF = numeric_limits<ld>::max();
5
#define cc(x) cout << fixed << setprecision(x);
6

7
ld fgcd(ld x, ld y) { // 实数域gcd
8
    return abs(y) < EPS ? abs(x) : fgcd(y, fmod(x, y));
9
}
10
template<class T, class S> bool equal(T x, S y) {
11
    return -EPS < x - y && x - y < EPS;
12
}
13
template<class T> int sign(T x) {
14
    if (-EPS < x && x < EPS) return 0;
15
    return x < 0 ? -1 : 1;
16
}

点线封装#

1
template<class T> struct Point { // 在C++17下使用 emplace_back 绑定可能会导致CE！
2
    T x, y;
3
    Point(T x_ = 0, T y_ = 0) : x(x_), y(y_) {} // 初始化
4
    template<class U> operator Point<U>() { // 自动类型匹配
5
        return Point<U>(U(x), U(y));
6
    }
7
    Point &operator+=(Point p) & { return x += p.x, y += p.y, *this; }
8
    Point &operator+=(T t) & { return x += t, y += t, *this; }
9
    Point &operator-=(Point p) & { return x -= p.x, y -= p.y, *this; }
10
    Point &operator-=(T t) & { return x -= t, y -= t, *this; }
11
    Point &operator*=(T t) & { return x *= t, y *= t, *this; }
12
    Point &operator/=(T t) & { return x /= t, y /= t, *this; }
13
    Point operator-() const { return Point(-x, -y); }
14
    friend Point operator+(Point a, Point b) { return a += b; }
15
    friend Point operator+(Point a, T b) { return a += b; }
16
    friend Point operator-(Point a, Point b) { return a -= b; }
17
    friend Point operator-(Point a, T b) { return a -= b; }
18
    friend Point operator*(Point a, T b) { return a *= b; }
19
    friend Point operator*(T a, Point b) { return b *= a; }
20
    friend Point operator/(Point a, T b) { return a /= b; }
21
    friend bool operator<(Point a, Point b) {
22
        return equal(a.x, b.x) ? a.y < b.y - EPS : a.x < b.x - EPS;
23
    }
24
    friend bool operator>(Point a, Point b) { return b < a; }
25
    friend bool operator==(Point a, Point b) { return !(a < b) && !(b < a); }
26
    friend bool operator!=(Point a, Point b) { return a < b || b < a; }
27
    friend auto &operator>>(istream &is, Point &p) {
28
        return is >> p.x >> p.y;
29
    }
30
    friend auto &operator<<(ostream &os, Point p) {
31
        return os << "(" << p.x << ", " << p.y << ")";
32
    }
33
};
34
template<class T> struct Line {
35
    Point<T> a, b;
36
    Line(Point<T> a_ = Point<T>(), Point<T> b_ = Point<T>()) : a(a_), b(b_) {}
37
    template<class U> operator Line<U>() { // 自动类型匹配
38
        return Line<U>(Point<U>(a), Point<U>(b));
39
    }
40
    friend auto &operator<<(ostream &os, Line l) {
41
        return os << "<" << l.a << ", " << l.b << ">";
42
    }
43
};

叉乘#

定义公式 $a\times b=|a||b|\sin \theta$ 。

1
template<class T> T cross(Point<T> a, Point<T> b) { // 叉乘
2
    return a.x * b.y - a.y * b.x;
3
}
4
template<class T> T cross(Point<T> p1, Point<T> p2, Point<T> p0) { // 叉乘 (p1 - p0) x (p2 - p0);
5
    return cross(p1 - p0, p2 - p0);
6
}

点乘#

定义公式 $a\times b=|a||b|\cos \theta$ 。

1
template<class T> T dot(Point<T> a, Point<T> b) { // 点乘
2
    return a.x * b.x + a.y * b.y;
3
}
4
template<class T> T dot(Point<T> p1, Point<T> p2, Point<T> p0) { // 点乘 (p1 - p0) * (p2 - p0);
5
    return dot(p1 - p0, p2 - p0);
6
}

欧几里得距离公式#

最常用的距离公式。需要注意，开根号会丢失精度，如无强制要求，先不要开根号，留到最后一步一起开。

1
template <class T> ld dis(T x1, T y1, T x2, T y2) {
2
    ld val = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2);
3
    return sqrt(val);
4
}
5
template <class T> ld dis(Point<T> a, Point<T> b) {
6
    return dis(a.x, a.y, b.x, b.y);
7
}

曼哈顿距离公式#

1
template <class T> T dis1(Point<T> p1, Point<T> p2) { // 曼哈顿距离公式
2
    return abs(p1.x - p2.x) + abs(p1.y - p2.y);
3
}

将向量转换为单位向量#

1
Point<ld> standardize(Point<ld> vec) { // 转换为单位向量
2
    return vec / sqrt(vec.x * vec.x + vec.y * vec.y);
3
}

向量旋转#

将当前向量移动至原点后顺时针旋转 $90^{\circ}$ ，即获取垂直于当前向量的、起点为原点的向量。在计算垂线时非常有用。例如，要想获取点 $a$ 绕点 $o$ 顺时针旋转 $90^{\circ}$ 后的点，可以这样书写代码：auto ans = o + rotate(o, a); ；如果是逆时针旋转，那么只需更改符号即可：auto ans = o - rotate(o, a); 。

1
template<class T> Point<T> rotate(Point<T> p1, Point<T> p2) { // 旋转
2
    Point<T> vec = p1 - p2;
3
    return {-vec.y, vec.x};
4
}

平面角度与弧度#

弧度角度相互转换#

1
ld toDeg(ld x) { // 弧度转角度
2
    return x * 180 / PI;
3
}
4
ld toArc(ld x) { // 角度转弧度
5
    return PI / 180 * x;
6
}

正弦定理#

$\dfrac{a}{\sin A}=\dfrac{b}{\sin B}=\dfrac{c}{\sin C}=2R$ ，其中 $R$ 为三角形外接圆半径；

余弦定理（已知三角形三边，求角）#

$\cos C=\dfrac{a^2+b^2-c^2}{2ab},\cos B=\dfrac{a^2+c^2-b^2}{2ac},\cos A=\dfrac{b^2+c^2-a^2}{2bc}$ 。可以借此推导出三角形面积公式 $S_{\triangle ABC}=\dfrac{ab\cdot\sin C}{2}=\dfrac{bc\cdot\sin A}{2}=\dfrac{ac\cdot\sin B}{2}$ 。

注意，计算格式是：由 $b,c,a$ 三边求 $\angle A$ ；由 $a, c, b$ 三边求 $\angle B$ ；由 $a, b, c$ 三边求 $\angle C$ 。

1
ld angle(ld a, ld b, ld c) { // 余弦定理
2
    ld val = acos((a * a + b * b - c * c) / (2.0 * a * b)); // 计算弧度
3
    return val;
4
}

求两向量的夹角#

能够计算 $[0^{\circ},180^{\circ}]$ 区间的角度。

1
ld angle(Point<ld> a, Point<ld> b) {
2
    ld val = abs(cross(a, b));
3
    return abs(atan2(val, a.x * b.x + a.y * b.y));
4
}

向量旋转任意角度#

逆时针旋转，转换公式： $\left\{\begin{matrix} x'=x\cos \theta-y\sin \theta \\ y'=x\sin \theta+y\cos \theta \end{matrix}\right.$

1
Point<ld> rotate(Point<ld> p, ld rad) {
2
    return {p.x * cos(rad) - p.y * sin(rad), p.x * sin(rad) + p.y * cos(rad)};
3
}

点绕点旋转任意角度#

逆时针旋转，转换公式： $\left\{\begin{matrix} x'=(x_0-x_1)\cos\theta+(y_0-y_1)\sin\theta+x_1 \\ y'=(x_1-x_0)\sin\theta+(y_0-y_1)\cos\theta+y_1 \end{matrix}\right.$

1
Point<ld> rotate(Point<ld> a, Point<ld> b, ld rad) {
2
    ld x = (a.x - b.x) * cos(rad) + (a.y - b.y) * sin(rad) + b.x;
3
    ld y = (b.x - a.x) * sin(rad) + (a.y - b.y) * cos(rad) + b.y;
4
    return {x, y};
5
}

平面点线相关#

点是否在直线上（三点是否共线）#

1
template<class T> bool onLine(Point<T> a, Point<T> b, Point<T> c) {
2
    return sign(cross(b, a, c)) == 0;
3
}
4
template<class T> bool onLine(Point<T> p, Line<T> l) {
5
    return onLine(p, l.a, l.b);
6
}

点是否在向量（直线）左侧#

需要注意，向量的方向会影响答案；点在向量上时不视为在左侧。

1
template<class T> bool pointOnLineLeft(Pt p, Lt l) {
2
    return cross(l.b, p, l.a) > 0;
3
}

两点是否在直线同侧/异侧#

1
template<class T> bool pointOnLineSide(Pt p1, Pt p2, Lt vec) {
2
    T val = cross(p1, vec.a, vec.b) * cross(p2, vec.a, vec.b);
3
    return sign(val) == 1;
4
}
5
template<class T> bool pointNotOnLineSide(Pt p1, Pt p2, Lt vec) {
6
    T val = cross(p1, vec.a, vec.b) * cross(p2, vec.a, vec.b);
7
    return sign(val) == -1;
8
}

两直线相交交点#

在使用前需要先判断直线是否平行。

1
Pd lineIntersection(Ld l1, Ld l2) {
2
    ld val = cross(l2.b - l2.a, l1.a - l2.a) / cross(l2.b - l2.a, l1.a - l1.b);
3
    return l1.a + (l1.b - l1.a) * val;
4
}

两直线是否平行/垂直/相同#

1
template<class T> bool lineParallel(Lt p1, Lt p2) {
2
    return sign(cross(p1.a - p1.b, p2.a - p2.b)) == 0;
3
}
4
template<class T> bool lineVertical(Lt p1, Lt p2) {
5
    return sign(dot(p1.a - p1.b, p2.a - p2.b)) == 0;
6
}
7
template<class T> bool same(Line<T> l1, Line<T> l2) {
8
    return lineParallel(Line{l1.a, l2.b}, {l1.b, l2.a}) &&
9
           lineParallel(Line{l1.a, l2.a}, {l1.b, l2.b}) && lineParallel(l1, l2);
10
}

点到直线的最近距离与最近点#

1
pair<Pd, ld> pointToLine(Pd p, Ld l) {
2
    Pd ans = lineIntersection({p, p + rotate(l.a, l.b)}, l);
3
    return {ans, dis(p, ans)};
4
}

如果只需要计算最近距离，下方的写法可以减少书写的代码量，效果一致。

1
template<class T> ld disPointToLine(Pt p, Lt l) {
2
    ld ans = cross(p, l.a, l.b);
3
    return abs(ans) / dis(l.a, l.b); // 面积除以底边长
4
}

点是否在线段上#

1
template<class T> bool pointOnSegment(Pt p, Lt l) { // 端点也算作在直线上
2
    return sign(cross(p, l.a, l.b)) == 0 && min(l.a.x, l.b.x) <= p.x && p.x <= max(l.a.x, l.b.x) &&
3
           min(l.a.y, l.b.y) <= p.y && p.y <= max(l.a.y, l.b.y);
4
}
5
template<class T> bool pointOnSegment(Pt p, Lt l) { // 端点不算
6
    return pointOnSegment(p, l) && min(l.a.x, l.b.x) < p.x && p.x < max(l.a.x, l.b.x) &&
7
           min(l.a.y, l.b.y) < p.y && p.y < max(l.a.y, l.b.y);
8
}

点到线段的最近距离与最近点#

1
pair<Pd, ld> pointToSegment(Pd p, Ld l) {
2
    if (sign(dot(p, l.b, l.a)) == -1) { // 特判到两端点的距离
3
        return {l.a, dis(p, l.a)};
4
    } else if (sign(dot(p, l.a, l.b)) == -1) {
5
        return {l.b, dis(p, l.b)};
6
    }
7
    return pointToLine(p, l);
8
}

点在直线上的投影点（垂足）#

1
Pd project(Pd p, Ld l) { // 投影
2
    Pd vec = l.b - l.a;
3
    ld r = dot(vec, p - l.a) / (vec.x * vec.x + vec.y * vec.y);
4
    return l.a + vec * r;
5
}

线段的中垂线#

1
template<class T> Lt midSegment(Lt l) {
2
    Pt mid = (l.a + l.b) / 2; // 线段中点
3
    return {mid, mid + rotate(l.a, l.b)};
4
}

两线段是否相交及交点#

该扩展版可以同时返回相交状态和交点，分为四种情况： $0$ 代表不相交； $1$ 代表普通相交； $2$ 代表重叠（交于两个点）； $3$ 代表相交于端点。需要注意，部分运算可能会使用到直线求交点，此时务必保证变量类型为浮点数！

1
template<class T> tuple<int, Pt, Pt> segmentIntersection(Lt l1, Lt l2) {
2
    auto [s1, e1] = l1;
3
    auto [s2, e2] = l2;
4
    auto A = max(s1.x, e1.x), AA = min(s1.x, e1.x);
5
    auto B = max(s1.y, e1.y), BB = min(s1.y, e1.y);
6
    auto C = max(s2.x, e2.x), CC = min(s2.x, e2.x);
7
    auto D = max(s2.y, e2.y), DD = min(s2.y, e2.y);
8
    if (A < CC || C < AA || B < DD || D < BB) {
9
        return {0, {}, {}};
10
    }
11
    if (sign(cross(e1 - s1, e2 - s2)) == 0) {
12
        if (sign(cross(s2, e1, s1)) != 0) {
13
            return {0, {}, {}};
14
        }
15
        Pt p1(max(AA, CC), max(BB, DD));
16
        Pt p2(min(A, C), min(B, D));
17
        if (!pointOnSegment(p1, l1)) {
18
            swap(p1.y, p2.y);
19
        }
20
        if (p1 == p2) {
21
            return {3, p1, p2};
22
        } else {
23
            return {2, p1, p2};
24
        }
25
    }
26
    auto cp1 = cross(s2 - s1, e2 - s1);
27
    auto cp2 = cross(s2 - e1, e2 - e1);
28
    auto cp3 = cross(s1 - s2, e1 - s2);
29
    auto cp4 = cross(s1 - e2, e1 - e2);
30
    if (sign(cp1 * cp2) == 1 || sign(cp3 * cp4) == 1) {
31
        return {0, {}, {}};
32
    }
33
    // 使用下方函数时请使用浮点数
34
    Pd p = lineIntersection(l1, l2);
35
    if (sign(cp1) != 0 && sign(cp2) != 0 && sign(cp3) != 0 && sign(cp4) != 0) {
36
        return {1, p, p};
37
    } else {
38
        return {3, p, p};
39
    }
40
}

如果不需要求交点，那么使用快速排斥+跨立实验即可，其中重叠、相交于端点均视为相交。

1
template<class T> bool segmentIntersection(Lt l1, Lt l2) {
2
    auto [s1, e1] = l1;
3
    auto [s2, e2] = l2;
4
    auto A = max(s1.x, e1.x), AA = min(s1.x, e1.x);
5
    auto B = max(s1.y, e1.y), BB = min(s1.y, e1.y);
6
    auto C = max(s2.x, e2.x), CC = min(s2.x, e2.x);
7
    auto D = max(s2.y, e2.y), DD = min(s2.y, e2.y);
8
    return A >= CC && B >= DD && C >= AA && D >= BB &&
9
           sign(cross(s1, s2, e1) * cross(s1, e1, e2)) == 1 &&
10
           sign(cross(s2, s1, e2) * cross(s2, e2, e1)) == 1;
11
}

极角排序#

从 $x$ 轴正半轴开始.

1
template <class S>
2
int halfPlane(Point<S> a) {
3
    if (cmp(a.y, 0) == 0) {
4
        return cmp(a.x, 0);
5
    }
6
    return cmp(a.y, 0); // 如果需要从负半轴开始则把这里改成负的
7
}
8
template <class S>
9
bool polarCmp(Point<S> a, Point<S> b) {
10
    if (halfPlane(a) != halfPlane(b)) {
11
        return halfPlane(a) > halfPlane(b);
12
    }
13
    if (sgn(crs(a, b)) == 0) {
14
        return len2(a) < len2(b);
15
    }
16
    return sgn(crs(a, b)) > 0;
17
}

平面圆相关（浮点数处理）#

点到圆的最近点#

同时返回最近点与最近距离。需要注意，当点为圆心时，这样的点有无数个，此时我们视作输入错误，直接返回圆心。

1
pair<Pd, ld> pointToCircle(Pd p, Pd o, ld r) {
2
    Pd U = o, V = o;
3
    ld d = dis(p, o);
4
    if (sign(d) == 0) { // p 为圆心时返回圆心本身
5
        return {o, 0};
6
    }
7
    ld val1 = r * abs(o.x - p.x) / d;
8
    ld val2 = r * abs(o.y - p.y) / d * ((o.x - p.x) * (o.y - p.y) < 0 ? -1 : 1);
9
    U.x += val1, U.y += val2;
10
    V.x -= val1, V.y -= val2;
11
    if (dis(U, p) < dis(V, p)) {
12
        return {U, dis(U, p)};
13
    } else {
14
        return {V, dis(V, p)};
15
    }
16
}

根据圆心角获取圆上某点#

将圆上最右侧的点以圆心为旋转中心，逆时针旋转 rad 度。

1
Point<ld> getPoint(Point<ld> p, ld r, ld rad) {
2
    return {p.x + cos(rad) * r, p.y + sin(rad) * r};
3
}

直线是否与圆相交及交点#

$0$ 代表不相交； $1$ 代表相切； $2$ 代表相交。

1
tuple<int, Pd, Pd> lineCircleCross(Ld l, Pd o, ld r) {
2
    Pd P = project(o, l);
3
    ld d = dis(P, o), tmp = r * r - d * d;
4
    if (sign(tmp) == -1) {
5
        return {0, {}, {}};
6
    } else if (sign(tmp) == 0) {
7
        return {1, P, {}};
8
    }
9
    Pd vec = standardize(l.b - l.a) * sqrt(tmp);
10
    return {2, P + vec, P - vec};
11
}

线段是否与圆相交及交点#

$0$ 代表不相交； $1$ 代表相切； $2$ 代表相交于一个点； $3$ 代表相交于两个点。

1
tuple<int, Pd, Pd> segmentCircleCross(Ld l, Pd o, ld r) {
2
    auto [type, U, V] = lineCircleCross(l, o, r);
3
    bool f1 = pointOnSegment(U, l), f2 = pointOnSegment(V, l);
4
    if (type == 1 && f1) {
5
        return {1, U, {}};
6
    } else if (type == 2 && f1 && f2) {
7
        return {3, U, V};
8
    } else if (type == 2 && f1) {
9
        return {2, U, {}};
10
    } else if (type == 2 && f2) {
11
        return {2, V, {}};
12
    } else {
13
        return {0, {}, {}};
14
    }
15
}

两圆是否相交及交点#

$0$ 代表内含； $1$ 代表相离； $2$ 代表相切； $3$ 代表相交。

1
tuple<int, Pd, Pd> circleIntersection(Pd p1, ld r1, Pd p2, ld r2) {
2
    ld x1 = p1.x, x2 = p2.x, y1 = p1.y, y2 = p2.y, d = dis(p1, p2);
3
    if (sign(abs(r1 - r2) - d) == 1) {
4
        return {0, {}, {}};
5
    } else if (sign(r1 + r2 - d) == -1) {
6
        return {1, {}, {}};
7
    }
8
    ld a = r1 * (x1 - x2) * 2, b = r1 * (y1 - y2) * 2, c = r2 * r2 - r1 * r1 - d * d;
9
    ld p = a * a + b * b, q = -a * c * 2, r = c * c - b * b;
10
    ld cosa, sina, cosb, sinb;
11
    if (sign(d - (r1 + r2)) == 0 || sign(d - abs(r1 - r2)) == 0) {
12
        cosa = -q / p / 2;
13
        sina = sqrt(1 - cosa * cosa);
14
        Point<ld> p0 = {x1 + r1 * cosa, y1 + r1 * sina};
15
        if (sign(dis(p0, p2) - r2)) {
16
            p0.y = y1 - r1 * sina;
17
        }
18
        return {2, p0, p0};
19
    } else {
20
        ld delta = sqrt(q * q - p * r * 4);
21
        cosa = (delta - q) / p / 2;
22
        cosb = (-delta - q) / p / 2;
23
        sina = sqrt(1 - cosa * cosa);
24
        sinb = sqrt(1 - cosb * cosb);
25
        Pd ans1 = {x1 + r1 * cosa, y1 + r1 * sina};
26
        Pd ans2 = {x1 + r1 * cosb, y1 + r1 * sinb};
27
        if (sign(dis(ans1, p1) - r2)) ans1.y = y1 - r1 * sina;
28
        if (sign(dis(ans2, p2) - r2)) ans2.y = y1 - r1 * sinb;
29
        if (ans1 == ans2) ans1.y = y1 - r1 * sina;
30
        return {3, ans1, ans2};
31
    }
32
}

两圆相交面积#

上述所言四种相交情况均可计算，之所以不使用三角形面积计算公式是因为在计算过程中会出现“负数”面积（扇形面积与三角形面积的符号关系会随圆的位置关系发生变化），故公式全部重新推导，这里采用的是扇形面积减去扇形内部的那个三角形的面积。

1
ld circleIntersectionArea(Pd p1, ld r1, Pd p2, ld r2) {
2
    ld x1 = p1.x, x2 = p2.x, y1 = p1.y, y2 = p2.y, d = dis(p1, p2);
3
    if (sign(abs(r1 - r2) - d) >= 0) {
4
        return PI * min(r1 * r1, r2 * r2);
5
    } else if (sign(r1 + r2 - d) == -1) {
6
        return 0;
7
    }
8
    ld theta1 = angle(r1, dis(p1, p2), r2);
9
    ld area1 = r1 * r1 * (theta1 - sin(theta1 * 2) / 2);
10
    ld theta2 = angle(r2, dis(p1, p2), r1);
11
    ld area2 = r2 * r2 * (theta2 - sin(theta2 * 2) / 2);
12
    return area1 + area2;
13
}

三点确定一圆#

1
tuple<int, Pd, ld> getCircle(Pd A, Pd B, Pd C) {
2
    if (onLine(A, B, C)) { // 特判三点共线
3
        return {0, {}, 0};
4
    }
5
    Ld l1 = midSegment(Line{A, B});
6
    Ld l2 = midSegment(Line{A, C});
7
    Pd O = lineIntersection(l1, l2);
8
    return {1, O, dis(A, O)};
9
}

求解点到圆的切线数量与切点#

1
pair<int, vector<Point<ld>>> tangent(Point<ld> p, Point<ld> A, ld r) {
2
    vector<Point<ld>> ans; // 储存切点
3
    Point<ld> u = A - p;
4
    ld d = sqrt(dot(u, u));
5
    if (d < r) {
6
        return {0, {}};
7
    } else if (sign(d - r) == 0) { // 点在圆上
8
        ans.push_back(u);
9
        return {1, ans};
10
    } else {
11
        ld ang = asin(r / d);
12
        ans.push_back(getPoint(A, r, -ang));
13
        ans.push_back(getPoint(A, r, ang));
14
        return {2, ans};
15
    }
16
}

求解两圆的内公、外公切线数量与切点#

同时返回公切线数量以及每个圆的切点。

1
tuple<int, vector<Point<ld>>, vector<Point<ld>>> tangent(Point<ld> A, ld Ar, Point<ld> B, ld Br) {
2
    vector<Point<ld>> a, b; // 储存切点
3
    if (Ar < Br) {
4
        swap(Ar, Br);
5
        swap(A, B);
6
        swap(a, b);
7
    }
8
    int d = disEx(A, B), dif = Ar - Br, sum = Ar + Br;
9
    if (d < dif * dif) { // 内含，无
10
        return {0, {}, {}};
11
    }
12
    ld base = atan2(B.y - A.y, B.x - A.x);
13
    if (d == 0 && Ar == Br) { // 完全重合，无数条外公切线
14
        return {-1, {}, {}};
15
    }
16
    if (d == dif * dif) { // 内切，1条外公切线
17
        a.push_back(getPoint(A, Ar, base));
18
        b.push_back(getPoint(B, Br, base));
19
        return {1, a, b};
20
    }
21
    ld ang = acos(dif / sqrt(d));
22
    a.push_back(getPoint(A, Ar, base + ang)); // 保底2条外公切线
23
    a.push_back(getPoint(A, Ar, base - ang));
24
    b.push_back(getPoint(B, Br, base + ang));
25
    b.push_back(getPoint(B, Br, base - ang));
26
    if (d == sum * sum) { // 外切，多1条内公切线
27
        a.push_back(getPoint(A, Ar, base));
28
        b.push_back(getPoint(B, Br, base + PI));
29
    } else if (d > sum * sum) { // 相离，多2条内公切线
30
        ang = acos(sum / sqrt(d));
31
        a.push_back(getPoint(A, Ar, base + ang));
32
        a.push_back(getPoint(A, Ar, base - ang));
33
        b.push_back(getPoint(B, Br, base + ang + PI));
34
        b.push_back(getPoint(B, Br, base - ang + PI));
35
    }
36
    return {a.size(), a, b};
37
}

平面三角形相关（浮点数处理）#

三角形面积#

1
ld area(Point<ld> a, Point<ld> b, Point<ld> c) {
2
    return abs(cross(b, c, a)) / 2;
3
}

三角形外心#

三角形外接圆的圆心，即三角形三边垂直平分线的交点。

1
template<class T> Pt center1(Pt p1, Pt p2, Pt p3) { // 外心
2
    return lineIntersection(midSegment({p1, p2}), midSegment({p2, p3}));
3
}

三角形内心#

三角形内切圆的圆心，也是三角形三个内角的角平分线的交点。其到三角形三边的距离相等。

1
Pd center2(Pd p1, Pd p2, Pd p3) { // 内心
2
    #define atan2(p) atan2(p.y, p.x) // 注意先后顺序
3
    Line<ld> U = {p1, {}}, V = {p2, {}};
4
    ld m, n, alpha;
5
    m = atan2((p2 - p1));
6
    n = atan2((p3 - p1));
7
    alpha = (m + n) / 2;
8
    U.b = {p1.x + cos(alpha), p1.y + sin(alpha)};
9
    m = atan2((p1 - p2));
10
    n = atan2((p3 - p2));
11
    alpha = (m + n) / 2;
12
    V.b = {p2.x + cos(alpha), p2.y + sin(alpha)};
13
    return lineIntersection(U, V);
14
}

三角形垂心#

三角形的三条高线所在直线的交点。锐角三角形的垂心在三角形内；直角三角形的垂心在直角顶点上；钝角三角形的垂心在三角形外。

1
Pd center3(Pd p1, Pd p2, Pd p3) { // 垂心
2
    Ld U = {p1, p1 + rotate(p2, p3)}; // 垂线
3
    Ld V = {p2, p2 + rotate(p1, p3)};
4
    return lineIntersection(U, V);
5
}

平面直线方程转换#

浮点数计算直线的斜率#

一般很少使用到这个函数，因为斜率的取值不可控（例如接近平行于 $x,y$ 轴时）。需要注意，当直线平行于 $y$ 轴时斜率为 inf 。

1
template <class T> ld slope(Pt p1, Pt p2) { // 斜率，注意 inf 的情况
2
    return (p1.y - p2.y) / (p1.x - p2.x);
3
}
4
template <class T> ld slope(Lt l) {
5
    return slope(l.a, l.b);
6
}

分数精确计算直线的斜率#

调用分数四则运算精确计算斜率，返回最简分数，只适用于整数计算。

1
template<class T> Frac<T> slopeEx(Pt p1, Pt p2) {
2
    Frac<T> U = p1.y - p2.y;
3
    Frac<T> V = p1.x - p2.x;
4
    return U / V; // 调用分数精确计算
5
}

两点式转一般式#

返回由三个整数构成的方程，在输入较大时可能找不到较小的满足题意的一组整数解。可以处理平行于 $x,y$ 轴、两点共点的情况。

1
template<class T> tuple<T, T, T> getfun(Lt p) {
2
    T A = p.a.y - p.b.y, B = p.b.x - p.a.x, C = p.a.x * A + p.a.y * B;
3
    if (A < 0) { // 符号调整
4
        A = -A, B = -B, C = -C;
5
    } else if (A == 0) {
6
        if (B < 0) {
7
            B = -B, C = -C;
8
        } else if (B == 0 && C < 0) {
9
            C = -C;
10
        }
11
    }
12
    if (A == 0) { // 数值计算
13
        if (B == 0) {
14
            C = 0; // 共点特判
15
        } else {
16
            T g = fgcd(abs(B), abs(C));
17
            B /= g, C /= g;
18
        }
19
    } else if (B == 0) {
20
        T g = fgcd(abs(A), abs(C));
21
        A /= g, C /= g;
22
    } else {
23
        T g = fgcd(fgcd(abs(A), abs(B)), abs(C));
24
        A /= g, B /= g, C /= g;
25
    }
26
    return tuple{A, B, C}; // Ax + By = C
27
}

一般式转两点式#

由于整数点可能很大或者不存在，故直接采用浮点数；如果与 $x,y$ 轴有交点则取交点。可以处理平行于 $x,y$ 轴的情况。

1
Line<ld> getfun(int A, int B, int C) { // Ax + By = C
2
    ld x1 = 0, y1 = 0, x2 = 0, y2 = 0;
3
    if (A && B) { // 正常
4
        if (C) {
5
            x1 = 0, y1 = 1. * C / B;
6
            y2 = 0, x2 = 1. * C / A;
7
        } else { // 过原点
8
            x1 = 1, y1 = 1. * -A / B;
9
            x2 = 0, y2 = 0;
10
        }
11
    } else if (A && !B) { // 垂直
12
        if (C) {
13
            y1 = 0, x1 = 1. * C / A;
14
            y2 = 1, x2 = x1;
15
        } else {
16
            x1 = 0, y1 = 1;
17
            x2 = 0, y2 = 0;
18
        }
19
    } else if (!A && B) { // 水平
20
        if (C) {
21
            x1 = 0, y1 = 1. * C / B;
22
            x2 = 1, y2 = y1;
23
        } else {
24
            x1 = 1, y1 = 0;
25
            x2 = 0, y2 = 0;
26
        }
27
    } else { // 不合法，请特判
28
        assert(false);
29
    }
30
    return {{x1, y1}, {x2, y2}};
31
}

抛物线与 x 轴是否相交及交点#

$0$ 代表没有交点； $1$ 代表相切； $2$ 代表有两个交点。

1
tuple<int, ld, ld> getAns(ld a, ld b, ld c) {
2
    ld delta = b * b - a * c * 4;
3
    if (delta < 0.) {
4
        return {0, 0, 0};
5
    }
6
    delta = sqrt(delta);
7
    ld ans1 = -(delta + b) / 2 / a;
8
    ld ans2 = (delta - b) / 2 / a;
9
    if (ans1 > ans2) {
10
        swap(ans1, ans2);
11
    }
12
    if (sign(delta) == 0) {
13
        return {1, ans2, 0};
14
    }
15
    return {2, ans1, ans2};
16
}

平面多边形#

两向量构成的平面四边形有向面积#

1
template<class T> T areaEx(Point<T> p1, Point<T> p2, Point<T> p3) {
2
    return cross(b, c, a);
3
}

判断四个点能否组成矩形/正方形#

可以处理浮点数、共点的情况。返回分为三种情况： $2$ 代表构成正方形； $1$ 代表构成矩形； $0$ 代表其他情况。

1
template<class T> int isSquare(vector<Pt> x) {
2
    sort(x.begin(), x.end());
3
    if (equal(dis(x[0], x[1]), dis(x[2], x[3])) && sign(dis(x[0], x[1])) &&
4
        equal(dis(x[0], x[2]), dis(x[1], x[3])) && sign(dis(x[0], x[2])) &&
5
        lineParallel(Lt{x[0], x[1]}, Lt{x[2], x[3]}) &&
6
        lineParallel(Lt{x[0], x[2]}, Lt{x[1], x[3]}) &&
7
        lineVertical(Lt{x[0], x[1]}, Lt{x[0], x[2]})) {
8
        return equal(dis(x[0], x[1]), dis(x[0], x[2])) ? 2 : 1;
9
    }
10
    return 0;
11
}

点是否在任意多边形内#

射线法判定， $t$ 为穿越次数，当其为奇数时即代表点在多边形内部；返回 $2$ 代表点在多边形边界上。

1
template<class T> int pointInPolygon(Point<T> a, vector<Point<T>> p) {
2
    int n = p.size();
3
    for (int i = 0; i < n; i++) {
4
        if (pointOnSegment(a, Line{p[i], p[(i + 1) % n]})) {
5
            return 2;
6
        }
7
    }
8
    int t = 0;
9
    for (int i = 0; i < n; i++) {
10
        auto u = p[i], v = p[(i + 1) % n];
11
        if (u.x < a.x && v.x >= a.x && pointOnLineLeft(a, Line{v, u})) {
12
            t ^= 1;
13
        }
14
        if (u.x >= a.x && v.x < a.x && pointOnLineLeft(a, Line{u, v})) {
15
            t ^= 1;
16
        }
17
    }
18
    return t == 1;
19
}

线段是否在任意多边形内部#

1
template<class T>
2
bool segmentInPolygon(Line<T> l, vector<Point<T>> p) {
3
// 线段与多边形边界不相交且两端点都在多边形内部
4
#define L(x, y) pointOnLineLeft(x, y)
5
    int n = p.size();
6
    if (!pointInPolygon(l.a, p)) return false;
7
    if (!pointInPolygon(l.b, p)) return false;
8
    for (int i = 0; i < n; i++) {
9
        auto u = p[i];
10
        auto v = p[(i + 1) % n];
11
        auto w = p[(i + 2) % n];
12
        auto [t, p1, p2] = segmentIntersection(l, Line(u, v));
13
        if (t == 1) return false;
14
        if (t == 0) continue;
15
        if (t == 2) {
16
            if (pointOnSegment(v, l) && v != l.a && v != l.b) {
17
                if (cross(v - u, w - v) > 0) {
18
                    return false;
19
                }
20
            }
21
        } else {
22
            if (p1 != u && p1 != v) {
23
                if (L(l.a, Line(v, u)) || L(l.b, Line(v, u))) {
24
                    return false;
25
                }
26
            } else if (p1 == v) {
27
                if (l.a == v) {
28
                    if (L(u, l)) {
29
                        if (L(w, l) && L(w, Line(u, v))) {
30
                            return false;
31
                        }
32
                    } else {
33
                        if (L(w, l) || L(w, Line(u, v))) {
34
                            return false;
35
                        }
36
                    }
37
                } else if (l.b == v) {
38
                    if (L(u, Line(l.b, l.a))) {
39
                        if (L(w, Line(l.b, l.a)) && L(w, Line(u, v))) {
40
                            return false;
41
                        }
42
                    } else {
43
                        if (L(w, Line(l.b, l.a)) || L(w, Line(u, v))) {
44
                            return false;
45
                        }
46
                    }
47
                } else {
48
                    if (L(u, l)) {
49
                        if (L(w, Line(l.b, l.a)) || L(w, Line(u, v))) {
50
                            return false;
51
                        }
52
                    } else {
53
                        if (L(w, l) || L(w, Line(u, v))) {
54
                            return false;
55
                        }
56
                    }
57
                }
58
            }
59
        }
60
    }
61
    return true;
62
}

任意多边形的面积#

1
template<class T> ld area(vector<Point<T>> P) {
2
    int n = P.size();
3
    ld ans = 0;
4
    for (int i = 0; i < n; i++) {
5
        ans += cross(P[i], P[(i + 1) % n]);
6
    }
7
    return ans / 2.0;
8
}

皮克定理#

绘制在方格纸上的多边形面积公式可以表示为 $S=n+\dfrac{s}{2}-1$ ，其中 $n$ 表示多边形内部的点数、 $s$ 表示多边形边界上的点数。一条线段上的点数为 $\gcd(|x_1-x_2|,|y_1-y_2|)+1$ 。

任意多边形上/内的网格点个数（仅能处理整数）#

皮克定理用。

1
int onPolygonGrid(vector<Point<int>> p) { // 多边形上
2
    int n = p.size(), ans = 0;
3
    for (int i = 0; i < n; i++) {
4
        auto a = p[i], b = p[(i + 1) % n];
5
        ans += gcd(abs(a.x - b.x), abs(a.y - b.y));
6
    }
7
    return ans;
8
}
9
int inPolygonGrid(vector<Point<int>> p) { // 多边形内
10
    int n = p.size(), ans = 0;
11
    for (int i = 0; i < n; i++) {
12
        auto a = p[i], b = p[(i + 1) % n], c = p[(i + 2) % n];
13
        ans += b.y * (a.x - c.x);
14
    }
15
    ans = abs(ans);
16
    return (ans - onPolygonGrid(p)) / 2 + 1;
17
}

二维凸包#

获取二维静态凸包（Andrew算法）#

flag 用于判定凸包边上的点、重复的顶点是否要加入到凸包中，为 $0$ 时代表加入凸包（不严格）；为 $1$ 时不加入凸包（严格）。时间复杂度为 $\mathcal O(N\log N)$ 。

1
template<class T> vector<Point<T>> staticConvexHull(vector<Point<T>> A, int flag = 1) {
2
    int n = A.size();
3
    if (n <= 2) { // 特判
4
        return A;
5
    }
6
    vector<Point<T>> ans(n * 2);
7
    sort(A.begin(), A.end());
8
    int now = -1;
9
    for (int i = 0; i < n; i++) { // 维护下凸包
10
        while (now > 0 && cross(A[i], ans[now], ans[now - 1]) <= 0) {
11
            now--;
12
        }
13
        ans[++now] = A[i];
14
    }
15
    int pre = now;
16
    for (int i = n - 2; i >= 0; i--) { // 维护上凸包
17
        while (now > pre && cross(A[i], ans[now], ans[now - 1]) <= 0) {
18
            now--;
19
        }
20
        ans[++now] = A[i];
21
    }
22
    ans.resize(now);
23
    return ans;
24
}

二维动态凸包#

固定为 int 型，需要重新书写 Line 函数，cmp 用于判定边界情况。可以处理如下两个要求：

动态插入点 $(x,y)$ 到当前凸包中；
判断点 $(x,y)$ 是否在凸包上或是在内部（包括边界）。

1
template<class T> bool turnRight(Pt a, Pt b) {
2
    return cross(a, b) < 0 || (cross(a, b) == 0 && dot(a, b) < 0);
3
}
4
struct Line {
5
    static int cmp;
6
    mutable Point<int> a, b;
7
    friend bool operator<(Line x, Line y) {
8
        return cmp ? x.a < y.a : turnRight(x.b, y.b);
9
    }
10
    friend auto &operator<<(ostream &os, Line l) {
11
        return os << "<" << l.a << ", " << l.b << ">";
12
    }
13
};
14

15
int Line::cmp = 1;
16
struct UpperConvexHull : set<Line> {
17
    bool contains(const Point<int> &p) const {
18
        auto it = lower_bound({p, 0});
19
        if (it != end() && it->a == p) return true;
20
        if (it != begin() && it != end() && cross(prev(it)->b, p - prev(it)->a) <= 0) {
21
            return true;
22
        }
23
        return false;
24
    }
25
    void add(const Point<int> &p) {
26
        if (contains(p)) return;
27
        auto it = lower_bound({p, 0});
28
        for (; it != end(); it = erase(it)) {
29
            if (turnRight(it->a - p, it->b)) {
30
                break;
31
            }
32
        }
33
        for (; it != begin() && prev(it) != begin(); erase(prev(it))) {
34
            if (turnRight(prev(prev(it))->b, p - prev(prev(it))->a)) {
35
                break;
36
            }
37
        }
38
        if (it != begin()) {
39
            prev(it)->b = p - prev(it)->a;
40
        }
41
        if (it == end()) {
42
            insert({p, {0, -1}});
43
        } else {
44
            insert({p, it->a - p});
45
        }
46
    }
47
};
48
struct ConvexHull {
49
    UpperConvexHull up, low;
50
    bool empty() const {
51
        return up.empty();
52
    }
53
    bool contains(const Point<int> &p) const {
54
        Line::cmp = 1;
55
        return up.contains(p) && low.contains(-p);
56
    }
57
    void add(const Point<int> &p) {
58
        Line::cmp = 1;
59
        up.add(p);
60
        low.add(-p);
61
    }
62
    bool isIntersect(int A, int B, int C) const {
63
        Line::cmp = 0;
64
        if (empty()) return false;
65
        Point<int> k = {-B, A};
66
        if (k.x < 0) k = -k;
67
        if (k.x == 0 && k.y < 0) k.y = -k.y;
68
        Point<int> P = up.upper_bound({{0, 0}, k})->a;
69
        Point<int> Q = -low.upper_bound({{0, 0}, k})->a;
70
        return sign(A * P.x + B * P.y - C) * sign(A * Q.x + B * Q.y - C) > 0;
71
    }
72
    friend ostream &operator<<(ostream &out, const ConvexHull &ch) {
73
        for (const auto &line : ch.up) out << "(" << line.a.x << "," << line.a.y << ")";
74
        cout << "/";
75
        for (const auto &line : ch.low) out << "(" << -line.a.x << "," << -line.a.y << ")";
76
        return out;
77
    }
78
};

点与凸包的位置关系#

$0$ 代表点在凸包外面； $1$ 代表在凸壳上； $2$ 代表在凸包内部。

1
template<class T> int contains(Point<T> p, vector<Point<T>> A) {
2
    int n = A.size();
3
    bool in = false;
4
    for (int i = 0; i < n; i++) {
5
        Point<T> a = A[i] - p, b = A[(i + 1) % n] - p;
6
        if (a.y > b.y) {
7
            swap(a, b);
8
        }
9
        if (a.y <= 0 && 0 < b.y && cross(a, b) < 0) {
10
            in = !in;
11
        }
12
        if (cross(a, b) == 0 && dot(a, b) <= 0) {
13
            return 1;
14
        }
15
    }
16
    return in ? 2 : 0;
17
}

闵可夫斯基和#

计算两个凸包合成的大凸包。

1
template<class T> vector<Point<T>> mincowski(vector<Point<T>> P1, vector<Point<T>> P2) {
2
    int n = P1.size(), m = P2.size();
3
    vector<Point<T>> V1(n), V2(m);
4
    for (int i = 0; i < n; i++) {
5
        V1[i] = P1[(i + 1) % n] - P1[i];
6
    }
7
    for (int i = 0; i < m; i++) {
8
        V2[i] = P2[(i + 1) % m] - P2[i];
9
    }
10
    vector<Point<T>> ans = {P1.front() + P2.front()};
11
    int t = 0, i = 0, j = 0;
12
    while (i < n && j < m) {
13
        Point<T> val = sign(cross(V1[i], V2[j])) > 0 ? V1[i++] : V2[j++];
14
        ans.push_back(ans.back() + val);
15
    }
16
    while (i < n) ans.push_back(ans.back() + V1[i++]);
17
    while (j < m) ans.push_back(ans.back() + V2[j++]);
18
    return ans;
19
}

半平面交#

计算多条直线左边平面部分的交集。

1
template<class T> vector<Point<T>> halfcut(vector<Line<T>> lines) {
2
    sort(lines.begin(), lines.end(), [&](auto l1, auto l2) {
3
        auto d1 = l1.b - l1.a;
4
        auto d2 = l2.b - l2.a;
5
        if (sign(d1) != sign(d2)) {
6
            return sign(d1) == 1;
7
        }
8
        return cross(d1, d2) > 0;
9
    });
10
    deque<Line<T>> ls;
11
    deque<Point<T>> ps;
12
    for (auto l : lines) {
13
        if (ls.empty()) {
14
            ls.push_back(l);
15
            continue;
16
        }
17
        while (!ps.empty() && !pointOnLineLeft(ps.back(), l)) {
18
            ps.pop_back();
19
            ls.pop_back();
20
        }
21
        while (!ps.empty() && !pointOnLineLeft(ps[0], l)) {
22
            ps.pop_front();
23
            ls.pop_front();
24
        }
25
        if (cross(l.b - l.a, ls.back().b - ls.back().a) == 0) {
26
            if (dot(l.b - l.a, ls.back().b - ls.back().a) > 0) {
27
                if (!pointOnLineLeft(ls.back().a, l)) {
28
                    assert(ls.size() == 1);
29
                    ls[0] = l;
30
                }
31
                continue;
32
            }
33
            return {};
34
        }
35
        ps.push_back(lineIntersection(ls.back(), l));
36
        ls.push_back(l);
37
    }
38
    while (!ps.empty() && !pointOnLineLeft(ps.back(), ls[0])) {
39
        ps.pop_back();
40
        ls.pop_back();
41
    }
42
    if (ls.size() <= 2) {
43
        return {};
44
    }
45
    ps.push_back(lineIntersection(ls[0], ls.back()));
46
    return vector(ps.begin(), ps.end());
47
}

三维几何必要初始化#

点线面封装#

1
struct Point3 {
2
    ld x, y, z;
3
    Point3(ld x_ = 0, ld y_ = 0, ld z_ = 0) : x(x_), y(y_), z(z_) {}
4
    Point3 &operator+=(Point3 p) & {
5
        return x += p.x, y += p.y, z += p.z, *this;
6
    }
7
    Point3 &operator-=(Point3 p) & {
8
        return x -= p.x, y -= p.y, z -= p.z, *this;
9
    }
10
    Point3 &operator*=(Point3 p) & {
11
        return x *= p.x, y *= p.y, z *= p.z, *this;
12
    }
13
    Point3 &operator*=(ld t) & {
14
        return x *= t, y *= t, z *= t, *this;
15
    }
16
    Point3 &operator/=(ld t) & {
17
        return x /= t, y /= t, z /= t, *this;
18
    }
19
    friend Point3 operator+(Point3 a, Point3 b) { return a += b; }
20
    friend Point3 operator-(Point3 a, Point3 b) { return a -= b; }
21
    friend Point3 operator*(Point3 a, Point3 b) { return a *= b; }
22
    friend Point3 operator*(Point3 a, ld b) { return a *= b; }
23
    friend Point3 operator*(ld a, Point3 b) { return b *= a; }
24
    friend Point3 operator/(Point3 a, ld b) { return a /= b; }
25
    friend auto &operator>>(istream &is, Point3 &p) {
26
        return is >> p.x >> p.y >> p.z;
27
    }
28
    friend auto &operator<<(ostream &os, Point3 p) {
29
        return os << "(" << p.x << ", " << p.y << ", " << p.z << ")";
30
    }
31
};
32
struct Line3 {
33
    Point3 a, b;
34
};
35
struct Plane {
36
    Point3 u, v, w;
37
};

其他函数#

1
ld len(P3 p) { // 原点到当前点的距离计算
2
    return sqrt(p.x * p.x + p.y * p.y + p.z * p.z);
3
}
4
P3 crossEx(P3 a, P3 b) { // 叉乘
5
    P3 ans;
6
    ans.x = a.y * b.z - a.z * b.y;
7
    ans.y = a.z * b.x - a.x * b.z;
8
    ans.z = a.x * b.y - a.y * b.x;
9
    return ans;
10
}
11
ld cross(P3 a, P3 b) {
12
    return len(crossEx(a, b));
13
}
14
ld dot(P3 a, P3 b) { // 点乘
15
    return a.x * b.x + a.y * b.y + a.z * b.z;
16
}
17
P3 getVec(Plane s) { // 获取平面法向量
18
    return crossEx(s.u - s.v, s.v - s.w);
19
}
20
ld dis(P3 a, P3 b) { // 三维欧几里得距离公式
21
    ld val = (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) + (a.z - b.z) * (a.z - b.z);
22
    return sqrt(val);
23
}
24
P3 standardize(P3 vec) { // 将三维向量转换为单位向量
25
    return vec / len(vec);
26
}

三维点线面相关#

空间三点是否共线#

其中第二个函数是专门用来判断给定的三个点能否构成平面的，因为不共线的三点才能构成平面。

1
bool onLine(P3 p1, P3 p2, P3 p3) { // 三点是否共线
2
    return sign(cross(p1 - p2, p3 - p2)) == 0;
3
}
4
bool onLine(Plane s) {
5
    return onLine(s.u, s.v, s.w);
6
}

四点是否共面#

1
bool onPlane(P3 p1, P3 p2, P3 p3, P3 p4) { // 四点是否共面
2
    ld val = dot(getVec({p1, p2, p3}), p4 - p1);
3
    return sign(val) == 0;
4
}

空间点是否在线段上#

1
bool pointOnSegment(P3 p, L3 l) {
2
    return sign(cross(p - l.a, p - l.b)) == 0 && min(l.a.x, l.b.x) <= p.x &&
3
           p.x <= max(l.a.x, l.b.x) && min(l.a.y, l.b.y) <= p.y && p.y <= max(l.a.y, l.b.y) &&
4
           min(l.a.z, l.b.z) <= p.z && p.z <= max(l.a.z, l.b.z);
5
}
6
bool pointOnSegmentEx(P3 p, L3 l) { // pointOnSegment去除端点版
7
    return sign(cross(p - l.a, p - l.b)) == 0 && min(l.a.x, l.b.x) < p.x &&
8
           p.x < max(l.a.x, l.b.x) && min(l.a.y, l.b.y) < p.y && p.y < max(l.a.y, l.b.y) &&
9
           min(l.a.z, l.b.z) < p.z && p.z < max(l.a.z, l.b.z);
10
}

空间两点是否在线段同侧#

当给定的两点与线段不共面、点在线段上时返回 $false$ 。

1
bool pointOnSegmentSide(P3 p1, P3 p2, L3 l) {
2
    if (!onPlane(p1, p2, l.a, l.b)) { // 特判不共面
3
        return 0;
4
    }
5
    ld val = dot(crossEx(l.a - l.b, p1 - l.b), crossEx(l.a - l.b, p2 - l.b));
6
    return sign(val) == 1;
7
}

两点是否在平面同侧#

点在平面上时返回 $false$ 。

1
bool pointOnPlaneSide(P3 p1, P3 p2, Plane s) {
2
    ld val = dot(getVec(s), p1 - s.u) * dot(getVec(s), p2 - s.u);
3
    return sign(val) == 1;
4
}

空间两直线是否平行/垂直#

1
bool lineParallel(L3 l1, L3 l2) {
2
    return sign(cross(l1.a - l1.b, l2.a - l2.b)) == 0;
3
}
4
bool lineVertical(L3 l1, L3 l2) {
5
    return sign(dot(l1.a - l1.b, l2.a - l2.b)) == 0;
6
}

两平面是否平行/垂直#

1
bool planeParallel(Plane s1, Plane s2) {
2
    ld val = cross(getVec(s1), getVec(s2));
3
    return sign(val) == 0;
4
}
5
bool planeVertical(Plane s1, Plane s2) {
6
    ld val = dot(getVec(s1), getVec(s2));
7
    return sign(val) == 0;
8
}

空间两直线是否是同一条#

1
bool same(L3 l1, L3 l2) {
2
    return lineParallel(l1, l2) && lineParallel({l1.a, l2.b}, {l1.b, l2.a});
3
}

两平面是否是同一个#

1
bool same(Plane s1, Plane s2) {
2
    return onPlane(s1.u, s2.u, s2.v, s2.w) && onPlane(s1.v, s2.u, s2.v, s2.w) &&
3
           onPlane(s1.w, s2.u, s2.v, s2.w);
4
}

直线是否与平面平行#

1
bool linePlaneParallel(L3 l, Plane s) {
2
    ld val = dot(l.a - l.b, getVec(s));
3
    return sign(val) == 0;
4
}

空间两线段是否相交#

1
bool segmentIntersection(L3 l1, L3 l2) { // 重叠、相交于端点均视为相交
2
    if (!onPlane(l1.a, l1.b, l2.a, l2.b)) { // 特判不共面
3
        return 0;
4
    }
5
    if (!onLine(l1.a, l1.b, l2.a) || !onLine(l1.a, l1.b, l2.b)) {
6
        return !pointOnSegmentSide(l1.a, l1.b, l2) && !pointOnSegmentSide(l2.a, l2.b, l1);
7
    }
8
    return pointOnSegment(l1.a, l2) || pointOnSegment(l1.b, l2) || pointOnSegment(l2.a, l1) ||
9
           pointOnSegment(l2.b, l2);
10
}
11
bool segmentIntersection1(L3 l1, L3 l2) { // 重叠、相交于端点不视为相交
12
    return onPlane(l1.a, l1.b, l2.a, l2.b) && !pointOnSegmentSide(l1.a, l1.b, l2) &&
13
           !pointOnSegmentSide(l2.a, l2.b, l1);
14
}

空间两直线是否相交及交点#

当两直线不共面、两直线平行时返回 $false$ 。

1
pair<bool, P3> lineIntersection(L3 l1, L3 l2) {
2
    if (!onPlane(l1.a, l1.b, l2.a, l2.b) || lineParallel(l1, l2)) {
3
        return {0, {}};
4
    }
5
    auto [s1, e1] = l1;
6
    auto [s2, e2] = l2;
7
    ld val = 0;
8
    if (!onPlane(l1.a, l1.b, {0, 0, 0}, {0, 0, 1})) {
9
        val = ((s1.x - s2.x) * (s2.y - e2.y) - (s1.y - s2.y) * (s2.x - e2.x)) /
10
              ((s1.x - e1.x) * (s2.y - e2.y) - (s1.y - e1.y) * (s2.x - e2.x));
11
    } else if (!onPlane(l1.a, l1.b, {0, 0, 0}, {0, 1, 0})) {
12
        val = ((s1.x - s2.x) * (s2.z - e2.z) - (s1.z - s2.z) * (s2.x - e2.x)) /
13
              ((s1.x - e1.x) * (s2.z - e2.z) - (s1.z - e1.z) * (s2.x - e2.x));
14
    } else {
15
        val = ((s1.y - s2.y) * (s2.z - e2.z) - (s1.z - s2.z) * (s2.y - e2.y)) /
16
              ((s1.y - e1.y) * (s2.z - e2.z) - (s1.z - e1.z) * (s2.y - e2.y));
17
    }
18
    return {1, s1 + (e1 - s1) * val};
19
}

直线与平面是否相交及交点#

当直线与平面平行、给定的点构不成平面时返回 $false$ 。

1
pair<bool, P3> linePlaneCross(L3 l, Plane s) {
2
    if (linePlaneParallel(l, s)) {
3
        return {0, {}};
4
    }
5
    P3 vec = getVec(s);
6
    P3 U = vec * (s.u - l.a), V = vec * (l.b - l.a);
7
    ld val = (U.x + U.y + U.z) / (V.x + V.y + V.z);
8
    return {1, l.a + (l.b - l.a) * val};
9
}

两平面是否相交及交线#

当两平面平行、两平面为同一个时返回 $false$ 。

1
pair<bool, L3> planeIntersection(Plane s1, Plane s2) {
2
    if (planeParallel(s1, s2) || same(s1, s2)) {
3
        return {0, {}};
4
    }
5
    P3 U = linePlaneParallel({s2.u, s2.v}, s1) ? linePlaneCross({s2.v, s2.w}, s1).second
6
                                               : linePlaneCross({s2.u, s2.v}, s1).second;
7
    P3 V = linePlaneParallel({s2.w, s2.u}, s1) ? linePlaneCross({s2.v, s2.w}, s1).second
8
                                               : linePlaneCross({s2.w, s2.u}, s1).second;
9
    return {1, {U, V}};
10
}

点到直线的最近点与最近距离#

1
pair<ld, P3> pointToLine(P3 p, L3 l) {
2
    ld val = cross(p - l.a, l.a - l.b) / dis(l.a, l.b); // 面积除以底边长
3
    ld val1 = dot(p - l.a, l.a - l.b) / dis(l.a, l.b);
4
    return {val, l.a + val1 * standardize(l.a - l.b)};
5
}

点到平面的最近点与最近距离#

1
pair<ld, P3> pointToPlane(P3 p, Plane s) {
2
    P3 vec = getVec(s);
3
    ld val = dot(vec, p - s.u);
4
    val = abs(val) / len(vec); // 面积除以底边长
5
    return {val, p - val * standardize(vec)};
6
}

空间两直线的最近距离与最近点对#

1
tuple<ld, P3, P3> lineToLine(L3 l1, L3 l2) {
2
    P3 vec = crossEx(l1.a - l1.b, l2.a - l2.b); // 计算同时垂直于两直线的向量
3
    ld val = abs(dot(l1.a - l2.a, vec)) / len(vec);
4
    P3 U = l1.b - l1.a, V = l2.b - l2.a;
5
    vec = crossEx(U, V);
6
    ld p = dot(vec, vec);
7
    ld t1 = dot(crossEx(l2.a - l1.a, V), vec) / p;
8
    ld t2 = dot(crossEx(l2.a - l1.a, U), vec) / p;
9
    return {val, l1.a + (l1.b - l1.a) * t1, l2.a + (l2.b - l2.a) * t2};
10
}

三维角度与弧度#

空间两直线夹角的 cos 值#

任意位置的空间两直线。

1
ld lineCos(L3 l1, L3 l2) {
2
    return dot(l1.a - l1.b, l2.a - l2.b) / len(l1.a - l1.b) / len(l2.a - l2.b);
3
}

空间两平面夹角的 cos 值#

1
ld planeCos(Plane s1, Plane s2) {
2
    P3 U = getVec(s1), V = getVec(s2);
3
    return dot(U, V) / len(U) / len(V);
4
}

直线与平面夹角的 sin 值#

1
ld linePlaneSin(L3 l, Plane s) {
2
    P3 vec = getVec(s);
3
    return dot(l.a - l.b, vec) / len(l.a - l.b) / len(vec);
4
}

空间多边形#

正N棱锥体积公式#

棱锥通用体积公式 $V=\dfrac{1}{3}Sh$ ，当其恰好是棱长为 $l$ 的正 $n$ 棱锥时，有公式 $\displaystyle V=\frac{l^3\cdot n}{12\tan \frac{\pi}{n}}\cdot\sqrt{1-\frac{1}{4\cdot \sin^2\frac{\pi}{n}}}$ 。

1
ld V(ld l, int n) { // 正n棱锥体积公式
2
    return l * l * l * n / (12 * tan(PI / n)) * sqrt(1 - 1 / (4 * sin(PI / n) * sin(PI / n)));
3
}

四面体体积#

1
ld V(P3 a, P3 b, P3 c, P3 d) {
2
    return abs(dot(d - a, crossEx(b - a, c - a))) / 6;
3
}

点是否在空间三角形上#

点位于边界上时返回 $false$ 。

1
bool pointOnTriangle(P3 p, P3 p1, P3 p2, P3 p3) {
2
    return pointOnSegmentSide(p, p1, {p2, p3}) && pointOnSegmentSide(p, p2, {p1, p3}) &&
3
           pointOnSegmentSide(p, p3, {p1, p2});
4
}

线段是否与空间三角形相交及交点#

只有交点在空间三角形内部时才视作相交。

1
pair<bool, P3> segmentOnTriangle(P3 l, P3 r, P3 p1, P3 p2, P3 p3) {
2
    P3 x = crossEx(p2 - p1, p3 - p1);
3
    if (sign(dot(x, r - l)) == 0) {
4
        return {0, {}};
5
    }
6
    ld t = dot(x, p1 - l) / dot(x, r - l);
7
    if (t < 0 || t - 1 > 0) { // 不在线段上
8
        return {0, {}};
9
    }
10
    bool type = pointOnTriangle(l + (r - l) * t, p1, p2, p3);
11
    if (type) {
12
        return {1, l + (r - l) * t};
13
    } else {
14
        return {0, {}};
15
    }
16
}

空间三角形是否相交#

相交线段在空间三角形内部时才视作相交。

1
bool triangleIntersection(vector<P3> a, vector<P3> b) {
2
    for (int i = 0; i < 3; i++) {
3
        if (segmentOnTriangle(b[i], b[(i + 1) % 3], a[0], a[1], a[2]).first) {
4
            return 1;
5
        }
6
        if (segmentOnTriangle(a[i], a[(i + 1) % 3], b[0], b[1], b[2]).first) {
7
            return 1;
8
        }
9
    }
10
    return 0;
11
}

常用结论#

平面几何结论归档#

hypot 函数可以直接计算直角三角形的斜边长；
边心距是指正多边形的外接圆圆心到正多边形某一边的距离，边长为 $s$ 的正 $n$ 角形的边心距公式为 $\displaystyle a=\frac{t}{2\cdot\tan \frac{\pi}{n}}$ ，外接圆半径为 $R$ 的正 $n$ 角形的边心距公式为 $a=R\cdot \cos \dfrac{\pi}{n}$ ；
三角形外接圆半径为 $\dfrac{a}{2\sin A}=\dfrac{abc}{4S}$ ，其中 $S$ 为三角形面积，内切圆半径为 $\dfrac{2S}{a+b+c}$ ；
由小正三角形拼成的大正三角形，耗费的小三角形数量即为构成一条边的小三角形数量的平方。如下图，总数量即为 $4^2$ See。
正 $n$ 边形圆心角为 $\dfrac{360^{\circ}}{n}$ ，圆周角为 $\dfrac{180^{\circ}}{n}$ 。定义正 $n$ 边形上的三个顶点 $A,B$ 和 $C$ （可以不相邻），使得 $\angle ABC=\theta$ ，当 $n\le 360$ 时， $\theta$ 可以取 $1^{\circ}$ 到 $179^{\circ}$ 间的任何一个整数 See。
某一点 $B$ 到直线 $AC$ 的距离公式为 $\dfrac{|\vec{BA}\times \vec{BC}|}{|AC|}$ ，等价于 $\dfrac{|aX+bY+c|}{\sqrt{a^2+b^2}}$ 。
atan(y / x) 函数仅用于计算第一、四象限的值，而 atan2(y, x) 则允许计算所有四个象限的正反切，在使用这个函数时，需要尽量保证 $x$ 和 $y$ 的类型为整数型，如果使用浮点数，实测会慢十倍。
在平面上有奇数个点 $A_0,A_1,\dots,A_n$ 以及一个点 $X_0$ ，构造 $X_1$ 使得 $X_0,X_1$ 关于 $A_0$ 对称、构造 $X_2$ 使得 $X_1,X_2$ 关于 $A_1$ 对称、……、构造 $X_j$ 使得 $X_{j-1},X_j$ 关于 $A_{(j-1)\mod n}$ 对称。那么周期为 $2n$ ，即 $A_0$ 与 $A_{2n}$ 共点、 $A_1$ 与 $A_{2n+1}$ 共点 See 。
已知 $A\ (x_A, y_A)$ 和 $X\ (x_X,y_X)$ 两点及这两点的坐标，构造 $Y$ 使得 $X,Y$ 关于 $A$ 对称，那么 $Y$ 的坐标为 $(2\cdot x_A-x_X,2\cdot y_A-y_X)$ 。
海伦公式：已知三角形三边长 $a,b$ 和 $c$ ，定义 $p=\dfrac{a+b+c}{2}$ ，则 $S_{\triangle}=\sqrt{p(p-a)(p-b)(p-c)}$ ，在使用时需要注意越界问题，本质是铅锤定理，一般多使用叉乘计算三角形面积而不使用该公式。
棱台体积 $V=\frac{1}{3}(S_1+S_2+\sqrt{S_1S_2})\cdot h$ ，其中 $S_1,S_2$ 为上下底面积。
正棱台侧面积 $\frac{1}{2}(C_1+C_2)\cdot L$ ，其中 $C_1,C_2$ 为上下底周长， $L$ 为斜高（上下底对应的平行边的距离）。
球面积 $4\pi r^2$ ，体积 $\frac{4}{3}\pi r^3$ 。
正三角形面积 $\dfrac{\sqrt 3 a^2}{4}$ ，正四面体面积 $\dfrac{\sqrt 2 a^3}{12}$ 。
设扇形对应的圆心角弧度为 $\theta$ ，则面积为 $S=\frac{\theta}{2}\cdot R^2$ 。

立体几何结论归档#

已知向量 $\vec{r}=\{x,y,z\}$ ，则该向量的三个方向余弦为 $\cos \alpha =\dfrac{x}{|\vec r|}=\dfrac{x}{\sqrt{x^2+y^2+z^2}}; \ \cos \beta = \dfrac{y}{|\vec r|};\ \cos \gamma =\dfrac{z}{|\vec r|}$ 。其中 $\alpha,\beta,\gamma\in [0,\pi]$ ， $\cos^2\alpha+\cos^2\beta+\cos^2\gamma=1$ 。

常用例题#

将平面某点旋转任意角度#

题意：给定平面上一点 $(a,b)$ ，输出将其逆时针旋转 $d$ 度之后的坐标。

1
signed main() {
2
    int a, b, d;
3
    cin >> a >> b >> d;
4

5
    ld l = hypot(a, b); // 库函数，求直角三角形的斜边
6
    ld alpha = atan2(b, a) + toArc(d);
7

8
    cout << l * cos(alpha) << " " << l * sin(alpha) << endl;
9
}

平面最近点对（set解）#

借助 set ，在严格 $\mathcal O(N\log N)$ 复杂度内求解，比常见的分治法稍快。

1
template<class T> T sqr(T x) {
2
    return x * x;
3
}
4

5
using V = Point<int>;
6
signed main() {
7
    int n;
8
    cin >> n;
9

10
    vector<V> in(n);
11
    for (auto &it : in) {
12
        cin >> it;
13
    }
14

15
    int dis = disEx(in[0], in[1]); // 设定阈值
16
    sort(in.begin(), in.end());
17

18
    set<V> S;
19
    for (int i = 0, h = 0; i < n; i++) {
20
        V now = {in[i].y, in[i].x};
21
        while (dis && dis <= sqr(in[i].x - in[h].x)) { // 删除超过阈值的点
22
            S.erase({in[h].y, in[h].x});
23
            h++;
24
        }
25
        auto it = S.lower_bound(now);
26
        for (auto k = it; k != S.end() && sqr(k->x - now.x) < dis; k++) {
27
            dis = min(dis, disEx(*k, now));
28
        }
29
        if (it != S.begin()) {
30
            for (auto k = prev(it); sqr(k->x - now.x) < dis; k--) {
31
                dis = min(dis, disEx(*k, now));
32
                if (k == S.begin()) break;
33
            }
34
        }
35
        S.insert(now);
36
    }
37
    cout << sqrt(dis) << endl;
38
}

平面若干点能构成的最大四边形的面积（简单版，暴力枚举）#

题意：平面上存在若干个点，保证没有两点重合、没有三点共线，你需要从中选出四个点，使得它们构成的四边形面积是最大的，注意这里能组成的四边形可以不是凸四边形。

暴力枚举其中一条对角线后枚举剩余两个点， $\mathcal O(N^3)$ 。

1
signed main() {
2
    int n;
3
    cin >> n;
4
    vector<Pi> in(n);
5
    for (auto &it : in) {
6
        cin >> it;
7
    }
8
    ld ans = 0;
9
    for (int i = 0; i < n; i++) {
10
        for (int j = i + 1; j < n; j++) { // 枚举对角线
11
            ld l = 0, r = 0;
12
            for (int k = 0; k < n; k++) { // 枚举第三点
13
                if (k == i || k == j) continue;
14
                if (pointOnLineLeft(in[k], {in[i], in[j]})) {
15
                    l = max(l, triangleS(in[k], in[j], in[i]));
16
                } else {
17
                    r = max(r, triangleS(in[k], in[j], in[i]));
18
                }
19
            }
20
            if (l * r != 0) { // 确保构成的是四边形
21
                ans = max(ans, l + r);
22
            }
23
        }
24
    }
25
    cout << ans << endl;
26
}

平面若干点能构成的最大四边形的面积（困难版，分类讨论+旋转卡壳）#

题意：平面上存在若干个点，可能存在多点重合、共线的情况，你需要从中选出四个点，使得它们构成的四边形面积是最大的，注意这里能组成的四边形可以不是凸四边形、可以是退化的四边形。

当凸包大小 $\le 2$ 时，说明是退化的四边形，答案直接为 $0$ ；大小恰好为 $3$ 时，说明是凹四边形，我们枚举不在凸包上的那一点，将两个三角形面积相减既可得到答案；大小恰好为 $4$ 时，说明是凸四边形，使用旋转卡壳求解。

1
using V = Point<int>;
2
signed main() {
3
    int Task = 1;
4
    for (cin >> Task; Task; Task--) {
5
        int n;
6
        cin >> n;
7

8
        vector<V> in_(n);
9
        for (auto &it : in_) {
10
            cin >> it;
11
        }
12
        auto in = staticConvexHull(in_, 0);
13
        n = in.size();
14

15
        int ans = 0;
16
        if (n > 3) {
17
            ans = rotatingCalipers(in);
18
        } else if (n == 3) {
19
            int area = triangleAreaEx(in[0], in[1], in[2]);
20
            for (auto it : in_) {
21
                if (it == in[0] || it == in[1] || it == in[2]) continue;
22
                int Min = min({triangleAreaEx(it, in[0], in[1]), triangleAreaEx(it, in[0], in[2]), triangleAreaEx(it, in[1], in[2])});
23
                ans = max(ans, area - Min);
24
            }
25
        }
26

27
        cout << ans / 2;
28
        if (ans % 2) {
29
            cout << ".5";
30
        }
31
        cout << endl;
32
    }
33
}

线段将多边形切割为几个部分#

题意：给定平面上一线段与一个任意多边形，求解线段将多边形切割为几个部分；保证线段的端点不在多边形内、多边形边上，多边形顶点不位于线段上，多边形的边不与线段重叠；多边形端点按逆时针顺序给出。下方的几个样例均合法，答案均为 $3$ 。

当线段切割多边形时，本质是与多边形的边交于两个点、或者说是与多边形的两条边相交，设交点数目为 $x$ ，那么答案即为 $\frac{x}{2}+1$ 。于是，我们只需要计算交点数量即可，先判断某一条边是否与线段相交，再判断边的两个端点是否位于线段两侧。

1
signed main() {
2
    Pi s, e;
3
    cin >> s >> e; // 读入线段
4

5
    int n;
6
    cin >> n;
7
    vector<Pi> in(n);
8
    for (auto &it : in) {
9
        cin >> it; // 读入多边形端点
10
    }
11

12
    int cnt = 0;
13
    for (int i = 0; i < n; i++) {
14
        Pi x = in[i], y = in[(i + 1) % n];
15
        cnt += (pointNotOnLineSide(x, y, {s, e}) && segmentIntersection(Line{x, y}, {s, e}));
16
    }
17
    cout << cnt / 2 + 1 << endl;
18
}

平面若干点能否构成凸包（暴力枚举）#

题意：给定平面上若干个点，判断其是否构成凸包 See 。

可以直接使用凸包模板，但是代码较长；在这里我们使用暴力枚举试点，也能以 $\mathcal O(N)$ 的复杂度通过。当两个向量的叉乘 $\le0$ 时说明其夹角大于等于 $180\degree$ ，使用这一点即可判定。

1
signed main() {
2
    int n;
3
    cin >> n;
4

5
    vector<Point<ld>> in(n);
6
    for (auto &it : in) {
7
        cin >> it;
8
    }
9

10
    for (int i = 0; i < n; i++) {
11
        auto A = in[(i - 1 + n) % n];
12
        auto B = in[i];
13
        auto C = in[(i + 1) % n];
14
        if (cross(A - B, C - B) > 0) {
15
            cout << "No\n";
16
            return 0;
17
        }
18
    }
19
    cout << "Yes\n";
20
}

凸包上的点能构成的最大三角形（暴力枚举）#

可以直接使用凸包模板，但是代码较长；在这里我们使用暴力枚举试点，也能以 $\mathcal O(N)$ 的复杂度通过。

另外补充一点性质：所求三角形的反互补三角形一定包含了凸包上的所有点（可以在边界）。通俗的说，构成的三角形是这个反互补三角形的中点三角形。如下图所示，点 $A$ 不在 $\triangle BCE$ 的反互补三角形内部，故 $\triangle BCE$ 不是最大三角形； $\triangle ACE$ 才是。

1
signed main() {
2
    int n;
3
    cin >> n;
4

5
    vector<Point<int>> in(n);
6
    for (auto &it : in) {
7
        cin >> it;
8
    }
9

10
    #define S(x, y, z) triangleAreaEx(in[x], in[y], in[z])
11

12
    int i = 0, j = 1, k = 2;
13
    while (true) {
14
        int val = S(i, j, k);
15
        if (S((i + 1) % n, j, k) > val) {
16
            i = (i + 1) % n;
17
        } else if (S((i - 1 + n) % n, j, k) > val) {
18
            i = (i - 1 + n) % n;
19
        } else if (S(i, (j + 1) % n, k) > val) {
20
            j = (j + 1) % n;
21
        } else if (S(i, (j - 1 + n) % n, k) > val) {
22
            j = (j - 1 + n) % n;
23
        } else if (S(i, j, (k + 1) % n) > val) {
24
            k = (k + 1) % n;
25
        } else if (S(i, j, (k - 1 + n) % n) > val) {
26
            k = (k - 1 + n) % n;
27
        } else {
28
            break;
29
        }
30
    }
31
    cout << i + 1 << " " << j + 1 << " " << k + 1 << endl;
32
}

凸包上的点能构成的最大四角形的面积（旋转卡壳）#

由于是凸包上的点，所以保证了四边形一定是凸四边形，时间复杂度 $\mathcal O(N^2)$ 。

1
template<class T> T rotatingCalipers(vector<Point<T>> &p) {
2
    #define S(x, y, z) triangleAreaEx(p[x], p[y], p[z])
3
    int n = p.size();
4
    T ans = 0;
5
    auto nxt = [&](int i) -> int {
6
        return i == n - 1 ? 0 : i + 1;
7
    };
8
    for (int i = 0; i < n; i++) {
9
        int p1 = nxt(i), p2 = nxt(nxt(nxt(i)));
10
        for (int j = nxt(nxt(i)); nxt(j) != i; j = nxt(j)) {
11
            while (nxt(p1) != j && S(i, j, nxt(p1)) > S(i, j, p1)) {
12
                p1 = nxt(p1);
13
            }
14
            if (p2 == j) {
15
                p2 = nxt(p2);
16
            }
17
            while (nxt(p2) != i && S(i, j, nxt(p2)) > S(i, j, p2)) {
18
                p2 = nxt(p2);
19
            }
20
            ans = max(ans, S(i, j, p1) + S(i, j, p2));
21
        }
22
    }
23
    return ans;
24
    #undef S
25
}

判断一个凸包是否完全在另一个凸包内#

题意：给定一个凸多边形 $A$ 和一个凸多边形 $B$ ，询问 $B$ 是否被 $A$ 包含，分别判断严格/不严格包含。例题。

考虑严格包含，使用 $A$ 点集计算出凸包 $T_1$ ，使用 $A,B$ 两个点集计算出不严格凸包 $T_2$ ，如果包含，那么 $T_1$ 应该与 $T_2$ 完全相等；考虑不严格包含，在计算凸包 $T_2$ 时严格即可。最终以 $\mathcal O(N)$ 复杂度求解，且代码不算很长。

组合数学#

基础#

依赖ModInt, 否则需要手动取模.

1
template <int N, int MOD = 998244353>
2
struct Factorial {
3
    using mint = ModInt<MOD>;
4
    array<mint, N + 5> fact, factinv;
5
    Factorial() {
6
        fact[0] = 1;
7
        for (int i = 1; i <= N; i++) {
8
            fact[i] = fact[i - 1] * i;
9
        }
10
        factinv[N] = fact[N].inv();
11
        for (int i = N - 1; i >= 0; i--) {
12
            factinv[i] = factinv[i + 1] * (i + 1);
13
        }
14
    }
15
    mint nCr(int n, int r) {
16
        if (!(0 <= r && r <= n)) {
17
            return 0;
18
        }
19
        return fact[n] * factinv[r] * factinv[n - r];
20
    }
21
    mint nPr(int n, int r) {
22
        if (!(0 <= r && r <= n)) {
23
            return 0;
24
        }
25
        return fact[n] * factinv[n - r];
26
    }
27
};

模运算#

模运算类#

1
constexpr int MOD = 998'244'353;
2

3
template <int P>
4
struct ModInt {
5
    u32 val;
6
    constexpr ModInt(ll x = 0) : val((x %= P) < 0 ? x + P : x) {}
7
    // 注意i128等类型需要先转成ll
8
    using M = ModInt;
9
    M &operator+=(const M &p) {
10
        if ((val += p.val) >= P) {
11
            val -= P;
12
        }
13
        return *this;
14
    }
15
    M &operator-=(const M &p) {
16
        if ((val += P - p.val) >= P) {
17
            val -= P;
18
        }
19
        return *this;
20
    }
21
    M &operator*=(const M &p) {
22
        val = u64(val) * p.val % P;
23
        return *this;
24
    }
25
    M &operator/=(const M &p) { return (*this) *= p.inv(); }
26
    constexpr M operator-() const { return M(val ? P - val : u32(0)); }
27
    constexpr M operator+(const M &p) const { return M(*this) += p; }
28
    constexpr M operator-(const M &p) const { return M(*this) -= p; }
29
    constexpr M operator*(const M &p) const { return M(*this) *= p; }
30
    constexpr M operator/(const M &p) const { return M(*this) /= p; }
31
    constexpr bool operator==(const M &p) const { return val == p.val; }
32
    constexpr bool operator!=(const M &p) const { return val != p.val; }
33
    M pow(ll n) const {
34
        if (n < 0) return inv().pow(-n);
35
        M ret(1), mul(val);
36
        while (n > 0) {
37
            if (n & 1) ret *= mul;
38
            mul *= mul;
39
            n >>= 1;
40
        }
41
        return ret;
42
    }
43
    M inv() const {
44
        int a = val, b = P, u = 1, v = 0, t;
45
        while (b > 0) {
46
            t = a / b;
47
            a -= t * b;
48
            u -= t * v;
49
            swap(a, b);
50
            swap(u, v);
51
        }
52
        return M(u);
53
    }
54
    constexpr static M raw(u32 x) {
55
        M ret;
56
        ret.val = x;
57
        return ret;
58
    }
59
    constexpr static int mod() {
60
        return P;
61
    }
62
};
63
using mint = ModInt<MOD>;

数论#

exCRT#

1
array<i128, 2> excrt(const vector<ll>& b, const vector<ll>& m) {
2
    int n = b.size();
3
    i128 nowb = b[0], nowm = m[0];
4
    for (int i = 1; i < n; i++) {
5
        ll A = nowm % m[i], B = m[i];
6
        i128 C = b[i] - nowb;
7
        auto [y1, _, g] = exgcd(A, B);
8
        if (C % g != 0) {
9
            return { -1, -1 };
10
        }
11

12
        i128 M = m[i] / g;
13
        i128 y = (i128)y1 * (C / g) % M;
14
        y = (y < 0) ? y + M : y;
15
        nowb += y * nowm;
16
        nowm *= M;
17
        nowb %= nowm;
18
    }
19
    return { nowb, nowm };
20
}

exGCD#

1
i128 abs128(i128 x) {
2
    return x < 0 ? -x : x;
3
}
4
array<i128, 3> exgcd(i128 a, i128 b) {
5
    if (b == 0) {
6
        if (a >= 0) return {1, 0, a};
7
        else return {-1, 0, -a};
8
    }
9
    auto [x, y, g] = exgcd(b, a % b);
10
    return {y, x - a / b * y, g};
11
}
12
i128 norm_mod(i128 x, i128 mod) {
13
    mod = abs128(mod);
14
    x %= mod;
15
    if (x < 0) x += mod;
16
    return x;
17
}
18
tuple<bool, i128, i128> solveLDE(i128 a, i128 b, i128 c) {
19
    // solve ax + by = c
20
    if (a == 0 && b == 0) {
21
        if (c == 0) return {true, 0, 0};
22
        else return {false, 0, 0};
23
    }
24
    auto [x, y, g] = exgcd(a, b);
25
    if (c % g != 0) {
26
        return {false, 0, 0};
27
    }
28
    i128 k = c / g;
29
    x *= k;
30
    y *= k;
31
    // 把 x 规范到 [0, |b / g| - 1]
32
    if (b != 0) {
33
        i128 step_x = b / g;
34
        i128 nx = norm_mod(x, step_x);
35
        // x = x0 + (b / g) * t
36
        // y = y0 - (a / g) * t
37
        i128 t = (nx - x) / step_x;
38
        x = nx;
39
        y = y - (a / g) * t;
40
    }
41
    return {true, x, y};
42
}

线性筛#

1
constexpr int N = 1'000'000;
2
constexpr int PIN = 78'498;
3
bitset<N + 5> not_prime;
4
int primes[PIN + 5], pcnt = 0;
5
// 1e9 : 50'847'534
6
// 1e8 : 5'761'455
7
// 1e7 : 664'579
8
// 1e6 : 78'498
9
// 1e5 : 9'592
10
// 1e4 : 1'229
11
// 1e3 : 168
12
// 1e2 : 25
13
// 1e1 : 4
14
// 下界N / log(N), 上界(>61000时) N / (log(N) - 1.1)",
15

16
// int minp[N + 5]; // 最小质因子
17
// int phi[N + 5]; // 欧拉函数
18
// int mu[N + 5]; // 莫比乌斯函数
19
// int d[N + 5], times[N + 5]; // 约数计数函数
20
// ll sigma[N + 5]; int g[N + 5]; // 约数和函数
21
void sieve() {
22
    // phi[1] = 1;
23
    // mu[1] = 1;
24
    // d[1] = 1;
25
    // sigma[1] = 1;
26
    for (int i = 2; i <= N; i++) {
27
        if (!not_prime[i]) {
28
            primes[pcnt] = i, pcnt++;
29
            // minp[i] = i;
30
            // phi[i] = i - 1;
31
            // mu[i] = -1;
32
            // d[i] = 2, times[i] = 1;
33
            // sigma[i] = i + 1, g[i] = i + 1;
34
        }
35
        for (int j = 0; j < pcnt && i * primes[j] <= N; j++) {
36
            int pj = primes[j];
37
            not_prime[i * pj] = 1;
38
            // minp[i * pj] = pj;
39
            if (i % primes[j] == 0) {
40
                // phi[i * pj] = phi[i] * pj;
41
                // mu[i * pj] = 0;
42
                // d[i * pj] = d[i] / (times[i] + 1) * (times[i] + 2), times[i * pj] = times[i] + 1;
43
                // g[i * pj] = g[i] * pj + 1, sigma[i * pj] = sigma[i] / g[i] * g[i * pj];
44
                break;
45
            }
46
            // phi[i * pj] = phi[i] * (pj - 1);
47
            // mu[i * pj] = -mu[i];
48
            // d[i * pj] = d[i] * 2, times[i * pj] = 1;
49
            // sigma[i * pj] = sigma[i] * sigma[pj], g[i * pj] = pj + 1;
50
        }
51
    }
52
}

分段筛#

1
constexpr ll R = 1'000'000'300;
2
constexpr int SQR = 32'000;
3
constexpr int PIS = 4'000;
4
// sqrt(10) < 3.17
5
bitset<SQR + 5> not_prime;
6
int primes[PIS + 5], pcnt = 0;
7
void init_sieve() {
8
    not_prime[0] = not_prime[1] = 1;
9
    for (ll i = 2; i * i <= R; i++) {
10
        if (!not_prime[i]) {
11
            primes[pcnt] = i;
12
            pcnt++;
13
            for (ll j = i * i; j <= SQR; j += i) {
14
                not_prime[j] = 1;
15
            }
16
        }
17
    }
18
}
19

20
vector<ll> get_primes(ll l, ll r) {
21
    vector<bool> notp(r - l + 1);
22
    if (l == 0) {
23
        notp[0] = notp[1] = 1;
24
    } else if (l == 1) {
25
        notp[0] = 1;
26
    }
27

28
    for (int i = 0; i < pcnt; i++) {
29
        ll p = primes[i];
30
        if (p * p > r) {
31
            break;
32
        }
33
        for (ll j = max(2LL, (l + p - 1) / p); j * p <= r; j++) {
34
            notp[j * p - l] = 1;
35
        }
36
    }
37
    vector<ll> ret;
38
    for (ll i = l; i <= r; i++) {
39
        if (!notp[i - l]) {
40
            ret.emplace_back(i);
41
        }
42
    }
43
    return ret;
44
}

Miller-Rabin#

1
bool primetest(ll n) {
2
    if (n < 2 || n % 6 % 4 != 1) {
3
        return (n | 1) == 3;
4
    }
5
    ll s = __builtin_ctzll(n - 1), d = n >> s;
6
    for (ll a : {2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
7
        ll p = pow(a % n, d, n), i = s;
8
        while (p != 1 && p != n - 1 && a % n && i--) {
9
            p = mul(p, p, n);
10
        }
11
        if (p != n - 1 && i != s) {
12
            return false;
13
        }
14
    }
15
    return true;
16
}

Pollard-rho#

1
mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());
2

3
ll rho(ll n) {
4
    if (n % 2 == 0) return 2;
5
    if (n % 3 == 0) return 3;
6

7
    while (true) {
8
        ll c = uniform_int_distribution<ll>(1, n - 1)(rng);
9
        ll x = uniform_int_distribution<ll>(0, n - 1)(rng);
10
        ll y = x, d = 1;
11

12
        auto f = [&](ll x) -> ll {
13
            return (mul(x, x, n) + c) % n;
14
        };
15

16
        while (d == 1) {
17
            x = f(x);
18
            y = f(f(y));
19
            d = gcd(abs(x - y), n);
20
        }
21

22
        if (d != n) return d;
23
    }
24
}
25
vector<ll> factorize(ll x) {
26
    vector<ll> res;
27

28
    auto dfs = [&](auto self, ll x) -> void {
29
        if (x == 1) return;
30
        if (primetest(x)) {
31
            res.push_back(x);
32
            return;
33
        }
34

35
        ll d = rho(x);
36
        self(self, d);
37
        self(self, x / d);
38
    };
39

40
    dfs(dfs, x);
41
    sort(res.begin(), res.end());
42
    return res;
43
}

BSGS#

1
namespace BSGS {
2
ll mul(ll a, ll b, ll mod) {
3
    return (i128)a * b % mod;
4
}
5

6
ll power(ll a, ll b, ll mod) {
7
    ll ret = 1 % mod;
8
    while (b) {
9
        if (b & 1) ret = mul(ret, a, mod);
10
        a = mul(a, a, mod);
11
        b >>= 1;
12
    }
13
    return ret;
14
}
15

16
ll exgcd(ll a, ll b, ll &x, ll &y) {
17
    if (b == 0) {
18
        x = 1, y = 0;
19
        return a;
20
    }
21
    ll x1, y1;
22
    ll g = exgcd(b, a % b, x1, y1);
23
    x = y1;
24
    y = x1 - a / b * y1;
25
    return g;
26
}
27

28
ll inv(ll a, ll mod) {
29
    ll x, y;
30
    ll g = exgcd(a, mod, x, y);
31
    if (g != 1) return -1;
32
    x %= mod;
33
    if (x < 0) x += mod;
34
    return x;
35
}
36
ll bsgs(ll a, ll b, ll m) {
37
    a %= m;
38
    b %= m;
39

40
    if (m == 1) return 0;
41
    if (b == 1) return 0;
42

43
    ll n = sqrtl(m) + 1;
44
    HashMap<ll, ll> mp;
45
    // b * a^q
46
    ll cur = b;
47
    for (ll q = 0; q < n; q++) {
48
        mp[cur] = q;
49
        cur = mul(cur, a, m);
50
    }
51
    ll an = power(a, n, m);
52
    cur = 1;
53
    for (ll p = 1; p <= n + 1; p++) {
54
        cur = mul(cur, an, m);
55
        auto it = mp.find(cur);
56
        if (it != mp.end()) {
57
            ll x = p * n - it->second;
58
            if (x >= 0) return x;
59
        }
60
    }
61
    return -1;
62
}
63

64
ll exbsgs(ll a, ll b, ll m) {
65
    a %= m;
66
    b %= m;
67
    if (m == 1) return 0;
68
    if (b == 1) return 0;
69
    ll cnt = 0;
70
    ll cur = 1;
71
    while (true) {
72
        ll g = gcd(a, m);
73
        if (g == 1) break;
74
        if (b % g != 0) return -1;
75
        b /= g;
76
        m /= g;
77
        cur = mul(cur, a / g, m);
78
        cnt++;
79
        if (cur == b) return cnt;
80
    }
81
    ll iv = inv(cur, m);
82
    ll rhs = mul(b, iv, m);
83
    ll t = bsgs(a, rhs, m);
84
    if (t == -1) return -1;
85
    return cnt + t;
86
}
87
};

杜教筛#

1
constexpr int N = 1'000'000;
2
constexpr int PIN = 78'498;
3
bitset<N + 5> not_prime;
4
int primes[PIN + 5], pcnt = 0;
5
int phi[N + 5];
6
int mu[N + 5];
7
ll smu[N + 5];
8
ll sphi[N + 5];
9
unordered_map<ll, ll> mp_mu;
10
unordered_map<ll, i128> mp_phi;
11
void sieve() {
12
    phi[1] = 1;
13
    mu[1] = 1;
14
    for (int i = 2; i <= N; i++) {
15
        if (!not_prime[i]) {
16
            primes[pcnt] = i, pcnt++;
17
            phi[i] = i - 1;
18
            mu[i] = -1;
19
        }
20
        for (int j = 0; j < pcnt && 1LL * i * primes[j] <= N; j++) {
21
            int pj = primes[j];
22
            not_prime[i * pj] = 1;
23

24
            if (i % pj == 0) {
25
                phi[i * pj] = phi[i] * pj;
26
                mu[i * pj] = 0;
27
                break;
28
            }
29

30
            phi[i * pj] = phi[i] * (pj - 1);
31
            mu[i * pj] = -mu[i];
32
        }
33
    }
34

35
    for (int i = 1; i <= N; i++) {
36
        smu[i] = smu[i - 1] + mu[i];
37
        sphi[i] = sphi[i - 1] + phi[i];
38
    }
39
}
40

41
ll getsmu(ll n) {
42
    // 莫比乌斯函数的前缀和
43
    if (n <= N) {
44
        return smu[n];
45
    }
46
    if (mp_mu.count(n)) {
47
        return mp_mu[n];
48
    }
49

50
    ll ans = 1;
51

52
    for (ll l = 2, r; l <= n; l = r + 1) {
53
        ll q = n / l;
54
        r = n / q;
55
        ans -= (r - l + 1) * getsmu(q);
56
    }
57
    return mp_mu[n] = ans;
58
}
59

60
i128 getsphi(ll n) {
61
    // 欧拉函数的前缀和
62
    if (n <= N) {
63
        return sphi[n];
64
    }
65
    if (mp_phi.count(n)) {
66
        return mp_phi[n];
67
    }
68

69
    i128 ans = (i128)n * (n + 1) / 2;
70

71
    for (ll l = 2, r; l <= n; l = r + 1) {
72
        ll q = n / l;
73
        r = n / q;
74
        ans -= (i128)(r - l + 1) * getsphi(q);
75
    }
76

77
    return mp_phi[n] = ans;
78
}

杂项结论#

麦乐鸡定理#

对于互质整数 $n, m$ 来说, $an + bm(a \ge 0, b \ge 0)$ 不能表示出的最大整数为 $nm - n - m$ .

抽屉原理#

$n$ 个物体放到 $m$ 个抽屉中, 抽屉中物品的最大值大于等于 $ceil(n / m)$

Dilworth定理#

偏序集中, 最大链长度等于最多反链划分数, 对偶结论亦成立.

取模与整数除法#

x % i = x - i * (x / i)

欧拉函数性质#

$[1, n]$ 中与 $n$ 互质的数之和为 $\frac{n\varphi(n)}{2}$
若 $a, b$ 互质, 则 $\varphi(ab) = \varphi(a)\varphi(b)$
若 $p$ 是 $n$ 的质因子: 若 $p$ 的次数大于等于2, 则 $\varphi(n) = \varphi(\frac{n}{p})p$ ; 否则 $\varphi(n) = \varphi(\frac{n}{p}) (p - 1)$
$\sum_{d|n}\varphi(d) = n$
$\sum_{i = 1}^n\gcd(i, n) = \sum_{d|n}\varphi(d)\frac{n}{d}$ . (欧拉反演)

代数#

多项式#

FFT#

1
using cd = complex<double>;
2
using numbers::pi;
3
const cd I(0, 1);
4
vector<int> rev;
5
vector<cd> Wn;
6

7
void fft(vector<cd>& a, bool invert = 0) {
8
    int n = a.size();
9
    if ((int)rev.size() != n) {
10
        rev.assign(n, 0);
11
        for (int i = 1; i < n; i++) {
12
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (n >> 1));
13
        }
14
    }
15
    if ((int)Wn.size() != n) {
16
        Wn.assign(n, cd(0, 0));
17
        for (int i = 0; i < n; i++) {
18
            Wn[i] = cd(cos(pi / n * i), sin(pi / n * i));
19
        }
20
    }
21
    if (invert) {
22
        for (int i = 1; i < n; i++) {
23
            if (i < n - i) {
24
                swap(a[i], a[n - i]);
25
            }
26
        }
27
    }
28
    for (int i = 0; i < n; i++) {
29
        if (i < rev[i]) {
30
            swap(a[i], a[rev[i]]);
31
        }
32
    }
33
    for (int m = 1; m < n; m <<= 1) {
34
        for (int i = 0; i < n; i += m << 1) {
35
            for (int j = 0; j < m; j++) {
36
                cd w = Wn[1LL * j * n / m];
37
                cd x = a[i + j];
38
                cd y = a[i + j + m] * w;
39

40
                a[i + j] = x + y;
41
                a[i + j + m] = x - y;
42
            }
43
        }
44
    }
45
    if (invert) {
46
        for (int i = 0; i < n; i++) {
47
            a[i] /= n;
48
        }
49
    }
50
}
51

52
vector<ll> convolution(const vector<ll>& a, const vector<ll>& b) {
53
    if (a.empty() || b.empty()) {
54
        return {};
55
    }
56
    int sz = (int)a.size() + b.size() - 1;
57
    int n = 1;
58
    while (n < sz) {
59
        n <<= 1;
60
    }
61
    vector<cd> A(n), B(n);
62
    for (int i = 0; i < (int)a.size(); i++) {
63
        A[i] = a[i];
64
    }
65
    for (int i = 0; i < (int)b.size(); i++) {
66
        B[i] = b[i];
67
    }
68
    for (int i = 0; i < n; i++) {
69
        A[i] += I * B[i];
70
    }
71
    fft(A);
72
    for (int i = 0; i < n; i++) {
73
        B[i] = conj(A[i ? n - i : 0]);
74
    }
75
    for (int i = 0; i < n; i++) {
76
        cd p = A[i];
77
        cd q = B[i];
78
        A[i] = (p + q) * 0.5;
79
        B[i] = (q - p) * 0.5 * I;
80
    }
81
    for (int i = 0; i < n; i++) {
82
        A[i] *= B[i];
83
    }
84
    fft(A, 1);
85
    vector<ll> c(sz);
86
    for (int i = 0; i < sz; i++) {
87
        double x = A[i].real();
88
        c[i] = x < 0 ? (ll)(x - 0.5) : (ll)(x + 0.5);
89
    }
90
    return c;
91
}

NTT#

1
template<class Z>
2
void ntt(vector<Z>& a, bool invert = 0) {
3
    int n = a.size();
4
    vector<int> rev(n);
5
    for (int i = 0; i < n; i++) {
6
        rev[i] = ((i & 1) * (n / 2)) | (rev[i / 2] / 2);
7
        if (i < rev[i]) {
8
            swap(a[i], a[rev[i]]);
9
        }
10
    }
11
    Z g = 3;
12

13
    for (int len = 2; len <= n; len *= 2) {
14
        Z wlen = g.pow((Z::mod() - 1) / len);
15
        if (invert) {
16
            wlen = wlen.inv();
17
        }
18

19
        for (int i = 0; i < n; i += len) {
20
            Z w = 1;
21
            for (int j = 0; j < len / 2; j++) {
22
                Z u = a[i + j];
23
                Z v = a[i + j + len / 2] * w;
24
                a[i + j] = u + v;
25
                a[i + j + len / 2] = u - v;
26
                w *= wlen;
27
            }
28
        }
29
    }
30

31
    if (invert) {
32
        Z ni = Z(n).inv();
33
        for (int i = 0; i < n; i++) {
34
            a[i] *= ni;
35
        }
36
    }
37
}
38

39
template<class Z>
40
vector<Z> convolution(vector<Z> a, vector<Z> b) {
41
    if (a.empty() || b.empty()) {
42
        return {};
43
    }
44
    int sz = (int)a.size() + b.size() - 1;
45
    int n = 1;
46
    while (n < sz) {
47
        n *= 2;
48
    }
49

50
    a.resize(n), b.resize(n);
51
    ntt(a), ntt(b);
52
    for (int i = 0; i < n; i++) {
53
        a[i] *= b[i];
54
    }
55
    ntt(a, 1);
56

57
    a.resize(sz);
58
    return a;
59
}

三模NTT#

1
// NTT...
2
ll mul(ll a, ll b, ll p) {
3
    return (i128)a * b % p;
4
}
5
ll qkp(ll a, ll b, ll p) {
6
    ll ret = 1 % p;
7
    while (b > 0) {
8
        if (b & 1) {
9
            ret = mul(ret, a, p);
10
        }
11
        b >>= 1;
12
        a = mul(a, a, p);
13
    }
14
    return ret;
15
}
16
ll inv(ll x, ll p) {
17
    return qkp(x, p - 2, p);
18
}
19
vector<ll> mtt(const vector<ll>& a, const vector<ll>& b, ll mod) {
20
    if (a.empty() || b.empty()) {
21
        return {};
22
    }
23
    constexpr int M1 = 998244353;
24
    constexpr int M2 = 1004535809;
25
    constexpr int M3 = 469762049;
26
    using Z1 = ModInt<M1>;
27
    using Z2 = ModInt<M2>;
28
    using Z3 = ModInt<M3>;
29
    vector<Z1> a1(a.begin(), a.end()), b1(b.begin(), b.end());
30
    vector<Z2> a2(a.begin(), a.end()), b2(b.begin(), b.end());
31
    vector<Z3> a3(a.begin(), a.end()), b3(b.begin(), b.end());
32

33
    auto c1 = convolution(a1, b1);
34
    auto c2 = convolution(a2, b2);
35
    auto c3 = convolution(a3, b3);
36

37
    ll inv12 = inv(M1, M2);
38
    ll modm3 = mul(M1, M2, M3);
39
    ll inv123 = inv(modm3, M3);
40

41
    int sz = c1.size();
42
    vector<ll> ans(sz);
43
    for (int i = 0; i < sz; i++) {
44
        ll v1 = c1[i].val, v2 = c2[i].val, v3 = c3[i].val;
45

46
        ll k1 = (v2 - v1 % M2 + M2) % M2 * inv12 % M2;
47
        ll x12 = v1 + k1 * M1;
48

49
        ll k2 = (v3 - x12 % M3 + M3) % M3 * inv123 % M3;
50

51
        ll m12 = mul(M1, M2, mod);
52
        ll val = (x12 % mod + k2 % mod * m12 % mod) % mod;
53

54
        ans[i] = (val + mod) % mod;
55
    }
56
    return ans;
57
}

MTT#

1
vector<ll> mtt(vector<ll> a, vector<ll> b, ll mod) {
2
    if (a.empty() || b.empty()) {
3
        return {};
4
    }
5
    int sz = (int)a.size() + b.size() - 1;
6
    if ((ll)a.size() * b.size() <= 30000) {
7
        vector<ll> c(sz);
8
        for (int i = 0; i < (int)a.size(); i++) {
9
            a[i] %= mod;
10
            if (a[i] < 0) {
11
                a[i] += mod;
12
            }
13
        }
14
        for (int i = 0; i < (int)b.size(); i++) {
15
            b[i] %= mod;
16
            if (b[i] < 0) {
17
                b[i] += mod;
18
            }
19
        }
20
        for (int i = 0; i < (int)a.size(); i++) {
21
            for (int j = 0; j < (int)b.size(); j++) {
22
                c[i + j] = (c[i + j] + (i128)a[i] * b[j]) % mod;
23
            }
24
        }
25
        return c;
26
    }
27
    int n = 1;
28
    while (n < sz) {
29
        n <<= 1;
30
    }
31
    ll M = sqrt((long double)mod) + 1;
32
    vector<cd> a0(n), a1(n), b0(n), b1(n);
33
    for (int i = 0; i < (int)a.size(); i++) {
34
        a[i] %= mod;
35
        if (a[i] < 0) {
36
            a[i] += mod;
37
        }
38
        a0[i] = a[i] / M;
39
        a1[i] = a[i] % M;
40
    }
41
    for (int i = 0; i < (int)b.size(); i++) {
42
        b[i] %= mod;
43
        if (b[i] < 0) {
44
            b[i] += mod;
45
        }
46
        b0[i] = b[i] / M;
47
        b1[i] = b[i] % M;
48
    }
49
    auto fft2 = [&](vector<cd>& x, vector<cd>& y) {
50
        for (int i = 0; i < n; i++) {
51
            x[i] += I * y[i];
52
        }
53
        fft(x);
54
        for (int i = 0; i < n; i++) {
55
            y[i] = conj(x[i ? n - i : 0]);
56
        }
57
        for (int i = 0; i < n; i++) {
58
            cd p = x[i];
59
            cd q = y[i];
60
            x[i] = (p + q) * 0.5;
61
            y[i] = (q - p) * 0.5 * I;
62
        }
63
    };
64
    fft2(a0, a1);
65
    fft2(b0, b1);
66
    vector<cd> p(n), q(n);
67
    for (int i = 0; i < n; i++) {
68
        p[i] = a0[i] * b0[i] + I * a1[i] * b0[i];
69
        q[i] = a0[i] * b1[i] + I * a1[i] * b1[i];
70
    }
71
    fft(p, 1);
72
    fft(q, 1);
73
    auto num = [&](double x) -> ll {
74
        ll v = x < 0 ? (ll)(x - 0.5) : (ll)(x + 0.5);
75
        v %= mod;
76
        if (v < 0) {
77
            v += mod;
78
        }
79
        return v;
80
    };
81
    vector<ll> c(sz);
82
    ll M1 = M % mod;
83
    ll M2 = (i128)M1 * M1 % mod;
84
    for (int i = 0; i < sz; i++) {
85
        ll c00 = num(p[i].real());
86
        ll c10 = num(p[i].imag());
87
        ll c01 = num(q[i].real());
88
        ll c11 = num(q[i].imag());
89
        c[i] = (
90
            (i128)M2 * c00
91
            + (i128)M1 * ((c10 + c01) % mod)
92
            + c11
93
        ) % mod;
94
    }
95

96
    return c;
97
}

FWT/FMT#

1
template <typename Z>
2
void fwt_or(vector<Z>& a, bool invert = 0) {
3
    int n = a.size();
4
    for (int len = 1; len < n; len <<= 1) {
5
        for (int i = 0; i < n; i += len << 1) {
6
            for (int j = 0; j < len; j++) {
7
                if (!invert) {
8
                    a[i + j + len] += a[i + j];
9
                } else {
10
                    a[i + j + len] -= a[i + j];
11
                }
12
            }
13
        }
14
    }
15
}
16
template <typename Z>
17
void fwt_and(vector<Z>& a, bool invert = 0) {
18
    int n = a.size();
19
    for (int len = 1; len < n; len <<= 1) {
20
        for (int i = 0; i < n; i += len << 1) {
21
            for (int j = 0; j < len; j++) {
22
                if (!invert) {
23
                    a[i + j] += a[i + j + len];
24
                } else {
25
                    a[i + j] -= a[i + j + len];
26
                }
27
            }
28
        }
29
    }
30
}
31

32
template <typename Z>
33
void fwt_xor(vector<Z>& a, bool invert = 0) {
34
    int n = a.size();
35
    for (int len = 1; len < n; len <<= 1) {
36
        for (int i = 0; i < n; i += len << 1) {
37
            for (int j = 0; j < len; j++) {
38
                Z x = a[i + j];
39
                Z y = a[i + j + len];
40
                a[i + j] = x + y;
41
                a[i + j + len] = x - y;
42
            }
43
        }
44
    }
45

46
    if (invert) {
47
        Z invn = Z(n).inv();
48
        for (int i = 0; i < n; i++) {
49
            a[i] *= invn;
50
        }
51
    }
52
}
53
template <typename Z>
54
vector<Z> convolution_or(vector<Z> a, vector<Z> b) {
55
    int n = a.size();
56
    assert((int)b.size() == n);
57

58
    fwt_or(a);
59
    fwt_or(b);
60

61
    for (int i = 0; i < n; i++) {
62
        a[i] *= b[i];
63
    }
64

65
    fwt_or(a, 1);
66
    return a;
67
}
68

69
template <typename Z>
70
vector<Z> convolution_and(vector<Z> a, vector<Z> b) {
71
    int n = a.size();
72
    assert((int)b.size() == n);
73

74
    fwt_and(a);
75
    fwt_and(b);
76

77
    for (int i = 0; i < n; i++) {
78
        a[i] *= b[i];
79
    }
80

81
    fwt_and(a, 1);
82
    return a;
83
}
84

85
template <typename Z>
86
vector<Z> convolution_xor(vector<Z> a, vector<Z> b) {
87
    int n = a.size();
88
    assert((int)b.size() == n);
89

90
    fwt_xor(a);
91
    fwt_xor(b);
92

93
    for (int i = 0; i < n; i++) {
94
        a[i] *= b[i];
95
    }
96

97
    fwt_xor(a, 1);
98
    return a;
99
}

求导积分#

1
template <typename Z>
2
vector<Z> deriv(const vector<Z>& a) {
3
    int n = a.size();
4
    vector<Z> ret(max(0, n - 1));
5
    for (int i = 0; i < n - 1; i++) {
6
        ret[i] = a[i + 1] * Z(i + 1);
7
    }
8
    return ret;
9
}
10
template <typename Z>
11
vector<Z> intgr(const vector<Z>& a) {
12
    int n = a.size();
13
    vector<Z> ret(n + 1);
14
    for (int i = 1; i <= n; i++) {
15
        ret[i] = a[i - 1] / Z(i);
16
    }
17
    return ret;
18
}

求逆#

1
// 卷积...
2
template<typename Z>
3
vector<Z> inv(vector<Z> a, int n) {
4
    // find B that AB === 1 (mod x^n)
5
    int m = 1;
6
    vector<Z> b(1);
7
    b[0] = a[0].inv();
8
    while (m < n) {
9
        // B = 2B' - A(B'^2)
10
        int len = min(m * 2, n);
11
        vector<Z> A(min(len, (int)a.size()));
12
        for (int i = 0; i < (int)A.size(); i++) {
13
            A[i] = a[i];
14
        }
15
        vector<Z> B(b);
16
        B.resize(len);
17
        auto C = convolution(A, convolution(B, B));
18
        C.resize(len);
19
        b.resize(len);
20
        for (int i = 0; i < len; i++) {
21
            b[i] = B[i] * 2 - C[i];
22
        }
23
        m *= 2;
24
    }
25
    b.resize(n);
26
    return b;
27
}

除法#

1
// 卷积, 求逆...
2
template <typename Z>
3
pair<vector<Z>, vector<Z>> divmod(vector<Z> a, vector<Z> b) {
4
    int n = a.size();
5
    int m = b.size();
6
    if (n < m) {
7
        return { vector<Z>{ Z(0) }, a };
8
    }
9
    int len = n - m + 1;
10
    vector<Z> ra = a;
11
    vector<Z> rb = b;
12
    reverse(ra.begin(), ra.end());
13
    reverse(rb.begin(), rb.end());
14
    ra.resize(len);
15
    rb.resize(len);
16
    auto q = convolution(ra, inv(rb, len));
17
    q.resize(len);
18
    reverse(q.begin(), q.end());
19
    auto c = convolution(q, b);
20
    for (int i = 0; i < min((int)a.size(), (int)c.size()); i++) {
21
        a[i] -= c[i];
22
    }
23
    a.resize(m - 1);
24
    return { q, a };
25
}

开根#

1
// 卷积, 求逆...
2
template <typename Z>
3
vector<Z> sqrt(const vector<Z>& a, int n) {
4
    // assert(a[0] == 1);
5
    int m = 1;
6
    vector<Z> b(1);
7
    b[0] = 1;
8

9
    Z inv2 = Z(2).inv();
10
    while (m < n) {
11
        // B = (B' + A / B') / 2
12
        int len = min(m * 2, n);
13
        vector<Z> A(len);
14
        for (int i = 0; i < min((int)a.size(), len); i++) {
15
            A[i] = a[i];
16
        }
17
        b.resize(len);
18
        vector<Z> B(b);
19
        b = convolution(A, inv(B, len));
20
        b.resize(len);
21
        for (int i = 0; i < len; i++) {
22
            b[i] += B[i];
23
        }
24
        for (int i = 0; i < len; i++) {
25
            b[i] *= inv2;
26
        }
27
        m *= 2;
28
    }
29
    b.resize(n);
30
    return b;
31
}

指数/对数#

1
// 卷积, 求导, 积分...
2
template <typename Z>
3
vector<Z> ln(const vector<Z>& a, int n) {
4
    // B = ln(A) = intgr(deriv(A) * inv(A))
5
    // assert(a[0] == 1);
6
    vector<Z> A(n);
7
    for (int i = 0; i < n; i++) {
8
        if (i < a.size()) {
9
            A[i] = a[i];
10
        }
11
    }
12
    auto B = intgr(convolution(deriv(A), inv(A, n)));
13
    B.resize(n);
14
    return B;
15
}
16

17
template <typename Z>
18
vector<Z> exp(vector<Z> a, int n) {
19
    // assert(a[0] == 0);
20
    int m = 1;
21
    vector<Z> b(1);
22
    b[0] = 1;
23
    while (m < n) {
24
        // B = B' * (1 - ln(B') + A)
25
        int len = min(m * 2, n);
26
        vector<Z> A(len);
27
        for (int i = 0; i < min((int)a.size(), len); i++) {
28
            A[i] = a[i];
29
        }
30
        vector<Z> B(b);
31
        B.resize(len);
32
        auto C = ln(B, len);
33
        for (int i = 0; i < len; i++) {
34
            C[i] = A[i] - C[i];
35
        }
36
        C[0] += 1;
37
        b = convolution(B, C);
38
        b.resize(len);
39
        m *= 2;
40
    }
41
    b.resize(n);
42
    return b;
43
}

幂#

1
template <typename Z>
2
vector<Z> pow(vector<Z> a, ll b, int n) {
3
    if (n == 0) return {};
4
    if (b == 0) {
5
        vector<Z> ret(n);
6
        ret[0] = 1;
7
        return ret;
8
    }
9
    a.resize(n);
10
    int t = n;
11
    for (int i = 0; i < n; i++) {
12
        if (a[i] != Z(0)) {
13
            t = i;
14
            break;
15
        }
16
    }
17
    vector<Z> ans(n);
18
    if (t == n) {
19
        return ans;
20
    }
21
    if ((__int128)t * b >= n) {
22
        return ans;
23
    }
24
    int s = t * b;
25
    int m = n - s;
26
    vector<Z> c(m);
27
    for (int i = 0; i < m && i + t < n; i++) {
28
        c[i] = a[i + t] / a[t];
29
    }
30
    c = ln(c, m);
31
    for (int i = 0; i < m; i++) {
32
        c[i] *= Z(b);
33
    }
34
    c = exp(c, m);
35
    Z l = a[t].pow(b);
36
    for (int i = 0; i < m; i++) {
37
        ans[i + s] = c[i] * l;
38
    }
39
    return ans;
40
}

多点求值#

1
// 卷积, 求逆...
2
template <typename Z>
3
vector<Z> mulT(vector<Z> a, vector<Z> b, int need) {
4
    vector<Z> ret(need);
5
    if (a.empty() || b.empty() || need == 0) {
6
        return ret;
7
    }
8
    int n = b.size();
9
    reverse(b.begin(), b.end());
10
    auto c = convolution(a, b);
11
    for (int i = 0; i < need; i++) {
12
        int p = i + n - 1;
13
        if (p < (int)c.size()) {
14
            ret[i] = c[p];
15
        }
16
    }
17
    return ret;
18
}
19

20
template <typename Z>
21
vector<Z> multipoint_eval(vector<Z> f, vector<Z> x) {
22
    int m = x.size();
23
    vector<Z> ans(m);
24
    if (m == 0) {
25
        return ans;
26
    }
27
    if (f.empty()) {
28
        return ans;
29
    }
30
    int n = max((int)f.size(), m);
31
    x.resize(n);
32
    vector<vector<Z>> q(4 * n);
33
    auto build = [&](auto self, int p, int l, int r) -> void {
34
        if (r - l == 1) {
35
            q[p] = {Z(1), -x[l]};
36
            return;
37
        }
38
        int mid = (l + r) / 2;
39
        self(self, p << 1, l, mid);
40
        self(self, p << 1 | 1, mid, r);
41
        q[p] = convolution(q[p << 1], q[p << 1 | 1]);
42
    };
43
    auto work = [&](auto self, int p, int l, int r, const vector<Z>& num) -> void {
44
        if (r - l == 1) {
45
            if (l < m) {
46
                ans[l] = num.empty() ? Z(0) : num[0];
47
            }
48
            return;
49
        }
50
        int mid = (l + r) / 2;
51
        auto vl = mulT(num, q[p << 1 | 1], mid - l);
52
        auto vr = mulT(num, q[p << 1], r - mid);
53
        self(self, p << 1, l, mid, vl);
54
        self(self, p << 1 | 1, mid, r, vr);
55
    };
56
    build(build, 1, 0, n);
57
    auto root = mulT(f, inv(q[1], n), n);
58
    work(work, 1, 0, n, root);
59
    return ans;
60
}

快速插值#

1
template <typename Z>
2
vector<Z> interpolate(vector<Z> x, vector<Z> y) {
3
    int n = x.size();
4
    assert((int)y.size() == n);
5
    if (n == 0) {
6
        return {};
7
    }
8
    vector<vector<Z>> q(4 * n);
9
    auto build = [&](auto self, int p, int l, int r) -> void {
10
        if (r - l == 1) {
11
            // x - x_l
12
            q[p] = {-x[l], Z(1)};
13
            return;
14
        }
15

16
        int mid = (l + r) / 2;
17

18
        self(self, p << 1, l, mid);
19
        self(self, p << 1 | 1, mid, r);
20

21
        q[p] = convolution(q[p << 1], q[p << 1 | 1]);
22
    };
23

24
    build(build, 1, 0, n);
25
    auto d = multipoint_eval(deriv(q[1]), x);
26
    vector<Z> w(n);
27
    for (int i = 0; i < n; i++) {
28
        w[i] = y[i] / d[i];
29
    }
30

31
    auto work = [&](auto self, int p, int l, int r) -> vector<Z> {
32
        if (r - l == 1) {
33
            return vector<Z>{w[l]};
34
        }
35

36
        int mid = (l + r) / 2;
37

38
        auto L = self(self, p << 1, l, mid);
39
        auto R = self(self, p << 1 | 1, mid, r);
40

41
        auto A = convolution(L, q[p << 1 | 1]);
42
        auto B = convolution(R, q[p << 1]);
43

44
        int len = max(A.size(), B.size());
45
        A.resize(len);
46
        B.resize(len);
47

48
        for (int i = 0; i < len; i++) {
49
            A[i] += B[i];
50
        }
51

52
        return A;
53
    };
54
    auto ans = work(work, 1, 0, n);
55
    ans.resize(n);
56
    return ans;
57
}

快速阶乘#

1
struct FastFactorial {
2
    using ld = long double;
3
    static constexpr int MAXD = 18;
4
    static constexpr int LIM = (1 << MAXD) + 5;
5
    static constexpr int BASE = 1 << 16;
6
    static constexpr int SF = 16;
7
    static constexpr int MASK = BASE - 1;
8
    struct Complex {
9
        ld r, i;
10
        friend Complex operator+(Complex a, Complex b) {
11
            return {a.r + b.r, a.i + b.i};
12
        }
13
        friend Complex operator-(Complex a, Complex b) {
14
            return {a.r - b.r, a.i - b.i};
15
        }
16
        friend Complex operator*(Complex a, Complex b) {
17
            return {a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r};
18
        }
19
        Complex& operator/=(int x) {
20
            r /= x;
21
            i /= x;
22
            return *this;
23
        }
24
    };
25
    ll mod, base2;
26
    vector<vector<int>> rev;
27
    vector<Complex> rt[2][MAXD + 1];
28
    vector<Complex> tr, tr1, tr2, tr3, tr4, tr5, tr6;
29
    vector<ll> m13, m14, m23, m24;
30
    vector<ll> f, g, h, ifac;
31
    vector<ll> val, tmp1, tmp2;
32
    FastFactorial() {
33
        rev.resize(MAXD + 1);
34
        tr.resize(LIM);
35
        tr1.resize(LIM);
36
        tr2.resize(LIM);
37
        tr3.resize(LIM);
38
        tr4.resize(LIM);
39
        tr5.resize(LIM);
40
        tr6.resize(LIM);
41
        m13.resize(LIM);
42
        m14.resize(LIM);
43
        m23.resize(LIM);
44
        m24.resize(LIM);
45
        f.resize(LIM);
46
        g.resize(LIM);
47
        h.resize(LIM);
48
        ifac.resize(LIM);
49
        val.resize(LIM);
50
        tmp1.resize(LIM);
51
        tmp2.resize(LIM);
52
        pre_fft();
53
    }
54
    ll norm(ll x) {
55
        x %= mod;
56
        if (x < 0) x += mod;
57
        return x;
58
    }
59
    ll qpow(ll a, ll b) {
60
        ll ret = 1 % mod;
61
        a = norm(a);
62
        while (b > 0) {
63
            if (b & 1) ret = (i128)ret * a % mod;
64
            a = (i128)a * a % mod;
65
            b >>= 1;
66
        }
67
        return ret;
68
    }
69
    void pre_fft() {
70
        const ld PI = acosl(-1.0L);
71
        for (int d = 1; d <= MAXD; d++) {
72
            int len = 1 << d;
73
            rev[d].assign(len, 0);
74
            for (int i = 1; i < len; i++) {
75
                rev[d][i] = (rev[d][i >> 1] >> 1) | ((i & 1) << (d - 1));
76
            }
77
        }
78
        for (int d = 1; d <= MAXD; d++) {
79
            int m = 1 << (d - 1);
80
            rt[0][d].resize(m);
81
            rt[1][d].resize(m);
82
            for (int i = 0; i < m; i++) {
83
                ld ang = PI * i / m;
84
                rt[0][d][i] = {cosl(ang), sinl(ang)};
85
                rt[1][d][i] = {cosl(ang), -sinl(ang)};
86
            }
87
        }
88
    }
89
    void fft(Complex* a, int len, int d, int inv) {
90
        for (int i = 1; i < len; i++) {
91
            if (i < rev[d][i]) {
92
                swap(a[i], a[rev[d][i]]);
93
            }
94
        }
95
        for (int k = 1, dep = 1; k < len; k <<= 1, dep++) {
96
            for (int s = 0; s < len; s += k << 1) {
97
                auto* w = rt[inv][dep].data();
98
                for (int i = s; i < s + k; i++, w++) {
99
                    Complex x = a[i];
100
                    Complex y = a[i + k] * (*w);
101
                    a[i] = x + y;
102
                    a[i + k] = x - y;
103
                }
104
            }
105
        }
106
        if (inv) {
107
            for (int i = 0; i < len; i++) {
108
                a[i] /= len;
109
            }
110
        }
111
    }
112
    void dbdft(ll* a, int len, int d, Complex* op1, Complex* op2) {
113
        for (int i = 0; i < len; i++) {
114
            tr[i] = {(ld)(a[i] >> SF), (ld)(a[i] & MASK)};
115
        }
116
        fft(tr.data(), len, d, 0);
117
        tr[len] = tr[0];
118
        for (int i = 0; i < len; i++) {
119
            Complex p = tr[i];
120
            Complex q = tr[len - i];
121

122
            op1[i] = Complex{p.r + q.r, p.i - q.i} * Complex{0.5, 0};
123
            op2[i] = Complex{p.r - q.r, p.i + q.i} * Complex{0, -0.5};
124
        }
125
    }
126
    ll round_mod(ld x) {
127
        ll v = x < 0 ? (ll)(x - 0.5) : (ll)(x + 0.5);
128
        v %= mod;
129
        if (v < 0) v += mod;
130
        return v;
131
    }
132
    void dbidft(Complex* a, int len, int d, ll* op1, ll* op2) {
133
        fft(a, len, d, 1);
134
        for (int i = 0; i < len; i++) {
135
            op1[i] = round_mod(a[i].r);
136
            op2[i] = round_mod(a[i].i);
137
        }
138
    }
139
    void poly_mul(ll* a, ll* b, ll* c, int len, int d) {
140
        dbdft(a, len, d, tr1.data(), tr2.data());
141
        dbdft(b, len, d, tr3.data(), tr4.data());
142
        for (int i = 0; i < len; i++) {
143
            tr5[i] = tr1[i] * tr3[i] + Complex{0, 1} * (tr2[i] * tr4[i]);
144
            tr6[i] = tr2[i] * tr3[i] + Complex{0, 1} * (tr1[i] * tr4[i]);
145
        }
146
        dbidft(tr5.data(), len, d, m13.data(), m24.data());
147
        dbidft(tr6.data(), len, d, m23.data(), m14.data());
148
        for (int i = 0; i < len; i++) {
149
            c[i] = (i128)m13[i] * base2 % mod;
150
            c[i] = (c[i] + (i128)(m23[i] + m14[i]) % mod * BASE + m24[i]) % mod;
151
        }
152
    }
153
    void init_ifac(int lim) {
154
        lim = min(lim, LIM - 1);
155
        ifac[0] = 1;
156
        if (lim >= 1) ifac[1] = 1;
157
        for (int i = 2; i <= lim; i++) {
158
            ifac[i] = (mod - mod / i) * ifac[mod % i] % mod;
159
        }
160
        for (int i = 1; i <= lim; i++) {
161
            ifac[i] = ifac[i] * ifac[i - 1] % mod;
162
        }
163
    }
164
    void shift_eval(ll del, int cur, ll* ip, ll* op) {
165
        int len = 1;
166
        int d = 0;
167
        while (len <= cur + cur + cur) {
168
            len <<= 1;
169
            d++;
170
        }
171
        for (int i = 0; i <= cur; i++) {
172
            f[i] = ip[i] * ifac[i] % mod * ifac[cur - i] % mod;
173
        }
174
        for (int i = cur - 1; i >= 0; i -= 2) {
175
            if (f[i]) f[i] = mod - f[i];
176
        }
177
        int total = cur + cur + 1;
178
        ll prod = 1;
179
        for (int i = 0; i < total; i++) {
180
            g[i] = norm(del - cur + i);
181
            prod = (i <= cur ? (i128)prod * g[i] % mod : prod);
182
        }
183
        h[0] = 1;
184
        for (int i = 0; i < total; i++) {
185
            h[i + 1] = (i128)h[i] * g[i] % mod;
186
        }
187
        ll inv_all = qpow(h[total], mod - 2);
188
        for (int i = total - 1; i >= 0; i--) {
189
            ll x = g[i];
190
            g[i] = (i128)inv_all * h[i] % mod;
191
            inv_all = (i128)inv_all * x % mod;
192
        }
193
        for (int i = cur + 1; i < len; i++) {
194
            f[i] = 0;
195
        }
196
        for (int i = total; i < len; i++) {
197
            g[i] = 0;
198
        }
199
        poly_mul(f.data(), g.data(), h.data(), len, d);
200
        ll cur_prod = prod;
201
        for (int i = 0; i <= cur; i++) {
202
            op[i] = h[i + cur] * cur_prod % mod;
203
            cur_prod = (i128)cur_prod * g[i] % mod;
204
            cur_prod = (i128)cur_prod * norm(del + i + 1) % mod;
205
        }
206
    }
207
    void build_block_values(int B) {
208
        int hb = 0;
209
        for (int x = B; x; x >>= 1) {
210
            hb++;
211
        }
212
        val[0] = 1;
213
        int cur = 0;
214
        ll invB = qpow(B, mod - 2);
215
        for (int z = hb; z >= 0; z--) {
216
            if (cur != 0) {
217
                shift_eval(cur + 1, cur, val.data(), tmp1.data());
218
                for (int i = 0; i <= cur; i++) {
219
                    val[cur + i + 1] = tmp1[i];
220
                }
221
                val[cur * 2 + 1] = 0;
222
                shift_eval((ll)cur * invB % mod, cur << 1, val.data(), tmp2.data());
223
                cur <<= 1;
224
                for (int i = 0; i <= cur; i++) {
225
                    val[i] = val[i] * tmp2[i] % mod;
226
                }
227
            }
228
            if ((B >> z) & 1) {
229
                for (int i = 0; i <= cur; i++) {
230
                    val[i] = val[i] * ((ll)B * i + cur + 1) % mod;
231
                }
232
                cur |= 1;
233
                val[cur] = 1;
234
                for (int i = 1; i <= cur; i++) {
235
                    val[cur] = val[cur] * ((ll)cur * B + i) % mod;
236
                }
237
            }
238
        }
239
    }
240
    ll brute_fact(ll n) {
241
        ll ans = 1 % mod;
242
        for (ll i = 1; i <= n; i++) {
243
            ans = (i128)ans * i % mod;
244
        }
245
        return ans;
246
    }
247
    ll fact_raw(ll n) {
248
        if (n == 0) return 1 % mod;
249
        // 这个阈值可以调整, 不要太小
250
        if (n <= 500000) {
251
            return brute_fact(n);
252
        }
253
        int B = sqrt((long double)n);
254
        while ((ll)B * B < n) B++;
255
        while (B > 1 && (ll)(B - 1) * (B - 1) >= n) B--;
256
        init_ifac(2 * B + 10);
257
        build_block_values(B);
258
        ll ans = 1 % mod;
259
        ll i = 0;
260
        int id = 0;
261
        while (i + B <= n) {
262
            ans = ans * val[id] % mod;
263
            i += B;
264
            id++;
265
        }
266
        for (ll j = i + 1; j <= n; j++) {
267
            ans = (i128)ans * j % mod;
268
        }
269
        return ans;
270
    }
271

272
    ll operator()(ll n, ll p) {
273
        mod = p;
274
        base2 = (ll)BASE * BASE % mod;
275
        if (n >= mod) return 0;
276
        if (n > mod - 1 - n) {
277
            ll r = mod - 1 - n;
278
            ll ans = fact_raw(r);
279
            ans = qpow(ans, mod - 2);
280
            if ((r + 1) & 1) {
281
                ans = (mod - ans) % mod;
282
            }
283
            return ans;
284
        }
285
        return fact_raw(n);
286
    }
287
};
288
FastFactorial f;

多点快速阶乘#

1
template <typename Z>
2
vector<Z> batch_inv(const vector<Z>& a) {
3
    int n = a.size();
4
    vector<Z> pre(n + 1), suf(n + 1), ret(n);
5
    pre[0] = 1;
6
    for (int i = 0; i < n; i++) {
7
        pre[i + 1] = pre[i] * a[i];
8
    }
9
    Z inv_all = pre[n].inv();
10
    for (int i = n - 1; i >= 0; i--) {
11
        ret[i] = inv_all * pre[i];
12
        inv_all *= a[i];
13
    }
14
    return ret;
15
}
16

17
template <typename Z>
18
vector<Z> many_factorials(const vector<int>& q) {
19
    constexpr int LIM = 100000;
20
    int n = q.size();
21
    int mod = Z::mod();
22
    vector<Z> ans(n, Z(1));
23
    vector<int> need_inv(n);
24
    vector<pair<int, int>> reg;
25
    vector<pair<int, int>> odd;
26
    reg.reserve(n);
27
    odd.reserve(n * 8);
28
    for (int id = 0; id < n; id++) {
29
        int t = q[id];
30
        if (t >= mod) {
31
            ans[id] = 0;
32
            continue;
33
        }
34
        if (t >= mod / 2) {
35
            int r = mod - 1 - t;
36
            if ((r & 1) == 0) {
37
                ans[id] = -ans[id];
38
            }
39
            need_inv[id] = 1;
40
            t = r;
41
        }
42
        long long pow2 = 0;
43
        while (t > LIM) {
44
            int x = (t - 1) / 2;
45
            odd.emplace_back(x, id);
46
            pow2 += t / 2;
47
            t >>= 1;
48
        }
49
        ans[id] *= Z(2).pow(pow2 % (mod - 1));
50
        reg.emplace_back(t, id);
51
    }
52
    sort(reg.begin(), reg.end());
53
    int cur = 0;
54
    Z prod = 1;
55
    for (auto [t, id] : reg) {
56
        while (cur < t) {
57
            cur++;
58
            prod *= cur;
59
        }
60
        ans[id] *= prod;
61
    }
62
    sort(odd.begin(), odd.end());
63
    int ocur = -1;
64
    Z oprod = 1;
65
    for (auto [x, id] : odd) {
66
        while (ocur < x) {
67
            ocur++;
68
            oprod *= Z(2LL * ocur + 1);
69
        }
70
        ans[id] *= oprod;
71
    }
72
    vector<Z> need;
73
    vector<int> pos;
74
    for (int i = 0; i < n; i++) {
75
        if (need_inv[i]) {
76
            pos.emplace_back(i);
77
            need.emplace_back(ans[i]);
78
        }
79
    }
80
    auto invs = batch_inv(need);
81
    for (int i = 0; i < (int)pos.size(); i++) {
82
        ans[pos[i]] = invs[i];
83
    }
84
    return ans;
85
}

线性代数#

矩阵和列向量#

1
template <class S, int N>
2
struct ColVector {
3
    array<S, N> a;
4
    ColVector(S def = S()) {
5
        fill(a.begin(), a.end(), def);
6
    }
7

8
    S& operator[](int i) {
9
        return a[i];
10
    }
11
    const S& operator[](int i) const {
12
        return a[i];
13
    }
14
};
15

16
template <typename S, int N>
17
struct Matrix {
18
    array<array<S, N>, N> a;
19
    Matrix(S def = S()) {
20
        for (int i = 0; i < N; i++) {
21
            fill(a[i].begin(), a[i].end(), def);
22
        }
23
    }
24
    using M = Matrix;
25

26
    S* operator[](int i) {
27
        return a[i].data();
28
    }
29
    const S* operator[](int i) const {
30
        return a[i].data();
31
    }
32

33
    // 单位矩阵
34
    constexpr static M e() {
35
        M ret;
36
        for (int i = 0; i < N; i++) {
37
            ret[i][i] = 1;
38
        }
39
        return ret;
40
    }
41

42
    // 加法
43
    friend M operator+(M lt, const M& rt) {
44
        for (int i = 0; i < N; i++) {
45
            for (int j = 0; j < N; j++) {
46
                lt[i][j] += rt[i][j];
47
            }
48
        }
49
        return lt;
50
    }
51

52
    // 乘法
53
    friend M operator*(const M& lt, const M& rt) {
54
        M ret;
55
        for (int i = 0; i < N; i++) {
56
            for (int k = 0; k < N; k++) {
57
                if (lt[i][k] == S()) {
58
                    continue;
59
                }
60

61
                for (int j = 0; j < N; j++) {
62
                    ret[i][j] += lt[i][k] * rt[k][j];
63
                }
64
            }
65
        }
66
        return ret;
67
    }
68

69
    M pow(ll b) const {
70
        M ret = e();
71
        M t = *this;
72
        while (b > 0) {
73
            if (b & 1) {
74
                ret = ret * t;
75
            }
76

77
            t = t * t;
78
            b >>= 1;
79
        }
80
        return ret;
81
    }
82

83
    ColVector<S, N> operator*(const ColVector<S, N>& vec) {
84
        ColVector<S, N> ret;
85
        for (int i = 0; i < N; i++) {
86
            for (int k = 0; k < N; k++) {
87
                if (a[i][k] == S()) {
88
                    continue;
89
                }
90

91
                ret[i] = ret[i] + (a[i][k] * vec[k]);
92
            }
93
        }
94
        return ret;
95
    }
96
};
97

98
template <int N>
99
struct Matrix_mp {
100
    array<array<ll, N>, N> a;
101
    Matrix_mp() {
102
        for (int i = 0; i < N; i++) {
103
            fill(a[i].begin(), a[i].end(), INFLL);
104
        }
105
    }
106

107
    ll* operator[](int i) {
108
        return a[i].data();
109
    }
110
    const ll* operator[](int i) const {
111
        return a[i].data();
112
    }
113

114
    static Matrix_mp e() {
115
        Matrix_mp ret;
116
        for (int i = 0; i < N; i++) {
117
            ret[i][i] = 0;
118
        }
119
        return ret;
120
    }
121

122
    friend Matrix_mp operator+(Matrix_mp lt, const Matrix_mp& rt) {
123
        for (int i = 0; i < N; i++) {
124
            for (int j = 0; j < N; j++) {
125
                lt[i][j] = min(lt[i][j], rt[i][j]);
126
            }
127
        }
128
        return lt;
129
    }
130
    friend Matrix_mp operator*(const Matrix_mp& lt, const Matrix_mp& rt) {
131
        Matrix_mp ret;
132
        for (int i = 0; i < N; i++) {
133
            for (int k = 0; k < N; k++) {
134
                if (lt[i][k] >= INFLL) {
135
                    continue;
136
                }
137

138
                for (int j = 0; j < N; j++) {
139
                    ll now;
140
                    if (rt[k][j] >= INFLL) {
141
                        now = INFLL;
142
                    } else {
143
                        now = lt[i][k] + rt[k][j];
144
                    }
145
                    ret[i][j] = min(ret[i][j], now);
146
                }
147
            }
148
        }
149
        return ret;
150
    }
151

152
    Matrix_mp pow(ll b) const {
153
        Matrix_mp ret = e();
154
        Matrix_mp t = *this;
155
        while (b > 0) {
156
            if (b & 1) {
157
                ret = ret * t;
158
            }
159

160
            t = t * t;
161
            b >>= 1;
162
        }
163
        return ret;
164
    }
165

166
    ColVector<ll, N> operator*(const ColVector<ll, N>& vec) {
167
        ColVector<ll, N> ret(INFLL);
168
        for (int i = 0; i < N; i++) {
169
            for (int k = 0; k < N; k++) {
170
                if (a[i][k] >= INFLL) {
171
                    continue;
172
                }
173
                ll now;
174
                if (a[i][k] >= INFLL || vec[k] >= INFLL) {
175
                    now = INFLL;
176
                } else {
177
                    now = a[i][k] + vec[k];
178
                }
179
                ret[i] = min(ret[i], now);
180
            }
181
        }
182
        return ret;
183
    }
184
};

高斯消元(整数)#

1
template <typename S, int N>
2
struct GaussJordan {
3
    Matrix<S, N> A;
4
    ColVector<S, N> x;
5
    ColVector<S, N> b;
6

7
    int solve(int n, int m) {
8
        int r = 0;
9
        for (int c = 0; r < n && c < m; c++) {
10
            int pivot = r;
11
            while (pivot < n && A[pivot][c] == S()) {
12
                pivot++;
13
            }
14
            if (pivot == n) {
15
                continue;
16
            }
17
            if (pivot != r) {
18
                swap(A.a[pivot], A.a[r]);
19
                swap(b.a[pivot], b.a[r]);
20
            }
21
            S inv = S(1) / A[r][c];
22
            for (int j = c; j < m; j++) {
23
                A[r][j] = A[r][j] * inv;
24
            }
25
            b[r] = b[r] * inv;
26

27
            for (int i = 0; i < n; i++) {
28
                if (i != r && A[i][c] != S()) {
29
                    S factor = A[i][c];
30
                    for (int j = c; j < m; j++) {
31
                        A[i][j] = A[i][j] - factor * A[r][j];
32
                    }
33
                    b[i] = b[i] - factor * b[r];
34
                }
35
            }
36
            r++;
37
        }
38

39
        for (int i = r; i < n; i++) {
40
            if (b[i] != S()) {
41
                return -1;
42
            }
43
        }
44

45
        for (int i = 0; i < n; i++) {
46
            x[i] = S();
47
        }
48

49
        for (int i = 0; i < r; i++) {
50
            int maj = -1;
51
            for (int j = 0; j < m; j++) {
52
                if (A[i][j] != S()) {
53
                    maj = j;
54
                    break;
55
                }
56
            }
57
            if (maj != -1) {
58
                x[maj] = b[i];
59
            }
60
        }
61
        if (r < n) {
62
            return n - r;
63
        } else {
64
            return 0;
65
        }
66
    }
67
};

高斯消元(实数)#

1
// for real
2
template <typename S, int N>
3
struct GaussJordan {
4
    Matrix<S, N> A;
5
    ColVector<S, N> x;
6
    ColVector<S, N> b;
7

8
    int solve(int n, int m) {
9
        int r = 0;
10
        for (int c = 0; r < n && c < m; c++) {
11
            int pivot = r;
12
            for (int i = r + 1; i < n; i++) {
13
                if (fabs(A[i][c]) > fabs(A[pivot][c])) {
14
                    pivot = i;
15
                }
16
            }
17
            if (fabs(A[pivot][c] < EPS)) {
18
                continue;
19
            }
20

21
            if (pivot != r) {
22
                swap(A.a[pivot], A.a[r]);
23
                swap(b.a[pivot], b.a[r]);
24
            }
25
            S inv = S(1) / A[r][c];
26
            for (int j = c; j < m; j++) {
27
                A[r][j] = A[r][j] * inv;
28
            }
29
            b[r] = b[r] * inv;
30

31
            for (int i = 0; i < n; i++) {
32
                if (i != r && A[i][c] != S()) {
33
                    S factor = A[i][c];
34
                    for (int j = c; j < m; j++) {
35
                        A[i][j] = A[i][j] - factor * A[r][j];
36
                    }
37
                    b[i] = b[i] - factor * b[r];
38
                }
39
            }
40
            r++;
41
        }
42

43
        for (int i = r; i < n; i++) {
44
            if (fabs(b[i]) > EPS) {
45
                return -1;
46
            }
47
        }
48

49
        for (int i = 0; i < n; i++) {
50
            x[i] = S();
51
        }
52

53
        for (int i = 0; i < r; i++) {
54
            int maj = -1;
55
            for (int j = 0; j < m; j++) {
56
                if (fabs(A[i][j]) > EPS) {
57
                    maj = j;
58
                    break;
59
                }
60
            }
61
            if (maj != -1) {
62
                x[maj] = b[i];
63
            }
64
        }
65
        if (r < n) {
66
            return n - r;
67
        } else {
68
            return 0;
69
        }
70
    }
71
};

高斯消元(异或)#

1
template <int DIM>
2
struct GaussXor {
3
    array<bitset<DIM>, DIM> mat; // 存增广矩阵, 行从 0 到 n - 1, 列从 0 到 m
4

5
    array<int, DIM> ans;
6
    // 解异或线性方程组, 无解返回-1, 唯一解返回0, 无穷解返回自由变元数
7
    int solve(int n, int m) {
8
        // n 个方程, m 个未知量(即系数矩阵为 n 行 m 列)
9
        int r = 0, c = 0;
10
        for (; r < n && c < m; c++) {
11
            int pivot = r;
12
            while (pivot < n && !mat[pivot][c]) {
13
                pivot++;
14
            }
15
            if (pivot == n) {
16
                continue;
17
            }
18
            if (pivot != r) {
19
                swap(mat[pivot], mat[r]);
20
            }
21
            for (int i = 0; i < n; i++) {
22
                if (i != r && mat[i][c]) {
23
                    mat[i] ^= mat[r];
24
                }
25
            }
26
            r++;
27
        }
28

29
        for (int i = r; i < n; i++) {
30
            if (mat[i][m]) {
31
                return -1;
32
            }
33
        }
34
        if (r < m) {
35
            return m - r;
36
        }
37
        for (int i = 0; i < m; i++) {
38
            if (i < n) {
39
                ans[i] = mat[i][m];
40
            }
41
        }
42
        return 0;
43
    }
44
};

异或线性基#

1
template <int DIM = 64>
2
struct LinearBasis {
3
    array<ull, DIM> p = {};
4
    int cnt = 0;
5
    bool zero = 0;
6

7
    bool insert(ull x) {
8
        for (int i = DIM - 1; i >= 0; i--) {
9
            if ((x & (1ULL << i)) == 0) {
10
                continue;
11
            }
12
            if (!p[i]) {
13
                p[i] = x;
14
                cnt++;
15
                return true;
16
            }
17
            x ^= p[i];
18
        }
19
        zero = true;
20
        return false;
21
    }
22

23
    bool check(ull x) {
24
        for (int i = DIM - 1; i >= 0; i--) {
25
            if ((x & (1ULL << i)) == 0) {
26
                continue;
27
            }
28
            if (!p[i]) {
29
                return false;
30
            }
31
            x ^= p[i];
32
        }
33
        return true;
34
    }
35

36
    ull max() {
37
        ull ret = 0;
38
        for (int i = DIM - 1; i >= 0; i--) {
39
            if ((ret ^ p[i]) > ret) {
40
                ret ^= p[i];
41
            }
42
        }
43
        return ret;
44
    }
45

46
    ull min() {
47
        if (zero) {
48
            return 0;
49
        }
50
        for (int i = 0; i < DIM; i++) {
51
            if (p[i]) {
52
                return p[i];
53
            }
54
        }
55
    }
56

57
    // 合并
58
    friend bool operator+(const LinearBasis& a, const LinearBasis& b) {
59
        LinearBasis ret = a;
60
        for (int i = 0; i < DIM; i++) {
61
            if (b.p[i]) {
62
                ret.insert(b.p[i]);
63
            }
64
        }
65
        return ret;
66
    }
67
};

字符串#

滚动哈希#

1
// 随机数生成器
2

3
// mod = (1ULL << 61) - 1, 便于rolling hash使用
4
struct modint61 {
5
    static constexpr u64 mod = (1ULL << 61) - 1;
6
    u64 val;
7

8
    constexpr modint61(ll x = 0) : val(((x %= static_cast<ll>(mod)) < 0) ? (x + static_cast<ll>(mod)) : x) {}
9

10
    modint61 &operator+=(const modint61 &a) {
11
        val = ((val += a.val) >= mod) ? (val - mod) : val;
12
        return *this;
13
    }
14
    modint61 &operator-=(const modint61 &a) {
15
        val = ((val -= a.val) >= mod) ? (val + mod) : val;
16
        return *this;
17
    }
18
    modint61 &operator*=(const modint61 &a) {
19
        const unsigned __int128 y = static_cast<unsigned __int128>(val) * a.val;
20
        val = (y >> 61) + (y & mod);
21
        val = (val >= mod) ? (val - mod) : val;
22
        return *this;
23
    }
24
    modint61 operator-() const { return modint61(val ? mod - val : u64(0)); }
25
    modint61 &operator/=(const modint61 &a) { return (*this *= a.inv()); }
26
    modint61 operator+(const modint61 &p) const { return modint61(*this) += p; }
27
    modint61 operator-(const modint61 &p) const { return modint61(*this) -= p; }
28
    modint61 operator*(const modint61 &p) const { return modint61(*this) *= p; }
29
    modint61 operator/(const modint61 &p) const { return modint61(*this) /= p; }
30
    bool operator==(const modint61 &p) const { return val == p.val; }
31
    bool operator!=(const modint61 &p) const { return val != p.val; }
32
    modint61 inv() const {
33
        ll a = val, b = mod, u = 1, v = 0, t;
34
        while (b > 0) {
35
            t = a / b;
36
            a -= t * b;
37
            u -= t * v;
38
            swap(a, b);
39
            swap(u, v);
40
        }
41
        return modint61(u);
42
    }
43
    modint61 pow(ll n) const {
44
        assert(n >= 0);
45
        modint61 ret(1), mul(val);
46
        while (n > 0) {
47
            if (n & 1) ret *= mul;
48
            mul *= mul, n >>= 1;
49
        }
50
        return ret;
51
    }
52
};
53
struct StringHash {
54
    using m61 = modint61;
55
    static m61 base1, base2;
56
    static vector<m61> pow1, pow2;
57

58
    vector<m61> h1, h2;
59

60
    // 初始化base
61
    static void init() {
62
        if (base1.val != 0) return;
63
        base1 = m61(RNG(1000, m61::mod - 1));
64
        base2 = m61(RNG(1000, m61::mod - 1));
65
        while (base2 == base1) base2 = m61(RNG(1000, modint61::mod - 1));
66
        pow1 = {1};
67
        pow2 = {1};
68
    }
69

70
    // 预处理幂次, 按需扩展
71
    static void expand(int n) {
72
        if (pow1.empty()) init();
73
        while ((int)pow1.size() <= n) {
74
            pow1.push_back(pow1.back() * base1);
75
            pow2.push_back(pow2.back() * base2);
76
        }
77
    }
78

79
    StringHash() {}
80

81
    // O(|s|)
82
    StringHash(const string& s) {
83
        int n = s.size();
84
        expand(n);
85
        h1.resize(n + 1);
86
        h2.resize(n + 1);
87
        for (int i = 0; i < n; i++) {
88
            h1[i + 1] = h1[i] * base1 + s[i];
89
            h2[i + 1] = h2[i] * base2 + s[i];
90
        }
91
    }
92

93
    // 查询子串 s[l...r-1] 的哈希值 (左闭右开)
94
    // 返回 pair<u64, u64>
95
    pair<u64, u64> query(int l, int r) const {
96
        // H[l...r-1] = H[r] - H[l] * B^(r-l)
97
        m61 res1 = h1[r] - h1[l] * pow1[r - l];
98
        m61 res2 = h2[r] - h2[l] * pow2[r - l];
99
        return {res1.val, res2.val};
100
    }
101

102
    // 合并两个 Hash 值 (对应 s1 + s2)
103
    // len_right: 右边子串的长度
104
    static pair<u64, u64> merge(pair<u64, u64> h_left, pair<u64, u64> h_right, int len_right) {
105
        expand(len_right);
106
        m61 res1 = m61(h_left.first) * pow1[len_right] + modint61(h_right.first);
107
        m61 res2 = m61(h_left.second) * pow2[len_right] + modint61(h_right.second);
108
        return {res1.val, res2.val};
109
    }
110
};
111

112
// 全局成员定义
113
inline modint61 StringHash::base1;
114
inline modint61 StringHash::base2;
115
inline vector<modint61> StringHash::pow1;
116
inline vector<modint61> StringHash::pow2;

杂#

随机数生成#

1
// random number giver
2
inline u64 RNG_64() {
3
    static u64 x_ = u64(chrono::duration_cast<chrono::nanoseconds>(chrono::high_resolution_clock::now().time_since_epoch()).count()) * 10150724397891781847ULL;
4
    x_ ^= x_ << 7;
5
    return x_ ^= x_ >> 9;
6
}
7
// gives random number of [0, lim)
8
inline u64 RNG(u64 lim) { return RNG_64() % lim; }
9
// gives random number of [l, r)
10
inline ll RNG(ll l, ll r) { return l + RNG_64() % (r - l); }