CC BY 4.0 (除特别声明或转载文章外)
如果这篇博客帮助到你,可以请我喝一杯咖啡~
本页是一个汇总页。
我的算法竞赛模板
本页是一个汇总页。
{% for post in site.tags[‘算法竞赛模板’] %}
{{ post.title }}
{{post.content}}
{% endfor %}
数据结构
以下数据结构均采用 ll 作为值类型,应用时根据需求调整。
typedef long long ll;
const ll INF = 1e9; //表示(值)正无穷,且两个正无穷相加不会溢出
const int NPOS = -1; //表示(下标)不存在
离散化
在 vector 基础上的离散化,使用 push_back()向其中插值,init()排序并离散化,ask 查询离散化之后的值,at/[]运算符查离散前的值。
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)]);
}
};
树状数组
模板中 Base 是对应的基础版本,支持单点修改区间查询。
一维
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); }
};
动态开点线段树
使用示例,支持区间线性变换、区间查询(最大值最小值区间和)。
这样写改可持久化也很方便,只要改down
函数为每次都新建节点即可,示例。
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);
}
}
};
树上莫队
按照欧拉序分块,使用 Tarjan 在生成欧拉序的同时预处理所有询问的 lca,预处理时间复杂度$O(n+q)$。 h 为查询图,即如果有一个询问(u,v),即在 h 上连$u\to v,v\to u$。多个询问边有序插入 h。
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
使用示例,如果要修改模数或者直接使用unsigned long long
的自然溢出的话直接修改 Mod 即可。
使用unsigned long long
的自然溢出快了 5 倍,但是容易被卡。
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;
}
};
暴力回文
时间复杂度$O(n^2)$,常数低,但会被ababababa
这样的数据卡。
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 为开头的回文串长度。
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 算法
使用示例,适用于边权为正的情况(无论有向图还是无向图),用于求单源最短路。
直接给出其优先队列优化的版本。另外,由于priority_queue
并不提供修改优先级的操作,为避免重复扩展,这里选择将新元素直接插入队列并在运行时判断该点是否被处理,并不影响结果的正确性。
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。稀疏图上性能优秀。
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。
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)$。
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 求费用流
使用示例,定义一条边的费用为流量*边长,求总费用最小的最大流。BellmanFord 算法找增广路,可能被卡但是可以跑负费用流(最大费用流)。
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 算法了。
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(作为终点)、其他顶点出度=入度。
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 的代码。
对于一个无向图的子图,当删除其中任意一个点后,不改变图内点的连通性,这样的子图叫做点的双连通子图。而当子图的边数达到最大时,叫做点的双连通分量。下面给出求点双连通分量的代码。
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
次,把所有点收缩到了一起。
无向图求边双连通分量&有向图求强连通分量
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
。
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)$。
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
即为最大匹配数。
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 )$,稀疏图上效果明显。
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
。
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;
}
}
}
}
};
树形图
最小生成树
无向图
同时给出 Prim 算法(生成新树)、Kruskal 算法(消耗小)。
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)$。
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\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$。
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减去这个值
}
}
}
};
欧拉筛
欧拉函数$\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
中。
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)$
前缀和
欧拉函数前缀和$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,看情况选用快速乘法。
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 时交换的对应位置(即保存的是置换)。
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;
}
};
快速数论变换
原理和 FFT 相同,解决特殊情况下 FFT 的浮点误差,并且可以在同余系进行变换。
对于形如$m=2^nc+1$的费马素数,记其原根为 g,则旋转因子为$g^{(m-1)/n}$,满足$g^{m-1}=1$且$2^n\mid m-1$。
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 |
快速沃尔什变换
如果要在同余系中进行运算,则下面代码需要修改。
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; }
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}$。
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项的系数$ |
置换
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]$中。
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}$
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,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
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;
}
};
异或线性基
若要查询第 k 小子集异或和,则把 k 写成二进制,对于是 1 的第 i 位,把从低位到高位第 i 个不为 0 的数异或进答案。若要判断是否有非空子集的异或和为 0,如果不存在自由基,那么说明只有空集的异或值为 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);
}
堆的操作
和 STL 的函数push_heap
、pop_heap
用法完全相同,实现的是大根堆。
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]$上单调/凹凸性唯一。
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);
}
};
插值法
拉格朗日插值法:插值多项式和插值基函数的形式对称,容易编程。但是,增加节点时,需要重新计算每一个插值基函数。要在$\pmod p$意义下进行的话,那么 p 只能是质数。 牛顿插值法:当插值节点增加时,之前已计算的结果仍然能用,每增加一个节点,只要再增加一项即可,从而避免了重复性计算。如果要 mod 非质数的话,那么就要用牛顿插值法。
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。
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 自带大小比较为字典序比较, !=
、 ==
运算,其余需要时一定记得重载!
减法,当被减数小于减数时减为 0。
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 高精度
代替整型进行位运算,更方便并且可以处理超过最大整形范围大小的位集合。 你可以把 bitset 看作可以位运算的 bool 数组,换言之,bitset 的大小是固定的。因此,用 bitset 做状态压缩是很方便的,也可以方便的实现集合的交并补操作。 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 ÷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 内置位运算
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);