CC BY 4.0 （除特别声明或转载文章外）
如果这篇博客帮助到你，可以请我喝一杯咖啡~
本页是一个汇总页。
我的算法竞赛模板
本页是一个汇总页。
{% for post in site.tags[‘算法竞赛模板’] %}
{{ post.title }}
{{post.content}}
{% endfor %}
数据结构
以下数据结构均采用 ll 作为值类型，应用时根据需求调整。
1
2
3
typedef long long ll;
const ll INF = 1e9; //表示（值）正无穷，且两个正无穷相加不会溢出
const int NPOS = -1; //表示（下标）不存在
离散化
在 vector 基础上的离散化，使用 push_back()向其中插值，init()排序并离散化，ask 查询离散化之后的值，at/[]运算符查离散前的值。
1
2
3
4
5
struct Ranker : vector<ll>
{
void init() { sort(begin(), end()), resize(unique(begin(), end()) - begin()); }
int ask(ll x) const { return lower_bound(begin(), end(), x) - begin(); }
};
并查集
1
2
3
4
5
6
7
8
9
10
11
12
13
14
struct UnionfindSet : vector<int>
{
UnionfindSet(int n) : vector<int>(n)
{
for (int i = 0; i < n; ++i)
at(i) = i;
}
void merge(int u, int w)
{
if (w = ask(w), u = ask(u), w != u)
at(w) = u;
}
int ask(int u) { return at(u) != u ? at(u) = ask(at(u)) : u; }
};
单调队列和单调栈
每次取队尾就是单调栈，取队头就是单调队列。
1
2
3
4
5
6
7
8
9
10
11
typedef pair<ll, int> pli;
struct Monotone : deque<pli>
{
void push(const pli &p, int k)
{
while (!empty() && back().first >= p.first)
pop_back();
for (push_back(p); p.second - front().second >= k;)
pop_front();
}
};
ST 表
$O(n\log n)$预处理，$O(1)$求静态区间最小值。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
/*
//可选优化
#define log2(n) LOG2[n]
struct Log : vector<ll>
{
Log(int N, ll E) : vector<ll>(N, -1)
{
for (int i = 1; i < N; ++i)
at(i) = at(i / E) + 1;
}
} LOG2(N, 2);
*/
struct SparseTable
{
vector<vector<ll>> f;
SparseTable(const vector<ll> &a) : f(log2(a.size()) + 1, a)
{
for (int k = 0; k + 1 < f.size(); ++k)
for (int i = 0; i + (1 << k) < a.size(); ++i)
f[k + 1][i] = min(f[k][i], f[k][i + (1 << k)]);
}
ll ask(int l, int r)
{
int k = log2(r - l + 1);
return min(f[k][l], f[k][r + 1 - (1 << k)]);
}
};
树状数组
模板中 Base 是对应的基础版本，支持单点修改区间查询。
一维
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
struct Fenwick
{
struct BaseFenwick
{
vector<ll> v;
BaseFenwick(int n) : v(n, 0) {}
void add(int x, ll w)
{
for (; x < v.size(); x += x & -x)
v[x] += w;
}
ll ask(int x)
{
ll ans = 0;
for (; x; x -= x & -x)
ans += v[x];
return ans;
}
};
pair<BaseFenwick, BaseFenwick> p;
Fenwick(int n) : p(n, n) {}
void add(int x, ll w) { p.first.add(x, w), p.second.add(x, x * w); }
void add(int l, int r, ll w) { add(l, w), add(r + 1, -w); }
ll ask(int x) { return (x + 1) * p.first.ask(x) - p.second.ask(x); }
ll ask(int l, int r) { return ask(r) - ask(l - 1); }
};
二维
高维的数据结构只要每一维维护低一维的数据（树套树）即可。其余数据结构亦同理。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
struct Fenwick2
{
struct BaseFenwick2
{
vector<Fenwick> v;
BaseFenwick2(int r, int c) : v(r, c) {}
void add(int x, int b, int t, ll w)
{
for (; x < v.size(); x += x & -x)
v[x].add(b, t, w);
}
ll ask(int x, int b, int t)
{
ll ans = 0;
for (; x; x -= x & -x)
ans += v[x].ask(b, t);
return ans;
}
};
pair<BaseFenwick2, BaseFenwick2> p;
Fenwick2(int r, int c) : p(BaseFenwick2(r, c), BaseFenwick2(r, c)) {}
void add(int x, int b, int t, ll w) { p.first.add(x, b, t, w), p.second.add(x, b, t, x * w); }
void add(int l, int b, int r, int t, ll w) { add(l, b, t, w), add(r + 1, b, t, -w); } //(l,b)~(r,t)
ll ask(int x, int b, int t) { return (x + 1) * p.first.ask(x, b, t) - p.second.ask(x, b, t); }
ll ask(int l, int b, int r, int t) { return ask(r, b, t) - ask(l - 1, b, t); }
};
动态开点线段树
使用示例，支持区间线性变换、区间查询（最大值最小值区间和）。
这样写改可持久化也很方便，只要改down函数为每次都新建节点即可，示例。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
struct SegmentTree
{
struct Seg
{
int l, r;
ll min, max, sum;
void upd(ll mul, ll add) { min = min * mul + add, max = max * mul + add, sum = sum * mul + add * (r - l + 1); }
friend Seg operator+(const Seg &lc, const Seg &rc) { return {lc.l, rc.r, std::min(lc.min, rc.min), std::max(lc.max, rc.max), lc.sum + rc.sum}; }
};
struct Node : Seg
{
int lc, rc;
ll mul, add;
};
vector<Node> v;
SegmentTree(int l, int r) { build(l, r); }
void build(int l, int r)
{
int rt = v.size();
v.push_back({});
v[rt].Seg::operator=({l, r, 0, 0, 0});
v[rt].lc = v[rt].rc = NPOS;
v[rt].mul = 1, v[rt].add = 0;
//if (l < r) //动态开点的时候注释掉本行和下一行
//down(rt), v[rt].Seg::operator=(v[v[rt].lc] + v[v[rt].rc]);
}
void down(int rt)
{
int m = v[rt].l + (v[rt].r - v[rt].l >> 1);
if (v[rt].lc == NPOS)
v[rt].lc = v.size(), build(v[rt].l, m);
//else //非可持久化的时候注释掉本行和下一行
//v.push_back(v[v[rt].lc]), v[rt].lc = v.size() - 1;
if (v[rt].rc == NPOS)
v[rt].rc = v.size(), build(m + 1, v[rt].r);
//else //非可持久化的时候注释掉本行和下一行
//v.push_back(v[v[rt].rc]), v[rt].rc = v.size() - 1;
upd(v[v[rt].lc].l, v[v[rt].lc].r, v[rt].mul, v[rt].add, v[rt].lc);
upd(v[v[rt].rc].l, v[v[rt].rc].r, v[rt].mul, v[rt].add, v[rt].rc);
v[rt].mul = 1, v[rt].add = 0;
}
void upd(int l, int r, ll mul, ll add, int rt = 0)
{
if (l <= v[rt].l && v[rt].r <= r)
return v[rt].mul *= mul, v[rt].add = v[rt].add * mul + add, v[rt].upd(mul, add);
down(rt);
if (r <= v[v[rt].lc].r)
upd(l, r, mul, add, v[rt].lc);
else if (l >= v[v[rt].rc].l)
upd(l, r, mul, add, v[rt].rc);
else
upd(l, v[v[rt].lc].r, mul, add, v[rt].lc), upd(v[v[rt].rc].l, r, mul, add, v[rt].rc);
v[rt].Seg::operator=(v[v[rt].lc] + v[v[rt].rc]);
}
Seg ask(int l, int r, int rt = 0)
{
if (l <= v[rt].l && v[rt].r <= r)
return v[rt];
down(rt);
if (r <= v[v[rt].lc].r)
return ask(l, r, v[rt].lc);
if (l >= v[v[rt].rc].l)
return ask(l, r, v[rt].rc);
return ask(l, v[v[rt].lc].r, v[rt].lc) + ask(v[v[rt].rc].l, r, v[rt].rc);
}
};
无旋 Treap
按子树大小分裂
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
struct FhqTreap
{
struct Node
{
int ch[2], siz, rev;
ll key, val, min, add;
void upd(ll v, int r)
{
val += v, min += v, add += v;
if (r)
rev ^= 1, swap(ch[0], ch[1]);
}
};
vector<Node> v;
int root;
FhqTreap() : v(1), root(0) {}
void down(int k)
{
if (!k)
return;
for (int i = 0, *ch = v[k].ch; i < 2; ++i)
if (ch[i])
v[ch[i]].upd(v[k].add, v[k].rev);
v[k].add = v[k].rev = 0;
}
void up(int k)
{
if (!k)
return;
v[k].siz = 1, v[k].min = v[k].val;
for (int i = 0, *ch = v[k].ch; i < 2; ++i)
if (ch[i])
v[k].siz += v[ch[i]].siz, v[k].min = min(v[k].min, v[ch[i]].min);
}
int merge(int a, int b)
{
if (!a || !b)
return a + b;
if (v[a].key < v[b].key)
return down(a), v[a].ch[1] = merge(v[a].ch[1], b), up(a), a;
return down(b), v[b].ch[0] = merge(a, v[b].ch[0]), up(b), b;
}
void split(int a, int s, int &l, int &r)
{
if (!s)
l = 0, r = a;
else if (v[v[a].ch[0]].siz < s)
down(a), split(v[a].ch[1], s - v[v[a].ch[0]].siz - 1, v[a].ch[1], r), up(l = a);
else
down(a), split(v[a].ch[0], s, l, v[a].ch[0]), up(r = a);
}
void push_back(ll d) { v.push_back(Node{{0, 0}, 1, 0, rand(), d, d, d}), root = merge(root, v.size() - 1); }
void insert(int x, ll d)
{
v.push_back(Node{{0, 0}, 1, 0, rand(), d, d, d});
int a, b, c;
split(root, x - 1, a, b), root = merge(merge(a, v.size() - 1), b);
}
void erase(int x)
{
int a, b, c;
split(root, x, a, b), split(a, x - 1, a, c), root = merge(a, b);
}
Node ask(int l, int r)
{
int a, b, c;
split(root, r, b, c), split(b, l - 1, a, b);
Node ret = v[b];
return root = merge(merge(a, b), c), ret;
}
void upd(int l, int r, ll add, int rev)
{
int a, b, c;
split(root, r, b, c), split(b, l - 1, a, b), v[b].upd(add, rev), root = merge(merge(a, b), c);
}
void revolve(int l, int r, int d)
{
int a, b, c, e = r - l + 1;
split(root, r, b, c), split(b, l - 1, a, b), split(b, (e - d % e) % e, b, e), root = merge(merge(a, merge(e, b)), c);
}
};
按值大小分裂
使用示例，即排序树。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
struct FhqTreap
{
struct Node
{
int ch[2], siz;
ll key, val;
};
vector<Node> v;
int root;
FhqTreap() : v(1), root(0) {}
void up(int k) { v[k].siz = v[v[k].ch[0]].siz + v[v[k].ch[1]].siz + 1; }
int merge(int a, int b)
{
if (!a || !b)
return a + b;
if (v[a].key < v[b].key)
return v[a].ch[1] = merge(v[a].ch[1], b), up(a), a;
return v[b].ch[0] = merge(a, v[b].ch[0]), up(b), b;
}
void splitVal(int a, ll w, int &l, int &r) //按值将树划分，使得左子树上的值恰小于w
{
if (!a)
l = r = 0;
else if (v[a].val > w)
splitVal(v[a].ch[0], w, l, v[a].ch[0]), up(r = a);
else
splitVal(v[a].ch[1], w, v[a].ch[1], r), up(l = a);
}
void insert(ll x)
{
int a, b;
v.push_back(Node{{0, 0}, 1, rand(), x}), splitVal(root, x, a, b), root = merge(merge(a, v.size() - 1), b);
}
void erase(ll x)
{
int a, b, c;
splitVal(root, x, a, b), splitVal(a, x - 1, a, c), root = merge(merge(a, merge(v[c].ch[0], v[c].ch[1])), b);
}
ll kth(int k)
{
for (int u = root, ls;;)
{
if (ls = v[v[u].ch[0]].siz, ls + 1 == k)
return v[u].val;
if (ls < k)
k -= ls + 1, u = v[u].ch[1];
else
u = v[u].ch[0];
}
}
int lower_bound(ll x) { return upper_bound(x - 1); }
int upper_bound(ll x)
{
int a, b, ret;
return splitVal(root, x, a, b), ret = v[a].siz + 1, root = merge(a, b), ret;
}
};
莫队
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
struct Mo
{
struct Query
{
int l, r, id;
bool operator<(const Query &n) const
{
return l / BS != n.l / BS ? l < n.l : r < n.r;
}
};
vector<Query> q;
int L, R;
void query(int l, int r) { q.push_back(Query{l, r, q.size()}); }
void rev(int x) {}
void cal(int id) {}
void ask()
{
L = 0, R = -1;
sort(q.begin(), q.end());
for (int i = 0; i < q.size(); ++i)
{
while (L < q[i].l)
rev(L++);
while (L > q[i].l)
rev(--L);
while (R < q[i].r)
rev(++R);
while (R > q[i].r)
rev(R--);
cal(q[i].id);
}
}
};
带修莫队
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
struct Mo
{
struct Update
{
int pos, NEW, OLD;
};
struct Query
{
int t, l, r, id;
bool operator<(const Query &n) const
{
return l / BS != n.l / BS ? l < n.l : r / BS != n.r / BS ? r < n.r : t < n.t;
}
};
vector<Update> cq;
vector<Query> q;
int T, L, R;
Mo() : cq(1) {}
void query(int x, int y) { q.push_back(Query{cq.size() - 1, x, y, q.size()}); }
void update(int x, int y) { cq.push_back(Update{x, y, t[x]}), t[x] = y; }
void set(int x, int d)
{
if (vis[x])
return rev(x), a[x] = d, rev(x);
a[x] = d;
}
void rev(int x) {}
void cal(int id) {}
void ask()
{
T = L = 0, R = -1;
sort(q.begin(), q.end());
for (int i = 0; i < q.size(); ++i)
{
while (T < q[i].t)
++T, set(cq[T].pos, cq[T].NEW);
while (T > q[i].t)
set(cq[T].pos, cq[T].OLD), --T;
while (L < q[i].l)
rev(L++);
while (L > q[i].l)
rev(--L);
while (R < q[i].r)
rev(++R);
while (R > q[i].r)
rev(R--);
cal(q[i].id);
}
}
};
树上莫队
按照欧拉序分块，使用 Tarjan 在生成欧拉序的同时预处理所有询问的 lca，预处理时间复杂度$O(n+q)$。 h 为查询图，即如果有一个询问(u,v)，即在 h 上连$u\to v,v\to u$。多个询问边有序插入 h。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
struct TreeMo : Graph
{
struct Query
{
int l, r, lca, id;
bool operator<(const Query &b) const
{
return l / BS != b.l / BS ? l < b.l : r < b.r;
}
};
vector<Query> q;
vector<int> dfp, dfi, dfo;
UnionFindSet ufs;
Graph h;
int L, R;
TreeMo(int n) : Graph(n), h(n), dfp(n * 2 + 1), dfi(n), dfo(n), ufs(n) {}
void query(int x, int y)
{
h.add(Edge{x, y}), h.add(Edge{y, x});
q.push_back(Query{0, 0, 0, q.size()});
}
void rev(int x) {}
void cal(int id) {}
void dfs(int u, int &cnt)
{
dfp[dfi[u] = ++cnt] = u;
for (int i = 0, k, to; i < v[u].a.size(); ++i)
if (k = v[u].a[i], to = e[k].second, !dfi[to])
dfs(to, cnt), ufs.merge(u, to);
dfp[dfo[u] = ++cnt] = u;
for (int i = 0, k, to, id; i < h.v[u].a.size(); ++i)
if (k = h.v[u].a[i], id = k / 2, to = h.e[k].second, dfo[to])
{
q[id].lca = ufs.fa(to);
q[id].l = q[id].lca != u ? dfo[u] : dfi[u];
q[id].r = dfi[to];
}
}
void ask(int root = 1)
{
dfs(root, BS = 0), BS = sqrt(BS);
sort(q.begin(), q.end());
L = 0, R = -1;
for (int i = 0; i < q.size(); ++i)
{
while (L < q[i].l)
rev(dfp[L++]);
while (L > q[i].l)
rev(dfp[--L]);
while (R < q[i].r)
rev(dfp[++R]);
while (R > q[i].r)
rev(dfp[R--]);
if (q[i].lca != dfp[L])
rev(q[i].lca);
cal(q[i].id);
if (q[i].lca != dfp[L])
rev(q[i].lca);
}
}
};
树上带修莫队
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
struct CapitalTreeMo : Graph
{
struct Update
{
int pos, NEW, OLD;
};
struct Query
{
int t, l, r, lca, id;
bool operator<(const Query &b) const
{
return l / BS != b.l / BS ? l < b.l : r / BS != b.r / BS ? r < b.r : t < b.t; //在BZOJ4129上去掉r/BS还快100ms?
}
};
vector<Update> cq;
vector<Query> q;
vector<int> dfp, dfi, dfo;
UnionFindSet ufs;
Graph h;
int T, L, R;
CapitalTreeMo(int n) : cq(1), Graph(n), h(n), dfp(n * 2 + 1), dfi(n), dfo(n), ufs(n) {}
void query(int x, int y)
{
h.add(Edge{x, y}), h.add(Edge{y, x});
q.push_back(Query{cq.size() - 1, 0, 0, 0, q.size()});
}
void update(int x, int y)
{
cq.push_back(Update{x, y, t[x]}), t[x] = y;
}
void dfs(int u, int &cnt)
{
dfp[dfi[u] = ++cnt] = u;
for (int i = 0, k, to; i < v[u].a.size(); ++i)
if (k = v[u].a[i], to = e[k].second, !dfi[to])
dfs(to, cnt), ufs.merge(u, to);
dfp[dfo[u] = ++cnt] = u;
for (int i = 0, k, to, id; i < h.v[u].a.size(); ++i)
if (k = h.v[u].a[i], id = k / 2, to = h.e[k].second, dfo[to])
{
q[id].lca = ufs.fa(to);
q[id].l = q[id].lca != u ? dfo[u] : dfi[u];
q[id].r = dfi[to];
}
}
void set(int u, int d)
{
if (vis[u])
return rev(u), a[u] = d, rev(u);
a[u] = d;
}
void rev(int u) {}
void cal(int id) {}
void ask(int root = 1)
{
dfs(root, BS = 0), BS = sqrt(BS);
sort(q.begin(), q.end());
T = L = 0, R = -1;
for (int i = 0; i < q.size(); ++i)
{
while (T < q[i].t)
++T, set(cq[T].pos, cq[T].NEW);
while (T > q[i].t)
set(cq[T].pos, cq[T].OLD), --T;
while (L < q[i].l)
rev(dfp[L++]);
while (L > q[i].l)
rev(dfp[--L]);
while (R < q[i].r)
rev(dfp[++R]);
while (R > q[i].r)
rev(dfp[R--]);
if (q[i].lca != dfp[L])
rev(q[i].lca);
cal(q[i].id);
if (q[i].lca != dfp[L])
rev(q[i].lca);
}
}
};
字符串/模式匹配
HashString
使用示例，如果要修改模数或者直接使用unsigned long long的自然溢出的话直接修改 Mod 即可。
使用unsigned long long的自然溢出快了 5 倍，但是容易被卡。
1
2
3
4
5
6
7
8
9
10
11
12
13
struct HashString : Mod
{
vector<ll> f, p;
HashString(const string &s, ll M = 1e9 + 7, ll P = 131) : Mod(M), f(s.size() + 1), p(s.size() + 1, 1)
{
for (int i = 0; i < s.size(); ++i)
{
f[i + 1] = add(mul(f[i], P), s[i]);
p[i + 1] = mul(p[i], P);
}
}
ll ask(int pos, int len) { return add(f[pos + len], -mul(f[pos], p[len])); } //从pos位置开始的长度为len的子串的hash值
};
KMP
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
struct KMP
{
vector<int> next;
string pattern;
KMP(const string &pattern) : pattern(pattern), next(pattern.size() + 1, -1)
{
for (int i = 0, j = -1; i < pattern.size();)
{
if (j == -1 || pattern[i] == pattern[j])
next[++i] = ++j;
else
j = next[j];
}
}
int find_in(const string &text)
{
for (int i = 0, j = 0;;)
{
if (j == pattern.size())
return i - j; //不return可得到所有匹配地址
if (i == text.size())
return -1;
if (j == -1 || text[i] == pattern[j])
++i, ++j;
else
j = next[j];
}
}
};
AC 自动机
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
struct AhoCorasick
{
struct Node
{
int ch[26], val, f, last;
int &to(char c)
{
return ch[c - 'a'];
} //如果不确定c的范围，使用map
};
vector<Node> v;
AhoCorasick() : v(1) {}
void getFail()
{
for (deque<int> q(1, v[0].last = v[0].f = 0); !q.empty(); q.pop_front())
for (char c = 'a'; c <= 'z'; ++c)
{
int r = q.front(), u = v[r].to(c), w = v[r].f;
if (!r && u)
{
q.push_back(u);
v[u].f = v[u].last = 0;
continue;
}
if (!u)
{
v[r].to(c) = v[w].to(c);
continue;
}
q.push_back(u);
while (w && !v[w].to(c))
w = v[w].f;
v[u].f = v[w].to(c);
v[u].last = v[v[u].f].val ? v[u].f : v[v[u].f].last;
}
}
void add(const string &s, int val, int u = 0)
{
for (int i = 0; i < s.size(); u = v[u].to(s[i++]))
if (!v[u].to(s[i]))
{
v[u].to(s[i]) = v.size();
v.push_back(Node());
}
v[u].val = val;
}
bool find_in(const string &s, int u = 0) //调用需要调用`getFail()`生成失配函数。
{
for (int i = 0; i < s.size(); ++i)
if (u = v[u].to(s[i]),
v[u].val || v[u].last)
return 1;
return 0;
}
};
暴力回文
时间复杂度$O(n^2)$，常数低，但会被ababababa这样的数据卡。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
int palindrome(const char *s)
{
int ans = 0;
for (int i = 0, b, e; s[i]; ++i)
{
for (b = i; s[i] == s[i + 1];)
++i;
for (e = i + 1; b && s[b - 1] == s[e];)
--b, ++e;
if (ans < e - b)
ans = e - b; //此时[b,e)为最大回文区间
}
return ans;
}
线性回文
对于一个位置 i，[i−f[i]+1,i+f[i]−1]是最长的以 i 为中心的奇回文串，g[i]−i 是最长的以 i 为开头的回文串长度。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
struct Manacher
{
vector<int> t, f, g;
Manacher(const string &s) : t(s.size() + 1 << 1, 0), f(t), g(t) //t初始值为s中没有出现过的值，g开始为0
{
for (int i = 0; i < s.size(); ++i)
t[i + 1 << 1] = s[i];
for (int i = 1, p = 0, m = 0; i < t.size(); ++i)
{
for (f[i] = i < m ? min(f[2 * p - i], m - i) : 1;
0 < i - f[i] && i + f[i] < t.size() &&
t[i - f[i]] == t[i + f[i]];)
++f[i];
if (m < i + f[i])
m = i + f[p = i];
}
for (int i = 2; i < t.size(); ++i)
if (g[i - f[i] + 1] < i + 1)
g[i - f[i] + 1] = i + 1;
for (int i = 1; i < t.size(); ++i)
if (g[i] < g[i - 1])
g[i] = g[i - 1];
}
int ask(int l, int r) //多次询问可做一个ST表
{
int ans = 0;
for (int i = l + 1 << 1, e = r + 1 << 1; i <= e; i += 2)
if (ans < g[i] - i)
ans = g[i] - i;
return ans;
}
};
后缀数组
m：字符集大小。
s：字符串，其中最后一位为加入的 0。
sa[i]：字典序第 i 小的是哪个后缀。
rk[i]：后缀 i 的排名。
h[i]：lcp(sa[i],sa[i−1])。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
struct SufArr
{
vector<int> sa, rk, h;
SufArr(const vector<int> &s, int m) : sa(s.size(), 0), rk(s), h(s.size(), 0)
{
vector<int> cnt(s.size() + m, 0);
for (int i = 0; i < s.size(); ++i)
++cnt[rk[i]];
for (int i = 1; i < m; ++i)
cnt[i] += cnt[i - 1];
for (int i = 0; i < s.size(); ++i)
sa[--cnt[rk[i]]] = i;
for (int k = 1, j = 0; k <= s.size() && j < s.size() - 1; k <<= 1)
{
for (int i = 0; i < s.size(); ++i)
{
if (j = sa[i] - k, j < 0)
j += s.size();
h[cnt[rk[j]]++] = j;
}
cnt[0] = sa[h[0]] = j = 0;
for (int i = 1; i < s.size(); ++i)
{
if (rk[h[i]] != rk[h[i - 1]] || rk[h[i] + k] != rk[h[i - 1] + k])
cnt[++j] = i;
sa[h[i]] = j;
}
swap(rk, sa), swap(sa, h);
}
for (int i = 0, k = 0, j = rk[0]; i < s.size() - 1; ++i, ++k)
for (; ~k && s[i] != s[sa[j - 1] + k]; j = rk[sa[j] + 1], --k)
h[j] = k;
}
};
图论
这里用类似邻接表的方法存图。有的算法可能需要邻接矩阵，详见线性代数部分。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
struct Graph
{
struct Vertex
{
vector<int> o, i; //相关出边和入边编号
int siz, dep, top, dfn; //树链剖分中使用，依次代表子树节点数、深度、所在链的顶端节点、dfs序
};
struct Edge : pair<int, int>
{
ll len, cap; //边长、容量，图论算法使用
};
vector<Vertex> v; //点集
vector<Edge> e; //边集
Graph(int n) : v(n) {}
void add(const Edge &ed)
{
if (ed.first == ed.second)
return; //如果有需要请拆点
v[ed.first].o.push_back(e.size());
v[ed.second].i.push_back(e.size());
e.push_back(ed);
}
int ch(int u, int i = 0) { return e[v[u].o[i]].second; } //u的第i个孩子节点
int fa(int u, int i = 0) { return e[v[u].i[i]].first; } //u的第i个父节点
};
最短路
Dijkstra 算法
使用示例，适用于边权为正的情况（无论有向图还是无向图），用于求单源最短路。
直接给出其优先队列优化的版本。另外，由于priority_queue并不提供修改优先级的操作，为避免重复扩展，这里选择将新元素直接插入队列并在运行时判断该点是否被处理，并不影响结果的正确性。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
struct Dijkstra : Graph
{
vector<ll> d;
vector<int> p;
Dijkstra(int n) : Graph(n) {}
void ask(int s)
{
d.assign(v.size(), INF);
p.assign(v.size(), e.size());
priority_queue<pair<ll, int>> q;
for (q.push(make_pair(d[s] = 0, s)); !q.empty();)
{
ll dis = -q.top().first;
int u = q.top().second;
if (q.pop(), d[u] < dis)
continue;
for (int i = 0, k, to; i != v[u].o.size(); ++i)
if (k = v[u].o[i], to = e[k].second,
d[to] > d[u] + e[k].len)
{
d[to] = d[u] + e[k].len, p[to] = k;
q.push(make_pair(-d[to], to));
}
}
}
};
BellmanFord 算法
使用示例，直接给出其队列优化、国内称之为 SPFA 算法的版本。较之 Dijkstra 算法，此算法不够快速稳定但是可以允许负边存在，当 s 到达负权回路时会直接返回 0。稀疏图上性能优秀。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
struct BellmanFord : Graph
{
vector<ll> d;
vector<int> p;
BellmanFord(int n) : Graph(n) {}
bool ask(int s)
{
d.assign(v.size(), INF);
p.assign(v.size(), e.size());
vector<int> cnt(v.size(), 0), flag(v.size(), d[s] = 0);
for (deque<int> q(cnt[s] = flag[s] = 1, s); !q.empty(); q.pop_front())
for (int u = q.front(), i = flag[u] = 0, k, to; i < v[u].o.size(); ++i)
if (k = v[u].o[i], to = e[k].second,
d[to] > d[u] + e[k].len)
{
d[to] = d[u] + e[k].len, p[to] = k;
if (!flag[to])
{
if (v.size() == ++cnt[to])
return 0;
flag[to] = 1, q.push_back(to);
}
}
return 1;
}
};
差分约束系统
按如下方式建图、跑 SPFA：
对每个不等式$x_i−x_j\leq d$，从$j$向$i$连一条边，边长为$d$。
若不等号的方向相反，即$x_i−x_j\geq d$，则在不等式两边同时乘以$-1$，变成$x_j−x_i\leq -d$，即从$i$到$j$连一条边，边长为$d$。
Floyed 求多源最短路
不连通的权置 INF。
1
2
3
4
5
6
7
8
void floyed(Mat &a, int n) // 其实一般手打即可，没必要套这个鬼模板…
{
for (int k = 0; k < n; ++k)
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
if (a[i][j] > a[i][k] + a[k][j])
a[i][j] = a[i][k] + a[k][j];
}
Astar 求 k 短路
使用示例，在这个比较坑的例子中需要在调用前补一句if(s==t)++k;。k 下标从 0 开始，即最短路==第 0 短路。
朴素的想法是使用priority_queue从原点出发向外探索，当取出终点 t 第 k 次时就得到第 k 短路，类似 bfs 的思想，缺陷是越往后状态数越多。
我们在反向图上从$t\to s$跑 Astar 算法，通过优先展开到 s 近的状态，使搜索方向靠近答案，而不是一层一层全都展开，估价函数$f\approx g+h$，f 是估计的 s 到 t 的距离，g 是到达当前点已经点的花费，h 是预计剩下的花费。这里 h 取当前点的距离到 s 距离，可通过从 s 跑一遍 Dijkstra 可以预处理得出。
Astar 算法是只有到达终点的时候才能统计答案，这导致可能拓展很多个状态才能得到一个用来更新答案的有效状态。例如一个 n 元环，当我们到达终点之后,可能还要拓展 n 次才能得到下一个状态，于是时间复杂度就被卡到$O(nk)$。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
struct Astar : Dijkstra
{
vector<ll> ans;
Astar(int n) : Dijkstra(n) {}
void ask(int s, int t, int k)
{
Dijkstra::ask(s);
ans.assign(k, INF);
if (d[t] == INF)
return;
vector<int> cnt(v.size(), 0);
priority_queue<pair<ll, int>> q;
for (q.push(make_pair(-d[t], t)); cnt[s] < k && !q.empty();)
{
ll dis = -q.top().first;
int u = q.top().second;
if (u == s)
ans[cnt[s]] = dis;
if (q.pop(), ++cnt[u] > k)
continue;
for (int i = 0, k; i < v[u].i.size(); ++i)
k = v[u].i[i], q.push(make_pair(d[u] - d[e[k].first] - e[k].len - dis, e[k].first));
}
}
};
网络流
ISAP 求最大流
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
struct ISAP : Graph
{
ll flow;
vector<ll> f;
vector<int> h, cur, gap;
ISAP(int n) : Graph(n) {}
void add(Edge ed)
{
Graph::add(ed);
swap(ed.first, ed.second), ed.cap = 0;
Graph::add(ed);
}
ll dfs(int s, int u, int t, ll r)
{
if (r == 0 || u == t)
return r;
ll _f, _r = 0;
for (int &i = cur[u], k; i < v[u].o.size(); ++i)
if (k = v[u].o[i], h[u] == h[e[k].second] + 1)
{
_f = dfs(s, e[k].second, t, min(r - _r, e[k].cap - f[k]));
f[k] += _f, f[k ^ 1] -= _f, _r += _f;
if (_r == r || h[s] >= v.size())
return _r;
}
if (!--gap[h[u]])
h[s] = v.size();
return ++gap[++h[u]], cur[u] = 0, _r;
}
void ask(int s, int t)
{
h.assign(v.size(), 0);
cur.assign(v.size(), 0);
gap.assign(v.size() + 2, 0);
/*
for (deque<int> q(h[t] = gap[t] = 1, t); !q.empty(); q.pop_front()) //优化，加了能快一点
for (int i = 0, u = q.front(), k, to; i < v[u].o.size(); ++i)
if (to = e[v[u].o[i]].second, !h[to])
++gap[h[to] = h[u] + 1], q.push_back(to);
*/
for (f.assign(e.size(), flow = 0); h[s] < v.size();)
flow += dfs(s, s, t, INF);
}
};
PrimalDual 求费用流
使用示例，定义一条边的费用为流量*边长，求总费用最小的最大流。性能优秀，只能跑非负权图。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
struct PrimalDual : Graph
{
ll flow, cost;
vector<ll> f;
PrimalDual(int n) : Graph(n) {}
void add(Edge ed)
{
Graph::add(ed);
swap(ed.first, ed.second), ed.cap = 0, ed.len *= -1;
Graph::add(ed);
}
void ask(int s, int t) //询问s到t的最小费用最大流，答案存在flow、cost中
{
vector<int> p(v.size(), e.size());
vector<ll> d(v.size(), INF), h(v.size(), 0);
for (f.assign(e.size(), flow = cost = 0);;)
{
priority_queue<pair<ll, int>> q;
for (q.push(make_pair(d[s] = 0, s)); !q.empty();)
{
ll dis = -q.top().first;
int u = q.top().second;
if (q.pop(), d[u] < dis)
continue;
for (int i = 0, k, to; i < v[u].o.size(); ++i)
if (k = v[u].o[i], to = e[k].second,
e[k].cap > f[k] && d[to] > d[u] + e[k].len + h[u] - h[to])
{
d[to] = d[u] + e[k].len + h[u] - h[to], p[to] = k;
q.push(make_pair(-d[to], to));
}
}
if (d[t] == INF)
return;
for (int i = 0; i < d.size(); ++i)
if (d[i] != INF)
h[i] += d[i], d[i] = INF;
ll _f = INF;
for (int u = t; u != s; u = e[p[u]].first)
_f = min(_f, e[p[u]].cap - f[p[u]]);
for (int u = t; u != s; u = e[p[u]].first)
cost += _f * e[p[u]].len, f[p[u]] += _f, f[p[u] ^ 1] -= _f;
flow += _f;
}
}
};
EK 求费用流
使用示例，定义一条边的费用为流量*边长，求总费用最小的最大流。BellmanFord 算法找增广路，可能被卡但是可以跑负费用流（最大费用流）。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
struct EdmondKarp : Graph
{
ll flow, cost;
vector<ll> f;
EdmondKarp(int n) : Graph(n) {}
void add(Edge ed)
{
Graph::add(ed);
swap(ed.first, ed.second), ed.cap = 0, ed.len *= -1;
Graph::add(ed);
}
void ask(int s, int t)
{
vector<int> p(v.size(), e.size());
for (f.assign(e.size(), flow = cost = 0);;)
{
vector<ll> d(v.size(), INF);
vector<int> flag(v.size(), d[s] = 0);
for (deque<int> q(flag[s] = 1, s); !q.empty(); q.pop_front())
for (int u = q.front(), i = flag[u] = 0, k, to; i < v[u].o.size(); ++i)
if (k = v[u].o[i], to = e[k].second,
e[k].cap > f[k] && d[to] > d[u] + e[k].len)
{
d[to] = d[u] + e[k].len, p[to] = k;
if (!flag[to])
q.push_back(to), flag[to] = 1;
}
if (d[t] == INF)
return;
ll _f = INF;
for (int u = t; u != s; u = e[p[u]].first)
_f = min(_f, e[p[u]].cap - f[p[u]]);
for (int u = t; u != s; u = e[p[u]].first)
cost += _f * e[p[u]].len, f[p[u]] += _f, f[p[u] ^ 1] -= _f;
flow += _f;
}
}
};
ZKW 求费用流
使用示例，定义一条边的费用为流量*边长，求总费用最小的最大流。
对于最终流量较大，而费用取值范围不大的图，或者是增广路径比较短的图（如二分图），zkw 算法都会比较快，原因是充分发挥优势。比如流多说明可以同一费用反复增广，费用窄说明不用改太多距离标号就会有新增广路，增广路径短可以显著改善最坏情况，因为即使每次就只增加一条边也可以很快凑成最短路。如果恰恰相反，流量不大，费用不小，增广路还较长，就不适合 zkw 算法了。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
struct ZKW : Graph
{
ll flow, cost;
vector<ll> h, f;
vector<int> vis;
ZKW(int n) : Graph(n) {}
void add(Edge ed)
{
Graph::add(ed);
swap(ed.first, ed.second), ed.cap = 0, ed.len *= -1;
Graph::add(ed);
}
ll dfs(int u, int t, ll r)
{
if (r == 0 || u == t)
return r;
if (vis[u])
return 0;
ll _f = vis[u] = 1, _r = 0;
for (int i = 0, k; r > _r && i < v[u].o.size(); ++i)
if (k = v[u].o[i], h[e[k].second] + e[k].len == h[u])
_f = dfs(e[k].second, t, min(r - _r, e[k].cap - f[k])), f[k] += _f, f[k ^ 1] -= _f, _r += _f;
return _r;
}
void ask(int s, int t)
{
h.assign(v.size(), 0);
vis.assign(v.size(), 0);
for (f.assign(e.size(), flow = cost = 0);;)
{
ll _f = dfs(s, t, INF), d = INF;
flow += _f, cost += _f * h[s];
for (int u = 0; u < v.size(); ++u)
for (int i = 0, k; vis[u] && i < v[u].o.size(); ++i)
if (k = v[u].o[i], !vis[e[k].second] && e[k].cap > f[k])
d = min(d, e[k].len + h[e[k].second] - h[e[k].first]);
if (d == INF)
return;
for (int i = 0; i < v.size(); ++i)
if (vis[i])
h[i] += d, vis[i] = 0;
}
}
};
上下界有源汇网络流
T 向 S 连容量正无穷的边，将有源汇转化为无源汇。每条边容量减去下界，设$in[i]$表示流入 i 的下界之和减去流出 i 的下界之和。新建超级源汇 SS,TT，对于$in[i]>0$的点，SS 向 i 连容量为$in[i]$的边。对于$in[i]<0$的点，i 向 TT 连容量为$−in[i]$的边。
求出以 SS,TT 为源汇的最大流，如果等于$\sum ini$则存在可行流。再求出以 S,T 为源汇的最大流即为最大流。
费用流：建完图后等价于求以 SS,TT 为源汇的的费用流。
上下界费用流：先把下界的费用加入答案。
判断边是否属于某一割集
在残余网络 (还有流量的边) 中求强连通分量，顶点不在同一 SCC 且满流的边。
判断边是否为全部最小割集的边： 在上一条的基础上，还要满足起点与 S 在同一 SCC，且终点与 T 在同一 SCC。
线性规划转费用流
首先添加松弛变量，将不等号都变为等号。分别用下一个式子减去上一个式子，如果每个变量只出现了两次且符号一正一负，那么可以转化为费用流。对于每个式子建立一个点，那么每个变量对应一条边，从一个点流出，向另一个点流入。这样，对于等式右边的常数，如果是正的，对应从源点向该点连一条流量 C，费用 0 的边；如果是负的对应从该点向汇点连一条流量 −C，费用 0 的边。对于每个变量，从它系数为正的式子向系数为负的式子连一条容量为 inf，费用为它在目标函数里系数的边。这样网络流模型就构造完毕了。
欧拉路
使用示例，给定无孤立结点图 G，若存在一条路，经过图中每边一次且仅一次，该条路称为欧拉路。
- 无向图：当仅当该图所有顶点的度数为偶数（此时为回路），或除两个度数为奇数外（作为路径的起点和终点）、其余全为偶数。
- 有向图：当仅当该图所有顶点出度=入度（此时为回路），或一个顶点出度=入度+1（作为起点）、另一个顶点入度=出度+1（作为终点）、其他顶点出度=入度。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
struct Fleury : Graph
{
vector<int> vis, cur, p;
Fleury(int n) : Graph(n) {}
void dfs(int u)
{
for (int &i = cur[u], k; i < v[u].i.size(); ++i) //遍历原图的反向图，这里加了一个「当前弧」优化
if (k = v[u].i[i], !vis[k] && !vis[k ^ 1]) //无向图需要同时检查反向边未被访问过
{
vis[k] = 1;
dfs(e[k].first);
p.push_back(k);
}
}
void ask() //查询欧拉回路，路径上边的序号按顺序存在p中
{
vis.assign(e.size(), 0), cur.assign(v.size(), 0), p.clear();
for (int i = 0; i < v.size(); ++i)
if (v[i].i.size() % 2)
return dfs(i);
dfs(0);
}
};
混合图欧拉回路判定
首先给无向边随便定一个方向，设$\deg x$为 x 连出去的边数 − 连入 x 的边数。若存在$\deg x$为奇数，或者图不连通，则无解。否则建立源点 S，汇点 T。对于一个点 x，若$\deg x>0$，则 S 向 x 连边，容量$\frac{\deg x}{2}$；若$\deg x<0$，则 x 向 T 连边，容量$-\frac{\deg x}{2}$。 对于一条定了向的无向边$x\to y$，x 向 y 连边，容量 1，求出最大流，若与 S 和 T 连的每条边都满流，则有解。
连通性
无向图求割和双连通分量
- 割边：在连通图中，删除了连通图的某条边后，图不再连通。这样的边被称为割边，也叫做桥。
- 割点：在连通图中，删除了连通图的某个点以及与这个点相连的边后，图不再连通。这样的点被称为割点。构造 dfs 搜索树，在树上有两类节点可以成为割点：
- 对根节点 u，若其有两棵或两棵以上的子树，则该根结点 u 为割点；
- 对非根非叶节点 u，若其中的某棵子树的节点均没有指向 u 的祖先节点的回边，说明删除 u 之后，根结点与该棵子树的节点不再连通；则节点 u 为割点。
对于一个无向图的子图，当删除其中任意一条边后，不改变图内点的连通性，这样的子图叫做边的双连通子图。而当子图的边数达到最大时，叫做边的双连通分量。原理是图中所有割边再求一次 SCC，可直接使用求 SCC 的代码。
对于一个无向图的子图，当删除其中任意一个点后，不改变图内点的连通性，这样的子图叫做点的双连通子图。而当子图的边数达到最大时，叫做点的双连通分量。下面给出求点双连通分量的代码。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
struct BiconnectedConnectedComponenet : Graph
{
vector<int> dep, bid, stak, cutPoint, cutEdge; //bid边的端点所属双连通块
BiconnectedConnectedComponenet(int n) : Graph(n) {}
void ask()
{
dep.assign(v.size(), v.size());
bid.assign(e.size(), e.size());
cutPoint.assign(v.size(), 0);
cutEdge.assign(e.size(), 0);
for (int i = 0; i < v.size(); ++i)
if (dep[i] == v.size())
dfs(i, v.size());
}
int dfs(int u, int fa)
{
int low = dep[u] = fa != v.size() ? dep[fa] + 1 : 0;
for (int i = 0, k, to, lowto; i < v[u].o.size(); ++i)
if (k = v[u].o[i], to = e[k].second, to != fa)
{
if (dep[to] == v.size())
{
stak.push_back(k);
low = min(low, lowto = dfs(to, u));
if (lowto > dep[u])
cutEdge[k] = cutEdge[k ^ 1] = 1;
if (lowto >= dep[u])
for (cutPoint[u] = fa != v.size() || i;;)
{
int x = stak.back();
stak.pop_back();
bid[x] = bid[x ^ 1] = k;
if (x == k)
break;
}
}
else if (dep[to] < dep[u])
{
stak.push_back(k);
low = min(low, dep[to]);
}
}
return low;
}
};
双连通图的构造
先求出所有的桥，然后删除这些桥边，剩下的每个连通块都是一个双连通子图。把每个双连通子图收缩为一个顶点，再把桥边加回来，最后的这个图一定是一棵树，边连通度为 1。统计出树中度为 1 的节点的个数，即为叶节点的个数，记为leaf。至少在树上添加(leaf+1)/2条边，就能使树达到边双连通：先把两个最近公共祖先最远的两个叶节点之间连接一条边，这样可以把这两个点到祖先的路径上所有点收缩到一起，因为一个形成的环一定是双连通的；然后再找两个最近公共祖先最远的两个叶节点，这样一对一对找完，恰好是(leaf+1)/2次，把所有点收缩到了一起。
无向图求边双连通分量&有向图求强连通分量
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
struct StronglyConnectedComponenet : Graph
{
vector<int> dep, sid, stak; //sid=点所属连通块内一点
StronglyConnectedComponenet(int n) : Graph(n) {}
void ask()
{
dep.assign(v.size(), v.size());
sid.assign(v.size(), v.size());
for (int i = 0; i < v.size(); ++i)
if (dep[i] == v.size())
dfs(i, v.size());
}
int dfs(int u, int fa)
{
int low = dep[u] = fa != v.size() ? dep[fa] + 1 : 0;
stak.push_back(u);
for (int i = 0, k, to; i < v[u].o.size(); ++i)
if (k = v[u].o[i], to = e[k].second, to != fa /*, 1*/) //求强连通分量把注释去掉，即允许走回边
{
if (dep[to] == v.size())
low = min(low, dfs(to, u));
else if (sid[to] == v.size())
low = min(low, dep[to]);
}
if (low == dep[u])
for (;;)
{
int x = stak.back();
stak.pop_back();
sid[x] = u;
if (x == u)
break;
}
return low;
}
};
2-SAT
使用示例，n 个布尔变量$x_0\dots x_{n-1}$，逻辑表达式$Y=(A_0+B_0)(A_1+B_1)\dots(A_{m-1}+B_{m-1})$，其中$A_i,B_i\in\lbrace x_j,\overline{x_j}\rbrace $，判断是否存在$x_0\dots x_{n-1}$的取值使得 Y 值为 1。因为$A+B=(\overline A\to B)(\overline B\to A)$，所以对于一个要求$A+B$，我们连$\overline A\to B,\overline B\to A$两条边。如果有一条边$A\to B$，意味着如果 A 成立那么 B 必然成立。若$\exists i,x_i,\overline{x_i}\in$同一 SCC，则不存在。
模板要求顶点数v.size()是偶数，i与i ^ 1互为相反变量且奇数下标代表布尔变量为真的情况。最后输出结果在ok中。不满足对任意 iok[i] ^ ok[i ^ 1] == 1时，ask函数返回false。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
struct TwoSat : StronglyConnectedComponenet
{
vector<int> ok;
TwoSat(int n) : StronglyConnectedComponenet(n) {}
void addOR(int a, int b) { add({a ^ 1, b}), add({b ^ 1, a}); }
void addXOR(int a, int b) { addOR(a, b), addOR(a ^ 1, b ^ 1); }
void addXNOR(int a, int b) { addOR(a, b ^ 1), addOR(a ^ 1, b); }
int ask()
{
StronglyConnectedComponenet::ask();
for (int i = 0; i < v.size(); i += 2)
if (sid[i] == sid[i ^ 1])
return 0;
vector<vector<int>> g(v.size());
vector<int> ind(v.size(), 0), cf(v.size(), 0), stak;
for (int i = 0; i < v.size(); ++i)
{
cf[sid[i]] = sid[i ^ 1];
for (int j = 0, k; j < v[i].o.size(); ++j)
if (sid[e[k = v[i].o[j]].second] != sid[i])
g[sid[e[k].second]].push_back(sid[i]), ++ind[sid[i]];
}
for (int i = 0; i < v.size(); ++i)
if (sid[i] == i && !ind[i])
stak.push_back(i);
for (ok.assign(v.size(), -1); !stak.empty();)
{
int u = stak.back();
stak.pop_back();
if (ok[u] < 0)
ok[u] = 1, ok[cf[u]] = 0;
for (int i = 0; i < g[u].size(); ++i)
if (!--ind[g[u][i]])
stak.push_back(g[u][i]);
}
for (int i = 0; i < v.size(); ++i)
if (i != sid[i])
ok[i] = ok[sid[i]];
return 1;
}
};
2-SAT 暴力搜索求字典序最小的解
使用示例，时间复杂度最坏为为$O(VE)$。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
struct TwoSat : Graph
{
vector<int> ok;
TwoSat(int n) : Graph(n) {}
void addOR(int a, int b) { add({a ^ 1, b}), add({b ^ 1, a}); }
void addXOR(int a, int b) { addOR(a, b), addOR(a ^ 1, b ^ 1); }
void addXNOR(int a, int b) { addOR(a, b ^ 1), addOR(a ^ 1, b); }
int dfs(int u, vector<int> &stak)
{
if (ok[u ^ 1])
return 0;
if (ok[u])
return 1;
ok[u] = 1;
stak.push_back(u);
for (int i = 0; i < v[u].o.size(); ++i)
if (!dfs(e[v[u].o[i]].second, stak))
return 0;
return 1;
}
int ask()
{
ok.assign(v.size(), 0);
for (int i = 0; i < v.size(); i += 2)
if (!ok[i] && !ok[i ^ 1])
{
vector<int> stak;
if (!dfs(i, stak))
{
for (int j = 0; j < stak.size(); ++j)
ok[stak[j]] = 0;
if (!dfs(i ^ 1, stak))
return 0;
}
}
return 1;
}
};
二分图
二分图的一个等价定义是：不含有含奇数条边的环的图。
完美匹配图中所有的顶点都是匹配点。
二分图的最小点覆盖（最小割）是指最少的顶点数使得二分图 G 中的每条边都至少与其中一个点相关联。二分图中，最小割=最大匹配。
二分图的最小边覆盖（最大独立集）是指用尽量少的不相交简单路径覆盖二分图中的所有顶点。二分图中，最小点覆盖+最小边覆盖=总点数。
Hall 定理：二分图中的两部分顶点组成的集合分别为 X,Y ，则有一组无公共点的边，一端恰好为组成 X 的点的充分必要条件是：X 中的任意 k 个点至少与 Y 中的 k 个点相邻。对于区间图只需要考虑极端情况，线段树维护。
关键点：一定在最大匹配中的点。 求出任意一个最大匹配，那么只需要考虑哪些匹配点不一定在。 假设是考虑左侧的点，右侧类似： 将匹配边从右往左，非匹配边从左到右，从左侧每个未匹配点开始 DFS 到的匹配点都不是关键点。
关键边：求出任意一个最大匹配，将匹配边从右到左，剩余边从左到右，求出 SCC。 对于一条边：若它位于当前匹配中，那么若两端点位于同一 SCC，则是可能在，否则必定在；若它不位于当前匹配中，那么若两端点位于同一 SCC，则是可能在，否则必定不在。
Hungary 求最大匹配
左边 nl 个点$0\dots nl-1$，右边 nr 个点$0\dots nr-1$，取n=max(nl,nr)，左 i 右 j 间相连则g[i].push_back(j);。
生成fl[i]表示左边第 i 个匹配右边第fl[i]个，fr同理。时间复杂度$O(n^3)$。未匹配的值为fl[i] == N。匹配是一个边集，其中任意两条边都没有公共顶点；扫一遍fl或fr判断有多少不等于N即为最大匹配数。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
struct Hungary
{
vector<int> g[N];
int n, fl[N], fr[N], vr[N];
int dfs(int x, int rt)
{
for (auto y : g[x])
if (vr[y] != rt)
if (vr[y] = rt, fr[y] == N || dfs(fr[y], rt))
return fl[fr[y] = x] = y, 1;
return 0;
}
void ask()
{
fill(fl, fl + n, N), fill(fr, fr + n, N), fill(vr, vr + n, N);
for (int i = 0; i < n; ++i)
if (fl[i] == N)
dfs(i, i); //找完全匹配时如果返回值为0可直接return
}
};
HopcroftKarp 求最大匹配
使用示例，时间复杂度$O(\sqrt{\vert V\vert }\vert E\vert )$，稀疏图上效果明显。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
struct HopcroftKarp
{
vector<int> g[N];
int n, fl[N], fr[N], hl[N], hr[N], q[N];
bool dfs(int x)
{
int c = hl[x] + 1, y = hl[x] = N;
for (auto y : g[x])
if (hr[y] == c)
if (hr[y] = N, fr[y] == N || dfs(fr[y]))
return fl[fr[y] = x] = y, 1;
return 0;
}
void ask()
{
fill(fl, fl + n, N), fill(fr, fr + n, N);
for (int x = 0; x < n; ++x)
for (auto y : g[x])
if (fr[y] == N)
{
fl[fr[y] = x] = y;
break;
}
for (int ql, qr, ok;;)
{
fill(hl, hl + n, N), fill(hr, hr + n, N);
for (int x = ql = qr = ok = 0; x < n; ++x)
if (fl[x] == N)
hl[q[qr++] = x] = 0;
while (ql < qr)
for (auto y : g[q[ql++]])
if (hr[y] == N)
{
hr[y] = hl[q[ql - 1]] + 1;
if (fr[y] == N)
ok = 1;
else if (hl[fr[y]] == N)
hl[q[qr++] = fr[y]] = hr[y] + 1;
}
if (!ok)
return;
for (int x = 0; x < n; ++x)
if (fl[x] == N)
dfs(x);
}
}
};
KuhnMunkres 求最优完备匹配
使用示例。
左 x 右 y 代价a[x][y]。最大费用流时，a 初始化 0；最大费用最大流时，a 初始化-N*INF。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
struct KuhnMunkres
{
ll a[N][N], hl[N], hr[N], slk[N];
int n, fl[N], fr[N], vl[N], vr[N], pre[N], q[N], ql, qr;
int check(int i)
{
if (vl[i] = 1, fl[i] != N)
return vr[q[qr++] = fl[i]] = 1;
while (i != N)
swap(i, fr[fl[i] = pre[i]]);
return 0;
}
void bfs(int s)
{
fill(slk, slk + n, INF), fill(vl, vl + n, 0), fill(vr, vr + n, 0);
for (vr[q[ql = 0] = s] = qr = 1;;)
{
for (ll d; ql < qr;)
for (int i = 0, j = q[ql++]; i < n; ++i)
if (!vl[i] && slk[i] >= (d = hl[i] + hr[j] - a[i][j]))
if (pre[i] = j, d)
slk[i] = d;
else if (!check(i))
return;
ll d = INF;
for (int i = 0; i < n; ++i)
if (!vl[i] && d > slk[i])
d = slk[i];
for (int i = 0; i < n; ++i)
{
if (vl[i])
hl[i] += d;
else
slk[i] -= d;
if (vr[i])
hr[i] -= d;
}
for (int i = 0; i < n; ++i)
if (!vl[i] && !slk[i] && !check(i))
return;
}
}
void ask()
{
fill(fl, fl + n, N), fill(fr, fr + n, N), fill(hr, hr + n, 0);
for (int i = 0; i < n; ++i)
hl[i] = *max_element(a[i], a[i] + n);
for (int j = 0; j < n; ++j)
bfs(j);
}
};
带花树
一般图最大匹配
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
struct Blossom : Graph
{
vector<int> f;
Blossom(int n) : Graph(n) {}
void ask()
{
f.assign(v.size(), v.size());
vector<int> vis(v.size());
for (int s = 0, cnt = vis.back(); s < v.size(); ++s)
if (f[s] == v.size())
{
UnionfindSet ufs(v.size());
vector<int> pre(v.size(), v.size()), flag(v.size(), 2);
for (deque<int> q(flag[s] = 1, s); !q.empty() && f[s] == v.size(); q.pop_front())
for (int i = 0, x = q.front(), y, a, b; i < v[x].o.size(); ++i)
if (y = e[v[x].o[i]].second, y != f[x] && flag[y] && ufs.ask(x) != ufs.ask(y))
{
if (flag[y] == 1)
{
for (a = x, b = y, ++cnt;; swap(a, b))
if (a != v.size())
{
if (vis[a = ufs.ask(a)] == cnt)
break;
vis[a] = cnt, a = f[a] != v.size() ? pre[f[a]] : v.size();
}
if (ufs.ask(x) != a)
pre[x] = y;
if (ufs.ask(y) != a)
pre[y] = x;
for (int p[2] = {x, y}, j = 0; j < 2; ++j)
for (int x = p[j], y, z; x != a; ufs.merge(y, x), ufs.merge(x = z, y))
{
if (ufs.ask(z = pre[y = f[x]]) != a)
pre[z] = y;
if (!flag[y])
flag[y] = 1, q.push_back(y);
if (!flag[z])
flag[z] = 1, q.push_back(z);
}
}
else if (f[y] != v.size())
pre[y] = x, flag[y] = 0, flag[f[y]] = 1, q.push_back(f[y]);
else
{
for (pre[y] = x; y != v.size();)
swap(y, f[f[y] = pre[y]]);
break;
}
}
}
}
};
树形图
最小生成树
无向图
同时给出 Prim 算法（生成新树）、Kruskal 算法（消耗小）。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
struct Prim : Graph
{
struct DisGreater
{
bool operator()(const Edge &e1, const Edge &e2)
{
return e1.len > e2.len;
}
};
ll ans;
vector<int> vis;
priority_queue<Edge, vector<Edge>, DisGreater> q;
Prim(const Graph &g, int root) : Tree(n), ans(0), vis(g.v.size(), 0) //生成新树，每条边都要有等长反向边
{
for (insert(root, g); !q.empty();)
{
Edge ed = q.top();
if (q.pop(), !vis[ed.second])
insert(ed.second, g), ans += ed.len, add(ed);
}
}
void insert(int u, const Graph &g) //把点和对应的相连的边加入集合
{
vis[u] = 1;
for (int i = 0, k; i < g.v[u].o.size(); ++i)
if (k = g.v[u].o[i], !vis[g.e[k].second])
q.push(g.e[k]);
}
};
ll kruskal(vector<Edge> &e, int n) //会清空边集e，每条边被认作无向边
{
ll ret = 0;
UnionFindSet ufs(n);
for (sort(e.begin(), e.end(), DisGreater()); !e.empty(); e.pop_back())
if (ufs.fa(e.back().from) != ufs.fa(e.back().to))
{
ufs.merge(e.back().from, e.back().to);
ret += e.back().len;
}
return /*ufs.siz>1?INF:*/ ret; //视情况选择去注释
}
有向图
指定以 root 为根，如果没有限定根那么新建一个虚拟点作为根，向所有边连边长最大边长+1的边，在最后生成的图中去掉此边。时间复杂度$O(VE)$。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
ll zhuLiu(vector<Edge> &e, int root, int n) //不存在返回INF
{
for (ll ret = 0;;)
{
vector<ll> in(n, INF);
vector<int> pre(n, n);
for (int i = 0, to; i < e.size(); ++i)
{
if (e[i].first == (to = e[i].second))
swap(e[i--], e.back()), e.pop_back();
else if (in[to] > e[i].len)
in[to] = e[i].len, pre[to] = e[i].first;
}
for (int i = in [root] = 0; i < n; ++i)
if (in[i] == INF)
return INF;
vector<int> id(n, n), vis(n, n);
int tn = 0;
for (int i = 0, v; i < n; ++i)
{
for (ret += in [v = i]; vis[v] != i && id[v] == n && v != root; v = pre[v])
vis[v] = i;
if (v != root && id[v] == n)
{
for (int u = pre[v]; u != v; u = pre[u])
id[u] = tn;
id[v] = tn++;
}
}
if (!tn)
return ret;
for (int i = 0; i < n; ++i)
if (id[i] == n)
id[i] = tn++;
for (int i = 0, v; i < e.size(); ++i)
if ((e[i].first = id[e[i].first]) != (e[i].second = id[v = e[i].second]))
e[i].len -= in[v];
n = tn, root = id[root];
}
}
树链剖分与 LCA
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
struct Diagram : Graph
{
Fenwick data; //暂用树状数组作为默认数据结构
Diagram(const Graph &g, int root) : Graph(g.v.size()), data(g.v.size())
{
build(root, g);
int cnt = v[root].dfn = v[root].dep = 1;
dfs(v[root].top = root, cnt);
}
void build(int u, const Graph &g) //无向图dfs建树，且重边在最前，u为根节点
{
v[u].siz = 1;
for (int i = 0, k, to; i < g.v[u].o.size(); ++i)
if (k = g.v[u].o[i], to = g.e[k].second, !v[to].siz) //没访问过的点siz默认0
{
build(to, g);
v[u].siz += v[to].siz;
Graph::add(g.e[k]);
if (v[ch(u)].siz < v[to].siz) //重边移到最前
swap(v[u].o.front(), v[u].o.back());
}
}
void dfs(int u, int &cnt)
{
for (int i = 0, to; i < v[u].o.size(); ++i)
{
v[to = ch(u, i)].dfn = ++cnt;
v[to].top = i ? to : v[u].top;
v[to].dep = v[u].dep + 1;
dfs(to, cnt);
}
}
int lca(int x, int y)
{
for (; v[x].top != v[y].top; x = fa(v[x].top))
if (v[v[x].top].dep < v[v[y].top].dep)
swap(x, y);
if (v[x].dep < v[y].dep)
swap(x, y);
return y;
}
ll ask(int x, int y)
{
ll ans = 0;
for (; v[x].top != v[y].top; x = fa(v[x].top))
{
if (v[v[x].top].dep < v[v[y].top].dep)
swap(x, y);
ans += data.ask(v[v[x].top].dfn, v[x].dfn);
}
if (v[x].dep < v[y].dep)
swap(x, y);
return ans += data.ask(v[y].dfn, v[x].dfn);
}
void add(int x, int y, ll pv)
{
for (; v[x].top != v[y].top; x = fa(v[x].top))
{
if (v[v[x].top].dep < v[v[y].top].dep)
swap(x, y);
data.add(v[v[x].top].dfn, v[x].dfn, pv);
}
if (v[x].dep < v[y].dep)
swap(x, y);
data.add(v[y].dfn, v[x].dfn, pv);
}
};
点剖（点分治）
使用示例，零号点为虚节点。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
struct TreeDiv : Graph
{
int root;
vector<int> mx, siz, vis;
TreeDiv(int n) : Graph(n), mx(n, n), siz(n), vis(n) {}
void dfsRoot(int u, int fa)
{
for (int i = mx[u] = siz[u] = 0, k, to; i < v[u].o.size(); ++i)
if (k = v[u].o[i], to = e[k].second, to != fa && !vis[to])
if (dfsRoot(to, u), siz[u] += siz[to], mx[u] < siz[to])
mx[u] = siz[to];
if (mx[u] < mx[0] - ++siz[u])
mx[u] = mx[0] - siz[u];
if (mx[root] >= mx[u])
root = u;
}
void dfsDis(int u, int fa, ll d)
{
//用d更新答案
for (int i = 0, k, to; i < v[u].o.size(); ++i)
if (k = v[u].o[i], to = e[k].second, to != fa && !vis[to])
dfsDis(to, u, d + e[k].len);
}
int cal(int u, ll d) //返回符合要求的点对数
{
return dfsDis(u, 0, d), /*得到答案*/;
}
void dfs(int u = 1)
{
dfsRoot(u, root = 0), ans += cal(u = root, 0), vis[u] = 1;
for (int i = 0, k, to; i < v[u].o.size(); ++i)
if (k = v[u].o[i], to = e[k].second, !vis[to])
ans -= cal(to, e[k].len), mx[0] = siz[to], dfs(to);
}
};
数论
辗转相除法
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
ll gcd(ll a, ll b) { return b ? gcd(b, a % b) : a; } //a、b的最大公约数
ll lcm(ll a, ll b) { return a / gcd(a, b) * b; } //a、b的最小公倍数
ll gcd(ll a, ll b, ll &x, ll &y) //扩展欧几里得，引用返回a*x+b*y=gcd(a,b)绝对值之和最小的解
{
if (!a)
return x = 0, y = 1, b;
ll d = gcd(b % a, a, y, x);
return x -= b / a * y, d;
}
ll gcd(ll a, ll b) //无取模求gcd
{
for (ll t = 1, c, d;;)
{
if (a == b)
return t * a;
if (a < b)
swap(a, b);
if (a & 1)
c = 0;
else
a >>= 1, c = 1;
if (b & 1)
d = 0;
else
b >>= 1, d = 1;
if (c && d)
t <<= 1;
else if (!c && !d)
a -= b;
}
}
裴蜀定理
对任何$a,b\in Z$和它们的最大公约数$d$，关于未知数$x,y$的线性不定方程（称为裴蜀等式）：$ax+by=c$当仅当$d\vert c$，可知有无穷多解。特别地，$ax+by=d$一定有解。
推论
$a,b$互质的充要条件是$ax+by=1$有整数解。
同余系运算
三种乘法的比较，结论是如果可以用__int128就用，否则就用第一种：
求乘法逆元的另外一种方法是用欧拉定理$x^{\phi(m)}\equiv1\pmod m$，x 的逆是$x^{\phi(m)-1}$。特别地，m 为素数时$\phi(m)=m-1$，此时 x 的逆就是pow(x,m-2,m)。
log 函数：m 为素数时求解模方程$a^x\equiv b\pmod m$。设 P 为质数，G 为 P 的原根，则$x^y\equiv b\pmod P$等价于$y\ ind\ x\equiv b\pmod{P-1}$，其中$G\ ind\ x\equiv x\pmod P$。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
struct Mod
{
const ll M, SM;
Mod(ll M) : M(M), SM(sqrt(M) + 0.5) {}
ll qadd(ll &a, ll b) const { return a += b, a >= M ? a -= M : a; } //假如a+b<2*M，就不必取模了，取模运算耗时很高
ll add(ll a, ll b) const { return qadd(a = (a + b) % M, M); } //考虑a和b不在同余系内甚至为负数的情况
ll mul(ll a, ll b) const { return add(a * b, M); } //ll类型是int的时候小心 a*b 爆精度
ll inv(ll a) const { return pow(a, M - 2); } //要求M为素数，否则return pow(a, phi(M) - 1);
ll pow(ll a, ll b) const
{
ll r = 1;
/*
if (b < 0)
b = -b, a = inv(a);
*/
for (a = add(a, M); b; b >>= 1, a = mul(a, a))
if (b & 1)
r = mul(r, a);
return r;
}
/*
ll mul(ll a, ll b) const { return add(a * b, -M * (ll)((long double)a / M * b)); } //long double 有效位数16~18位，模数过大时慎用！
ll mul(ll a, ll b) const //Head算法，无循环快速计算同余乘法，根据a*b是否爆ll替换a*b%M，需要a<M且b<M，可以调用时手动取模
{
ll c = a / SM, d = b / SM;
a %= SM, b %= SM;
ll e = add(add(a * d, b * c), c * d / SM * (SM * SM - M));
return add(add(a * b, e % SM * SM), add(c * d % SM, e / SM) * (SM * SM - M));
}
ll mul(ll a, ll b) const //龟速乘
{
ll r = 0;
for (a %= M; b; b >>= 1, qadd(a, a))
if (b & 1)
qadd(r, a);
return r;
}
ll inv(ll a) const //模m下a的乘法逆元，不存在返回-1（m为素数时a不为0必有逆元）
{
ll x, y, d = gcd(a, M, x, y);
return d == 1 ? add(x, M) : -1;
}
vector<ll> sol(ll a, ll b) const //解同余方程，返回ax=b(mod M)循环节内所有解
{
vector<ll> ans;
ll x, y, d = gcd(a, M, x, y);
if (b % d)
return ans;
ans.push_back(mul((b / d) % (M / d), x));
for (ll i = 1; i < d; ++i)
ans.push_back(add(ans.back(), M / d));
return ans;
}
ll log(ll a, ll b) const
{
unordered_map<ll, ll> x;
for (ll i = 0, e = 1; i <= SM; ++i, e = mul(e, a))
if (!x.count(e))
x[e] = i;
for (ll i = 0, v = inv(pow(a, SM)); i <= SM; ++i, b = mul(b, v))
if (x.count(b))
return i * SM + x[b];
return -1;
}
*/
};
中国剩余定理解同余方程组
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
ll crt(const vector<pair<ll, ll>> &v) //同余方程组，x%v[i].first==v[i].second，不存在返回-1
{
ll m = v[0].first, r = v[0].second, c, d, x, y, z;
for (int i = 1; i < v.size(); ++i)
{
if (c = v[i].second - r, d = gcd(m, v[i].first, x, y), c % d)
return -1;
gcd(m / d, z = v[i].first / d, x, y), r += c / d * x % z * m, r %= m *= z;
}
return r < 0 ? r + m : r;
}
struct Excrt
{
vector<ll> b2, a2;
int n;
void init()
{
scanf("%d", &n);
ll q, qq;
for (int i = 0; i < n; i++)
scanf("%lld%lld", &q, &qq), b2.push_back(q), a2.push_back(qq);
}
ll get()
{
ll M, ans, x = 0, y = 0;
M = 1, ans = 0;
for (int i = 0; i < n; i++)
{
ll a = M, b = b2[i], c = (a2[i] - ans % b + b) % b; //ax=c(mod b)
ll gcd = exgcd(a, b, x, y), bg = b / gcd;
if (c % gcd != 0)
{
return -1;
} //返回-1表示无解
x = mul(x, c / gcd, bg);
y = M / gcd * b;
ans = (x * M + ans) % y; //???k???????
if (ans < 0)
ans += y;
M = y; //M??k?m?lcm
}
return ans;
}
};
二次剩余
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
struct QuadraticRes
{
struct numb
{
ll x, y;
};
numb mul(const numb &a, const numb &b, ll p)
{
return (numb){(a.x * b.x % p + a.y * b.y % p * w % p) % p, (a.x * b.y + a.y * b.x) % p};
}
ll pow2(numb a, ll b, ll p)
{
numb ret = {1, 0}, tem = a;
while (b)
{
if (b & 1)
ret = mul(ret, tem, p);
tem = mul(tem, tem, p), b >>= 1;
}
return ret.x % p;
}
ll get(ll n, ll p)
{
if (Mod(p).pow(n, p / 2) == p - 1)
{
return -1;
} //此时无解
while (1)
{
ll a = rand();
w = (a * a - n + p) % p;
if (Mod(p).pow(w, p / 2) == p - 1)
{
return pow2((numb){a, 1}, p / 2 + 1, p); //还有一个解是p减去这个值
}
}
}
};
欧拉筛
欧拉函数$\phi(n)$是小于 n 的正整数中与 n 互素的数的数目。特别地，规定$\phi(1)=1$，易知 n>2 时都为偶数。
欧拉函数是积性函数，即对任意素数$p,q$满足下列关系：$\phi(pq)=\phi(p)\phi(q)=(p-1)(q-1)$对任何两个互质的正整数$x, m(m\geq2)$有欧拉定理：$x^{\phi(m)}\equiv1\pmod m$当 m 为素数 p 时，此式变为费马小定理：$x^{p-1}\equiv1\pmod p$利用欧拉函数和它本身不同质因数的关系，用筛法$O(N)$预处理某个范围内所有数的欧拉函数值，并求出素数表。同时，利用计算欧拉函数过程中求出的最小素因子 m，可以实现$O(log N)$的素因数分解。
同时求莫比乌斯函数$\mu(n)$，存在mu中。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
struct EulerSieve
{
vector<int> p, m, phi, mu; //素数序列，最小素因子，欧拉函数，莫比乌斯函数
EulerSieve(int N) : m(N, 0), phi(N, 0), mu(N, 0)
{
phi[1] = mu[1] = 1; //m[1]=0,m[i]==i可判断i是素数
for (long long i = 2, k; i < N; ++i) //防i*p[j]爆int
{
if (!m[i])
p.push_back(m[i] = i), phi[i] = i - 1, mu[i] = -1; //i是素数
for (int j = 0; j < p.size() && (k = i * p[j]) < N; ++j)
{
phi[k] = phi[i] * p[j];
if ((m[k] = p[j]) == m[i])
{
mu[k] = 0;
break;
}
phi[k] -= phi[i];
mu[k] = -mu[i];
}
}
}
};
直接求欧拉函数
1
2
3
4
5
6
7
8
9
10
11
ll phi(ll n)
{
ll phi = n;
for (ll i = 2; i * i <= n; ++i)
if (!(n % i))
for (phi = phi / i * (i - 1); !(n % i);)
n /= i;
if (n > 1)
phi = phi / n * (n - 1);
return phi;
}
常见数论函数变换
$\sum_{d\vert n}\mu(d)=[n=1]$
$\phi(n)=\sum_{i=1}^n[\gcd(i,n)=1]=\sum_{i=1}^n\sum_{k\mid i,k\mid n}\mu(k)=\sum_{k\mid n}\frac nk\mu(k)$
前缀和
欧拉函数前缀和$S_\phi(n)=\frac{(n+1)n}2-\sum_{d=1}^nS_\phi(\frac{n}{d})$
莫比乌斯函数前缀和$S_\mu(n)=1-\sum_{d=1}^nS_\mu(\frac{n}{d})$
莫比乌斯反演
若$f(n)=\sum_{d\vert n}g(d)$，则$g(n)=\sum_{d\vert n}\mu(d)f(\frac{n}{d})$
若$f(n)=\sum_{i=1}^nt(i)g(\frac{n}{i})$，则$g(n)=\sum_{i=1}^n\mu(i)t(i)f(\frac{n}{i})$（此时代$t(i)=[gcd(n,i)>1]$可得上式）
举例（其中$T=kd$）：
[\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)
=\sum_d d\sum_{i=1}^n\sum_{j=1}^m[\gcd (i,j)=d]
=\sum_{d}d\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}[\gcd (i,j)=1]
=\sum_{d}d\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}\sum_{k\mid i,k\mid j}\mu(k)
=\sum_d d\sum_k\mu(k)\sum_{k\mid i}^{\lfloor\frac nd\rfloor}\sum_{k\mid j}^{\lfloor\frac md\rfloor}
=\sum_{d}d\sum_k\mu(k)\lfloor\frac n{kd}\rfloor\lfloor\frac m{kd}\rfloor
=\sum_{T}\lfloor\frac nT\rfloor\lfloor\frac mT\rfloor\sum_{k\mid T}\frac Tk\mu(k)
=\sum_{T}\lfloor\frac nT\rfloor\lfloor\frac mT\rfloor\varphi(T)]
$\varphi(T)$可以使用线性筛预处理处理，我们就可以枚举$T$求上式了，时间复杂度$O(n)$。多组数据$n,m$询问上式，时间复杂度就变成了$O(Tn)$。事实上，$\lfloor\frac{n}{T}\rfloor$是不会轻易变化的，是过了连续的一段后才发生变化的，那么我们就可以计算出这一段的结束位置，对$\varphi$函数作前缀和，就可以直接分块了，这样的时间复杂度是$O(T\sqrt{n})$的。
PollardRho 大数素因子分解
时间复杂度$O(N^{1/4})$，数据多的时候可考虑欧拉筛优化。
注意这里模乘法很容易爆 long long，看情况选用快速乘法。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
struct PollardRho
{
bool isPrime(ll n, int S = 12) //MillerRabin素数测试，S为测试次数，用前S个素数测一遍，S=12可保证unsigned long long范围内无错；n<2请特判
{
static ll d, u, t, p[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
for (d = n - 1; !(d & 1);)
d >>= 1; //未对0，1做特判！
Mod mo(n);
for (int i = 0; i < S; ++i)
{
if (!(n % p[i]))
return n == p[i];
for (t = mo.pow(p[i], u = d); t != n - 1 && t != 1 && u != n - 1;)
t = mo.mul(t, t), u <<= 1;
if (t != n - 1 && !(u & 1))
return 0;
}
return 1;
}
void fac(ll n, vector<ll> &factor)
{
/*
if (n < e.m.size()) //欧拉筛预处理优化，可以防止很多小素数因子的情况
{
for (; n > 1; n /= e.m[n])
factor.push_back(e.m[n]);
return;
}
*/
if (isPrime(n))
return factor.push_back(n);
Mod mo(n);
for (ll c = 1;; ++c)
for (ll i = 0, k = 1, x = rand() % (n - 1) + 1, y, p;;)
{
if (++i == k)
y = x, k <<= 1;
if (x = mo.add(mo.mul(x, x), c), p = __gcd(abs(x - y), n), p == n)
break;
if (p > 1)
return fac(p, factor), fac(n / p, factor);
}
}
};
快速变换
蝴蝶变换（雷德变换）
保存 FTT 和 FNTT 时交换的对应位置（即保存的是置换）。
1
2
3
4
5
6
7
8
9
struct Rader : vector<int>
{
Rader(int n) : vector<int>(1 << int(ceil(log2(n))))
{
for (int i = at(0) = 0; i < size(); ++i)
if (at(i) = at(i >> 1) >> 1, i & 1)
at(i) += size() >> 1;
}
};
快速傅里叶变换
使用示例，高精度乘法。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
struct FFT : Rader
{
vector<complex<lf>> w;
FFT(int n) : Rader(n), w(size())
{
for (int i = 0; i < size(); ++i)
w[i] = polar(1.0, 2 * M_PI * i / size());
}
void fft(vector<complex<lf>> &x) const
{
for (int i = 0; i < x.size(); ++i)
if (i < at(i))
std::swap(x[i], x[at(i)]);
for (int i = 1; i < size(); i <<= 1)
for (int j = 0; j < i; ++j)
for (int k = j; k < size(); k += i << 1)
{
complex<lf> t = w[size() / (i << 1) * j] * x[k + i];
x[k + i] = x[k] - t, x[k] += t;
}
}
vector<ll> ask(const vector<ll> &a, const vector<ll> &b) const
{
vector<complex<lf>> xa(a.begin(), a.end()), xb(b.begin(), b.end());
xa.resize(size()), xb.resize(size()), fft(xa), fft(xb);
for (int i = 0; i < size(); ++i)
xa[i] *= xb[i];
fft(xa);
vector<ll> ans(size());
for (int i = 0; i < size(); ++i)
ans[i] = xa[(size() - i) & (size() - 1)].real() / size() + 0.5;
return ans;
}
vector<ll> askMod(const vector<ll> &a, const vector<ll> &b, ll M) const //任意模数
{
vector<complex<lf>> e(size()), c(size());
for (int i = 0; i < a.size(); ++i)
e[i].real(a[i] & 0x7fff), e[i].imag(a[i] >> 15);
for (int i = 0; i < b.size(); ++i)
c[i].real(b[i] & 0x7fff), c[i].imag(b[i] >> 15);
fft(e), fft(c);
vector<complex<lf>> d(c);
for (int i = 0; i < size(); ++i)
{
int fr = (size() - i) & (size() - 1);
c[i] *= 0.5 * complex<lf>(e[i].X + e[fr].X, e[i].Y - e[fr].Y);
d[i] *= 0.5 * complex<lf>(e[i].Y + e[fr].Y, e[fr].X - e[i].X);
}
fft(c), fft(d);
vector<ll> ans(size());
for (int i = 0; i < size(); ++i)
{
ll fr = (size() - i) & (size() - 1),
p = c[fr].X / size() + 0.5,
o = c[fr].Y / size() + 0.5,
x = d[fr].X / size() + 0.5,
u = d[fr].Y / size() + 0.5;
ans[i] = (p % M + ((o + x) % M << 15) + (u % M << 30)) % M;
}
return ans;
}
};
快速数论变换
原理和 FFT 相同，解决特殊情况下 FFT 的浮点误差，并且可以在同余系进行变换。
对于形如$m=2^nc+1$的费马素数，记其原根为 g，则旋转因子为$g^{(m-1)/n}$，满足$g^{m-1}=1$且$2^n\mid m-1$。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
struct FNTT : Rader, Mod
{
vector<ll> w;
FNTT(int N, ll M, ll G) : Rader(N), Mod(M), w(size(), pow(G, (M - 1) / size()))
{
for (int i = w[0] = 1; i < size(); ++i)
w[i] = mul(w[i], w[i - 1]);
}
void fntt(vector<ll> &x) const
{
for (int i = 0; i < size(); ++i)
if (i < at(i))
std::swap(x[i], x[at(i)]);
for (int i = 1, j; i < size(); i <<= 1)
for (int j = 0; j < i; ++j)
for (int k = j; k < size(); k += i << 1)
{
ll t = mul(w[size() / (i << 1) * j], x[k + i]);
qadd(x[k + i] = x[k], M - t), qadd(x[k], t);
}
}
vector<ll> ask(vector<ll> a, vector<ll> b) const
{
a.resize(size()), b.resize(size()), fntt(a), fntt(b);
for (int i = 0; i < size(); ++i)
a[i] = mul(a[i], b[i]);
fntt(a), reverse(a.begin() + 1, a.end());
ll u = inv(size());
for (int i = 0; i < size(); ++i)
a[i] = mul(a[i], u);
return a;
}
};
常见素数的原根
| r * 2 ^ k + 1 | r | k | g，mod(r * 2 ^ k + 1)的原根 |
|---|---|---|---|
| 3 | 1 | 1 | 2 |
| 5 | 1 | 2 | 2 |
| 17 | 1 | 4 | 3 |
| 97 | 3 | 5 | 5 |
| 193 | 3 | 6 | 5 |
| 257 | 1 | 8 | 3 |
| 7681 | 15 | 9 | 17 |
| 12289 | 3 | 12 | 11 |
| 40961 | 5 | 13 | 3 |
| 65537 | 1 | 16 | 3 |
| 786433 | 3 | 18 | 10 |
| 5767169 | 11 | 19 | 3 |
| 7340033 | 7 | 20 | 3 |
| 23068673 | 11 | 21 | 3 |
| 104857601 | 25 | 22 | 3 |
| 167772161 | 5 | 25 | 3 |
| 469762049 | 7 | 26 | 3 |
| 998244353 | 119 | 23 | 3 |
| 1004535809 | 479 | 21 | 3 |
| 2013265921 | 15 | 27 | 31 |
| 2281701377 | 17 | 27 | 3 |
| 3221225473 | 3 | 30 | 5 |
| 75161927681 | 35 | 31 | 3 |
| 77309411329 | 9 | 33 | 7 |
| 206158430209 | 3 | 36 | 22 |
| 2061584302081 | 15 | 37 | 7 |
| 2748779069441 | 5 | 39 | 3 |
| 6597069766657 | 3 | 41 | 5 |
| 39582418599937 | 9 | 42 | 5 |
| 79164837199873 | 9 | 43 | 5 |
| 263882790666241 | 15 | 44 | 7 |
| 1231453023109121 | 35 | 45 | 3 |
| 1337006139375617 | 19 | 46 | 3 |
| 3799912185593857 | 27 | 47 | 5 |
| 4222124650659841 | 15 | 48 | 19 |
| 7881299347898369 | 7 | 50 | 6 |
| 31525197391593473 | 7 | 52 | 3 |
| 180143985094819841 | 5 | 55 | 6 |
| 1945555039024054273 | 27 | 56 | 5 |
| 4179340454199820289 | 29 | 57 | 3 |
快速沃尔什变换
如果要在同余系中进行运算，则下面代码需要修改。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void fwt(vector<ll> &x, void f(ll &l, ll &r))
{
for (int i = 1; i < x.size(); i <<= 1)
for (int j = 0; j < i; ++j)
for (int k = j; k < x.size(); k += i << 1)
f(x[k], x[k + i]);
}
void fwt(ll *b, ll *e, void f(ll &l, ll &r)) //再给一个递归二分的代码便于理解
{
if (e - b < 2)
return;
ll *m = b + (e - b) / 2;
fwt(b, m, f), fwt(m, e, f);
while (m < e)
f(*(b++), *(m++));
}
AND
1
2
void tf(ll &l, ll &r) { l += r; }
void utf(ll &l, ll &r) { l -= r; }
OR
1
2
void tf(ll &l, ll &r) { r += l; }
void utf(ll &l, ll &r) { r -= l; }
XOR
1
2
3
4
5
6
void tf(ll &l, ll &r)
{
ll tl = l + r, tr = l - r;
l = tl, r = tr;
}
void utf(ll &l, ll &r) { tf(l, r), l >>= 1, r >>= 1; }
XNOR、NAND、NOR
直接用异或运算、与运算、或运算的方法求出来，然后将互反的两位交换即可。
Pell 方程
形如$x^2-Dy^2=1$（D 为任意正整数）的方程称为佩尔方程，必有最小正整数解$(x_0,y_0)$，用
$x_n=x_0x_{n-1}+Dy_0y_{n-1},y_n=y_0x_{n-1}+x_0y_{n-1}$
可递推方程的第 n 小整数解（可用矩阵快速幂求），同时还有
$2x_0x_n=x_{n-1}+x_{n+1},2x_0y_n=y_{n-1}+y_{n+1}$
Bertrand 猜想
$\forall n>3,\exist n<p<n\times 2$其中 n 为整数，p 为质数。
威尔逊定理
当且仅当 p 为素数时：$(p-1)!\equiv -1\pmod p$，
Jacobi’s Four Square Theorem
设$a^2+b^2+c^2+d^2=n$的自然数解个数为$r4(n)$，$d(n)$为$n$的约数和，由 Jacobi’s Four Square Theorem 可知，若$n$是奇数，则$r4(n)=8d(n)$，否则$r4(n)=24d(k)$，$k$是$n$去除所有$2$后的结果。
组合数学
组合数取模
为方便，记$C(n,m)=C_n^m=\binom{n}{m}$。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
struct Factorial : Mod
{
vector<ll> fac, ifac;
Factorial(int N, ll M) : fac(N, 1), ifac(N, 1), Mod(M)
{
for (int i = 2; i < N; ++i)
fac[i] = mul(fac[i - 1], i), ifac[i] = mul(M - M / i, ifac[M % i]);
for (int i = 2; i < N; ++i)
ifac[i] = mul(ifac[i], ifac[i - 1]);
}
ll c(int n, int m) { return mul(mul(fac[n], ifac[m]), ifac[n - m]); }
ll lucas(ll n, ll m) //卢卡斯定理求C(n,m)%M，适用于模数M小于N的情况，或者m较小的时候也可以暴力求
{
if (!m)
return 1;
if (n < m || n % M < m % M)
return 0;
if (n < M && m < M)
return c(n, m);
return mul(lucas(n / M, m / M), lucas(n % M, m % M));
}
};
组合数 LCM
$(n + 1)lcm(C(n,0),C(n,1),\dots,C(n,k))=lcm(n+1,n,n−1,n−2,\dots,n−k+1)$
区间 lcm 的维护：对于一个数，将其分解质因数，若有因子$p^k$，那么拆分出 k 个数 $p^1,p^2,\dots,p^k$，权值都为 p，那么区间$[l,r]$内所有数的 lcm 的答案=所有在该区间中出现过的数的权值之积，可持久化线段树维护之。
Stirling 数
第一类斯特林数
第一类斯特林数$S(p,k)$的一个的组合学解释是：将 p 个物体排成 k 个非空循环排列的方法数。
递推公式：$S(p,k)=(p−1)S(p−1,k)+S(p−1,k−1),1\leq k\leq p−1;S(p,0)=0,p\ge1;S(p,p)=1,p\ge0$
第二类斯特林数
第二类斯特林数$S(p,k)$的一个的组合学解释是：将 p 个物体划分成 k 个非空不可辨别的（可以理解为盒子没有编号）集合的方法数。
递推公式：$S(p,k)=kS(p−1,k)+S(p−1,k−1),1\leq k\leq p−1;S(p,0)=0,p\ge 1;S(p,p)=1,p\ge0$
卷积形式：$S(n,m)=\frac{1}{m!}\sum_{k=0}^m(-1)^kC(m,k)(m-k)^n=\sum_{k=0}^m\frac{(-1)^k}{k!}\frac{(m-k)^n}{(m-k)!}$
同时有转化：$x^k=\sum_{i=1}^ki!C(x,i)S(k,i)$
斯特林近似公式
$n!\approx\sqrt{2\pi n}(\frac{n}{e})^n$
小球入盒模型通解
| k 个球 | m 个盒子 | 空盒子 | 方案数 |
|---|---|---|---|
| 各不相同 | 各不相同 | 允许 | $m^k$ |
| 各不相同 | 各不相同 | 无 | $m!Stirling2(k,m)$ |
| 各不相同 | 完全相同 | 允许 | $\sum_{i=1}^mStirling2(k,i)$ |
| 各不相同 | 完全相同 | 无 | $Stirling2(k,m)$ |
| 完全相同 | 各不相同 | 允许 | $C(m+k−1,k)$ |
| 完全相同 | 各不相同 | 无 | $C(k−1,m−1)$ |
| 完全相同 | 完全相同 | 允许 | $\frac{1}{(1−x)(1−x^2)…(1−x^m)}的x^k项的系数$ |
| 完全相同 | 完全相同 | 无 | $\frac{x^m}{(1−x)(1−x^2)…(1−x^m)}的x^k项的系数$ |
置换
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
struct Permutation : vector<int>
{
Permutation(int n = 0) : vector<int>(n) {}
friend Permutation operator*(const Permutation &f, const Permutation &g)
{
Permutation ans(f.size());
for (int i = 0; i < f.size(); ++i)
ans[i] = g[f[i]];
return ans;
}
friend Permutation inv(const Permutation &f)
{
Permutation ans(f.size());
for (int i = 0; i < f.size(); ++i)
ans[f[i]] = i;
return ans;
}
friend vector<vector<int>> cycle(const Permutation &f)
{
vector<int> vis(f.size(), 0);
vector<vector<int>> ans;
for (int i = 0; i < f.size(); ++i)
if (!vis[i])
{
ans.push_back(vector<int>());
for (int j = i; !vis[j]; j = f[j])
vis[j] = 1, ans.back().push_back(j);
}
return ans;
}
};
生成字典序
下一排列
对给定的排列$a_1a_2\dots a_n$，找到$a_j$使得$a_j<a_{j+1},a_{j+1}>a_{j+2}>\dots>a_n$即这列数中最后一个相邻递增数对，然后把$a_{j+1},a_{j+2},\dots,a_n$中大于$a_j$的最小数放到位置 j，然后$a_j\dots a_n$中剩余的数从小到大排序放到$[j+1,n]$中。
1
2
3
4
5
6
7
8
9
10
11
bool nextPermutation(ll *b, ll *e) //标准库有这个函数next_permutation
{
ll *i = e - 1, *j = e - 2;
while (j >= b && *j >= *(j + 1))
--j;
if (j < b)
return 0;
while (*i <= *j)
--i;
return swap(*i, *j), reverse(j + 1, e), 1;
}
二项式反演
$f(n)=\sum_{k=0}^nC(n,k)g(k),g(n)=\sum_{k=0}^n(−1)^{n−k}C(n,k)f(k)$
第 k 小期望
$f(n,k)$表示有 n 个变量，和为 1，第 k 小的期望。
$f(n,k)=\frac{1}{n^2}+(1-\frac{1}{n})f(n-1,k-1),f(n,0)=0$
错排数
考虑一个有 n 个元素的排列，若一个排列中所有的元素都不在自己原来的位置上，那么这样的排列就称为原排 列的一个错排。
n 个元素的错排数$D_n$满足递推公式：$D_1=0,D_2=1,D_n=(n−1)(D_{n−2}+D_{n−1})$
通项：$D(n)=n![\frac{(-1)^2}{2!}+\dots+\frac{(-1)^{n-1} }{(n-1)!}+\frac{(-1)^n}{n!}]=\lfloor\frac{n!}{e}+\frac{1}{2}\rfloor$
Bonuli 数
$B_n = -\frac{1}{C(n+1,n)}(C(n+1,0)B_0+C(n+1,1)B_1+\dots+C(n+1,n-1)B_{n-1})=-\frac{1}{n+1}(C(n+1,0)B_0+C(n+1,1)B_1+\dots+C(n+1,n-1)B_{n-1})$
可用于计算任意正整数次数的幂和：$\sum_{i=1}^ni^k=\frac{1}{k+1}\sum_{j=0}^kC(k+1,j)B_jn^{k+1-j}$
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
struct Bonuli : Factorial
{
vector<ll> b;
Bonuli(int N, ll M) : Factorial(N, M), b(N, 0)
{
for (int i = b[0] = 1; i < N; ++i)
{
for (int j = 0; j < i; ++j)
b[i] = qadd(b[i], mul(b[j], c(i + 1, j)));
b[i] = qadd(M, -mul(b[i], mul(fac[i], ifac[i + 1])));
}
}
ll ask(ll n, int k)
{
ll r = 0, w = 1, u = add(n, 1);
for (int i = 1; i <= k + 1; ++i)
r = qadd(r, mul(mul(b[k + 1 - i], c(k + 1, i)), w = mul(w, u)));
return mul(r, mul(fac[k], ifac[k + 1]));
}
};
Catalan 数
$h_1=1,h_n=\frac{4n−2}{n+1}h_{n−1}=\frac{C(2n,n)}{n+1}=C(2n,n)−C(2n,n−1)$。
在一个格点阵列中，从$(0,0)$点走到$(n,m)$点且不经过对角线$x=y$的方法数：$C(n+m−1,m)−C(n+m−1,m−1),x>y;C(n+m,m)−C(n+m,m−1),x\ge y$。
常见的 Catalan 数：括号序的个数、凸多边形三角剖分的方案数等。
Bell 数
把 n 个带标号的物品划分为若干不相交集合的方案数称为贝尔数，其递推公式：$B_n=\sum_{i=0}^{N-1}C_{n-1}^iB_i$
前几项贝尔数：
1
1,2,5,15,52,203,877,4140,21147,115975,678570,4213597,27644437,190899322,1382958545,...
等价类容斥
考虑容斥，Bell(p)枚举所有等价情况。对于一种情况，强制了一个等价类里面的数都要相同，其它的可以相同也可以不同。
容斥系数为：$(−1)^{p−等价类个数}(每个等价类大小−1)!之积$。
Grey 码
格雷序列第 i 个是i^(i>>1)。长为 n 的 01 序列共$2^n$个，下标从$0\dots 2^n-1$。
扩展 Cayley 公式
对于 n 个点，m 个连通块的图，假设每个连通块有 a[i]个点，那么用 s−1 条边把它连通的方案数为$n^{s−2}a[1]a[2]\dots a[m]$。
超立方体
n 维超立方体有$2^{n−i}C(n,i)$个 i 维元素。
枚举位集 I 的非空子集 J
1
for(J=I; J; J=I&J−1) {}
线性代数
矩阵
1
typedef array<array<ll, N>, N> Mat;
乘法和快速幂
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
Mat operator*(const Mat &a, const Mat &b)
{
Mat r;
for (int i = 0; i < r.size(); ++i)
for (int j = 0; j < r.size(); ++j)
for (int k = r[i][j] = 0; k < r.size(); ++k)
M.qadd(r[i][j], M.mul(a[i][k], b[k][j]));
return r;
}
Mat pow(Mat a, ll b)
{
Mat r;
for (int i = 0; i < r.size(); ++i)
for (int j = 0; j < r[i].size(); ++j)
r[i][j] = i == j;
for (; b; b >>= 1, a = a * a)
if (b & 1)
r = r * a;
return r;
}
行列式
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
ll det(Mat a, int n)
{
ll ans = 1;
for (int i = 0; i < n; ++i)
{
for (int j = i + 1; j < n; ++j)
while (fabs(a[j][i]) > EPS)
{
ll t = a[i][i] / a[j][i];
for (int k = i; k < n; ++k)
a[i][k] -= t * a[j][k], swap(a[i][k], a[j][k]);
}
if (fabs(ans *= a[i][i]) < EPS)
return 0;
}
return ans;
}
高斯消元
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
void GaussElimination(Mat &a, int n) //a为增广矩阵，要求n*n的系数矩阵可逆，运行结束后a[i][n]为第i个未知数的值
{
for (int i = 0, r; i < n; ++i)
{
for (int j = r = i; j < n; ++j)
if (fabs(a[r][i]) < fabs(a[j][i]))
r = j;
if (r != i)
swap_ranges(a[r].begin(), a[r].begin() + n + 1, a[i]);
for (int j = n; j >= i; --j)
for (int k = i + 1; k < n; ++k)
a[k][j] -= a[k][i] * a[i][j] / a[i][i];
}
for (int i = n - 1; ~i; --i)
{
for (int j = i + 1; j < n; ++j)
a[i][n] -= a[j][n] * a[i][j];
a[i][n] /= a[i][i];
}
}
线性基
向量线性基
add 返回要插入的向量 z 是否与已插入的线性无关。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
struct Base
{
vector<vector<double>> v;
Base(int N) : v(N, vector<double>(N, 0)) {} //R^N的子空间
bool add(vector<double> x)
{
for (int i = 0; i < x.size(); ++i)
if (fabs(x[i]) > EPS)
{
if (fabs(v[i][i]) < EPS)
return v[i] = x, 1;
double t = x[i] / v[i][i];
for (int j = 0; j < x.size(); ++j)
x[j] -= t * v[i][j];
}
return 0;
}
};
异或线性基
若要查询第 k 小子集异或和，则把 k 写成二进制，对于是 1 的第 i 位，把从低位到高位第 i 个不为 0 的数异或进答案。若要判断是否有非空子集的异或和为 0，如果不存在自由基，那么说明只有空集的异或值为 0，需要高斯消元来判断。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
struct BaseXOR
{
vector<ll> a;
BaseXOR() : a(64, 0) {}
ll ask() //查询最大子集异或和
{
ll t = 0;
for (int i = a.size() - 1; ~i; --i)
t = max(t, t ^ a[i]);
return t;
}
bool add(ll x)
{
for (int i = a.size() - 1; ~i; --i)
if (x >> i & 1)
{
if (a[i])
x ^= a[i];
else
return a[i] = x, 1;
}
return 0;
}
bool check(ll x) //判断一个数是否能够被异或出，0根据需要特判
{
for (int i = a.size() - 1; ~i; --i)
if (x >> i & 1)
if (x ^= a[i], !x)
return 1;
return 0;
}
};
离散数学
归并排序求逆序对
1
2
3
4
5
6
7
8
9
10
ll merge_sort(ll *b, ll *e) //int存答案会爆
{
if (e - b < 2)
return 0;
ll *m = b + (e - b) / 2, ans = merge_sort(b, m) + merge_sort(m, e);
for (ll *i = b, *j = m; i < m && j < e; ++i)
for (; j < e && *j < *i; ++j)
ans += m - i;
return inplace_merge(b, m, e), ans;
}
约瑟夫问题
1
2
3
4
5
6
7
8
9
ll josephus(ll n, ll k) //编号0~n-1，每k个出列，时间复杂度O(min(n,k))
{
if (n < 3)
return k % n;
if (n < k)
return (Josephus(n - 1, k) + k) % n;
ll ret = Josephus(n - n / k, k) - n % k;
return ret < 0 ? ret + n : ret + ret / (k - 1);
}
堆的操作
和 STL 的函数push_heap、pop_heap用法完全相同，实现的是大根堆。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
ll pushHeap(ll *b, ll *e)
{
ll tmp = *--e;
for (int i = e - b, p;; b[i] = b[p], i = p)
if (i <= 0 || b[p = (i - 1) / 2] >= tmp)
return b[i] = tmp;
}
ll popHeap(ll *b, ll *e)
{
ll tmp = b[0];
for (int i = 0, siz = --e - b, ch;; b[i] = b[ch], i = ch)
{
if (ch = i * 2 + 1, ch + 1 < siz && b[ch] < b[ch + 1])
++ch;
if (i >= siz / 2 || b[siz] >= b[ch])
return b[i] = b[siz], b[siz] = tmp;
}
}
蔡勒公式
$w=(\lfloor\frac{c}{4}\rfloor-2c+y+\lfloor\frac{y}{4}\rfloor+\lfloor\frac{13(m+1)}{5}\rfloor+d-1)\mod7$
w：$0,1,\dots,6$对应周日，周一，$\dots$，周六
c：世纪减 1（即年份前两位数）。
y：年份后两位数。
m：月（$3\leq m\leq14$，即在蔡勒公式中，1、2 月要看作上一年的 13、14 月来计算）。
d：日。
曼哈顿距离的变换
$\mid x_1−x_2\mid +\mid y_1−y_2\mid=\max (\mid (x_1 + y_1)−(x_2 + y_2)\mid ,\mid (x_1 −y_1)−(x_2 −y_2)\mid )$
皮克定理
顶点坐标均是整点（或说正方形格点）的简单多边形中，面积 S 和内部格点数目 n、边上格点数目 m 的满足关系$S=n+\frac{m}{2}-1$。
矩形面积并、矩形面积交、矩形周长并（线段树、扫描线总结）
详见这篇博文。
数学分析
增长趋势
[n\to+\infty,\forall p,q>0,a>1,{(\ln n)}^q\ll n^p\ll a^n\ll n!\ll n^n]
积分表
反读可得导数表，此处略。
[\int k\,\mathrm{d}x=kx+C
\int x^a\,dx=\frac{x^{a+1}}{a+1}+C
\int\frac{1}{x}\,dx=\ln\mid x\mid +C
\int e^x\,dx=e^x + C
\int a^x\,dx=\frac{a^x}{\ln a}+C
\int\cos x\,dx=\sin x+C
\int\sin x\,dx=-\cos x+C
\int\frac{1}{cos^2x}\,dx=\int\sec^2 x\,dx=\tan x+C
\int\frac{1}{sin^2x}\,dx=\int\csc^2 x\,dx=-\cot x+C
\int\frac{1}{\sqrt{1-x^2}}\,dx=\arcsin x+C=-\arccos x+C
\int\frac{1}{1+x^2}\,dx=\arctan x+C=-arccot\,x+C
\int\sec x\tan x\,dx=\sec x+C
\int\csc x\cot x\,dx=-\csc x+C
\int\tan x\,dx=-\ln\mid \cos x\mid +C
\int\cot x\,dx=\ln\mid \sin x\mid +C
\int\sec x\,dx=\ln\mid \sec x+\tan x\mid +C
\int\csc x\,dx=\ln\mid \csc x-\cot x\mid +C
\int sh\,x\,dx=ch\,x+C
\int ch\,x\,dx=sh\,x+C
\int\frac{1}{x^2+a^2}\,dx=\frac{1}{a}\arctan\frac{x}{a}+C
\int\frac{1}{x^2-a^2}\,dx=\frac{1}{2a}\ln\mid \frac{x-a}{x+a}\mid +C
\int\frac{1}{\sqrt{a^2-x^2}}\,dx=\arcsin\frac{x}{a}+C
\int\frac{1}{\sqrt{x^2-a^2}}\,dx=\ln\mid x+\sqrt{x^2-a^2}\mid +C
\int\frac{1}{\sqrt{x^2+a^2}}\,dx=\ln\mid x+\sqrt{x^2+a^2}\mid +C\]
积分求几何量
弧长
若简单闭曲线
[\begin{cases}
x=x(t),
y=y(t),
\end{cases}
t\in[\alpha,\beta]]
端点处重合（$x(\alpha)=x(\beta),y(\alpha)=y(\beta)$）且其他地方不自交，$x(t),y(t)$连续且满足
[[x’(t)]^2+[y’(t)]^2\ne0,\forall t\in[\alpha,\beta]]
此时称曲线光滑，其长度
[s=\int_\alpha^\beta\sqrt{[x’(t)]^2+[y’(t)]^2}\,dt]
此式可对称推广到高维空间曲线。 极坐标下，
[r=r(\theta),\theta\in[\alpha,\beta]]
的长度为
[s=\int_\alpha^\beta\sqrt{[r(\theta)]^2+[r’(\theta)]^2}\,d\theta]
面积
若简单闭曲线
[\begin{cases}
x=x(t),
y=y(t),
\end{cases}
t\in[\alpha,\beta]]
端点处连续（$x(\alpha)=x(\beta),y(\alpha)=y(\beta)$）且其他地方不自交，$x(t),y(t)$都逐段有连续微商，则此闭合曲线围起来的有界区域面积
[S=-\int_\alpha^\beta x’(t)y(t)\,dt=-\int_\alpha^\beta y(t)\,dx(t)=-\oint_\Gamma y\,dx=\oint_\Gamma x\,dy]
等式右边称为曲线$\Gamma$上的积分，其计算方法是带入参数方程到定积分计算式中，积分上下限为始点与终点对应的参数值。下限并不总是小于上限，参数从下限到上限变化时对应曲线的正向（沿正向观察时，曲线所围的区域永远在左侧）。 极坐标下，连续非负曲线$r=r(\theta)$与向径$\theta=\alpha,\theta=\beta$，其中$0\leq\beta-\alpha\leq2\pi$所围成的平面图形面积
[S=\frac{1}{2}\int_\alpha^\beta r^2(\theta)\,d\theta]
体积
记立体过 x 点且垂直于 x 轴的截面面积为$S(x)$，则其体积
[V=\int_a^bS(x)\,dx]
连续曲线$y=f(x)\ge 0,x\in[a,b]$绕 x 轴旋转一周产生的旋转体体积
[V=\pi\int_a^by^2\,dx]
旋转体侧面积
若曲线由参数方程
[\begin{cases}
x=x(t),
y=y(t),
\end{cases}
t\in[\alpha,\beta]]
给出，则其绕 x 轴旋转体的侧面积
[s=2\pi\int_\alpha^\beta y(t)\sqrt{[x’(t)]^2+[y’(t)]^2}\,dt]
方向导数
设三元函数$u=f(x,y,z)$在点$P_0(x_0,y_0,z_0)$的某邻域内有定义，任意给定始于点$P_0$的射线$l$，$P(x,y,z)$为 l 上且含于定义域内的点。若极限
[\lim_{r(p,p_0)\to0^+}\frac{f(P)-f(P_0)}{r(P,P_0)}=\lim_{r(p,p_0)\to0^+}\frac{\Delta_lf(P_0)}{r(P,P_0)}]
存在，则称该极限值为函数$f$在点$P_0$沿方向$l$的方向导数，记为 $\frac{\partial f}{\partial l}\mid _{P_0}$或$\frac{\partial f(P_0)}{\partial l}$，$\frac{\Delta_lf(P_0)}{r(P,P_0)}$称为函数在$P_0$点沿$l$方向的增量。
特别地，$\frac{\partial f(P_0)}{\partial x}$就是函数在$P_0$点沿$x$轴正向的方向导数，$y,z$轴上的方向导数同理。若函数在$P_0$点可微，则其在$P_0$沿任何方向$l$的方向导数都存在，则有以下公式
[\frac{\partial f(P_0)}{\partial l}=(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y},\frac{\partial f}{\partial z})\mid _{P_0}\cdot\vec{l_0}]
其中$\vec{l_0}=(\cos\alpha,\cos\beta,cos\gamma)=\frac{1}{\rho}(\Delta x,\Delta y,\Delta z)$为$l$的方向余弦。
曲率
若曲线由参数方程
[\begin{cases}
x=x(t),
y=y(t),
\end{cases}
t\in[\alpha,\beta]]
给出且有二阶微商，则其在一点的曲率
[K=\frac{\mid y’‘x’-y’x’‘\mid }{[x’^2+y’^2]^{\frac{3}{2}}}]
若$y=f(x)$，则
[K=\frac{\mid y’‘\mid }{(1+y’^2)^\frac{3}{2}}]
同时记$\frac{1}{K}$为曲率半径。
空间曲线的切线与法平面
若已知曲线上一点$P(x_0,y_0,z_0)$处的切向量为$\tau(x_0,y_0,z_0)=(A,B,C)$则曲线在该点的切线方程为
[\frac{x-x_0}A=\frac{y-y_0}B=\frac{z-z_0}C]
法平面方程为
[A(x-x_0)+B(y-y_0)+C(z-z_0)=0]
当曲线由参数方程
[\begin{cases}
x=x(t),
y=y(t),
z=z(t),
\end{cases}
t\in[\alpha,\beta]]
给出时，曲线在 P 点的切向量为
[\tau=\pm(x’(t_0),y’(t_0),z’(t_0))]
更一般地，若曲线用两曲面的交线给出
[\begin{cases}
F(x,y,z)=0,
G(x,y,z)=0,
\end{cases}]
且在 P 点的某邻域能确定函数组$y=y(x),z=z(x)$满足$y_0=y(x_0),z_0=z(x_0)$，且$y’(x),z’(x)$存在，则曲线在 P 点的切向量
[\tau=\pm(\frac{\partial(F,G)}{\partial(y,z)},\frac{\partial(F,G)}{\partial(z,x)},\frac{\partial(F,G)}{\partial(x,y)})]
空间曲面的切平面与法线
若已知曲面上一点$P(x_0,y_0,z_0)$处的切平面的法向量为$\vec n=(A’,B’,C’)$则曲线在该点的法线方程为
[\frac{x-x_0}{A’}=\frac{y-y_0}{B’}=\frac{z-z_0}{C’}]
切平面方程为
[A’(x-x_0)+B’(y-y_0)+C’(z-z_0)=0]
当曲面方程为$\pi:F(x,y,z)=0$在曲面上任取一条过 P 的曲线，设其方程为
[\begin{cases}
x=x(t),
y=y(t),
z=z(t),
\end{cases}
t\in[\alpha,\beta]]
此时有$F(x(t),y(t),z(t))=0$令$t=t_0$两边对 t 求导，并写成向量的内积式，得
[(F_x,F_y,F_z)_P\cdot(x’(t_0),y’(t_0),z’(t_0))=0]
则曲线在 P 点的法向量为
[\vec{n}=\pm(F_x,F_y,F_z)_P]
若曲线由参数方程给出
[\begin{cases}
x=x(u,v),
y=y(u,v),
z=z(u,v),
\end{cases}]
则曲线在 P 点的法向量
[\vec{n}=\pm(\frac{\partial(y,z)}{\partial(u,v)},\frac{\partial(z,x)}{\partial(u,v)},\frac{\partial(x,y)}{\partial(u,v)})]
高阶导数与泰勒公式
用$f^{(n)}(x)$表示 f(x)的 n 阶导数，只要让余项<EPS即可计算指定函数到任意精确度，特别取 a=0 时称为麦克劳林公式。
[f(x)=f(a)+f^{(1)}(a)(x-a)+\frac{f^{(2)}(a)}{2!}(x-a)^2+\dots+\frac{f^{(n)}(a)}{n!}(x-a)^n+R_n(x)]
佩亚诺余项
[R_n(x)=o((x-a)^n)]
积分余项
[R_n(x)=\frac{1}{n!}\int_a^x(x-t)^nf^{(n+1)}(t)\,dt]
拉格朗日余项
[R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}(x-a)^{n+1},a<\xi<x]
柯西余项
[R_n(x)=\frac{(x-a)^{n+1}}{n!}(1-\theta)^nf^{(n+1)}(a+\theta(x-a)),0<\theta<1]
对数函数
[[\ln(1+x)]^{(n)}=(-1)^{n-1}(n-1)!(1+x)^{-n} \ln(1+x)=x-\frac{x^2}{2}+\frac{x^3}{3}-\frac{x^4}{4}+\dots+(-1)^{n-1}\frac{x^n}{n}+R_n(x)]
幂函数
[[(1+x)^a]^{(n)}=a(a-1)\dots(a-n+1)(1+x)^{a-n}
(1+x)^a=1+ax+\frac{a(a-1)}{2!}x^2+\dots+\frac{a(a-1)\dots(a-n+1)}{n!}x^n+R_n(x)]
三角函数
[(\sin x)^{(n)}=\sin(x+\frac{n\pi}{2})
\sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\dots+(-1)^{k-1}\frac{x^{2k-1}}{(2k-1)!}+R{2k}(x)
\R{2k}(x)=(-1)^k\frac{\cos\theta x}{(2k+1)!}x^{2k+1}
(\cos x)^{(n)}=\cos(x+\frac{n\pi}{2})
\cos x=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\dots+(-1)^{k-1}\frac{x^{2k-2}}{(2k-2)!}+R{2k-1}(x)
R{2k-1}(x)=(-1)^k\frac{\cos\theta x}{(2k)!}x^{2k}]
指数函数
[(e^x)^{(n)}=e^x
e^x=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\dots+\frac{x^n}{n!}+R_n(x)
R_n(x)=\frac{e^{\theta x}}{(n+1)!}x^{n+1},\xi=\theta x,0<\theta<1]
二元函数
设$f(x,y)$在$P_0(x_0,y_0)$的某邻域$O(P_0)$内有直到$n+1$阶连续偏导数，则对$O(P_0)$内$\forall(x_0+\Delta x,y_0+\Delta y),\exists\theta\in(0,1)$，使得
[f(x0+\Delta x,y_0+\Delta y)=\sum{k=0}^n\frac{1}{k!}(\frac{\partial}{\partial x}\Delta x+\frac{\partial}{\partial y}\Delta y)^kf(x_0,y_0)+R_n]
其中
[R_n=\frac{1}{(n+1)!}(\frac{\partial}{\partial x}\Delta x+\frac{\partial}{\partial y}\Delta y)^{n+1}f(x_0+\theta\Delta x,y_0+\theta\Delta y)]
级数部分和
幂级数
快速计算幂级数的部分和$\sum_{i=1}^ni^k\mod M$可借助伯努利数，详见组合数学模板。
[\sum_{i=1}^ni^1=\frac 1 2n(n+1)
\sum_{i=1}^ni^2=\frac 1 6n(n+1)(2n+1)
\sum_{i=1}^ni^3=\frac 1 4[n(n+1)]^2
\sum_{i=1}^ni^4=\frac 1{30}n(n+1)(2n+1)(3n^2+3n-1)
\sum_{i=1}^ni^5=\frac 1{12}[n(n+1)]^2(2n^2+2n-1)
\sum_{i=1}^ni^6=\frac 1{42}n(n+1)(2n+1)(3n^4+6n^3-3n+1)]
调和级数
[n\to\infty,\sum_{i=1}^n\frac 1 i\to\ln n+r,r\approx0.5772156649015328\dots]
二分求零点、三分求极值点
需要$f(x)$在区间$[l,r]$上单调/凹凸性唯一。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
lf bs(lf l, lf r, lf f(lf x))
{
if (r - l < EPS)
return l;
lf m = (l + r) / 2;
return sgn(f(l) * f(m)) < 0 ? bs(l, m, f) : ts(m, r, f);
}
lf ts(lf l, lf r, lf f(lf x))
{
if (r - l < EPS)
return l;
lf d = (r - l) / 3, lm = l + d, rm = r - d;
return f(lm) < f(rm) ? ts(l, rm, f) : ts(lm, r, f); //极小值
}
自适应辛普森求积分
这篇论文论证了加一个十五分之一的偏移收敛会比较快…
1
2
3
4
5
6
7
8
9
10
11
12
struct Simpson
{
lf simpson(lf a, lf b, lf f(lf x))
{
return (f(a) + 4 * f((a + b) / 2) + f(b)) * (b - a) / 6;
}
lf ask(lf a, lf b, lf f(lf x), lf e = EPS)
{
lf c = (a + b) / 2, L = simpson(a, c, f), R = simpson(c, b, f), delta = (L + R - simpson(a, b, f)) / 15;
return fabs(delta) < e ? L + R + delta : ask(a, c, f, e / 2) + ask(c, b, f, e / 2);
}
};
插值法
拉格朗日插值法：插值多项式和插值基函数的形式对称，容易编程。但是，增加节点时，需要重新计算每一个插值基函数。要在$\pmod p$意义下进行的话，那么 p 只能是质数。 牛顿插值法：当插值节点增加时，之前已计算的结果仍然能用，每增加一个节点，只要再增加一项即可，从而避免了重复性计算。如果要 mod 非质数的话，那么就要用牛顿插值法。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
typedef complex<lf> Coord;
#define X real()
#define Y imag()
struct Lagrange
{
lf ask(const vector<Coord> &p, lf x) //返回p确定的多项式函数在x处的值
{
lf ret = 0;
for (int i = 0; i < p.size(); ++i)
{
lf tmp = p[i].Y;
for (int j = 0; j < p.size(); ++j)
if (i != j)
tmp *= (x - p[j].X) / (p[i].X - p[j].X);
ret += tmp;
}
return ret;
}
vector<lf> ask(vector<Coord> p) //返回p确定的多项式系数向量
{
vector<lf> ret(p.size()), sum(p.size());
ret[0] = p[0].Y, sum[0] = 1;
for (int i = 1; i < p.size(); ++i)
{
for (int j = p.size() - 1; j >= i; --j)
p[j].Y = (p[j].Y - p[j - 1].Y) / (p[j].X - p[j - i].X);
for (int j = i; ~j; --j)
sum[j] = (j ? sum[j - 1] : 0) - sum[j] * p[i - 1].X,
ret[j] += sum[j] * p[i].Y;
}
return ret;
}
};
struct Newton
{
lf differenceQuotient(const vector<Coord> &p, int k) //计算差商
{
lf ret = 0;
for (int i = 0; i <= k; ++i)
{
lf tmp = p[i].Y;
for (int j = 0; j <= k; ++j)
if (i != j)
tmp /= p[i].X - p[j].X;
ret += tmp;
}
return ret;
}
lf ask(const vector<Coord> &p, lf x)
{
lf ret = p[0].Y;
for (int i = 1; i < p.size(); ++i)
{
lf tmp = differenceQuotient(p, i); //多次求，可O(n^3)预处理优化
for (int j = 0; j < i; ++j)
tmp *= x - p[j].X;
ret += tmp;
}
return ret;
}
};
[]
计算几何
二维
点和向量化为坐标 Coord 进行运算，使用 stl 中的 complex 实现。
复数相乘的几何意义为长度相乘，极角相加。
用直线上的一点 p 和方向向量 v 表示一条经过 p 的直线，直线上的所有点 q 满足 q=p+t*v，其中 t 是参数；当限制 t≥0 时，该参数方程表示射线；限制 0≤t≤1 时，该参数方程表示线段。
此外，如果已知线段端点 a1 和 a2，可以通过 Line(a1,a2-a1)来得到对应的参数形式。
Morley 定理：三角形每个内角的三等分线相交成等边三角形。
欧拉定理：平面图的点数 V、边数 E 和面数 F 满足 V+F-E=2。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
typedef double lf;
typedef complex<lf> Coord;
const lf EPS = 1e-9;
#define X real()
#define Y imag()
struct Line
{
Coord p, v;
Line(Coord p = Coord(), Coord v = Coord()) : p(p), v(v) {}
Coord point(lf t) const { return p + v * t; }
};
struct Circle
{
Coord c;
lf r;
Circle(Coord c = Coord(), lf r = 0) : c(c), r(r) {}
Coord point(lf t) const { return c + polar(r, t); } //t为参数，幅角
};
/*
Coord(lf x=0,lf y=0);//构造函数
lf real(Coord a);//a的实部（复平面的横坐标）,也可写作a.real()
lf imag(Coord a);//a的虚部（复平面的纵坐标）,也可写作a.imag()
lf abs(Coord a);//向量a的模长，或是点a到原点的距离
lf norm(Coord a);//abs的平方，比abs快，但是要注意浮点数精度溢出
lf arg(Coord a);//a的幅角，与atan2(a.real(),a.imag())等价
Coord polar(lf r,lf t);//极坐标生成方式，r为幅值，t为幅角
//运算符重载+、-、*、/（以及对应的赋值运算，但是赋值运算不能写在表达式中，详见参考地址）、<<、>>（输出括号形式的坐标）
*/
int sgn(lf d) { return (d > EPS) - (d < -EPS); }
bool operator!=(const Coord &A, const Coord &B) { return sgn(A.X - B.X) || sgn(A.Y - B.Y); } //不等运算符，涉及到浮点数比较要重写
bool operator==(const Coord &A, const Coord &B) { return !(A != B); }
bool cmpCoord(const Coord &A, const Coord &B) { return sgn(A.X - B.X) ? sgn(A.X - B.X) < 0 : sgn(A.Y - B.Y) < 0; } //复数没有小于运算，只能这样定义一个比较函数
bool cmpLine(const Line &A, const Line &B) { return sgn(arg(A.v) - arg(B.v)) < 0; } //按极角排序，求凸包中使用
lf Dot(const Coord &A, const Coord &B) { return A.X * B.X + A.Y * B.Y; }
lf Cross(const Coord &A, const Coord &B) { return A.X * B.Y - B.X * A.Y; }
lf Angle(const Coord &A, const Coord &B) { return acos(Dot(A, B) / abs(A) / abs(B)); }
lf Area2(const Coord &A, const Coord &B, const Coord &C) { return Cross(B - A, C - A); } //三角形ABC有向面积的两倍
Coord Rotate(const Coord &A, lf rad) { return A * polar(1.0, rad);} //向量A逆时针旋转rad弧度
Coord Normal(const Coord &A) //A的法向量，把A逆时针旋转九十度并长度化为1
{
lf L = abs(A);
return Coord(-A.Y / L, A.X / L);
}
bool onLeft(const Coord &P, const Line &L) { return sgn(Cross(L.v, P - L.p)) > 0; } //p是否在有向直线L左侧，不含线上
lf DistanceToLine(const Coord &P, const Line &L) { return Cross(L.v, P - L.p) / abs(L.v); } //点到直线距离（有向）
lf DistanceToLine(const Coord &P, const Coord &A, const Coord &B) { return DistanceToLine(P, Line(A, B - A)); }
lf DistanceToSegment(const Coord &P, const Coord &A, const Coord &B) //点到线段的距离（无向）
{
if (A == B)
return abs(P - A);
Coord v1 = B - A, v2 = P - A, v3 = P - B;
if (sgn(Dot(v1, v2)) < 0)
return abs(v2);
if (sgn(Dot(v1, v3)) > 0)
return abs(v3);
return fabs(DistanceToLine(P, Line(A, B - A)));
}
Coord getLineProjection(const Coord &P, const Line &L) { return L.point(Dot(L.v, P - L.p) / norm(L.v)); } //点在直线上的投影
Coord getLineProjection(const Coord &P, const Coord &A, const Coord &B) { return getLineProjection(P, Line(A, B - A)); }
Coord getSymmetry(const Coord &P, const Coord &O) { return O + O - P; } //P关于O的对称点
Coord getSymmetry(const Coord &P, const Line &L) { return getSymmetry(P, getLineProjection(P, L)); } //P关于L的对称点
Coord getLineIntersection(const Line &L1, const Line &L2) { return L1.point(Cross(L2.v, L1.p - L2.p) / Cross(L1.v, L2.v)); } //直线交点,须确保两直线相交
Coord getLineIntersection(const Coord &A1, const Coord &A2, const Coord &B1, const Coord &B2) { return getLineIntersection(Line(A1, A2 - A1), Line(B1, B2 - B1)); }
bool SegmentProperIntersection(const Coord &A1, const Coord &A2, const Coord &B1, const Coord &B2) //线段相交判定，交点不在一条线段的端点
{
lf C1 = Cross(A2 - A1, B1 - A1), C2 = Cross(A2 - A1, B2 - A1),
C3 = Cross(B2 - B1, A1 - B1), C4 = Cross(B2 - B1, A2 - B1);
return sgn(C1) * sgn(C2) < 0 && sgn(C3) * sgn(C4) < 0;
}
bool onSegment(const Coord &P, const Coord &A1, const Coord &A2) { return sgn(Dot(A1 - P, A2 - P)) < 0 && !sgn(Cross(A1 - P, A2 - P)); } //判断点是否在线段上，不包含端点
lf PolygonArea(const vector<Coord> &p) //计算多边形的有向面积，凸多边形即为面积
{
lf s = 0;
for (int i = 2; i < p.size(); ++i)
s += Area2(p[0], p[i - 1], p[i]);
return s / 2;
}
int inPolygon(const Coord &p, const vector<Coord> &poly) //点在多边形内的判定，转角法，正值为内部，0为外部，-1在边界上
{
int ans = 0;
for (int i = 0, k, d1, d2, n = poly.size(); i != n; ++i)
{
if (onSegment(p, poly[i], poly[(i + 1) % n]))
return -1; //在边界上
k = sgn(Cross(poly[(i + 1) % n] - poly[i], p - poly[i]));
d1 = sgn(poly[i].Y - p.Y);
d2 = sgn(poly[(i + 1) % n].Y - p.Y);
if (k > 0 && d1 <= 0 && d2 > 0)
++ans;
if (k < 0 && d2 <= 0 && d1 > 0)
--ans;
}
return ans;
}
vector<Coord> ConvexHull(vector<Coord> p, int collineation = 1) //获得凸包，不希望凸包的边上有输入点第二个参数传0
{
vector<Coord> ans;
sort(p.begin(), p.end(), cmpCoord); //先比横坐标再比纵坐标
for (int i = 0; i < p.size(); ++i) //求出下凸包
{
while (ans.size() > 1 && sgn(Area2(ans[ans.size() - 2], ans[ans.size() - 1], p[i])) < collineation)
ans.pop_back();
ans.push_back(p[i]);
}
for (int i = p.size() - 2, k = ans.size(); i >= 0; --i) //求出上凸包
{
while (ans.size() > k && -sgn(Area2(ans[ans.size() - 1], ans[ans.size() - 2], p[i])) < collineation)
ans.pop_back();
ans.push_back(p[i]);
}
if (p.size() > 1)
ans.pop_back();
return ans;
}
vector<Coord> cutPolygon(const vector<Coord> &poly, const Coord &A, const Coord &B) //用有向直线A->B切割多边形poly， 返回「左侧」。 如果退化，可能会返回一个单点或者线段，复杂度O(n^2)
{
vector<Coord> newpoly;
for (int i = 0, n = poly.size(); i != n; ++i)
{
Coord C = poly[i], D = poly[(i + 1) % n];
if (sgn(Cross(B - A, C - A)) >= 0)
newpoly.push_back(C);
if (!sgn(Cross(B - A, C - D)))
{
Coord ip = getLineIntersection(Line(A, B - A), Line(C, D - C));
if (onSegment(ip, C, D))
newpoly.push_back(ip);
}
}
return newpoly;
}
vector<Coord> getHalfPlaneIntersection(vector<Line> L) //半平面交
{
sort(L.begin(), L.end(), cmpLine); //按极角排序
vector<Coord> p(L.size(), Coord()); //p[i]为q[i]和q[i+1]的交点
int first = 0, last = 0; //双端队列的第一个元素和最后一个元素
vector<Line> q(L.size(), Line()); //双端队列
q[0] = L[0]; //队列初始化为只有一个半平面L[0]
for (int i = 0, n = L.size(); i != n; ++i)
{
while (first < last && !onLeft(p[last - 1], L[i]))
--last;
while (first < last && !onLeft(p[first], L[i]))
++first;
q[++last] = L[i];
if (!sgn(Cross(q[last].v, q[last - 1].v)))
{
--last;
if (onLeft(L[i].p, q[last]))
q[last] = L[i];
}
if (first < last)
p[last - 1] = getLineIntersection(q[last - 1], q[last]);
}
while (first < last && !onLeft(p[last - 1], q[first]))
--last; //删除无用平面
if (last - first <= 1)
return vector<Coord>(); //空集
p[last] = getLineIntersection(q[last], q[first]);
return vector<Coord>(p.begin() + first, p.begin() + last + 1); //从deque复制到输出中
}
int getLineCircleIntersection(const Line &L, const Circle &C, vector<Coord> &sol)
{
lf a = L.v.X,
b = L.p.X - C.c.X,
c = L.v.Y,
d = L.p.Y - C.c.Y,
e = a * a + c * c,
f = 2 * (a * b + c * d),
g = b * b + d * d - C.r * C.r,
delta = f * f - 4 * e * g;
if (sgn(delta) < 0)
return 0;
if (!sgn(delta))
return sol.push_back(L.point(-f / (2 * e))), 1;
sol.push_back(L.point((-f - sqrt(delta)) / (2 * e)));
sol.push_back(L.point((-f + sqrt(delta)) / (2 * e)));
return 2;
}
int getCircleIntersection(const Circle &C1, const Circle &C2, vector<Coord> &sol)
{
lf d = abs(C1.c - C2.c);
if (!sgn(d))
return sgn(C1.r - C2.r) ? 0 : -1; //重合返回-1
if (sgn(C1.r + C2.r - d) < 0 || sgn(fabs(C1.r - C2.r) - d) > 0) //外离或内含
return 0;
lf a = arg(C2.c - C1.c),
da = acos((C1.r * C1.r + d * d - C2.r * C2.r) / (2 * C1.r * d));
Coord p1 = C1.point(a - da), p2 = C1.point(a + da);
sol.push_back(p1);
if (p1 == p2)
return 1; //相切
return sol.push_back(p2), 2;
}
Line getTangent(const Coord &C, const Coord &P) { return Line(P, Normal(C - P)); } //圆心C，圆上一点P处切线
int getTangents(const Coord &p, const Circle &C, vector<Coord> &sol) //点到圆的切点，返回个数
{
Coord u = p - C.c;
lf d = abs(u);
if (d < C.r)
return 0; //点在圆内
if (!sgn(d - C.r)) //点在圆上
return sol.push_back(p), 1;
lf base = arg(u), ang = acos(C.r / d);
sol.push_back(C.point(base + ang));
sol.push_back(C.point(base - ang));
return 2;
}
int getTangents(Circle A, Circle &B, vector<Coord> &a, vector<Coord> &b) //公共切线的切点
{
int cnt = 0;
if (A.r < B.r)
swap(A, B), swap(a, b); //有时需标记交换
lf d = abs(A.c - B.c),
rdiff = A.r - B.r,
rsum = A.r + B.r;
if (sgn(d - rdiff) < 0)
return 0; //内含
lf base = arg(B.c - A.c);
if (!sgn(d) && !sgn(rdiff))
return -1; //重合，无穷多条切线
if (!sgn(d - rdiff)) //内切，外公切线
{
a.push_back(A.point(base));
b.push_back(B.point(base));
return 1;
}
//有外公切线的情形
lf ang = acos(rdiff / d);
a.push_back(A.point(base + ang));
b.push_back(B.point(base + ang));
a.push_back(A.point(base - ang));
b.push_back(B.point(base - ang));
cnt += 2;
if (!sgn(d - rsum))
{
a.push_back(A.point(base));
b.push_back(B.point(base + M_PI));
++cnt;
}
else if (sgn(d - rsum) > 0)
{
lf ang_in = acos(rsum / d);
a.push_back(A.point(base + ang_in));
b.push_back(B.point(base + ang_in + M_PI));
a.push_back(A.point(base - ang_in));
b.push_back(B.point(base - ang_in + M_PI));
cnt += 2;
}
return cnt;
}
lf AreaCircleWithTriangle(const Circle &C, Coord A, Coord B) //C和三角形OAB的相交面积，如果三角形顶点不在O上则把圆和三角形同时平移，直到有一个顶点在O上
{
int sg = sgn(Cross(A, B));
if (!sg || A == C.c || B == C.c)
return 0;
lf OA = abs(A - C.c), OB = abs(B - C.c), angle = Angle(A, B),
d = DistanceToLine(Coord(), A, B);
if (sgn(OA - C.r) <= 0 && sgn(OB - C.r) <= 0)
return Cross(A, B) / 2;
if (sgn(OA - C.r) >= 0 && sgn(OB - C.r) >= 0 && sgn(d - C.r) >= 0)
return sg * C.r * C.r * angle / 2;
if (sgn(OA - C.r) >= 0 && sgn(OB - C.r) >= 0 && sgn(d - C.r) < 0)
{
Coord prj = getLineProjection(Coord(), A, B);
if (!onSegment(prj, A, B))
return sg * C.r * C.r * angle / 2;
vector<Coord> p;
Line L = Line(A, B - A);
getLineCircleIntersection(L, C, p);
lf s1 = C.r * C.r * angle / 2,
s2 = C.r * C.r * Angle(p[0], p[1]) / 2;
s2 -= fabs(Cross(p[0], p[1]) / 2);
s1 = s1 - s2;
return sg * s1;
}
if (sgn(OB - C.r) < 0)
swap(A, B);
Line L = Line(A, B - A);
vector<Coord> inter;
getLineCircleIntersection(L, C, inter);
Coord inter_point = inter[!onSegment(inter[0], A, B)];
lf s = fabs(Cross(inter_point, A) / 2);
s += C.r * C.r * Angle(inter_point, B) / 2;
return s * sg;
}
lf AreaCircleWithPolygon(const Circle &C, const vector<Coord> &p)
{
lf ans = 0;
for (int i = 0; i < p.size(); ++i)
ans += AreaCircleWithTriangle(C, p[i], p[(i + 1) % p.size()]);
return fabs(ans);
}
Coord getGravityCenter(const vector<Coord> &p) //多边形重心
{
Coord a(0, 0);
lf am = 0, mj;
for (int i = 0; i < p.size(); ++i)
{
mj = Cross(p[i], p[(i + 1) % p.size()]);
a += mj * (p[i] + p[(i + 1) % p.size()]);
am += mj;
}
return a / am / 3.0;
}
三维
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
typedef double lf;
const lf EPS = 1e-9, INF = 1e9;
int sgn(lf d) { return (d > EPS) - (d < -EPS); }
struct Coord3
{
lf X, Y, Z;
friend bool operator!=(const Coord3 &a, const Coord3 &b) { return sgn(a.X - b.X) || sgn(a.Y - b.Y) || sgn(a.Z - b.Z); }
friend bool operator==(const Coord3 &a, const Coord3 &b) { return !(a != b); }
Coord3 &operator+=(const Coord3 &b) { return X += b.X, Y += b.Y, Z += b.Z, *this; }
friend Coord3 operator+(Coord3 a, const Coord3 &b) { return a += b; }
Coord3 &operator-=(const Coord3 &b) { return X -= b.X, Y -= b.Y, Z -= b.Z, *this; }
friend Coord3 operator-(Coord3 a, const Coord3 &b) { return a -= b; }
Coord3 &operator*=(lf d) { return X *= d, Y *= d, Z *= d, *this; }
friend Coord3 operator*(Coord3 a, lf d) { return a *= d; }
friend Coord3 operator*(lf d, Coord3 a) { return a *= d; }
Coord3 &operator/=(lf d) { return X /= d, Y /= d, Z /= d, *this; }
friend Coord3 operator/(Coord3 a, lf d) { return a /= d; }
friend lf Dot(const Coord3 &A, const Coord3 &B) { return A.X * B.X + A.Y * B.Y + A.Z * B.Z; }
friend Coord3 Cross(const Coord3 &A, const Coord3 &B) { return {A.Y * B.Z - A.Z * B.Y, A.Z * B.X - A.X * B.Z, A.X * B.Y - A.Y * B.X}; }
friend lf norm(const Coord3 &A) { return Dot(A, A); }
friend lf abs(const Coord3 &A) { return sqrt(norm(A)); }
friend lf Angle(const Coord3 &A, const Coord3 &B) { return acos(Dot(A, B) / abs(A) / abs(B)); }
friend lf Area2(Coord3 A, Coord3 B, Coord3 C) { return abs(Cross(B - A, C - A)); }
friend lf Volume6(Coord3 A, Coord3 B, Coord3 C, Coord3 D) { return Dot(D - A, Cross(B - A, C - A)); } //四面体体积
friend Coord3 Centroid(Coord3 A, Coord3 B, Coord3 C, Coord3 D) { return (A + B + C + D) / 4.0; } //四面体的重心
friend lf DistanceToPlane(Coord3 p, Coord3 p0, const Coord3 &n) { return Dot(p - p0, n) / abs(n); } //点p到平面p0-n的有向距离
friend Coord3 getPlaneProjection(Coord3 p, Coord3 p0, const Coord3 &n) { return p - n * Dot(p - p0, n); } //点p在平面p0-n上的投影。n必须为单位向量
friend Coord3 LinePlaneIntersection(Coord3 p1, Coord3 p2, Coord3 p0, Coord3 n) //直线p1-p2 与平面p0-n的交点，假设交点唯一存在
{
Coord3 v = p2 - p1;
lf t = Dot(n, p0 - p1) / Dot(n, p2 - p1); //分母为0，直线与平面平行或在平面上
return p1 + v * t; //如果是线段 判断t是否在0~1之间
}
friend lf DistanceToLine(Coord3 P, Coord3 A, Coord3 B) //点P到直线AB的距离
{
Coord3 v1 = B - A, v2 = P - A;
return abs(Cross(v1, v2)) / abs(v1);
}
friend lf DistanceToSeg(Coord3 P, Coord3 A, Coord3 B) //点到线段的距离
{
if (A == B)
return abs(P - A);
Coord3 v1 = B - A, v2 = P - A, v3 = P - B;
if (sgn(Dot(v1, v2)) < 0)
return abs(v2);
if (sgn(Dot(v1, v3)) > 0)
return abs(v3);
return fabs(DistanceToLine(P, A, B));
}
friend bool LineDistance3D(Coord3 p1, Coord3 u, Coord3 p2, Coord3 v, lf &s) //求异面直线 p1+s*u与p2+t*v的公垂线对应的s，如果平行|重合，返回0
{
lf b = Dot(u, u) * Dot(v, v) - Dot(u, v) * Dot(u, v);
if (!sgn(b))
return 0;
lf a = Dot(u, v) * Dot(v, p1 - p2) - Dot(v, v) * Dot(u, p1 - p2);
return s = a / b, 1;
}
friend bool SameSide(Coord3 p1, Coord3 p2, Coord3 a, Coord3 b) { return sgn(Dot(Cross(b - a, p1 - a), Cross(b - a, p2 - a))) >= 0; } //p1和p2是否在线段a-b的同侧
friend bool PointInTri(Coord3 PP, Coord3 P[3]) //点P在三角形P0,P1,p中
{
return SameSide(PP, P[0], P[1], P[2]) &&
SameSide(PP, P[1], P[0], P[2]) &&
SameSide(PP, P[2], P[0], P[1]);
}
friend bool TriSegIntersection(Coord3 P[3], Coord3 A, Coord3 B, Coord3 &PP) //三角形P0P1p是否和线段AB相交，如有则为PP
{
Coord3 n = Cross(P[1] - P[0], P[2] - P[0]);
if (sgn(Dot(n, B - A)) == 0)
return false; //线段A-B和平面P0P1p平行或共面
lf t = Dot(n, P[0] - A) / Dot(n, B - A); //平面A和直线P1-p有惟一交点
if (sgn(t) < 0 || sgn(t - 1) > 0)
return false; //不在线段AB上
return PointInTri(PP = A + (B - A) * t, P);
}
friend bool TriTriIntersection(Coord3 T1[3], Coord3 T2[3]) //空间两三角形是否相交
{
Coord3 P;
for (int i = 0; i < 3; ++i)
if (TriSegIntersection(T1, T2[i], T2[(i + 1) % 3], P) ||
TriSegIntersection(T2, T1[i], T1[(i + 1) % 3], P))
return 1;
return 0;
}
friend lf SegSegDistance(Coord3 a1, Coord3 b1, Coord3 a2, Coord3 b2, Coord3 &ans1, Coord3 &ans2) //空间两直线上最近点对 返回最近距离 点对保存在ans1 ans2中
{
Coord3 v1 = (a1 - b1), v2 = (a2 - b2);
Coord3 N = Cross(v1, v2);
Coord3 ab = (a1 - a2);
lf ans = Dot(N, ab) / abs(N);
Coord3 d1 = b1 - a1, d2 = b2 - a2, cd = Cross(d1, d2);
lf nd = norm(cd), t1 = Dot(Cross(a2 - a1, d2), cd) / nd, t2 = Dot(Cross(a2 - a1, d1), cd) / nd;
return ans1 = a1 + (b1 - a1) * t1, ans2 = a2 + (b2 - a2) * t2, fabs(ans);
}
friend bool InsideWithMinDistance(Coord3 PP, Coord3 *P, lf dist) //判断PP是否在三角形P中，并且到三条边的距离都至少为dist。保证P,A,B,C共面
{
return PointInTri(PP, P) &&
DistanceToLine(PP, P[0], P[1]) >= dist ||
DistanceToLine(PP, P[1], P[2]) >= dist ||
DistanceToLine(PP, P[2], P[0]) >= dist;
}
friend lf minBall(const vector<Coord3> &data, const lf eps = EPS * 1e-3) //模拟退火求最小球覆盖，EPS玄学调整；返回半径
{
lf step = 1, ans = INF;
for (Coord3 z{0, 0, 0}; step > eps; step *= 0.99)
{
int s = 0;
for (int i = 0; i < data.size(); ++i)
if (abs(z - data[s]) < abs(z - data[i]))
s = i;
ans = min(ans, abs(z - data[s]));
z -= (z - data[s]) * step;
}
return ans;
}
};
struct Sphere
{
Coord3 o;
lf r;
Coord3 point(lf lat, lf lng) const //经纬度确定球面上一点
{
lat *= M_PI / 180;
lng *= M_PI / 180;
return o + r * Coord3{cos(lat) * cos(lng), cos(lat) * sin(lng), sin(lat)};
}
};
struct ConvexPolyhedron //空间多边形和凸包问题
{
struct Face
{
int v[3];
Face(int a, int b, int c) { v[0] = a, v[1] = b, v[2] = c; }
Coord3 Normal(const vector<Coord3> &P) const { return Cross(P[v[1]] - P[v[0]], P[v[2]] - P[v[0]]); }
bool CanSee(const vector<Coord3> &P, int i) const { return Dot(P[i] - P[v[0]], Normal(P)) > 0; } //f是否能看见P[i]
};
vector<Face> faces;
vector<Coord3> p;
ConvexPolyhedron(vector<Coord3> P) : p(P)
{
for (int i = 0; i < p.size(); ++i)
P[i] += Coord3{randEPS(), randEPS(), randEPS()};
vector<vector<int>> vis(P.size(), vector<int>(P.size()));
faces.push_back(Face(0, 1, 2)); //由于已经进行扰动，前三个点不共线
faces.push_back(Face(2, 1, 0));
for (int i = 3; i < P.size(); ++i)
{
vector<Face> next;
for (int j = 0; j < faces.size(); ++j) //计算每条边的「左面」的可见性
{
Face &f = faces[j];
int res = f.CanSee(P, i);
if (!res)
next.push_back(f);
for (int k = 0; k < 3; ++k)
vis[f.v[k]][f.v[(k + 1) % 3]] = res;
}
for (int j = 0; j < faces.size(); ++j)
for (int k = 0; k < 3; ++k)
{
int a = faces[j].v[k], b = faces[j].v[(k + 1) % 3];
if (vis[a][b] != vis[b][a] && vis[a][b]) //(a,b)是分界线，左边对P[i]可见
next.push_back(Face(a, b, i));
}
swap(faces, next);
}
}
lf randEPS() { return (rand() / lf(RAND_MAX) - 0.5) * EPS; }
Coord3 centroid() //三维凸包重心
{
Coord3 C = p[0], tot{0, 0, 0};
lf totv = 0;
for (int i = 0; i < faces.size(); ++i)
{
Coord3 p1 = p[faces[i].v[0]], p2 = p[faces[i].v[1]], p3 = p[faces[i].v[2]];
lf v = -Volume6(p1, p2, p3, C);
totv += v;
tot += v * Centroid(p1, p2, p3, C);
}
return tot / totv;
}
lf dist(Coord3 C) //凸包内一点到表面最近距离
{
lf ans = INF;
for (int i = 0; i < faces.size(); ++i)
{
Coord3 p1 = p[faces[i].v[0]], p2 = p[faces[i].v[1]], p3 = p[faces[i].v[2]];
ans = min(ans, fabs(-Volume6(p1, p2, p3, C) / Area2(p1, p2, p3)));
}
return ans;
}
};
高精度
无符号整数
vector 自带大小比较为字典序比较， != 、 == 运算，其余需要时一定记得重载！
减法，当被减数小于减数时减为 0。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
struct Wint : vector<int> //继承vector
{
static const int width = 9, base = 1e9;
Wint(unsigned long long n = 0) //普通初始化，当整型数和Wint同时运算时会提升至Wint
{
for (; n; n /= base)
push_back(n % base);
}
explicit Wint(const string &s) //字符串初始化函数，未判断字符串合法情况
{
for (int len = int(s.size() - 1) / width + 1, b, e, i = 0; i < len; ++i)
for (e = s.size() - i * width, b = max(0, e - width), push_back(0); b != e; ++b)
back() = back() * 10 + s[b] - '0';
trim(0);
}
Wint &trim(bool up = 1) //去前导0，是否需要进位，很常用的小函数，为方便返回自身
{
for (int i = 1; up && i < size(); ++i)
{
if ((*this)[i - 1] < 0)
--(*this)[i], (*this)[i - 1] += base;
if ((*this)[i - 1] >= base)
(*this)[i] += (*this)[i - 1] / base, (*this)[i - 1] %= base;
}
while (!empty() && back() <= 0)
pop_back();
for (; up && !empty() && back() >= base; (*this)[size() - 2] %= base)
push_back(back() / base);
return *this;
}
friend istream &operator>>(istream &is, Wint &n)
{
string s; //懒
return is >> s, n = Wint(s), is;
}
friend ostream &operator<<(ostream &os, const Wint &n)
{
if (n.empty())
return os.put('0');
os << n.back();
char ch = os.fill('0');
for (int i = n.size() - 2; ~i; --i)
os.width(n.width), os << n[i];
return os.fill(ch), os;
}
friend bool operator<(const Wint &a, const Wint &b)
{
if (a.size() != b.size())
return a.size() < b.size();
for (int i = a.size() - 1; ~i; --i)
if (a[i] != b[i])
return a[i] < b[i];
return 0;
}
friend bool operator>(const Wint &a, const Wint &b) { return b < a; }
friend bool operator<=(const Wint &a, const Wint &b) { return !(a > b); }
friend bool operator>=(const Wint &a, const Wint &b) { return !(a < b); }
Wint &operator+=(const Wint &b)
{
if (size() < b.size())
resize(b.size()); //保证有足够的位数
for (int i = 0; i < b.size(); ++i)
(*this)[i] += b[i];
return trim(); //单独进位防自运算
}
friend Wint operator+(Wint a, const Wint &b) { return a += b; }
Wint &operator++() { return *this += 1; } //前置版本，懒
Wint operator++(int) //后置版本
{
Wint b(*this);
return ++*this, b;
}
Wint &operator-=(const Wint &b) //a<b会使a变为0
{
if (size() < b.size())
resize(b.size()); //保证有足够的位数
for (int i = 0; i < b.size(); ++i)
(*this)[i] -= b[i];
return trim(); //单独进位防自运算
}
friend Wint operator-(Wint a, const Wint &b) { return a -= b; }
Wint &operator--() { return *this -= 1; } //前置版本，懒
Wint operator--(int) //后置版本
{
Wint b(*this);
return --*this, b;
}
Wint &operator*=(const Wint &b) //高精度乘法，常规写法
{
Wint c;
c.assign(size() + b.size(), 0);
for (int j = 0, k, l; j < b.size(); ++j)
if (b[j]) //稀疏优化，特殊情况很有效
for (int i = 0; i < size(); ++i)
{
unsigned long long n = (*this)[i];
for (n *= b[j], k = i + j; n; n /= base)
c[k++] += n % base;
for (l = i + j; c[l] >= base || l + 1 < k; c[l++] %= base)
c[l + 1] += c[l] / base;
}
return swap(c), trim(0);
}
/*
Wint &operator*=(const Wint &b) //一种效率略高但对位宽有限制的写法
{
vector<unsigned long long> n(size() + b.size(), 0); //防爆int
//乘法算完后统一进位效率高，防止乘法溢出（unsigned long long范围0~1.8e19）
//位宽为9时size()不能超过18（十进制162位），位宽为8时size()不能超过1800（十进制14400位）等等。
for (int j = 0; j != b.size(); ++j)
if (b[j]) //稀疏优化，特殊情况很有效
for (int i = 0; i != size(); ++i)
n[i + j] += (unsigned long long)(*this)[i] * b[j];
for (int i = 1; i < n.size(); ++i) //这里用<防止位数0，单独进位防自运算
n[i] += n[i - 1] / base, n[i - 1] %= base;
return assign(n.begin(), n.end()), trim(0);
}
Wint &operator*=(const Wint &b) //fft优化乘法，注意double仅15位有效数字，调小Wint::width不超过2，计算自2*log2(base)+2*log2(len)<53
{
vector<ll> ax(begin(), end()), bx(b.begin(), b.end());
ax = FFT(size() + b.size()).ask(ax, bx);
for (int i = 1; i < ax.size(); ++i)
ax[i] += ax[i - 1] / base, ax[i - 1] %= base;
return assign(ax.begin(), ax.end()), trim(0);
}
Wint &operator*=(const Wint &b) //ntt优化，Wint::width不超过2
{
vector<ll> ax(begin(), end()), bx(b.begin(), b.end());
ax = FNTT(size() + b.size(), (7 << 26) + 1, 3).ask(ax, bx);
for (int i = 1; i < ax.size(); ++i)
ax[i] += ax[i - 1] / base, ax[i - 1] %= base;
return assign(ax.begin(), ax.end()), trim(0);
}
*/
friend Wint operator*(Wint a, const Wint &b) { return a *= b; }
Wint &operator/=(Wint b)
{
Wint r, c, d = b.base / (b.back() + 1);
*this *= d, b *= d, c.assign(size(), 0);
for (int i = size() - 1; ~i; --i)
{
r.insert(r.begin(), (*this)[i]);
unsigned long long s = 0;
for (int j = b.size(); j + 1 >= b.size(); --j) //b.size()==0肯定第一行就出问题的
s = s * b.base + (j < r.size() ? r[j] : 0);
for (d = c[i] = s / b.back(), d *= b; r < d; r += b)
--c[i];
r -= d;
}
return swap(c), trim(0); //r为加倍后的余数，可通过高精度除低精度得到真正余数，此处略
}
friend Wint operator/(Wint a, const Wint &b) { return a /= b; }
Wint &operator%=(const Wint &b) { return *this -= *this / b * b; }
friend Wint operator%(Wint a, const Wint &b) { return a %= b; }
//开平方，改自ZJU模板
bool cmp(long long c, int d, const Wint &b) const
{
if ((int)b.size() - (int)size() < d + 1 && c)
return 1;
long long t = 0;
for (int i = b.size() - 1, lo = -(base << 1); lo <= t && t <= 0 && ~i; --i)
if (t = t * base - b[i], 0 <= i - d - 1 && i - d - 1 < size())
t += (*this)[i - d - 1] * c;
return t > 0;
}
Wint &sub(const Wint &b, long long k, int d)
{
int l = b.size() + d;
for (int i = d + 1; i <= l; ++i)
{
long long tmp = (*this)[i] - k * b[i - d - 1];
if (tmp < 0)
{
(*this)[i + 1] += (tmp - base + 1) / base;
(*this)[i] = tmp - (tmp - base + 1) / base * base;
}
else
(*this)[i] = tmp;
}
for (int i = l + 1; i < size() && (*this)[i] < 0; ++i)
{
(*this)[i + 1] += ((*this)[i] - base + 1) / base;
(*this)[i] -= ((*this)[i] - base + 1) / base * base;
}
return trim(0);
}
friend Wint sqrt(Wint a)
{
Wint n;
n.assign(a.size() + 1 >> 1, 0);
for (int i = n.size() - 1, l, r; ~i; --i)
{
for (l = 0, r = a.base, n[i] = l + r >> 1; r - l > 1; n[i] = l + r >> 1)
{
if (n.cmp(n[i], i - 1, a))
r = n[i];
else
l = n[i];
}
a.sub(n, n[i], i - 1), n[i] += l + r >> 1;
}
for (int i = 0; i < n.size(); ++i)
n[i] >>= 1;
return n.trim(0);
}
/*
friend Wint sqrt(const Wint &a) //常规牛顿迭代实现的开平方算法，慢但是好敲
{
Wint b = a, c = (b + 1) / 2;
while (b != c)
swap(b, c), c = (b + a / b) / 2;
return c;
}
friend Wint sqrt(const Wint &a)
{
Wint ret, t;
ret.assign((a.size() + 1) >> 1, 0);
for (int i = ret.size() - 1, l, r; ~i; --i)
{
for (l = 0, r = a.base; r - l > 1;)
{
ret[i] = l + (r - l) / 2;
t = ret * ret;
if (a < t)
r = ret[i];
else
l = ret[i];
}
if (!l && i == ret.size() - 1)
ret.pop_back();
else
ret[i] = l;
}
return ret;
}
*/
};
分数
使用示范。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
struct Fraction
{
ll num, den;
Fraction(ll n = 0, ll d = 1) : num(n), den(d)
{
d = __gcd(num, den), num /= d, den /= d;
if (den < 0)
num = -num, den = -den;
}
friend Fraction operator+(const Fraction &A, const Fraction &B)
{
ll d = __gcd(A.den, B.den);
return Fraction(B.den / d * A.num + A.den / d * B.num, A.den / d * B.den);
}
Fraction &operator+=(const Fraction &c) { return *this = *this + c; }
Fraction operator-() const
{
Fraction r(*this);
return r.num = -r.num, r;
}
friend Fraction operator-(const Fraction &a, const Fraction &c) { return -c + a; }
Fraction &operator-=(const Fraction &c) { return *this = *this - c; }
friend Fraction operator*(const Fraction &A, const Fraction &B) { return Fraction(A.num * B.num, A.den * B.den); }
Fraction &operator*=(const Fraction &c) { return *this = *this * c; }
friend Fraction operator/(const Fraction &A, const Fraction &B) { return Fraction(A.num * B.den, A.den * B.num); }
Fraction &operator/=(const Fraction &c) { return *this = *this / c; }
friend Fraction operator%(const Fraction &a, const Fraction &c) { return Fraction(a.num * c.den % (c.num * a.den), a.den * c.den); }
Fraction &operator%=(const Fraction &c) { return *this = *this % c; }
friend bool operator==(const Fraction &a, const Fraction &b) { return a.num * b.den == a.den * b.num; }
friend bool operator!=(const Fraction &a, const Fraction &b) { return !(a == b); }
friend bool operator<(const Fraction &a, const Fraction &b) { return a.num * b.den < a.den * b.num; }
friend bool operator>(const Fraction &a, const Fraction &b) { return b < a; }
friend bool operator<=(const Fraction &a, const Fraction &b) { return !(a > b); }
friend bool operator>=(const Fraction &a, const Fraction &b) { return !(a < b); }
friend Fraction abs(Fraction f)
{
if (f.num < 0)
f.num = -f.num;
return f;
}
friend ostream &operator<<(ostream &os, const Fraction &f) { return !f.num ? os << 0 : f.den == 1 ? os << f.num : os << f.num << '/' << f.den; }
};
bitset 高精度
代替整型进行位运算，更方便并且可以处理超过最大整形范围大小的位集合。 你可以把 bitset 看作可以位运算的 bool 数组，换言之，bitset 的大小是固定的。因此，用 bitset 做状态压缩是很方便的，也可以方便的实现集合的交并补操作。 bitset 仅重载了相等不等和位运算符，原生不支持四则运算和大小比较，所以很少代替高精度数。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
typedef bitset<127> Bint;
/*
Bint b(unsigned long long u=0);//用u的低N位初始化b，N是一个可转成ULL类型的常量表达式，高位补0
Bint bs(string s,int pos,int m=string::npos,char zero='0',char one='1');//用s从pos位开始的m位初始化b，s中只含zero和one
b.size();//b的大小，即N
b.count();//b中1的个数
b[pos];//b中pos位的引用
b.set();//b全赋1
b.reset();//b全赋0
b.flip();//b全反转
b.to_ull();//b转成unsigned long long，b.size()>64时出错
b.to_string(char zero='0',char one='1');//按参数输出字符串
os<<b;//按'0'和'1'打印b
is>>b;//按'0'和'1'读入b，当下一个字符不是'0'或'1'或读到b.size()个数后停止
==、!=、<<、>>、&、|、^//保持它们在整型运算中的含义
*/
//大小比较，其他运算符类似
bool operator<(const Bint &a, const Bint &b)
{
for (int i = a.size() - 1; ~i; --i)
if (a[i] != b[i])
return a[i] < b[i];
return 0;
}
bool operator!=(const Bint &a, const Bint &b)
{
for (int i = a.size() - 1; ~i; --i)
if (a[i] != b[i])
return 1;
return 0;
}
//加法
Bint operator+(const Bint &a, const Bint &b) { return b.any() ? (a ^ b) + ((a & b) << 1) : a; }
Bint &operator+=(Bint &a, const Bint &b) { return a = a + b; }
//减法
Bint operator-(const Bint &a) { return Bint(1) + ~a; }
Bint &operator-=(Bint &a, const Bint &b) { return a += -b; }
Bint operator-(Bint a, const Bint &b) { return a -= b; }
//乘法
Bint operator*(Bint a, Bint b)
{
Bint r(0);
for (; b.any(); b >>= 1, a += a)
if (b[0])
r += a;
return r;
}
Bint &operator*=(Bint &a, const Bint &b) { return a = a * b; }
//整除，取模
Bint operator%=(Bint ÷nd, const Bint &divisor)
{
if (dividend < divisor || divisor.none())
return dividend;
Bint c, res = 0;
for (int k; divisor < dividend;)
{
for (k = 0, c = divisor; !(dividend < c); c <<= 1, ++k)
if (dividend < divisor + c)
{
res += 1 << k;
break;
}
if (dividend < divisor + c)
break;
res += 1 << (k - 1);
dividend -= c >> 1;
}
return dividend; //res是商
}
//输入输出，bitset已经原生重载了输入输出运算符，避免歧义。
istream &getb(istream &is, Bint &val)
{
int sign = 1, ch = is.get();
for (; !isdigit(ch) && ch != EOF; ch = is.get())
if (ch == '-')
sign = -sign;
for (val = 0; isdigit(ch); ch = is.get())
val = (val << 3) + (val << 1) + (ch ^ '0');
if (sign == -1)
val = -val;
return is.putback(ch);
}
ostream &putb(ostream &os, const Bint &val)
{
if (Bint(9) < val)
putb(os, val / 10);
return os.put(val.to_ulong() % 10 + '0');
}
C++语言相关
GCC 内置位运算
1
2
3
4
5
int __builtin_clz(unsigned int x); //求前缀0的个数
int __builtin_ctz(unsigned int x); //求后缀0的个数
int __builtin_ffs(unsigned int x); //x的二进制末尾最后一个1的位置，从1开始
int __builtin_popcount(unsigned int x); //x二进制中1的个数，相当于bitset::count()
int __builtin_parity(unsigned int x); //判断x的二进制中1的个数的奇（1）偶（0）性，这些函数都有相应的usigned long和usigned long long版本，只需要在函数名后面加上l或ll就可以了，比如__builtin_clzll
开栈
For C++
1
#pragma comment(linker, "/STACK:102400000,102400000") //For C++
For G++
1
2
3
4
5
6
int main() //For G++
{
int size = 256 << 20; //256MB
char *p = (char *)malloc(size) + size;
__asm__ __volatile__("movq %0, %%rsp\n" ::"r"(p)); //64bit，一定要最后写一句`exit(0);`退出程序，否则会得到非零退出的错误，可能RE。
}
读入优化
C 文件指针版
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
ll getll(FILE *in = stdin)
{
ll val = 0, sgn = 1, ch = getc(in);
for (; !isdigit(ch) && ch != EOF; ch = getc(in))
if (ch == '-')
sgn = -sgn;
for (; isdigit(ch); ch = getc(in))
val = val * 10 + ch - '0';
return ungetc(ch, in), sgn * val;
}
lf getlf(FILE *in = stdin)
{
lf val = getll(in), p = val < 0 ? -1 : 1;
ll ch = getc(in);
if (ch == '.')
for (ch = getc(in); isdigit(ch); ch = getc(in))
val += (p /= 10) * (ch - '0');
return ungetc(ch, in), val;
}
仿 C++IO 流沙雕版
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#define cin kin
struct Istream
{
char b[20 << 20], *i, *e; //20MB
Istream(FILE *in) : i(b), e(b + fread(b, sizeof(*b), sizeof(b), in)) {}
bool eof() const
{
return i == e;
}
Istream &operator>>(long long &val)
{
return val = strtoll(i, &i, 10 /*进制，取值2~36*/), *this;
}
Istream &operator>>(ll &val) //极限快
{
while (*i < '0')
++i; //无符号
for (val = 0; *i >= '0'; ++i)
val = (val << 3) + (val << 1) + *i - '0';
return *this;
}
Istream &operator>>(double &val)
{
return val = strtod(i, &i), *this;
}
Istream &operator>>(string &s)
{
while (!eof() && isspace(*i))
++i;
for (s.clear(); !eof() && !isspace(*i); ++i)
s += *i;
return *this;
}
} kin(stdin);