# [复习笔记] 模板整理

#### 最大公约数

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


#### 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];
inline void _push(int x, int y, int f) {
e[++ cnt] = (edge) {
};
}
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];
void Add(int x, int y, int f, int w) {
}
void add(int x, int y, int f, int 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));
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];
}


#### 闵可夫斯基和

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


#### 辛普森积分

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


#### 平面图转对偶图

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


#### 处理圆和多边形包含关系

posted @ 2022-01-30 20:57  HN-wrp  阅读(56)  评论(2编辑  收藏  举报