板子
目录
基本：缺省源，超快读入/输出，模大质数常用函数，浮点数常用函数，随机数据生成器，对拍器
网络流：最大流，最小费用最大流，有源汇有上下界最大流，有源汇有上下界最小流，无源汇有上下界最小费用可行流
字符串：扩展 KMP，后缀排序（SA-IS 算法），后缀自动机，后缀树
数据结构：并查集，配对堆，树状数组，动态树
图论：重链剖分，长链剖分求 \(k\) 级祖先，虚树，稳定婚姻（Gale-Shapley 算法），二分图最大匹配，一般图最大匹配
数论：数论基本操作，原根，离散对数，欧拉筛，万能欧几里得算法
多项式：多项式基本操作，多项式求导/积分，\(O\left(n\log n\right)\) 多项式乘法，\(O\left(n^2\right)\) 多项式乘法，\(O\left(n\log n\right)\) 多项式求逆，\(O\left(n^2\right)\) 多项式取模，\(O\left(n\log n\right)\) 多项式 \(\ln\)，\(O\left(n^2\right)\) 多项式 \(\ln\)，\(O\left(n\log n\right)\) 多项式 \(\exp\)，\(O\left(n^2\right)\) 多项式 \(\exp\)
线性代数：矩阵基本操作，拟阵交，带删除线性基
二维计算几何：二维计算几何基本操作，圆和直线交点，圆和一顶点在圆心的三角形的面积交，圆和多边形的面积交，二维凸包，多边形面积，凸包和半平面的面积交，三角形外接圆
三维计算几何：三维计算几何基本操作，\(O\left(n^2\right)\) 三维凸包
基本
缺省源
#include <bits/stdc++.h>
#define mp std::make_pair
#define pb push_back
#define eb emplace_back
#define fi first
#define se second
#define gc getchar
#define pc putchar
typedef unsigned char u8;
typedef unsigned short u16;
typedef unsigned int u32;
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
typedef std::pair <int, int> pii;
typedef std::pair <ll, ll> pll;
int read_int() {
char c = gc(); int ans = 0; bool neg = false;
while (!isdigit(c)) neg |= (c == '-'), c = gc();
while (isdigit(c)) ans = 10 * ans + c - '0', c = gc();
return neg ? -ans : ans;
}
void write_int(int x) {
if (x < 0) pc('-'), x = -x;
if (x > 9) write_int(x / 10);
pc(x % 10 + '0');
}
int min(int x, int y) {return x < y ? x : y;}
int max(int x, int y) {return x > y ? x : y;}
void _min(int &x, int y) {if (x > y) x = y;}
void _max(int &x, int y) {if (x < y) x = y;}
超快读入/输出
#undef gc
#undef pc
namespace io {
const int SIZE = (1 << 21) + 1;
char ibuf[SIZE], *iS, *iT, obuf[SIZE], *oS = obuf, *oT = oS + SIZE - 1;
inline char gc() {return (iS == iT ? (iT = (iS = ibuf) + fread(ibuf, 1, SIZE, stdin), (iS == iT ? EOF : *iS++)) : *iS++);}
inline void flush() {fwrite(obuf, 1, oS - obuf, stdout); oS = obuf;}
inline void pc(char x) {*oS++ = x; if (oS == oT) flush();}
struct _flush {~_flush() {flush();}} __flush;
}
using io::gc; using io::pc;
模大质数常用函数
ll fact[?], finv[?];
int fact_size = -1;
int plus(int x, int y) {return (x += y) >= mod ? x - mod : x;}
int minus(int x, int y) {return (x -= y) < 0 ? x + mod : x;}
int mul(int x, int y) {return (ll)x * y % mod;}
int fma(int x, int y, int z) {return (x + (ll)y * z) % mod;}
int fms(int x, int y, int z) {return (x + (ll)(mod - y) * z) % mod;}
int neg(int x) {return x == 0 ? 0 : mod - x;}
void _plus(int &x, int y) {if ((x += y) >= mod) x -= mod;}
void _minus(int &x, int y) {if ((x -= y) < 0) x += mod;}
void _mul(int &x, int y) {x = (ll)x * y % mod;}
void _fma(int &x, int y, int z) {x = (x + (ll)y * z) % mod;}
void _fms(int &x, int y, int z) {x = (x + (ll)(mod - y) * z) % mod;}
void _neg(int &x) {if (x) x = mod - x;}
int normal(int x) {return x < 0 ? x + mod : x;}
ll inv(int x) {
if (x <= fact_size) return fact[x - 1] * finv[x] % mod;
return x == 1 ? 1 : (mod - mod / x) * inv(mod % x) % mod;
}
ll sqr(int x) {return (ll)x * x % mod;}
ll sign(int x) {return x & 1 ? mod - 1 : 1;}
ll binom(int x, int y) {return x < y || y < 0 ? 0 : fact[x] * finv[y] % mod * finv[x - y] % mod;}
ll inv_binom(int x, int y) {return finv[x] * fact[y] % mod * fact[x - y] % mod;}
ll merge2(int x, int y) {return fact[x + y] * finv[x] % mod * finv[y] % mod;}
ll merge3(int x, int y, int z) {return fact[x + y + z] * finv[x] % mod * finv[y] % mod * finv[z] % mod;}
ll power(ll x, int y, ll ans = 1) {
for (; y; y >>= 1, x = x * x % mod)
if (y & 1) ans = ans * x % mod;
return ans;
}
void init_fact(int n) {
fact[0] = finv[0] = 1;
for (int i = 1; i <= n; i++) fact[i] = fact[i - 1] * i % mod;
finv[n] = inv(fact[n]);
for (int i = n; i >= 1; i--) finv[i - 1] = finv[i] * i % mod;
fact_size = n;
}
浮点数常用函数
bool le(ld x, ld y) {return x <= y + eps;}
bool ge(ld x, ld y) {return x >= y - eps;}
bool lt(ld x, ld y) {return x < y - eps;}
bool gt(ld x, ld y) {return x > y + eps;}
bool eq(ld x, ld y) {return le(x, y) && ge(x, y);}
ld sqr(ld x) {return x * x;}
随机数据生成器
std::mt19937_64 rnd(std::chrono::steady_clock::now().time_since_epoch().count());
int Random(int l, int r) {return rnd() % (r - l + 1) + l;}
对拍器
#include <bits/stdc++.h>
int test, gap = ?;
int main() {
while (true) {
system("?.exe"); double time_S = clock();
system("?.exe"); double time_T = clock();
system("?.exe");
if (system("fc ?.out ?.out > nul")) break;
if (++test % gap == 0)
printf("test %6d: Accepted! use time : %.5lf second\n", test, (time_T - time_S) / CLOCKS_PER_SEC);
}
printf("test %6d: Wrong Answer!\n", ++test);
return 0;
}
网络流
最大流
namespace maxflow {
const int max_N = ?, max_E = ?, INF = 2e9;
int nxt[2 * max_E], head[max_N], to[2 * max_E], cap[2 * max_E], tot;
int dist[max_N], _head[max_N], n, S, T;
void init(int _n, int _S, int _T) {
n = _n, S = _S, T = _T, tot = 1;
std::fill(head + 1, head + n + 1, 0);
}
void add_edge(int u, int v, int w) {
if (!w) return;
nxt[++tot] = head[u], head[u] = tot, to[tot] = v, cap[tot] = w;
nxt[++tot] = head[v], head[v] = tot, to[tot] = u, cap[tot] = 0;
}
bool BFS() {
std::fill(dist + 1, dist + n + 1, INF);
std::queue <int> Q; Q.push(S); dist[S] = 0;
while (!Q.empty()) {
int cur = Q.front(); Q.pop();
for (int i = head[cur]; i; i = nxt[i])
if (cap[i] > 0 && dist[to[i]] == INF)
dist[to[i]] = dist[cur] + 1, Q.push(to[i]);
}
return dist[T] < INF;
}
int dfs(int x, int low) {
if (x == T) return low; int used = 0, tmp;
for (int &i = _head[x]; i; i = nxt[i])
if (cap[i] > 0 && dist[to[i]] == dist[x] + 1) {
tmp = dfs(to[i], min(low - used, cap[i]));
if (tmp > 0) used += tmp, cap[i] -= tmp, cap[i ^ 1] += tmp;
if (used == low) break;
}
return used;
}
int solve() {
int ans = 0;
while (BFS()) std::copy(head + 1, head + n + 1, _head + 1), ans += dfs(S, INF);
return ans;
}
}
最小费用最大流
namespace costflow {
const int max_N = ?, max_E = ?, INF = 2e9;
int nxt[2 * max_E], head[max_N], to[2 * max_E], cap[2 * max_E], cost[2 * max_E], dist[max_N], _head[max_N];
int n, S, T, tot, maxflow, mincost;
bool inq[max_N], vis[max_N];
void init(int _n, int _S, int _T) {
n = _n, S = _S, T = _T, tot = 1;
std::fill(head + 1, head + n + 1, 0);
}
void add_edge(int u, int v, int w, int _cost) {
if (!w) return;
nxt[++tot] = head[u], head[u] = tot, to[tot] = v, cap[tot] = w, cost[tot] = _cost;
nxt[++tot] = head[v], head[v] = tot, to[tot] = u, cap[tot] = 0, cost[tot] = -_cost;
}
bool SPFA() {
std::fill(dist + 1, dist + n + 1, INF);
std::queue <int> Q; Q.push(T); dist[T] = 0, inq[T] = true;
while (!Q.empty()) {
int cur = Q.front(); Q.pop(); inq[cur] = false;
for (int i = head[cur]; i; i = nxt[i])
if (cap[i ^ 1] > 0 && dist[to[i]] > dist[cur] - cost[i]) {
dist[to[i]] = dist[cur] - cost[i];
if (!inq[to[i]]) Q.push(to[i]), inq[to[i]] = true;
}
}
return dist[S] < INF;
}
int dfs(int x, int low) {
if (x == T) return low; vis[x] = true; int used = 0, tmp;
for (int &i = _head[x]; i; i = nxt[i])
if (cap[i] > 0 && dist[to[i]] == dist[x] - cost[i] && !vis[to[i]]) {
tmp = dfs(to[i], min(low - used, cap[i]));
if (tmp > 0) cap[i] -= tmp, cap[i ^ 1] += tmp, used += tmp, mincost += tmp * cost[i];
if (used == low) break;
}
return used;
}
void solve() {
maxflow = mincost = 0;
while (SPFA())
std::copy(head + 1, head + n + 1, _head + 1),
std::fill(vis + 1, vis + n + 1, false),
maxflow += dfs(S, INF);
}
}
有源汇有上下界最大流
namespace bounded_maxflow {
int n, S, T, deg[maxflow::max_N];
void init(int _n, int _S, int _T) {
n = _n, S = _S, T = _T;
std::fill(deg + 1, deg + n + 1, 0);
maxflow::init(n + 2, n + 1, n + 2);
}
void add_edge(int u, int v, int l, int r) {
deg[u] -= l, deg[v] += l, maxflow::add_edge(u, v, r - l);
}
int solve() {
int sum = 0;
for (int i = 1; i <= n; i++) {
if (deg[i] > 0) maxflow::add_edge(n + 1, i, deg[i]), sum += deg[i];
if (deg[i] < 0) maxflow::add_edge(i, n + 2, -deg[i]);
}
maxflow::add_edge(T, S, maxflow::INF);
if (maxflow::solve() != sum) return -1;
for (int i = 1; i <= n; i++)
for (int j = maxflow::head[i]; j; j = maxflow::nxt[j])
if (maxflow::to[j] > n) maxflow::cap[j] = maxflow::cap[j ^ 1] = 0;
return maxflow::S = S, maxflow::T = T, maxflow::solve();
}
}
有源汇有上下界最小流
namespace bounded_minflow {
int n, S, T, deg[maxflow::max_N];
void init(int _n, int _S, int _T) {
n = _n, S = _S, T = _T;
std::fill(deg + 1, deg + n + 1, 0);
maxflow::init(n + 2, n + 1, n + 2);
}
void add_edge(int u, int v, int l, int r) {
deg[u] -= l, deg[v] += l, maxflow::add_edge(u, v, r - l);
}
int solve() {
int sum = 0;
for (int i = 1; i <= n; i++) {
if (deg[i] > 0) maxflow::add_edge(n + 1, i, deg[i]), sum += deg[i];
if (deg[i] < 0) maxflow::add_edge(i, n + 2, -deg[i]);
}
int tmp = maxflow::solve();
maxflow::add_edge(T, S, maxflow::INF);
if (tmp + maxflow::solve() != sum) return -1;
return sum - tmp;
}
}
无源汇有上下界最小费用可行流
namespace bounded_costflow {
int n, S, T, deg[costflow::max_N], ans;
bool ok;
void init(int _n, int _S, int _T) {
n = _n, S = _S, T = _T;
std::fill(deg + 1, deg + n + 1, 0);
costflow::init(n + 2, n + 1, n + 2);
}
void add_edge(int u, int v, int l, int r, int _cost) {
ans += l * _cost, deg[u] -= l, deg[v] += l;
costflow::add_edge(u, v, r - l, _cost);
}
void solve() {
int sum = 0;
for (int i = 1; i <= n; i++) {
if (deg[i] > 0) costflow::add_edge(n + 1, i, deg[i], 0), sum += deg[i];
if (deg[i] < 0) costflow::add_edge(i, n + 2, -deg[i], 0);
}
costflow::solve();
if (ok = (costflow::maxflow == sum)) ans += costflow::mincost;
}
}
字符串
扩展 KMP
namespace exKMP {
const int max_N = ?;
char A[max_N];
int n, nxt[max_N], extend[max_N];
void get_nxt() {
n = strlen(A + 1); nxt[1] = n;
if (n > 1) nxt[2] = std::mismatch(A + 2, A + n + 1, A + 1).fi - (A + 2);
for (int i = 3, tmp = 2; i <= n; i++) {
nxt[i] = max(0, min(nxt[i - tmp + 1], nxt[tmp] + tmp - i));
nxt[i] = std::mismatch(A + i + nxt[i], A + n + 1, A + nxt[i] + 1).fi - (A + i);
if (i + nxt[i] > tmp + nxt[tmp]) tmp = i;
}
}
void get_extend(char *B, int m) {
extend[1] = std::mismatch(A + 1, A + n + 1, B + 1).fi - (A + 1);
for (int i = 2, tmp = 1; i <= m; i++) {
extend[i] = max(0, min(nxt[i - tmp + 1], tmp + extend[tmp] - i));
extend[i] = std::mismatch(B + i + extend[i], B + m + 1, A + extend[i] + 1).fi - (B + i);
if (i + extend[i] > tmp + extend[tmp]) tmp = i;
}
}
}
后缀排序（SA-IS 算法）
namespace suffix_sort {
const int max_N = 2 * ::max_N;
int n, X[max_N], sa[max_N], rank[max_N], height[max_N], s[max_N];
bool type[max_N];
void induce_sort(int *s, bool *type, int n, int m, int *X, int X_size) {
static int cnt[max_N], pos[max_N];
std::fill(cnt + 1, cnt + m + 1, 0);
std::fill( sa + 1, sa + n + 1, 0);
for (int i = 1; i <= n; i++) cnt[s[i]]++;
std::partial_sum(cnt + 1, cnt + m + 1, cnt + 1);
std::copy(cnt + 1, cnt + m + 1, pos + 1);
for (int i = X_size; i; i--) sa[pos[s[X[i]]]--] = X[i];
std::copy(cnt, cnt + m, pos + 1);
for (int i = 1; i <= n; i++)
if (sa[i] > 1 && type[sa[i] - 1])
sa[++pos[s[sa[i] - 1]]] = sa[i] - 1;
for (int i = n; i >= 1; i--)
if (sa[i] > 1 && !type[sa[i] - 1])
sa[cnt[s[sa[i] - 1]]--] = sa[i] - 1;
}
void SAIS(int *s, bool *type, int *X, int n, int m) {
int X_size = 0; type[n] = false;
for (int i = n - 1; i >= 1; i--)
type[i] = (s[i] == s[i + 1] ? type[i + 1] : s[i] > s[i + 1]);
if (std::find(type + 1, type + n + 1, true) == type + n + 1)
return std::iota(sa + 1, sa + n + 1, 1);
std::fill(rank + 1, rank + n + 1, 0);
for (int i = 2; i <= n; i++)
if (type[i - 1] && !type[i]) X[rank[i] = ++X_size] = i;
induce_sort(s, type, n, m, X, X_size);
int *_s = s + n, next_m = 0;
for (int i = 1, j = 0, k = 0; i <= n; i++) {
if (!(k = rank[sa[i]])) continue;
if (!j || j == X_size || k == X_size || X[j + 1] - X[j] != X[k + 1] - X[k] ||
!std::equal(s + X[j], s + X[j + 1] + 1, s + X[k])) next_m++;
_s[k] = next_m, j = k;
}
if (next_m < X_size) SAIS(_s, type + n + 1, X + n + 1, X_size, next_m);
else for (int i = 1; i <= X_size; i++) sa[_s[i]] = i;
for (int i = 1; i <= X_size; i++) _s[i] = X[sa[i]];
induce_sort(s, type, n, m, _s, X_size);
}
void solve() {
for (int i = 1; i <= n; i++) s[i]++; s[++n] = 1;
SAIS(s, type, X, n, *std::max_element(s + 1, s + n + 1));
for (int i = 1; i <= n; i++) s[i]--; n--;
for (int i = 1; i <= n; i++) sa[i] = sa[i + 1];
for (int i = 1; i <= n; i++) rank[sa[i]] = i;
}
void get_height() {
for (int i = 1, j = 0, k; i <= n; i++) {
if (j > 0) j--;
if (rank[i] > 1) {
k = sa[rank[i] - 1];
while (i + j <= n && k + j <= n && s[i + j] == s[k + j]) j++;
height[i] = j;
}
}
}
}
后缀自动机
namespace SAM {
const int max_N = 2 * ?, max_S = ?;
int fa[max_N], trans[max_N][max_S], len[max_N];
int p, np, q, nq, cnt;
void init() {
std::fill(fa + 1, fa + cnt + 1, 0);
std::fill(len + 1, len + cnt + 1, 0);
for (int i = 1; i <= cnt; i++)
std::fill(trans[i], trans[i] + max_S, 0);
cnt = np = 1;
}
void extend(int x) {
p = np, np = ++cnt; len[np] = len[p] + 1;
while (p && !trans[p][x]) trans[p][x] = np, p = fa[p];
if (!p) {fa[np] = 1; return;} q = trans[p][x];
if (len[q] == len[p] + 1) {fa[np] = q; return;}
nq = ++cnt; std::copy(trans[q], trans[q] + max_S, trans[nq]);
len[nq] = len[p] + 1, fa[nq] = fa[q], fa[q] = fa[np] = nq;
while (p && trans[p][x] == q) trans[p][x] = nq, p = fa[p];
}
}
后缀树
namespace suffix_tree {
int pos[SAM::max_N];
int son[SAM::max_N][SAM::max_S], right[SAM::max_N];
bool vis[SAM::max_N];
void init(int n) {
std::fill(pos + 1, pos + 2 * n + 1, 0);
std::fill(right + 1, right + 2 * n + 1, 0);
std::fill(vis + 1, vis + 2 * n + 1, false);
for (int i = 1; i <= 2 * n; i++)
std::fill(son[i], son[i] + SAM::max_S, 0);
SAM::init(n);
}
void dfs(int x) {
for (int i = 0; i < SAM::max_S; i++) if (son[x][i])
dfs(son[x][i]), right[x] += right[son[x][i]];
}
void build(char *s, int n) {
init(n);
for (int i = n; i; i--)
SAM::extend(s[i]), pos[SAM::np] = i, right[SAM::np] = 1;
for (int i = 1; i <= SAM::cnt; i++) if (pos[i])
for (int j = SAM::fa[i], k = i; j && !vis[k]; k = j, j = SAM::fa[j])
pos[j] = pos[i], son[j][s[pos[i] + SAM::len[j]]] = k, vis[k] = true;
dfs(1);
}
}
数据结构
并查集
namespace DSU {
const int max_N = ?;
int fa[max_N], size[max_N];
void init(int n) {
std::iota(fa + 1, fa + n + 1, 1);
std::fill(size + 1, size + n + 1, 1);
}
int anc(int x) {return fa[x] == x ? x : fa[x] = anc(fa[x]);}
void link(int x, int y) {
if ((x = anc(x)) != (y = anc(y))) {
if (size[x] > size[y]) std::swap(x, y);
fa[x] = y, size[y] += size[x];
}
}
bool connected(int x, int y) {return anc(x) == anc(y);}
int num() {
int ans = 0;
for (int i = 1; i <= n; i++) ans += (fa[i] == i);
return ans;
}
}
配对堆
namespace pairing_heap {
const int max_N = ?;
int son[max_N][2], root[max_N], val[max_N];
void change(int rt, int x) {?}
void pushdown(int rt) {?}
int merge(int x, int y) {
if (!x || !y) return x | y;
pushdown(x), pushdown(y);
if (val[x] > val[y]) std::swap(x, y);
return son[y][1] = son[x][0], son[x][0] = y, x;
}
void link(int x, int y) {root[y] = merge(x, root[y]);}
int del(int rt) {
pushdown(rt);
int x = son[rt][0], y = son[x][1];
if (!x) return 0; if (!y) return x;
pushdown(x), pushdown(y);
son[rt][0] = son[y][1], son[x][1] = son[y][1] = 0;
return merge(merge(x, y), del(rt));
}
void pop(int x) {root[x] = del(root[x]);}
}
树状数组
namespace BIT {
const int max_N = ?;
int n, A[max_N];
void init(int _n) {n = _n, std::fill(A + 1, A + n + 1, 0);}
void modify(int x, int y) {while (x <= n) A[x] += y, x += (x & -x);}
int query(int x) {int ans = 0; while (x) ans += A[x], x &= (x - 1); return ans;}
}
动态树
namespace LCT {
const int max_N = ?;
int son[max_N][2], fa[max_N];
bool rev[max_N];
int dir(int x) {
if (son[DSU::anc(fa[x])][0] == x) return 0;
if (son[DSU::anc(fa[x])][1] == x) return 1;
return -1;
}
void push_up(int x) {?}
void add_tag(int x) {if (x) rev[x] = !rev[x], std::swap(son[x][0], son[x][1]);}
void push_down(int x) {
if (rev[x]) add_tag(son[x][0]), add_tag(son[x][1]), rev[x] = false;
}
void pull_down(int x) {if (dir(x) != -1) pull_down(DSU::anc(fa[x])); push_down(x);}
void rotate(int x) {
int f = fa[x], ff = fa[f], k = dir(x);
son[f][k] = son[x][k ^ 1], fa[son[x][k ^ 1]] = f;
if (dir(f) != -1) son[ff][dir(f)] = x; fa[x] = ff;
son[x][k ^ 1] = f, fa[f] = x; push_up(f), push_up(x);
}
void splay(int x) {
pull_down(x);
for (int y = fa[x]; dir(x) != -1; rotate(x), y = fa[x])
if (dir(y) != -1) rotate(dir(x) == dir(y) ? y : x);
}
void access(int x) {
for (int y = 0; x; x = fa[y = x]) splay(x), son[x][1] = y, push_up(x);
}
void make_root(int x) {
access(x), splay(x), add_tag(x);
}
void link(int x, int y) {
make_root(x), make_root(y), fa[x] = y;
}
void cut(int x, int y) {
make_root(x), access(y), splay(x);
son[x][1] = fa[y] = 0, push_up(x);
}
int find_root(int x) {
access(x), splay(x);
while (son[x][0]) x = son[x][0], push_down(x);
return splay(x), x;
}
}
图论
重链剖分
void dfs1(int x = 1) {
size[x] = 1, dep[x] = dep[fa[x]] + 1;
for (int i : T[x]) {
if (i == fa[x]) continue;
fa[i] = x, dfs1(i), size[x] += size[i];
if (size[i] > size[son[x]]) son[x] = i;
}
}
void dfs2(int x = 1, int _top = 1) {
Idfn[dfn[x] = ++*dfn] = x, top[x] = _top;
if (!son[x]) return; dfs2(son[x], _top);
for (int i : T[x])
if (i != fa[x] && i != son[x]) dfs2(i, i);
}
int LCA(int x, int y) {
while (top[x] != top[y])
(dep[top[x]] < dep[top[y]]) ? (y = fa[top[y]]) : (x = fa[top[x]]);
return dep[x] < dep[y] ? x : y;
}
int dist(int x, int y) {return dep[x] + dep[y] - 2 * dep[LCA(x, y)];}
int kth(int x, int k) {
int t = dep[x] - k;
while (dep[top[x]] > t) x = fa[top[x]];
return Idfn[dfn[x] - (dep[x] - t)];
}
长链剖分求 \(k\) 级祖先
namespace kth_ancestor {
const int max_N = ?, max_logN = ?;
int fa[max_N][max_logN], maxd[max_N], log_2[max_N], dep[max_N], son[max_N], top[max_N];
std::vector <int> vec[max_N][2];
void dfs1(int x) {
dep[x] = dep[fa[x][0]] + 1;
for (int i = 1; i <= log_2[dep[x]]; i++)
fa[x][i] = fa[fa[x][i - 1]][i - 1];
for (int i : T[x]) {
fa[i][0] = x, dfs1(i);
if (maxd[x] <= maxd[i]) maxd[x] = maxd[i] + 1, son[x] = i;
}
}
void dfs2(int x, int _top) {
top[x] = _top;
if (x == _top) {
vec[x][0].reserve(maxd[x] + 1); vec[x][1].reserve(maxd[x] + 1);
for (int i = 0, j = x; i <= maxd[x]; i++, j = son[j]) vec[x][0].pb(j);
for (int i = 0, j = x; i <= maxd[x] && j; i++, j = fa[j][0]) vec[x][1].pb(j);
}
if (!son[x]) return; dfs2(son[x], _top);
for (int i : T[x]) if (i != son[x]) dfs2(i, i);
}
int kth(int x, int k) {
if (k >= dep[x]) return rt; if (k == 0) return x;
int y = top[fa[x][log_2[k]]], d = dep[y] - (dep[x] - k);
return d < 0 ? vec[y][0][-d] : vec[y][1][d];
}
void init() {
for (int i = 2; i <= n; i++)
log_2[i] = log_2[i - 1] + (i == (i & -i));
dfs1(rt), dfs2(rt, rt);
}
}
虚树
namespace virtual_tree {
const int max_N = ?;
int nxt[max_N], head[max_N], to[max_N], st[max_N], top, tot, rt;
std::vector <int> V, vis;
void init() {
top = tot = 0; for (int i : vis) head[i] = 0; vis = V;
std::sort(V.begin(), V.end(), [&] (int x, int y) {return dfn[x] < dfn[y];});
V.erase(std::unique(V.begin(), V.end()), V.end());
}
void add_edge(int u, int v) {
nxt[++tot] = head[u], head[u] = tot, to[tot] = v;
}
void build(std::vector <int> &_V) {
V = _V; init(); st[++top] = V[0];
for (int i = 1; i < V.size(); i++) {
int lca = LCA(st[top], V[i]); vis.pb(lca);
while (top > 1 && dep[st[top - 1]] >= dep[lca])
add_edge(st[top - 1], st[top]), top--;
if (dep[st[top]] > dep[lca]) add_edge(lca, st[top]), top--;
if (st[top] != lca) st[++top] = lca;
st[++top] = V[i];
}
while (top > 1) add_edge(st[top - 1], st[top]), top--;
rt = st[top];
}
}
稳定婚姻（Gale-Shapley 算法）
namespace stable_match {
const int max_N = ?, INF = 2e9;
int n, valx[max_N][max_N], valy[max_N][max_N], linkx[max_N], linky[max_N];
int order[max_N][max_N], pos[max_N];
void init(int _n) {
n = _n;
for (int i = 1; i <= n; i++)
std::fill(valx[i] + 1, valx[i] + n + 1, -INF),
std::fill(valy[i] + 1, valy[i] + n + 1, -INF);
std::fill(linkx + 1, linkx + n + 1, 0);
std::fill(linky + 1, linky + n + 1, 0);
}
void set(int x, int y, int v1, int v2) {
valx[x][y] = v1, valy[y][x] = v2;
}
void link(int x, int y) {linky[x] = y, linkx[y] = x;}
void solve() {
std::fill(pos + 1, pos + n + 1, 1);
for (int i = 1; i <= n; i++)
std::iota(order[i] + 1, order[i] + n + 1, 1),
std::sort(order[i] + 1, order[i] + n + 1, [&] (int x, int y) {return valx[i][x] > valx[i][y];});
while (std::find(linky + 1, linky + n + 1, 0) <= linky + n)
for (int i = 1; i <= n; i++) if (!linky[i])
while (true) {
int y = order[i][pos[i]], x = linkx[y]; pos[i]++;
if (!x || valy[y][x] < valy[y][i]) {linky[x] = 0; link(i, y); break;}
}
}
}
二分图最大匹配
namespace bipartite_graph_matching {
const int max_N = ?;
int nl, nr, levell[max_N], levelr[max_N], linkl[max_N], linkr[max_N];
std::vector <int> G[max_N];
bool vis[max_N];
void init(int _nl, int _nr) {
nl = _nl, nr = _nr;
for (int i = 1; i <= nl; i++) G[i].clear();
std::fill(linkl + 1, linkl + nl + 1, 0);
std::fill(linkr + 1, linkr + nr + 1, 0);
}
void link(int x, int y) {linkr[x] = y, linkl[y] = x;}
void add_edge(int x, int y) {G[x].pb(y);}
bool bfs() {
std::queue <int> Q; bool res = false;
std::fill(levell + 1, levell + nl + 1, 0);
std::fill(levelr + 1, levelr + nr + 1, 0);
for (int i = 1; i <= nl; i++)
if (!linkr[i]) Q.push(i), levell[i] = 1;
while (!Q.empty()) {
int cur = Q.front(); Q.pop();
for (int i : G[cur])
if (levelr[i] == 0) {
levelr[i] = levell[cur] + 1;
if (linkl[i] && levell[linkl[i]] == 0)
levell[linkl[i]] = levelr[i] + 1, Q.push(linkl[i]);
else res = true;
}
}
return res;
}
bool dfs(int x) {
for (int i : G[x])
if (levelr[i] == levell[x] + 1 && !vis[i])
if (vis[i] = true, !linkl[i] || dfs(linkl[i])) return link(x, i), true;
return false;
}
int solve() {
int ans = 0;
while (bfs()) {
std::fill(vis + 1, vis + nr + 1, false);
for (int i = 1; i <= nl; i++)
if (!linkr[i] && dfs(i)) ans++;
}
return ans;
}
}
一般图最大匹配
namespace general_graph_matching {
const int max_N = ?;
int n, vis[max_N], col[max_N], match[max_N], pre[max_N];
std::vector <int> G[max_N];
std::queue <int> Q;
namespace DSU {
int fa[max_N];
void init() {std::iota(fa + 1, fa + n + 1, 1);}
int anc(int x) {return fa[x] == x ? x : fa[x] = anc(fa[x]);}
bool connected(int x, int y) {return anc(x) == anc(y);}
}
void init(int _n) {
n = _n;
std::fill(match + 1, match + n + 1, 0);
for (int i = 1; i <= n; i++) G[i].clear();
}
void add_edge(int u, int v) {G[u].pb(v), G[v].pb(u);}
int LCA(int x, int y) {
++*vis; x = DSU::anc(x), y = DSU::anc(y);
while (x) vis[x] = *vis, x = DSU::anc(pre[match[x]]);
while (vis[y] != *vis) y = DSU::anc(pre[match[y]]);
return y;
}
void contract(int x, int y, int w) {
while (DSU::anc(x) != w) {
pre[x] = y, y = match[x];
if (col[y] == 1) col[y] = 0, Q.push(y);
if (DSU::fa[x] == x) DSU::fa[x] = w;
if (DSU::fa[y] == y) DSU::fa[y] = w;
x = pre[y];
}
}
void augment(int x) {
while (x) match[x] = pre[x], std::swap(x, match[match[x]]);
}
bool find_augmenting_path(int s) {
std::fill(vis, vis + n + 1, 0);
std::fill(col + 1, col + n + 1, -1);
DSU::init(); Q.push(s); col[s] = 0;
while (!Q.empty()) {
int cur = Q.front(), rt; Q.pop();
for (int i : G[cur]) {
if (DSU::connected(cur, i) || col[i] == 1) continue;
if (col[i] == 0) rt = LCA(cur, i), contract(i, cur, rt), contract(cur, i, rt);
else {
col[i] = 1, pre[i] = cur;
if (!match[i]) {
while (!Q.empty()) Q.pop();
return augment(i), true;
}
else col[match[i]] = 0, Q.push(match[i]);
}
}
}
return false;
}
int solve() {
int ans = 0;
for (int i = 1; i <= n; i++)
if (!match[i] && find_augmenting_path(i)) ans++;
return ans;
}
}
数论
数论基本操作
int gcd(int x, int y) {return y ? gcd(y, x % y) : x;}
void exgcd(int a, int b, int &x, int &y, int &d) {
if (!b) x = 1, y = 0, d = a;
else exgcd(b, a % b, y, x, d), y -= a / b * x;
}
int inv(int a, int m) {
int b = m, x, y, d; exgcd(a, b, x, y, d);
return x %= m, x < 0 ? x + m : x;
}
原根
int get_root() {
std::vector <int> vec; int tmp = phi;
for (int i = 2; i * i <= p; i++) {
if (tmp % i) continue; tmp /= i, vec.pb(i);
while (tmp % i == 0) tmp /= i;
}
if (tmp != 1) vec.pb(tmp);
for (int i = 1; i < p; i++) {
bool ok = true;
for (int j : vec)
if (power(i, phi / j) == 1) {ok = false; break;}
if (ok) return i;
}
return -1;
}
离散对数
namespace BSGS {
std::unordered_map <int, int> Map;
int g, p, H, p1, p2;
void init(int _g, int _p, int _H) {
g = _g, p = _p, H = _H; Map.clear(); p1 = p2 = 1;
for (int i = 1; i <= H; i++, p1 = p1 * g % p);
for (int i = 0; i <= p; i += H, p2 = p2 * p1 % p) Map.ep(p2, i);
}
int solve(int x) {
for (int i = 0; i <= H; i++, x = x * g % p) {
auto it = Map.find(x);
if (it != Map.end()) {
int ans = it -> se - i;
return ans < 0 ? ans + p - 1 : ans;
}
}
return -1;
}
}
欧拉筛
void sieve(int n) {
not_prime[1] = true;
for (int i = 2; i <= n; i++) {
if (!not_prime[i]) prime[++*prime] = i;
for (int j = 1; j <= *prime && i * prime[j] <= n; j++) {
not_prime[i * prime[j]] = true;
if (i % prime[j] == 0) break;
}
}
}
万能欧几里得算法
data data_power(data x, ll y, data ans = data()) {
for (; y; y >>= 1, x *= x) if (y & 1) ans *= x;
return ans;
}
data solve(int p, int q, int r, int l, data U, data R) {
if (l == 0) return data();
if (p >= q) return solve(p % q, q, r, l, U, data_power(U, p / q) * R);
int h = (l * p + r) / q; if (h == 0) return data_power(R, l);
return data_power(R, (q - r - 1) / p) * U * solve(q, p, (q - r - 1) % p, h - 1, R, U) * data_power(R, l - (h * q - r - 1) / p);
}
多项式
多项式基本操作
typedef std::vector <int> poly;
namespace poly_base {
const int g = ?, max_L = ?, max_len = 1 << max_L;
u32 W[max_len + 7], rev[max_len + 7];
}
void prepare() {
using namespace poly_base;
for (int i = 1; i < max_len; i++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << max_L - 1);
ll wn = power(g, (mod - 1) / max_len);
for (int i = max_len / 2, p = 1; i < max_len; i++, p = p * wn % mod) W[i] = p;
for (int i = max_len / 2 - 1; i; i--) W[i] = W[i << 1];
}
void up_power2(poly &A) {
int l = 1; while (l < A.size()) l *= 2; A.resize(l);
}
void NTT(poly &A, int L, bool flag) {
using namespace poly_base; static ull F[max_len + 7], tmp;
for (u32 i = 0, o = max_L - L; i < A.size(); i++) F[rev[i] >> o] = A[i];
for (u32 l = 1; l < A.size(); l <<= 1) {
for (u32 i = 0; i < A.size(); i += 2 * l) {
ull *p1 = F + i, *p2 = p1 + l; u32 *p3 = W + l;
for (u32 j = 0; j < l; j++)
tmp = *p2 * *p3 % mod, *p2 = *p1 + mod - tmp, *p1 += tmp, p1++, p2++, p3++;
}
if (l == 32768) for (u32 i = 0; i < A.size(); i++) F[i] %= mod;
}
if (flag) {
A[0] = F[0] * (tmp = inv(A.size())) % mod;
for (u32 i = 1; i < A.size(); i++) A[i] = F[A.size() - i] * tmp % mod;
}
else for (u32 i = 0; i < A.size(); i++) A[i] = F[i] % mod;
}
poly &operator += (poly &A, poly B) {
if (A.size() < B.size()) A.resize(B.size());
for (int i = 0; i < B.size(); i++) _plus(A[i], B[i]);
return A;
}
poly operator + (poly A, poly B) {return A += B;}
poly &operator -= (poly &A, poly B) {
if (A.size() < B.size()) A.resize(B.size());
for (int i = 0; i < B.size(); i++) _minus(A[i], B[i]);
return A;
}
poly operator - (poly A, poly B) {return A -= B;}
poly &operator <<= (poly &A, int B) {
A.resize(A.size() + B);
for (int i = A.size() - 1; i >= B; i--) A[i] = A[i - B];
return std::fill(A.begin(), A.begin() + B, 0), A;
}
poly operator << (poly A, int B) {return A <<= B;}
poly &operator >>= (poly &A, int B) {
for (int i = B; i < A.size(); i++) A[i - B] = A[i];
return A.resize(A.size() - B), A;
}
poly operator >> (poly A, int B) {return A >>= B;}
多项式求导/积分
poly poly_derivative(poly A) {
for (int i = 1; i < A.size(); i++) A[i - 1] = (ll)A[i] * i % mod;
return A.pop_back(), A;
}
poly poly_integral(poly A) {
A.pb(0);
for (int i = A.size() - 1; i; i--) A[i] = A[i - 1] * inv(i) % mod;
return A[0] = 0, A;
}
\(O\left(n\log n\right)\) 多项式乘法
poly &operator *= (poly &A, poly B) {
int deg = A.size() + B.size() - 1, L = 0;
while ((1 << L) < deg) L++;
A.resize(1 << L), B.resize(1 << L);
NTT(A, L, false), NTT(B, L, false);
for (int i = 0; i < A.size(); i++) A[i] = (ll)A[i] * B[i] % mod;
return NTT(A, L, true), A.resize(deg), A;
}
poly operator * (poly A, poly B) {return A *= B;}
\(O\left(n^2\right)\) 多项式乘法
poly &operator *= (poly &A, poly B) {
static ull C[?];
int deg = A.size() + B.size() - 1;
std::fill(C, C + deg, 0);
for (int i = 0; i < A.size(); i++) {
for (int j = 0; j < B.size(); j++)
C[i + j] += (ull)A[i] * B[j];
if (!(i & 15))
for (int j = 0; j < deg; j++) C[j] %= mod;
}
A.resize(deg);
for (int i = 0; i < deg; i++) A[i] = C[i] % mod;
return A;
}
poly operator * (poly A, poly B) {return A.size() < B.size() ? A *= B : B *= A;}
\(O\left(n\log n\right)\) 多项式求逆
poly poly_inv(poly A) {
int deg = A.size(); up_power2(A);
poly B = {(int)inv(A[0])};
for (int l = 1, logl = 0; l < A.size(); l *= 2, logl++) {
poly A0 = {A.begin(), A.begin() + 2 * l}, B0 = B;
B.resize(4 * l), A0.resize(4 * l), B0.resize(4 * l);
NTT(A0, logl + 2, false), NTT(B0, logl + 2, false);
for (int i = 0; i < 4 * l; i++) B[i] = mul(B0[i], fms(2, A0[i], B0[i]));
NTT(B, logl + 2, true); B.resize(2 * l);
}
return B.resize(deg), B;
}
\(O\left(n^2\right)\) 多项式取模
poly operator %= (poly &A, poly B) {
while (A.size() >= B.size()) {
int tmp = (mod - A.back()) * inv(B.back()) % mod;
for (int i = 0; i < B.size(); i++)
_fma(A[A.size() - B.size() + i], tmp, B[i]);
A.pop_back();
}
return A;
}
poly operator % (poly A, poly B) {return A %= B;}
\(O\left(n\log n\right)\) 多项式 \(\ln\)
poly poly_ln(poly A) {
poly B = poly_derivative(A) * poly_inv(A);
return B.resize(A.size() - 1), poly_integral(B);
}
\(O\left(n^2\right)\) 多项式 \(\ln\)
poly poly_ln(poly A) {
poly B(A.size()); static ull tmp;
for (int i = 0; i + 1 < A.size(); i++) {
tmp = (ull)(mod - A[i + 1]) * (i + 1);
for (int j = 1; j <= i; j++) {
tmp += (ull)A[j] * B[i - j + 1];
if (!(j & 15)) tmp %= mod;
}
B[i + 1] = mod - tmp % mod;
}
for (int i = 1; i < A.size(); i++) B[i] = B[i] * inv(i) % mod;
return B;
}
\(O\left(n\log n\right)\) 多项式 \(\exp\)
poly poly_exp(poly A) {
int deg = A.size(); up_power2(A);
poly B = {1};
for (int l = 1; l < A.size(); l *= 2) {
B.resize(2 * l); poly C = poly_ln(B);
for (int i = 0; i < B.size(); i++) C[i] = minus(A[i], C[i]);
_plus(C[0], 1); B *= C; B.resize(2 * l);
}
return B.resize(deg), B;
}
\(O\left(n^2\right)\) 多项式 \(\exp\)
poly poly_exp(poly A) {
poly B(A.size()); B[0] = 1; static ull tmp;
for (int i = 1; i < A.size(); i++) A[i] = (ll)A[i] * i % mod;
for (int i = 0; i + 1 < A.size(); i++) {
tmp = 0;
for (int j = 0; j <= i; j++) {
tmp += (ull)B[j] * A[i - j + 1];
if (!(j & 15)) tmp %= mod;
}
B[i + 1] = tmp % mod * inv(i + 1) % mod;
}
return B;
}
线性代数
矩阵基本操作
typedef std::vector <std::vector <int> > matrix;
namespace matrix_base {
int size;
matrix O, I;
}
void init_matrix(int n) {
using namespace matrix_base;
size = n; O.clear(), I.clear();
O.resize(n), I.resize(n);
for (int i = 0; i < n; i++)
O[i].resize(n), I[i].resize(n), I[i][i] = 1;
}
matrix operator * (matrix A, matrix B) {
matrix C = matrix_base::O;
for (int i = 0; i < matrix_base::size; i++)
for (int k = 0; k < matrix_base::size; k++)
for (int j = 0; j < matrix_base::size; j++)
_fma(C[i][j], A[i][k], B[k][j]);
return C;
}
matrix matrix_power(matrix x, int y, matrix ans = matrix_base::I) {
for (; y; y >>= 1, x = x * x)
if (y & 1) ans = ans * x;
return ans;
}
拟阵交
namespace matriod_cross {
const int max_N = ?, INF = 2e9;
bool flag[max_N];
std::vector <int> G[max_N];
int n, pre[max_N], dist[max_N];
namespace matriod_1 {
void build() {?}
bool ins(int x) {?}
}
namespace matriod_2 {
void build() {?}
bool ins(int x) {?}
}
bool find_augmenting_path() {
for (int i = 1; i <= n + 2; i++) G[i].clear();
matriod_1::build(), matriod_2::build();
for (int i = 1; i <= n; i++) {
if (flag[i]) continue; bool tmp = true;
if (matriod_1::ins(i)) G[n + 1].pb(i); else tmp = false;
if (matriod_2::ins(i)) G[i].pb(n + 2); else tmp = false;
if (tmp) return flag[i] = true;
}
for (int i = 1; i <= n; i++) {
if (!flag[i]) continue; flag[i] = false;
matriod_1::build(), matriod_2::build();
for (int j = 1; j <= n; j++) {
if (i == j || flag[j]) continue;
if (matriod_1::ins(j)) G[i].pb(j);
if (matriod_2::ins(j)) G[j].pb(i);
}
flag[i] = true;
}
std::fill(dist + 1, dist + n + 3, INF);
std::queue <int> Q; Q.push(n + 1); dist[n + 1] = 0;
while (!Q.empty()) {
int cur = Q.front(); Q.pop();
for (int i : G[cur]) if (dist[i] == INF)
dist[i] = dist[cur] + 1, pre[i] = cur, Q.push(i);
}
if (dist[n + 2] == INF) return false;
for (int i = pre[n + 2]; i != n + 1; i = pre[i]) flag[i] = !flag[i];
return true;
}
int solve(int _n) {
n = _n; while (find_augmenting_path());
return std::accumulate(flag + 1, flag + n + 1, 0);
}
}
带删除线性基
namespace basis {
const int max_N = ?;
bitset A[max_N], from[max_N];
int id[max_N];
bool insert(int x) {
if (A[x].none()) return false;
from[x].reset(), from[x].set(x);
for (int i = 0; i < m; i++) {
if (!A[x].test(i)) continue;
if (id[i]) A[x] ^= A[id[i]], from[x] ^= from[id[i]];
else return id[i] = x, true;
}
return false;
}
void erase(int x) {
if (from[x].none()) return;
if (A[x].none()) {from[x].reset(); return;}
int pos = std::find(id, id + m, x) - id, g = 0;
for (int i = 1; i <= n; i++)
if (A[i].none() && from[i].test(x)) {g = i; break;}
if (g) {
for (int i = 1; i <= n; i++)
if (i != g && from[i].test(x)) from[i] ^= from[g];
from[g] = from[x], A[g] = A[x], from[x].reset(), A[x].reset(), id[pos] = g;
}
else {
for (int i = m - 1; i >= 0; i--)
if (from[id[i]].test(x)) {g = id[i], id[i] = 0; break;}
for (int i = 0; i < m; i++)
if (id[i] != g && from[id[i]].test(x)) from[id[i]] ^= from[g], A[id[i]] ^= A[g];
if (g == x) from[x].reset(), A[x].reset();
else from[g] = from[x], A[g] = A[x], from[x].reset(), A[x].reset(), id[pos] = g;
}
}
}
二维计算几何
二维计算几何基本操作
struct vec2 {
ld x, y;
vec2 (ld _x = 0, ld _y = 0) : x(_x), y(_y) {}
vec2 &read() {scanf("%Lf%Lf", &x, &y); return *this;}
vec2 operator - () const {return vec2(-x, -y);}
vec2 operator + (const vec2 &b) const {return vec2(x + b.x, y + b.y);}
vec2 operator - (const vec2 &b) const {return vec2(x - b.x, y - b.y);}
vec2 operator * (const ld &b) const {return vec2(x * b, y * b);}
vec2 operator / (const ld &b) const {return *this * (1.0 / b);}
ld operator * (const vec2 &b) const {return x * b.x + y * b.y;}
ld operator ^ (const vec2 &b) const {return x * b.y - y * b.x;}
ld norm2() {return sqr(x) + sqr(y);}
ld norm() {return sqrtl(sqr(x) + sqr(y));}
bool operator < (const vec2 &b) const {
return lt(x, b.x) ? true : (gt(x, b.x) ? false : lt(y, b.y));
}
bool operator == (const vec2 &b) const {
return eq(x, b.x) && eq(y, b.y);
}
bool operator << (const vec2 &b) const {
if (lt(y, 0) != lt(b.y, 0)) return lt(y, 0);
if (!eq(*this ^ b, 0)) return gt(*this ^ b, 0);
if (gt(x, 0) != gt(b.x, 0)) return gt(x, 0);
return fabs(x) + fabs(y) < fabs(b.x) + fabs(b.y);
}
ld arg() {return atan2l(y, x);}
vec2 trans(ld A, ld B, ld C, ld D) {return vec2(A * x + B * y, C * x + D * y);}
void output() {printf("%.15Lf %.15Lf\n", x, y);}
};
struct line {
ld A, B, C;
line (ld _A = 0, ld _B = 0, ld _C = 0) : A(_A), B(_B), C(_C) {}
line (vec2 u, vec2 v) : A(u.y - v.y), B(v.x - u.x), C(u ^ v) {}
vec2 norm_vec() {return vec2(A, B);}
ld norm2() {return sqr(A) + sqr(B);}
ld norm() {return sqrtl(sqr(A) + sqr(B));}
ld operator () (const vec2 &b) const {return A * b.x + B * b.y + C;}
};
ld dot(vec2 x, vec2 u, vec2 v) {return (u - x) * (v - x);}
ld cross(vec2 x, vec2 u, vec2 v) {return (u - x) ^ (v - x);}
ld angle(vec2 A, vec2 B) {return atan2l(A ^ B, A * B);}
int collinear(vec2 A, vec2 B, vec2 C) {
return eq(cross(C, A, B), 0) ? le(dot(C, A, B), 0) + 1 : 0;
}
vec2 intersection(line u, line v) {
return vec2(u.B * v.C - u.C * v.B, u.C * v.A - u.A * v.C) / (u.A * v.B - u.B * v.A);
}
bool parallel(line u, line v) {return eq(u.norm_vec() ^ v.norm_vec(), 0);}
bool perpendicular(line u, line v) {return eq(u.norm_vec() * v.norm_vec(), 0);}
bool same_dir(line u, line v) {return parallel(u, v) && gt(u.norm_vec() * v.norm_vec(), 0);}
line bisector(vec2 u, vec2 v) {return line(v.x - u.x, v.y - u.y, 0.5 * (u.norm2() - v.norm2()));}
ld dist(vec2 u, vec2 v) {return (u - v).norm();}
ld dist2(vec2 u, vec2 v) {return (u - v).norm2();}
ld dist(vec2 u, line v) {return v(u) / v.norm();}
ld dist2(vec2 u, line v) {return sqr(v(u)) / v.norm2();}
vec2 projection(vec2 u, line v) {return u - v.norm_vec() * (v(u) / v.norm2());}
vec2 symmetry(vec2 u, line v) {return u - v.norm_vec() * (2 * v(u) / v.norm2());}
圆和直线交点
std::vector <vec2> intersection(vec2 O, ld r2, line l) {
ld d = dist2(O, l);
if (gt(d, r2)) return std::vector <vec2> {};
if (eq(d, r2)) return std::vector <vec2> {projection(O, l)};
if (eq(d, 0)) {
vec2 dir = l.norm_vec() * sqrtl(r2 / l.norm2());
return std::vector <vec2> {O + dir.trans(0, -1, 1, 0), O + dir.trans(0, 1, -1, 0)};
}
else {
ld T = copysign(sqrtl(r2 / d - 1), l(O));
vec2 dir = -l.norm_vec() * (l(O) / l.norm2());
return std::vector <vec2> {O + dir.trans(1, T, -T, 1), O + dir.trans(1, -T, T, 1)};
}
}
圆和一顶点在圆心的三角形的面积交
ld inter_area(ld r2, vec2 A, vec2 B) {
if (eq(A ^ B, 0)) return 0;
std::vector <vec2> vec = intersection(vec2(), r2, line(A, B));
if (vec.empty()) return 0.5 * r2 * angle(A, B);
if (vec.size() == 1) vec.pb(vec[0]);
switch (2 * gt(A.norm2(), r2) + gt(B.norm2(), r2)) {
case 0 : return 0.5 * (A ^ B);
case 1 : return 0.5 * ((A ^ vec[1]) + r2 * angle(vec[1], B));
case 2 : return 0.5 * ((vec[0] ^ B) + r2 * angle(A, vec[0]));
case 3 : {
if (lt((A - vec[0]) * (B - vec[0]), 0))
return 0.5 * (r2 * (angle(A, vec[0]) + angle(vec[1], B)) + (vec[0] ^ vec[1]));
else return 0.5 * r2 * angle(A, B);
}
}
}
圆和多边形的面积交
ld inter_area(vec2 O, ld r2, std::vector <vec2> poly) {
ld ans = inter_area(r2, poly.back() - O, poly.front() - O);
for (int i = 1; i < poly.size(); i++)
ans += inter_area(r2, poly[i - 1] - O, poly[i] - O);
return ans;
}
二维凸包
std::vector <vec2> CH2D(std::vector <vec2> poly) {
if (poly.size() <= 2) return poly;
std::sort(poly.begin(), poly.end());
std::vector <vec2> ans1, ans2;
ans1.reserve(poly.size()), ans2.reserve(poly.size());
for (vec2 i : poly) {
while (ans1.size() > 1 && ge((ans1.back() - ans1[ans1.size() - 2]) ^ (i - ans1.back()), 0)) ans1.pop_back();
while (ans2.size() > 1 && le((ans2.back() - ans2[ans2.size() - 2]) ^ (i - ans2.back()), 0)) ans2.pop_back();
ans1.pb(i), ans2.pb(i);
}
std::reverse_copy(++ans1.begin(), --ans1.end(), std::back_inserter(ans2));
return ans2;
}
多边形面积
ld area(std::vector <vec2> poly) {
ld ans = poly.back() ^ poly.front();
for (int i = 1; i < poly.size(); i++) ans += poly[i - 1] ^ poly[i];
return 0.5 * ans;
}
凸包和半平面的面积交
ld inter_area(std::vector <vec2> poly, line l) {
std::vector <vec2> vec; vec.reserve(poly.size());
for (int i = 0, j = poly.size() - 1; i < poly.size(); j = i++) {
if (ge(l(poly[j]), 0)) vec.pb(poly[j]);
if (lt(l(poly[i]) * l(poly[j]), 0))
vec.pb(intersection(line(poly[i], poly[j]), l));
}
return vec.empty() ? 0 : area(vec);
}
三角形外接圆
vec2 circum_center(vec2 A, vec2 B, vec2 C) {
vec2 a = B - A, b = C - A; ld na = a.norm2(), nb = b.norm2();
return A + vec2(na * b.y - nb * a.y, nb * a.x - na * b.x) * (0.5 / (a ^ b));
}
ld circum_radius2(vec2 A, vec2 B, vec2 C) {
vec2 a = B - A, b = C - A;
return (a.norm2() * b.norm2() * (a - b).norm2()) / sqr(2 * (a ^ b));
}
三维计算几何
三维计算几何基本操作
struct vec3 {
ld x, y, z;
vec3 (ld _x = 0, ld _y = 0, ld _z = 0) : x(_x), y(_y), z(_z) {}
vec3 &read() {scanf("%Lf%Lf%Lf", &x, &y, &z); return *this;}
vec3 operator - () const {return vec3(-x, -y, -z);}
vec3 operator + (const vec3 &b) const {return vec3(x + b.x, y + b.y, z + b.z);}
vec3 operator - (const vec3 &b) const {return vec3(x - b.x, y - b.y, z - b.z);}
vec3 operator * (const ld &b) const {return vec3(x * b, y * b, z * b);}
vec3 operator / (const ld &b) const {return *this * (1.0 / b);}
ld operator * (const vec3 &b) const {return x * b.x + y * b.y + z * b.z;}
vec3 operator ^ (const vec3 &b) const {return vec3(y * b.z - z * b.y, z * b.x - x * b.z, x * b.y - y * b.x);}
ld norm2() {return sqr(x) + sqr(y) + sqr(z);}
ld norm() {return sqrtl(sqr(x) + sqr(y) + sqr(z));}
bool operator < (const vec3 &b) const {
if (lt(x, b.x)) return true; if (lt(b.x, x)) return false;
if (lt(y, b.y)) return true; if (lt(b.y, y)) return false;
return lt(z, b.z);
}
bool operator == (const vec3 &b) const {
return eq(x, b.x) && eq(y, b.y) && eq(z, b.z);
}
};
ld triple(vec3 A, vec3 B, vec3 C) {return A * (B ^ C);}
ld volume(vec3 A, vec3 B, vec3 C, vec3 D) {return triple(B - A, C - A, D - A);}
\(O\left(n^2\right)\) 三维凸包
namespace CH3D {
const int max_N = ?;
int n, cntF, cntG;
vec3 P[max_N], tmp[max_N];
bool vis[max_N][max_N];
struct triangle {
int v[3];
triangle (int v0 = 0, int v1 = 0, int v2 = 0) {
v[0] = v0, v[1] = v1, v[2] = v2;
}
vec3 norm_vec() {return (P[v[1]] - P[v[0]]) ^ (P[v[2]] - P[v[0]]);}
ld area() {return 0.5 * norm_vec().norm();}
}F[2 * max_N], G[2 * max_N];
vec3 rand_vec3() {
return vec3((ld)rand() / RAND_MAX, (ld)rand() / RAND_MAX, (ld)rand() / RAND_MAX);
}
bool see(triangle A, vec3 B) {return ge(((B - P[A.v[0]]) * A.norm_vec()), 0);}
void solve() {
F[1] = triangle(1, 2, 3), F[2] = triangle(3, 2, 1), cntF = 2;
for (int i = 1; i <= n; i++) tmp[i] = P[i], P[i] = P[i] + rand_vec3() * 1e-8;
for (int i = 4; i <= n; i++) {
cntG = 0;
for (int j = 1, tmp; j <= cntF; j++) {
if (!(tmp = see(F[j], P[i]))) G[++cntG] = F[j];
for (int k = 0; k < 3; k++)
vis[F[j].v[k]][F[j].v[(k + 1) % 3]] = tmp;
}
for (int j = 1; j <= cntF; j++) for (int k = 0; k < 3; k++)
if (vis[F[j].v[k]][F[j].v[(k + 1) % 3]] && !vis[F[j].v[(k + 1) % 3]][F[j].v[k]])
G[++cntG] = triangle(F[j].v[k], F[j].v[(k + 1) % 3], i);
cntF = cntG; std::copy(G + 1, G + cntG + 1, F + 1);
}
std::copy(tmp + 1, tmp + n + 1, P + 1);
}
}
