# 我的算法竞赛模板

{% for post in site.tags[‘算法竞赛模板’] %}

{{post.content}}

{% endfor %}

# 数据结构

typedef long long ll;
const ll INF = 1e9;  //表示（值）正无穷，且两个正无穷相加不会溢出
const int NPOS = -1; //表示（下标）不存在


## 离散化

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


## 并查集

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


## 单调队列和单调栈

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)$求静态区间最小值。

/*
//可选优化
#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)]);
}
};


## 树状数组

### 一维

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


### 二维

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


## 动态开点线段树

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

### 按子树大小分裂

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


### 按值大小分裂

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


## 莫队

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


### 带修莫队

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


### 树上莫队

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


### 树上带修莫队

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

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

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 自动机

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


### 暴力回文

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


### 线性回文

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])。

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


# 图论

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 算法

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 算法

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


### Floyed 求多源最短路

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 短路

Astar 算法是只有到达终点的时候才能统计答案，这导致可能拓展很多个状态才能得到一个用来更新答案的有效状态。例如一个 n 元环，当我们到达终点之后,可能还要拓展 n 次才能得到下一个状态，于是时间复杂度就被卡到$O(nk)$。

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 求最大流

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 求费用流

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 求费用流

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 求费用流

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]$的边。

## 欧拉路

• 无向图：当仅当该图所有顶点的度数为偶数（此时为回路），或除两个度数为奇数外（作为路径的起点和终点）、其余全为偶数。
• 有向图：当仅当该图所有顶点出度=入度（此时为回路），或一个顶点出度=入度+1（作为起点）、另一个顶点入度=出度+1（作为终点）、其他顶点出度=入度。
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);
}
};


## 连通性

### 无向图求割和双连通分量

• 割边：在连通图中，删除了连通图的某条边后，图不再连通。这样的边被称为割边，也叫做桥。
• 割点：在连通图中，删除了连通图的某个点以及与这个点相连的边后，图不再连通。这样的点被称为割点。构造 dfs 搜索树，在树上有两类节点可以成为割点：
• 对根节点 u，若其有两棵或两棵以上的子树，则该根结点 u 为割点；
• 对非根非叶节点 u，若其中的某棵子树的节点均没有指向 u 的祖先节点的回边，说明删除 u 之后，根结点与该棵子树的节点不再连通；则节点 u 为割点。

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


### 无向图求边双连通分量&有向图求强连通分量

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

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 暴力搜索求字典序最小的解

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


## 二分图

Hall 定理：二分图中的两部分顶点组成的集合分别为 X,Y ，则有一组无公共点的边，一端恰好为组成 X 的点的充分必要条件是：X 中的任意 k 个点至少与 Y 中的 k 个点相邻。对于区间图只需要考虑极端情况，线段树维护。

### Hungary 求最大匹配

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 求最大匹配

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 求最优完备匹配

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


## 带花树

### 一般图最大匹配

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


## 树形图

### 最小生成树

#### 无向图

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; //视情况选择去注释
}


#### 有向图

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

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


### 点剖（点分治）

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


# 数论

## 辗转相除法

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$互质的充要条件是$ax+by=1$有整数解。

## 同余系运算

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$。

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


### 中国剩余定理解同余方程组

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


### 二次剩余

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减去这个值
}
}
}
};


## 欧拉筛

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


### 直接求欧拉函数

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)$

#### 莫比乌斯反演

[\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 大数素因子分解

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


## 快速变换

### 蝴蝶变换（雷德变换）

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


### 快速傅里叶变换

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


### 快速数论变换

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 + 1rkg，mod(r * 2 ^ k + 1)的原根
3112
5122
17143
97355
193365
257183
768115917
1228931211
409615133
655371163
78643331810
576716911193
73400337203
2306867311213
10485760125223
1677721615253
4697620497263
998244353119233
1004535809479213
2013265921152731
228170137717273
32212254733305
7516192768135313
773094113299337
20615843020933622
206158430208115377
27487790694415393
65970697666573415
395824185999379425
791648371998739435
26388279066624115447
123145302310912135453
133700613937561719463
379991218559385727475
4222124650659841154819
78812993478983697506
315251973915934737523
1801439850948198415556
194555503902405427327565
417934045419982028929573

### 快速沃尔什变换

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

void tf(ll &l, ll &r) { l += r; }
void utf(ll &l, ll &r) { l -= r; }


#### OR

void tf(ll &l, ll &r) { r += l; }
void utf(ll &l, ll &r) { r -= l; }


#### XOR

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


## Pell 方程

$x_n=x_0x_{n-1}+Dy_0y_{n-1},y_n=y_0x_{n-1}+x_0y_{n-1}$

$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 为质数。

# 组合数学

## 组合数取模

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)$

## Stirling 数

### 斯特林近似公式

$n!\approx\sqrt{2\pi n}(\frac{n}{e})^n$

k 个球m 个盒子空盒子方案数

## 置换

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


## 生成字典序

### 下一排列

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 个元素的错排数$D_n$满足递推公式：$D_1=0,D_2=1,D_n=(n−1)(D_{n−2}+D_{n−1})$

## 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})$

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)$。

## Bell 数

1,2,5,15,52,203,877,4140,21147,115975,678570,4213597,27644437,190899322,1382958545,...


## 超立方体

n 维超立方体有$2^{n−i}C(n,i)$个 i 维元素。

## 枚举位集 I 的非空子集 J

for(J=I; J; J=I&J−1) {}


# 线性代数

## 矩阵

typedef array<array<ll, N>, N> Mat;


### 乘法和快速幂

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


### 行列式

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


### 高斯消元

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 是否与已插入的线性无关。

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


### 异或线性基

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


# 离散数学

## 归并排序求逆序对

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


## 约瑟夫问题

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


## 堆的操作

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 )$

# 数学分析

## 增长趋势

[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’(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]]

[S=-\int_\alpha^\beta x’(t)y(t)\,dt=-\int_\alpha^\beta y(t)\,dx(t)=-\oint_\Gamma y\,dx=\oint_\Gamma x\,dy]

[S=\frac{1}{2}\int_\alpha^\beta r^2(\theta)\,d\theta]

### 体积

[V=\int_a^bS(x)\,dx]

[V=\pi\int_a^by^2\,dx]

### 旋转体侧面积

[\begin{cases} x=x(t),
y=y(t), \end{cases} t\in[\alpha,\beta]]

[s=2\pi\int_\alpha^\beta y(t)\sqrt{[x’(t)]^2+[y’(t)]^2}\,dt]

### 方向导数

[\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)}]

[\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}]

### 曲率

[\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}}}]

[K=\frac{\mid y’‘\mid }{(1+y’^2)^\frac{3}{2}}]

### 空间曲线的切线与法平面

[\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]]

[\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}]

[\tau=\pm(\frac{\partial(F,G)}{\partial(y,z)},\frac{\partial(F,G)}{\partial(z,x)},\frac{\partial(F,G)}{\partial(x,y)})]

### 空间曲面的切平面与法线

[\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]]

[(F_x,F_y,F_z)_P\cdot(x’(t_0),y’(t_0),z’(t_0))=0]

[\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}]

[\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(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(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^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]

## 二分求零点、三分求极值点

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); //极小值
}


## 自适应辛普森求积分

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


## 插值法

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


[]

# 计算几何

## 二维

Morley 定理：三角形每个内角的三等分线相交成等边三角形。

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


## 三维

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 自带大小比较为字典序比较， !=== 运算，其余需要时一定记得重载！

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


## 分数

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 高精度

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 &dividend, 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 内置位运算

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++

#pragma comment(linker, "/STACK:102400000,102400000") //For C++


### For G++

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 文件指针版

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 流沙雕版

#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);