一晚上写了17KB代码...
目录
最大公约数
i64 gcd(i64 x, i64 y) { return y ? gcd(y, x % y) : x; }
快速幂
int power(int x, int k) {
int res = 1;
while(k) {
if(k & 1) res = 1ll * res * x % P;
x = 1ll * x * x % P; k >>= 1;
}
return res;
}
exgcd
void exgcd(i64 a, i64 b, i64 &x, i64 &y) {
if(! b) return x = 1, y = 0, void();
exgcd(b, a % b, x, y);
i64 t = x; x = y; y = t - a / b * y;
}
线性求逆元
// inv[i] = (P - P / i) * inv[P % i]
inv[0] = inv[1] = 1;
for(int i = 2; i <= n; i ++) inv[i] = 1ll * (P - P / i) * inv[P % i] % P;
线性筛
void sieve(int n = 1e6) {
static int vis[N], p[N], tot;
for(int i = 2; i <= n; i ++) {
if(! vis[i]) p[++ tot] = i;
for(int j = 1; j <= tot && i * p[j] <= n; j ++) {
vis[i * p[j]] = 1;
if(i % p[j] == 0) break;
}
}
}
光速乘
// 时代变了
inline i64 MulMod(i64 x, i64 y, i64 mod) {
return (__int128) x * y % mod;
}
Miller-Rabin 素数测试
inline bool Miller_Rabbin(i64 n) {
static int pri[10] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 61};
lep (i, 0, 9) if(n % pri[i] == 0) return n == pri[i];
lep (i, 0, 9) {
i64 r = n - 1, b = 0;
while(! (r & 1)) r >>= 1, b ++;
i64 w = power(pri[i], r, n), p = w;
lep (j, 1, b) {
w = mul(w, w, n);
if(w == 1 && p != n - 1 && p != 1) return 0;
p = w;
}
if(w != 1) return 0;
}
return 1;
}
整除分块
i64 solve(int n) {
i64 ans = 0;
for(int l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
ans += (S[r] - S[l - 1]) * F(n / l);
ans %= P;
}
return ans;
}
PollardRho
inline i64 run(i64 x, i64 n, i64 c) {
return (mul(x, x, n) + c) % n;
}
inline i64 Pollard(i64 n) {
if(n == 4) return 2;
i64 c = Rand() % (n - 1) + 1, a = 0, b = 0;
a = run(a, n, c);
b = run(b, n, c); b = run(b, n, c);
while(a ^ b) {
i64 g = __gcd(abs(a - b), n);
if(g != 1) return g;
a = run(a, n, c);
b = run(b, n, c); b = run(b, n, c);
}
return n;
}
inline i64 Pollard_Rho(i64 n) {
if(Miller_Rabbin(n)) return -1;
i64 res;
while((res = Pollard(n)) == n);
return res;
}
KMP 字符串匹配
// get border
for(int i = 2, j = 0; i <= n; i ++) {
while(j && s[j + 1] != s[i]) j = p[j];
if(s[j + 1] == s[i]) j ++;
p[i] = j;
}
Aho-Corasick 自动机
void build() {
std :: queue<int> q;
for(int i = 0; i < 26; i ++) if(ch[0][i]) q.push(ch[0][i]);
while(q.size()) {
int x = q.front(); q.pop();
val[x] |= val[fail[x]];
for(int i = 0; i < 26; i ++)
if(ch[x][i]) {
fail[ch[x][i]] = ch[fail[x]][i];
q.push(ch[x][i]);
}
else ch[x][i] = ch[fail[x]][i];
}
}
后缀自动机
int node = 1, las = 1, ch[N][26], len[N], fa[N];
void extend(int c) {
int p = ++ node, now = las; las = p;
len[p] = len[now] + 1;
while(now && ! ch[now][c]) now = fa[now];
if(! now) return fa[p] = 1, void();
int x = ch[now][c];
if(len[x] == len[now] + 1) return fa[p] = x, void();
int y = ++ node; len[y] = len[now] + 1;
fa[y] = fa[x]; fa[x] = fa[p] = y;
memcpy(ch[y], ch[x], sizeof(ch[x]));
while(now && ch[now][c] == x) ch[now][c] = y, now = fa[now];
}
Z函数
n = strlen(_s + 1); m = strlen(_t + 1);
//std :: swap(n, m); std :: swap(_s, _t);
for(int i = 1; i <= n; i ++) s[i] = _s[i];
s[n + 1] = '*';
for(int i = 1; i <= m; i ++) s[i + n + 1] = _t[i];
z[1] = n;
int l = 0, r = 0;
for(int i = 2; i <= m + n + 1; i ++) {
if(i <= r) z[i] = std :: min(z[i - l + 1], r - i + 1);
else z[i] = 0;
while(s[z[i] + 1] == s[i + z[i]]) z[i] ++;
if(i + z[i] - 1 > r) r = i + z[i] - 1, l = i;
}
Manacher
s[++ len] = '#';
lep (i, 1, n) s[++ len] = str[i], s[++ len] = '#';
s[++ len] = '%';
int mx = 0, ps = 0, res = 0;
lep (i, 2, len) {
if(i < mx) p[i] = min(p[ps * 2 - i], mx - i + 1);
else p[i] = 1;
while(s[i - p[i]] == s[i + p[i]]) p[i] ++;
res = max(res, p[i] - 1);
if(i + p[i] - 1 > mx) mx = i + p[i] - 1, ps = i;
}
回文自动机
inline int getfail(int x, int dex) {
while(s[dex - len[x] - 1] != s[dex]) x = fail[x];
return x;
}
inline int New(int x) { len[++ tot] = x; return tot; }
void work() {
// 偶根指向奇根
s[0] = -1; fail[0] = 1; las = 0;
len[0] = 0; len[1] = -1; tot = 1;
for(int i = 1; i <= n; i ++) {
s[i] -= 'a';
int p = getfail(las, i);
if(! ch[p][s[i]]) {
int q = New(len[p] + 2);
fail[q] = ch[getfail(fail[p], i)][s[i]];
ch[p][s[i]] = q;
}
cnt[las = ch[p][s[i]]] ++;
}
i64 ans = 0;
for(int i = tot; i >= 1; i --)
cnt[fail[i]] += cnt[i],
ans = std :: max(ans, 1ll * cnt[i] * len[i]);
}
倍增 LCA
void dfs(int x, int fx) {
dep[x] = dep[fx] + 1;
fa[x][0] = fx;
for(int i = 1; i <= lg[dep[x]]; i ++)
fa[x][i] = fa[fa[x][i - 1]][i - 1];
for(int y : e[x]) if(y != fx) dfs(y, x);
}
int lca(int x, int y) {
if(dep[x] < dep[y]) std :: swap(x, y);
while(dep[x] > dep[y]) x = fa[x][lg[dep[x] - dep[y]]];
if(x == y) return x;
for(int i = lg[dep[x]]; ~ i; i --)
if(fa[x][i] ^ fa[y][i]) x = fa[x][i], y = fa[y][i];
return fa[x][0];
}
ST表 LCA
int dfn[N], st[N << 1][20], lg[N << 1], dep[N];
void dfs(int x, int fx) {
st[++ dfn[0]][0] = x; dep[x] = dep[fx] + 1;
for(int y : e[x]) if(y != fx) {
dfs(y, x); st[++ dfn[0]][0] = x;
}
}
int lca(int x, int y) {
x = dfn[x]; y = dfn[y];
if(x > y) std :: swap(x, y);
int d = lg[y - x + 1];
return dep[st[x][d]] < dep[st[y - (1 << d) + 1][d]] ? st[x][d] : st[y - (1 << d) + 1][d];
}
void build() {
for(int i = 2; i <= n; i ++) lg[i] = lg[i >> 1] + 1;
dfs(1, 0);
for(int i = 1; (1 << i) <= (n << 1); i ++)
for(int j = 1; j + (1 << i) - 1 <= (n << 1); j ++)
if(dep[st[j][i - 1]] < dep[st[j + (1 << (i - 1))][i - 1]])
st[j][i] = st[j][i - 1];
else
st[j][i] = st[j + (1 << (i - 1))][i - 1];
}
轻重链剖分
int dep[N], top[N], sz[N], son[N], fa[N];
void dfs1(int x, int fx) {
sz[x] = 1; dep[x] = dep[fx] + 1; fa[x] = fx;
for(int y : e[x]) if(y != fx) {
dfs1(y, x); sz[x] += sz[y];
if(sz[y] > sz[son[x]]) son[x] = y;
}
}
void dfs2(int x, int topx) {
top[x] = topx; dfn[x] = ++ tot;
if(son[x]) dfs2(son[x], topx);
for(int y : e[x]) if(y != fx && y != son[x]) dfs2(y, y);
}
点分治
int vis[N], rt, all, sz[N];
void findrt(int x, int fx) {
sz[x] = 1;
int mx = 0;
for(int y : e[x]) if(! vis[y] && y != fx) {
findrt(y, x);
sz[x] += sz[y];
mx = std :: max(mx, sz[y]);
}
mx = std :: max(mx, all - sz[y]);
if(mx * 2 <= all || rt == 0) rt = x;
}
void solve(int x) {
calc(x);
vis[x] = 1;
for(int y : e[x]) if(! vis[y]) {
all = sz[y];
findrt(y, 0);
findrt(rt, 0);
fa[rt] = x;
solve(rt);
}
}
虚树
inline int isfa(int x, int y) {
return dfn[x] <= dfn[y] && dfn[y] <= dfn[x] + sz[x] - 1;
}
void build(std :: vector<int> &vec) {
std :: sort(vec.begin(), vec.end(), [&] (int a, int b) {
return dfn[a] < dfn[b];
} );
std :: vector<int> node = vec;
node.push_back(1);
for(int i = 0; i < vec.size() - 1; i ++)
node.push_back(lca(vec[i], vec[i + 1]));
std :: sort(node.begin(), node.end());
node.erase(std :: unique(node.begin(), node.end()), node.end());
std :: sort(node.begin(), node.end(), [&] (int a, int b) {
return dfn[a] < dfn[b];
} );
std :: vector<int> stk;
stk.push_back(node[0]);
for(int i = 1; i < node.size(); i ++) {
int x = node[i];
while(! isfa(stk.back(), x)) stk.pop_back();
push(stk.back(), x), stk.push_back(x);
}
}
长链剖分
就硬剖。
缩点
void tarjan(int x) {
stk[++ top] = x;
dfn[x] = low[x] = ++ dfn[0];
for(int y : e[x]) {
if(! dfn[y]) tarjan(y), low[x] = std :: min(low[x], low[y]);
else if(! co[y]) low[x] = std :: min(low[x], dfn[y]);
}
if(low[x] == dfn[x]) {
co[x] = ++ col;
while(stk[top] != x) {
co[stk[top]] = col;
top --;
}
top --;
}
}
割点
void tarjan(int x, int fx) {
dfn[x] = low[x] = ++ num;
int sz = 0;
for(int i = head[x]; i; i = e[i].next) {
int y = e[i].to;
if(y == fx) continue;
if(! dfn[y]) {
tarjan(y, x), low[x] = min(low[x], low[y]), sz ++;
if(low[y] >= dfn[x] && x != rt) bj[x] = 1;
}
else {
low[x] = min(low[x], dfn[y]);
}
}
if(x == rt && sz >= 2) bj[x] = 1;
}
割边
void tarjan(int x, int lst) {
dfn[x] = low[x] = ++ dfn[0];
for(int i = head[x]; i; i = e[i].next) if(i != lst) {
int y = e[i].to;
if(! dfn[y]) {
tarjan(y, i ^ 1), low[x] = std :: min(low[x], low[y]);
if(low[y] > dfn[x]) vis[i] = vis[i ^ 1] = 1;
}
else low[x] = std :: min(low[x], dfn[y]);
}
}
圆方树
struct edge { int to, next; } e[M];
int cnt, head[N];
inline void add(int x, int y) { e[++ cnt] = {y, head[x]}; head[x] = cnt; }
int dfn[N], low[N], num, stk[N], top;
int w[N], sz[N];
vector<int> T[N];
inline void tarjan(int x) {
stk[++ top] = x;
dfn[x] = low[x] = ++ num;
for(R int i = head[x]; i; i = e[i].next) {
int y = e[i].to;
if(!dfn[y]) {
tarjan(y);
low[x] = min(low[x], low[y]);
if(low[y] == dfn[x]) {
tot ++; w[tot] = 0;
while(1){
T[tot].push_back(stk[top]);
T[stk[top]].push_back(tot);
top --; w[tot] ++;
if(stk[top + 1] == y) break;
}
w[tot] ++;
T[tot].push_back(x);
T[x].push_back(tot);
}
}
else low[x] = min(low[x], dfn[y]);
}
}
Dijskra
略
Floyed
for(int k = 1; k <= n; k ++)
for(int i = 1; i <= n; i ++)
for(int j = 1; j <= n; j ++)
d[i][j] = std :: min(d[i][j], d[i][k] + d[k][j]);
最大流
struct MaxFlow {
static const int N = 2e3 + 10;
static const int M = 5e7 + 10;
static const int INF = 0x3f3f3f3f;
struct edge {
int to, next, f;
} e[M];
int cnt = 1, head[N];
inline void _push(int x, int y, int f) {
e[++ cnt] = (edge) {
y, head[x], f
};
head[x] = cnt;
}
inline void push(int x, int y, int f) {
_push(x, y, f);
_push(y, x, 0);
}
int dep[N], cur[N];
int S, T;
bool bfs() {
queue<int> q;
lep(i, 1, T) dep[i] = INF, cur[i] = head[i];
q.push(S);
dep[S] = 0;
while (q.size()) {
int x = q.front();
q.pop();
for (int i = head[x]; i; i = e[i].next)
if (e[i].f) {
int y = e[i].to;
if (dep[y] == INF) {
dep[y] = dep[x] + 1;
q.push(y);
if (y == T)
return true;
}
}
}
return false;
}
int dfs(int x, int lim) {
if (x == T || lim == 0)
return lim;
int tmp, flow = 0;
for (int &i = cur[x]; i; i = e[i].next)
if (dep[e[i].to] == dep[x] + 1 && (tmp = dfs(e[i].to, min(lim, e[i].f)))) {
e[i].f -= tmp;
e[i ^ 1].f += tmp;
flow += tmp;
lim -= tmp;
if (lim == 0)
break;
}
return flow;
}
int dinic() {
int ans = 0;
while (bfs())
ans += dfs(S, INF);
return ans;
}
} Maxflow;
最小费用最大流
struct Net {
struct edge {
int to, next, f, w;
} e[M << 1];
int cnt, head[N];
void Add(int x, int y, int f, int w) {
e[++ cnt] = (edge) {y, head[x], f, w}; head[x] = cnt;
}
void add(int x, int y, int f, int w) {
Add(x, y, f, w); Add(y, x, 0, -w);
}
int S, T;
int flow, cost;
int vis[N], dis[N], cur[N];
bool spfa() {
queue<int> q;
memset(vis, 0, sizeof(vis));
memset(dis, 0x3f, sizeof(dis));
memcpy(cur, head, sizeof(cur));
q.push(S); dis[S] = 0;
while(q.size()) {
int x = q.front(); q.pop();
vis[x] = 0;
for(int i = head[x]; i; i = e[i].next) {
int y = e[i].to;
if(dis[y] > dis[x] + e[i].w && e[i].f) {
dis[y] = dis[x] + e[i].w;
if(! vis[y]) { vis[y] = 1; q.push(y); }
}
}
}
return dis[T] != 0x3f3f3f3f;
}
int dfs(int x, int f) {
if(x == T) {
flow += f;
cost += f * dis[T];
return f;
}
vis[x] = 1;
int res = 0, tmp;
for(int &i = cur[x]; i; i = e[i].next) {
int y = e[i].to;
if(vis[y]) continue;
if(e[i].f && dis[y] == dis[x] + e[i].w) {
tmp = dfs(y, min(f - res, e[i].f));
e[i].f -= tmp;
e[i ^ 1].f += tmp;
res += tmp;
if(res == f) break;
}
}
return res;
}
void MCMF() {
while(spfa()) dfs(S, 0x3f3f3f3f);
cout << flow << ' ' << cost << endl;
}
Net() {
cnt = 1;
}
} net;
二分图最大匹配
bool dfs(int x){
for(int i = 1; i <= n_g; i++)
if(e[x][i] && !used[i]){
used[i] = 1;
if(!boy[i] || dfs(boy[i])){
boy[i] = x; girl[x] = i; return true;
}
}
return false;
}
最小树形图
咕咕
FHQ treap
#define ls(x) t[x].ch[0]
#define rs(x) t[x].ch[1]
struct Node {
int ch[2];
int rnd, val;
} t[N];
....
LCT
int fa[N], ch[N][2];
int val[N], sum[N];
int rv[N];
#define ls(x) ch[x][0]
#define rs(x) ch[x][1]
int get(int x) { return x == rs(fa[x]); }
inline int nroot(int x) { return x == ls(fa[x]) || x == rs(fa[x]); }
void update(int x) {
sum[x] = sum[ls(x)] ^ sum[rs(x)] ^ val[x];
}
void reverse(int x) {
swap(ls(x), rs(x));
rv[x] ^= 1;
}
void pushdown(int x) {
if(rv[x]) {
if(ls(x)) reverse(ls(x));
if(rs(x)) reverse(rs(x));
rv[x] ^= 1;
}
}
void pushall(int x) {
if(nroot(x)) pushall(fa[x]); pushdown(x);
}
void rotate(int x) {
int y = fa[x], z = fa[y], k = get(x), w = ch[x][! k];
if(nroot(y)) ch[z][get(y)] = x; ch[x][! k] = y; ch[y][k] = w;
if(w) fa[w] = y; fa[y] = x; fa[x] = z; update(y);
}
void splay(int x) {
pushall(x);
while(nroot(x)) {
if(nroot(fa[x]))
rotate(get(x) ^ get(fa[x]) ? x : fa[x]);
rotate(x);
} update(x);
}
void access(int x) {
for(R int y = 0; x; x = fa[y = x]) {
splay(x); rs(x) = y; update(x);
}
}
void makeroot(int x) {
access(x); splay(x); reverse(x);
}
void split(int x, int y) {
makeroot(x); access(y); splay(y);
}
int findrt(int x) {
access(x); splay(x);
while(ls(x)) x = ls(x);
splay(x); return x;
}
void link(int x, int y) {
makeroot(x); access(y); splay(y);
if(fa[x] != y) fa[x] = y;
}
void cut(int x, int y) {
makeroot(x);
if(findrt(y) == x && fa[y] == x && ls(y) == 0) {
fa[y] = rs(x) = 0;
update(x);
}
}
ODT
struct Node {
int l, r;
mutable int val;
Node(int _l = 0, int _r = 0, int _val = 0) : l(_l), r(_r), val(_val) {}
inline bool operator < (const Node &t) const {
return l < t.l;
}
} ;
std :: set<Node> st;
std :: set<Node> :: iterator split(int pos) {
auto it = st.lower_bound(Node(pos));
if(it != st.end() && it -> l == pos) return it;
it --;
int ll = it -> l, lr = it -> r, lv = it -> val;
st.erase(it);
st.insert(Node(ll, pos - 1, lv));
return st.insert(Node(pos, lr, lv)).first;
}
void solve(int l, int r) {
auto itr = split(r + 1), itl = split(l);
for(auto it = itl; it != itr; ++ it)
// do something
}
哈希表
用邻接表写
NTT
const int G = 3;
const int Gi = power(3, P - 2);
int lim, bit, rev[N];
void init(int n) {
lim = 1; bit = 0;
while(lim <= n) lim <<= 1, bit ++;
Lep (i, 0, lim) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
}
void NTT(int *A, int type) {
Lep (i, 0, lim) if(i < rev[i]) swap(A[i], A[rev[i]]);
for(int dep = 1; dep < lim; dep <<= 1) {
int Wn = power(type == 1 ? G : Gi, (P - 1) / (dep << 1));
for(int j = 0, len = (dep << 1); j < lim; j += len) {
int w = 1;
for(int k = 0; k < dep; k ++, w = 1ll * w * Wn % P) {
int x = A[j + k], y = 1ll * A[j + k + dep] * w % P;
A[j + k] = add(x, y);
A[j + k + dep] = dec(x, y);
}
}
}
if(type == -1) {
int inv = power(lim, P - 2);
Lep (i, 0, lim) A[i] = 1ll * A[i] * inv % P;
}
}
行列式
int sgn = 1;
lep (i, 1, n) {
if(A[i][i] == 0) {
lep (j, i + 1, n) if(A[j][i] != 0) { swap(A[i], A[j]); sgn = P - sgn; break; }
}
if(A[i][i] == 0) return printf("%d\n", 0), 0;
lep (j, i + 1, n) {
if(A[i][i] < A[j][i]) swap(A[i], A[j]), sgn = P - sgn;
while(A[j][i]) {
swap(A[i], A[j]);
int d = A[j][i] / A[i][i];
lep (k, i, n) A[j][k] = (A[j][k] - 1ll * A[i][k] * d % P + P) % P;
sgn = P - sgn;
}
}
}
int ans = sgn;
lep (i, 1, n) ans = 1ll * ans * A[i][i] % P;
计算几何基础
const double eps = 1e-8;
#define lt(x, y) ((x) < (y) - eps)
#define gt(x, y) ((x) > (y) + eps)
#define le(x, y) ((x) <= (y) + eps)
#define ge(x, y) ((x) >= (y) - eps)
#define eq(x, y) (le(x, y) && ge(x, y))
#define dot(x, y, z) (((y) - (x)) * ((z) - (x)))
#define cross(x, y, z) (((y) - (x)) ^ ((z) - (x)))
struct vec2 {
double x, y;
vec2 (double _x = 0, double _y = 0) : x(_x), y(_y) {}
inline vec2 operator - () const { return vec2(- x, - y); }
inline vec2 operator + (const vec2 &t) const { return vec2(x + t.x, y + t.y); }
inline vec2 operator - (const vec2 &t) const { return vec2(x - t.x, y - t.y); }
inline vec2 operator * (double k) const { return vec2(x * k, y * k); }
inline vec2 operator / (double k) const { return * this * (1.0 / k); }
inline double operator * (const vec2 & t) const { return x * t.x + y * t.y; }
inline double operator ^ (const vec2 & t) const { return x * t.y - y * t.x; }
inline double norm() { return sqrt(x * x + y * y); }
inline double norm2() { return x * x + y * y; }
inline bool operator < (const vec2 &t) const { return lt(x, t.x) || le(x, t.x) && lt(y, t.y); }
inline bool operator == (const vec2 &t) const { return eq(x, t.x) && eq(y, t.y); }
inline void rotate(db t) { * this = vec2(cos(t) * x - sin(t) * y, sin(t) * x + cos(t) * y); }
} ;
struct line {
vec2 p1, p2;
line() {}
line(vec2 _p1, vec2 _p2) : p1(_p1), p2(_p2) {}
line(double k, double b) { p1 = vec2(0, b); p2 = vec2(1000, b + k * 1000); }
inline vec2 direct() const { return p2 - p1; }
} ;
inline bool parallel(const line &a, const line &b) {
return eq(a.direct() ^ b.direct(), 0);
}
凸包
采用求一次上凸包然后求一次下凸包在拼起来实现。
\(type = 0\) 是下凸包， 反之是上凸包。
inline void graham(vector<vec2> &vec, int type) {
if(vec.size() == 0) { cerr << "Hull is empty !\n" << endl; exit(0); }
sort(vec.begin(), vec.end());
vector<vec2> rec; rec.pb(* vec.begin()); int sz = 0;
for(int i = 1; i < vec.size(); i ++) {
if(type == 0) while(sz >= 1 && le(cross(rec[sz - 1], rec[sz], vec[i]), 0)) rec.pop_back(), sz --;
else while(sz >= 1 && ge(cross(rec[sz - 1], rec[sz], vec[i]), 0)) rec.pop_back(), sz --;
rec.push_back(vec[i]); sz ++;
}
swap(vec, rec);
}
inline void graham_full(vector<vec2> &vec) {
vector<vec2> v1 = vec, v2 = vec;
graham(v1, 0); graham(v2, 1);
v1.pop_back(); for(int i = v2.size() - 1; i >= 1; i --) v1.push_back(v2[i]); swap(vec, v1);
}
凸包直径
采用旋转卡壳实现。每次就看移动指针以后面积会不会变大就好了。
inline double convDiameter(vector<vec2> vec) {
if(vec.size() == 2) { return (vec[0] - vec[1]).norm2(); }
vec.pb(vec[0]);
int j = 2, n = vec.size() - 1;
double res = 0;
for(int i = 0; i < vec.size() - 1; i ++) {
while(abs(cross(vec[i], vec[i + 1], vec[j])) < abs(cross(vec[i], vec[i + 1], vec[j + 1]))) j = (j + 1) % n;
Max(res, max((vec[i] - vec[j]).norm2(), (vec[i + 1] - vec[j]).norm2()));
}
return res;
}
直线交点
用三角形面积然后定比分点实现。
inline vec2 intersection(const line &l1, const line &l2) {
double ls = cross(l1.p1, l1.p2, l2.p1);
double rs = cross(l1.p1, l1.p2, l2.p2);
return l2.p1 + (l2.p2 - l2.p1) * ls / (ls - rs);
}
半平面交
把直线按照极角排序，平行的保留内侧在前， 注意排序判断是否在内侧的时候使用严格大于，从头开始加， 如果队头或者队尾的两条直线交点在当前直线右侧则弹出， 如果最后队尾两条直线交点在队头直线右侧则弹出队尾。
直线是逆时针存的。
inline void HPI(vector<line> &lv) {
vector<pair<line, double> > sorted(lv.size());
for(int i = 0; i < lv.size(); i ++) sorted[i].fi = lv[i], sorted[i].se = atan2(lv[i].direct().y, lv[i].direct().x);
sort(sorted.begin(), sorted.end(), [] (auto a, auto b) -> bool {
if(eq(a.se, b.se)) {
if(lt(cross(a.fi.p1, a.fi.p2, b.fi.p2), 0)) return 1;
else return 0;
}
else return a.se < b.se;
} );
for(int i = 0; i < lv.size(); i ++) lv[i] = sorted[i].fi;
deque<line> q;
q.push_back(lv[0]);
for(int i = 1; i < lv.size(); i ++) if(! parallel(lv[i], lv[i - 1])) {
while(q.size() > 1) {
vec2 p = intersection(* --q.end(), * -- -- q.end());
if(lt(cross(lv[i].p1, lv[i].p2, p), 0)) q.pop_back();
else break;
}
while(q.size() > 1) {
vec2 p = intersection(* q.begin(), * ++ q.begin());
if(lt(cross(lv[i].p1, lv[i].p2, p), 0)) q.pop_front();
else break;
}
q.push_back(lv[i]);
}
while(q.size() > 1) {
vec2 p = intersection(* --q.end(), * -- -- q.end());
if(lt(cross(q.begin() -> p1, q.begin() -> p2, p), 0)) q.pop_back();
else break;
}
lv = vector<line> (q.size());
for(int i = 0; i < q.size(); i ++) lv[i] = q[i];
}
闵可夫斯基和
要求 \(v1, v2\) 是凸包且极角排序完毕。
把所有边拿下来归并即可。注意最后弹出一个初始点。
inline vector<vec2> Minkowski(vector<vec2> v1, vector<vec2> v2) {
// v1, v2 is sorted
vector<vec2> s1(v1.size()), s2(v2.size());
for(int i = 1; i < s1.size(); i ++) s1[i - 1] = v1[i] - v1[i - 1]; s1[s1.size() - 1] = v1[0] - v1[s1.size() - 1];
for(int i = 1; i < s2.size(); i ++) s2[i - 1] = v2[i] - v2[i - 1]; s2[s2.size() - 1] = v2[0] - v2[s2.size() - 1];
vector<vec2> hull(v1.size() + v2.size() + 1);
int p1 = 0, p2 = 0, cnt = 0;
hull[cnt ++] = v1[0] + v2[0];
while(p1 < s1.size() && p2 < s2.size()) {
hull[cnt] = hull[cnt - 1] + (ge(s1[p1] ^ s2[p2], 0) ? s1[p1 ++] : s2[p2 ++]);
cnt ++;
}
while(p1 < s1.size()) hull[cnt] = hull[cnt - 1] + s1[p1 ++], cnt ++;
while(p2 < s2.size()) hull[cnt] = hull[cnt - 1] + s2[p2 ++], cnt ++;
hull.pop_back();
return hull;
}
辛普森积分
积分公式 ： \((R - L) \times (f(L) + f(R) + 4f(M)) \times \frac{1}{6}\)
db calc(db L, db R) {
db mid = (L + R) / 2;
return (f(L) + f(R) + 4 * f(mid)) / 6 * (R - L);
}
db simpson(db L, db R, int dep) {
if(abs(L - R) <= 1e-9 && dep >= 4) return 0;
db mid = (L + R) / 2;
if(abs(calc(L, R) - calc(L, mid) - calc(mid, R)) <= 1e-9 && dep >= 4) return calc(L, R);
return simpson(L, mid, dep + 1) + simpson(mid, R, dep + 1);
}
圆圆交点
把一个圆放到圆心， 另外一个旋转到水平， 解方程以后转回去。
平面图转对偶图
使用最小左转法。把所有边拆为两条有向边， 然后在每个点上极角排序， 每次一个边在终点找反向边然后看极角序在反向边前面的一条边记为 \(nxt\) ， 然后跳 \(nxt\) 编号， 根据面积正负判断是否是无界域。
void build() {
for(int i = 1; i <= n; i ++) sort(e[i].begin(), e[i].end());
for(int i = 2; i <= tot; i ++) {
int y = E[i].y;
auto it = lower_bound(e[y].begin(), e[y].end(), E[i ^ 1]);
if(it == e[y].begin()) it = -- e[y].end(); else -- it;
nxt[i] = it -> id;
}
for(int i = 2; i <= tot; i ++) if(! pos[i]) {
pos[i] = pos[nxt[i]] = ++ cnt;
int u = E[i].x;
for(int j = nxt[i]; E[j].y != E[i].x; j = nxt[j], pos[j] = cnt)
s[cnt] += ( (nod[E[j].x] - nod[u]) ^ (nod[E[j].y] - nod[u]) );
if(s[cnt] < 1e-9) rt = cnt;
}
for(int i = 2; i <= tot; i ++) {
hv[pos[i]].pb(i);
g[pos[i]].pb( {E[i].w, pos[i ^ 1]} ) ;
}
}
处理圆和多边形包含关系
使用 \(set\) 直接维护， 具体见 [SCOI2012]Blinker 的噩梦 。