---
# System prepended metadata

title: 2021 高階競技程式設計 14組 魔法書

---

# 競程魔法書

## C++模板
**vimrc**
``` vim=
set encoding=utf-8
set nu
syntax on
set ai
set shiftwidth=4
set tabstop=4
set expandtab
set smartindent
set backspace=2
set ruler
set ic
set number
set mouse=a
set hlsearch

inoremap {<CR> {<CR>}<Esc>ko
inoremap " ""<Esc>i
inoremap ' ''<Esc>i
```
**模板**
```cpp=
#include <bits/stdc++.h>

using namespace std;
#define INITIO()    ios_base::sync_with_stdio(false);cin.tie(NULL);cout<<setprecision(numeric_limits<double>::max_digits10+1);
#define MEM(i,j)    memset(i,j,sizeof i);
#define X           first
#define Y           second
#define pb          push_back
#define mp          make_pair
#define ALL(a)      a.begin(),a.end()
typedef             long long ll;
typedef             pair<int ,int> pii;
const double Eps  = 1e-8;
const ll Inf      = 1e18+7;
const int Mod     = 1e9+7;

int main (){
    INITIO();
    int tt;
    cin >> tt;
    while (tt--) {
        
    }
    return 0;
}
```
## STL

### Map STL

```cpp=
map<T, T> m;
unordered_map<T, T> umap; //c++11

//放東西
m.insert();
m[] = ;

//找東西
auto iter = m.find();
if(iter != m.end())
    cout<<iter->second<<endl;
```

### Queue STL

```cpp=
queue<T> q;
priority_queue<T> pq;

q.push();
q.pop();
q.front();
q.back();
```

### 排列組合 C(n,r)
>stackoverflow 來的，所以好像有overflow可能
```cpp=
long long C(int n, int r) {
    if(r > n - r) r = n - r;
    long long ans = 1;
    int i;

    for(i = 1; i <= r; i++) {
        ans *= n - r + i;
        ans /= i;
    }

    return ans;
}
```

### Sort & Cmp

>**sort** 預設升冪 , 改降冪使用 **greater<T>()**

```cpp=
bool myfunction (int a,int b) {
    return (a<b);
}
sort(All(s), myfunction);

//用lambda，A是某種struct
sort( values.begin( ), values.end( ), [ ]( const A& l, const A& r )
{
   return l.key < r.key;
});
```

## 快速冪
>一般的
```cpp=
ll Power(ll n,ll time){
    ll a = 1;
    while(time){
        if(time & 1) a = a*n % Mod;
        n = n*n % Mod;
        time>>=1;
    }
    return a;
}
```
### 矩陣快速冪
>費式數列 , 比**dp**解更快
```cpp=
ll Pow(ll t){
    t--;
    vector<ll> a = {1,1,1,0},b = {1,1,1,0};

    while(t){
        vector<ll> c = a,d = b;

        if(t&(ll)1){
            for(int i=0;i<2;i++)
                for(int j=0;j<2;j++)
                    a[2*i+j] = c[2*i]*b[j]+c[2*i+1]*b[2+j];
        }

        for(int i=0;i<2;i++)
            for(int j=0;j<2;j++)
                b[2*i+j] = d[2*i]*d[j]+d[2*i+1]*d[2+j];
        t>>=(ll)1;
    }
    return a[0];
}
```

## 質數
>線性篩法
```cpp=
ll n;
cin >> n;
bool arr[n+1];
MEM(arr,true);
for(int i=2;i<=sqrt(n);++i)
    if(arr[i]) 
        for(int j = i*i ; j<=n ;j+=i) 
            arr[j] = false;
    
for(int k=1;k<=n;k++)
    if(arr[k]) cout << k <<endl;
```

## 最小生成樹 MST
>求長度應用例
>Kruskal演算法 , 路徑壓縮沒啟發式合併
```cpp=
const int N = 5e3+5,M = 2e5+5;
int boss[N]={0};

struct Edge{
    int a;
    int b;
    int val;
} e[M];

bool cmp(Edge a,Edge b){
    return a.val<b.val;
}

int find(int idx){
    if(boss[idx]==idx) return idx;
    return boss[idx] = find(boss[idx]);
}
```
``` cpp=
int n,m,length=0,en=0;
cin>>n>>m;
    
for(int i=1;i<n+1;i++){
    boss[i]=i;
}
for(auto &it : e){
    cin>>it.a>>it.b>>it.val;
}

sort(e,e+m,cmp);

for(int j=0;j<m;j++){
    if(find(e[j].a) == find(e[j].b)) continue;
    boss[find(e[j].a)] = find(e[j].b);
    length+=e[j].val;
    en++;
}

if(en<n-1) cout << "orz"<<endl;
else cout <<length<<endl;

```
### Prim’s algorithm 
```cpp=
#define MAX 10000

struct node {
    int id,cost;
};
vector<node> E[MAX];

struct edge {
    int u,v,cost;
    bool operator<(const edge& y) const { return cost > y.cost; }
};
priority_queue<edge> Q;

bool vis[MAX];

//use
{   
    // init
    for(int i = 0; i < n; i++)     
        vis[i] = false;
    while(m--) {
        cin >> u >> v >> cost;
        E[u].push_back({v,cost});
        E[v].push_back({u,cost});
    }

    // Prim's Algorithm
    Q.push({0,0,0}); 
    while(!Q.empty()) {
        edge e = Q.top();
        Q.pop();
        u = e.v;
        if(!vis[u]) {
            ans += e.cost;
            for(auto v : E[u]) {
                if(!vis[v.id])
                    Q.push({u,v.id,v.cost});
            }
        }
        vis[u] = true;
    }
    cout << ans << '\n';
}
```

## 最短路

**Floyd**
```cpp= 
ll dp[Maxn][Maxn];

ll n,k,q;
cin >> n >> k >> q;

for(int i = 0;i<n ;i++){
    for(int j =0;j<n ;j++){
        if(i == j) dp[i][j] = 0;
        else dp[i][j] = Inf;
    }
}

while(k--){
    ll a,b,w;
    cin >> a >> b >> w;
    dp[a][b] = w;
}

for(int i =0;i<n;i++)
    for(int k=0;k<n;k++)
        for(int j = k;j <n;j++)
            dp[i][j] = min(dp[i][j],dp[i][k]+dp[k][j]);

```
## Dijkstra
單點搜尋
$O((E+N)log(N))$
```
struct Node{
    ll n,w;
    bool operator<(const node &lhs) const {
        return lhs.w < w; //min heap
    }
};
```
初始:
```
MEM(dist,Inf); //不連通設Inf
MEM(vis,false); //沒有拜訪
priority_queue<Node> q;
q.push((Node){1,dist[1]=0}); //看起點更改
```
循環:
```
while(q.size()){
    node now = q.top();q.pop();
    vis[now.n] = true;
    
    for(auto nxt : e[now.n]){
        if(vis[nxt.n]) continue;
        const int &updater = dist[now.n] + nxt.w; 
        if(update < dist[nxt.n]) {
            dist[nxt.n] = update; //鬆弛
            q.push((Node){nxt.n,dist[nxt.n]});
        }
    }
}
```
## Bellman-Ford+檢查負環
單點搜尋
$O(NE)$
```
int w[N][N];
int d[N];
int parent[N];

void bellman_ford(int source){
    for(int i=0;i<n;i++) d[i]=Inf;
    d[source] = 0;
    parent[source] = source;
    
    for(int i=0;i<n-1;i++)
        for(int a=0;a<n;++a)
            for(int b=0;b<n;b++)
                if(d[a] != Inf && w[a][b] != Inf)
                    if(d[a]+w[a][b] < d[b]){
                        d[b] = d[a] + w[a][b];
                        parent[b] = a;
                    }
}
```
負環檢查
* 如果還能鬆弛代表有負環
```
bool is_negative(){
    for(int a=0;a<n;++a)
        for(int b=0;b<n;++b)
            if(d[a] != Inf && w[a][b] != Inf)
                if(d[a]+w[a][b] < d[b])
                    return true;
        
    
    return false;
}
```
遞迴版本 視情況修改 (一開始刻我以為是 dijkstra,但能處理處理負環)
```
ll bell(int start,int end,int val){
    ll ans = Inf;
    
    if(dis[start]<=val) return Inf;
    else dis[start] = val;
    if(start == end) return val;

    for(int i=0;i<E[start].size();i++) 
        ans = min(ans,bell(E[start][i].X,end,E[start][i].Y+val));

    return ans;
}
```
## Floyd-Warshell + 檢查負環
$O(N^3)$
多源搜尋
```
for(int i = 0;i<n ;i++){
    for(int j =0;j<n ;j++){
        if(i == j) dp[i][j] = 0;
        else dp[i][j] = Inf;
    }
}
...
for(int k =0;k<n;k++)
    for(int i=0;i<n;i++)
        for(int j = 0;j <n;j++)
            dp[i][j] = min(dp[i][j],dp[i][k]+dp[k][j]);
```
負環
```
bool is_negative(){
    for (int i=0; i<n; i++)
        if (d[i][i] < 0) //自己到自己必須是0否則有負環
            return true;
    return false;
}
```
## 關節點
```cpp=
#define MAX 105
int low[MAX],dep[MAX];
bool ap[MAX];
int d = 0;
vector<int> adj[MAX];
    
void dfs(int u, int w == -1) {
    int child{0};
    low[u] = dep[u] = ++d;
    for (auto& v : adj[u]) {
        if (v == w) continue; // 剛剛過來的點
        if (!dep[v]) { // 尚未拜訪的點
            ++child;
            dfs(v, u);
            low[u] = min(low[u], low[v]); // 更新 back edge (u 先往 child 走可能會有更好的 back edge)
            if (low[v] >= dep[u]) // v 走 back edge 無法到 u 的祖先
                ap[u] = true;
        } 
        else  // 已經拜訪過得的點
            low[u] = min(low[u], dep[v]); // 更新 back edge (u 可能可以直接走到祖先)
    }
    if (w == -1) 
        ap[u] = (child > 1); // 檢查 root 是否為關節點 (子樹>1)
    --d;
}
```
## 最低共同祖先 
```cpp=
// pointer jumping technique (binary lifting)
#define MAX 30001
vector<int> node[MAX];
int DPWay[MAX][15]; /*點往上走 2^i level 所到的點，2^15=32768 > 30001*/
int treeLevel[MAX]; /*點高度*/
int treePath[MAX]; /*紀錄當前各level的點*/

void dfs(int now, int pre, int level) {
    treeLevel[now] = level;
    treePath[level] = now;
    for(int i = 1; (1 << i) < level; i++)
        DPWay[now][i] = treePath[level - (1 << i)];     
    for(int i = 0; i < node[now].size(); i++) {
        if(node[now][i].first == pre)
            continue;
        DPWay[node[now][i].first][0] = now;
        dfs(node[now][i].first, now, level + 1);
    }
}

int LCA(int A,int B) {
    int levelA = treeLevel[A], levelB = treeLevel[B];
    for(int i = (int)log2(abs(levelA - levelB)); i >= 0; --i) {
        if(abs(levelA - levelB) >= (1 << i)) {
            if(levelA > levelB) {
                A = DPWay[A][i];
                levelA -= (1 << i);
            }
            else {
                B = DPWay[B][i];
                levelB -= (1 << i);
            }
        }
    }
    // A B 走到同一個高度 (假設樹高小於32，我們先看距離有沒有超過16，有的話就往上走16步；剩下的距離一定小於16，我們再看它有沒有超過8；以此類推，4,2,1，就到終點了。)
    
    if (A == B) 
        return A;
    for (int i = (int)log2(levelA); i >= 0; --i)
        if (DPWay[A][i] != DPWay[B][i]) {
            A = DPWay[A][i];
            B = DPWay[B][i];
        }
    // u, v 一起走到 lca(u, v) 的下方

    return DPWay[A][0];
    // 回傳 lca(u, v)
}
```
    
## 線段樹
```cpp=
template <typename value_t>
class SGT {
    public:
    int n;
    vector<value_t> t;
    vector<optional<value_t>> lz;
    int left(int tv) { return tv + 1; }
    int right(int tv,int tl, int tm) { return tv + (tm-tl) * 2 ; }
    void build(const vector<value_t>& v, int tv, int tl, int tr) {
        if(tr - tl == 1) {
            t[tv] = v[tl];
        }
        else {
            int tm = (tl+tr) / 2;
            build(v,left(tv),tl,tm);
            build(v,right(tv,tl,tm),tm,tr);
            t[tv] = merge(t[left(tv)],t[right(tv,tl,tm)]);
        }
    }
    value_t merge(const value_t& x, const value_t& y) {
        return x + y;
    }
    void update(int tv, int len, const value_t& x) {
        if(!lz[tv])
            lz[tv] = x;
        else
            lz[tv] = lz[tv].value() + x;
        t[tv] += len * x;
    }
    void push(int tv, int tl, int tr) {
        if(!lz[tv])
            return;
        int tm = (tl+tr) / 2;
        update(left(tv),tm-tl,lz[tv].value());
        update(right(tv,tl,tm),tr-tm,lz[tv].value());
        lz[tv].reset();
    }
    void modify(int l, int r, const value_t& x, int tv, int tl, int tr) {
        if(l == tl && r == tr)
            update(tv,tr-tl,x);
        else {
            push(tv,tl,tr);
            int tm = (tl+tr) / 2;
            if(r <= tm)
                modify(l,r,x,left(tv),tl,tm);
            else if(l >= tm)
                modify(l,r,x,right(tv,tl,tm),tm,tr);
            else {
                modify(l,tm,x,left(tv),tl,tm);
                modify(tm,r,x,right(tv,tl,tm),tm,tr);
            }
            t[tv] = merge(t[left(tv)],t[right(tv,tl,tm)]);
        }
    }
    value_t query(int l, int r, int tv, int tl, int tr) {
        if(l == tl && r == tr)
            return t[tv];
        push(tv,tl,tr);
        int tm = (tl+tr) / 2;
        if(r <= tm)
            return query(l,r,left(tv),tl,tm);
        else if(l >= tm)
            return query(l,r,right(tv,tl,tm),tm,tr);
        else
            return merge(query(l,tm,left(tv),tl,tm),query(tm,r,right(tv,tl,tm),tm,tr)); 
    }
public:
    SGT(const vector<value_t>& v) : n{v.size()}, t(2*n), lz(2*n) {
        build(v,1,0,n);
    }
    void modify(int l, int r, const value_t& x) { modify(l,r,x,1,0,n); }
    value_t query(int l,int r) { return query(l,r,1,0,n); }
};
    
// a[0] a[1] a[2]
//   1    2    3
// l = 1     r = 3   k = 5
// query(l-1,r)  modify(l-1,r,k)

```
### 懶人標記
```cpp=
//新值是value，修改區間是[l,r]，節點v的範圍是[L,R]
void modify(type value, int l, int r, int L, int R, vertex v){
    if(l == L && r == R){
        打懶標在v上;
        return;
    }
    int M = (L + R) / 2;
    if(r <= M) modify(value, l, r, L, M, v的左子節點);
    else if(l > M) modify(value, l, r, M + 1, R, v的右子節點);
    else{
        modify(value, l, M, L, M, v的左子節點);
        modify(value, M + 1, r, M + 1, R, v的右子節點);
    }
    用兩個子節點的答案更新v的答案;
}
```
修改
```cpp=
//查詢範圍是[l,r]，節點v的區間範圍是[L,R]
type query(int l, int r, int L, int R, vertex v){
    if(l == L && r == R)
        return ans of v;
    把節點v的懶標下推; //<----
    int M = (L + R) / 2;
    if(r <= M)
        return query(l, r, L, M, v的左子節點);
    else if(l > M)
        return query(l, r, M + 1, R, v的右子節點);
    else{
        return 用query(l, M, L, M, v的左子節點)
            和query(M + 1, r, M + 1, R, v的右子節點)
            得出區間[l,r]的答案;
    }
}
```
    
## Computation Geometry  
```cpp=
struct Point{
	double x, y;
	Point(double x = 0, double y = 0) : x(x), y(y){}
	Point operator+(const Point &b)const{  // vector addition
		return Point(x + b.x, y + b.y);
	}
	Point operator-(const Point &b)const{  // vector subtraction
		return Point(x - b.x, y - b.y);
	}
	Point operator*(double d)const{        // scale (multiplication)
		return Point(x * d, y * d);
	}
	Point operator/(double d)const{        // scale (division)
		return Point(x / d, y / d);
	}
	double dot(const Point &b)const{       // vector dot
		return x * b.x + y * b.y;
	}
	double cross(const Point &b)const{     // vector cross
		return x * b.y - y * b.x;
	}
};

// premise: three points are collinear
bool btw(const Point &p1, const Point &p2, const Point &p3){
	return (p1 - p3).dot(p2 - p3) <= 0;
}

bool collinearity(const Point &p1, const Point &p2, const Point &p3){
	return (p1 - p3).cross(p2 - p3) == 0;
}

bool pointOnSegment(const Point &p1, const Point &p2, const Point &p3){
	return collinearity(p1, p2, p3) && btw(p1, p2, p3);
}

int ori(const Point &p1, const Point &p2, const Point &p3){
	double d = (p2 - p1).cross(p3 - p1);
	if(d == 0)
		return 0;
	return d > 0 ? 1 : -1;
}

bool seg_intersect(const Point &p1, const Point &p2, const Point &p3, const Point &p4){
	int a123 = ori(p1, p2, p3);
	int a124 = ori(p1, p2, p4);
	int a341 = ori(p3, p4, p1);
	int a342 = ori(p3, p4, p2);
	if(a123 == 0 && a124 == 0)
		return btw(p1, p2, p3) || btw(p1, p2, p4) || btw(p3, p4, p1) || btw(p3, p4, p2);
	return a123 * a124 <= 0 && a341 * a342 <= 0;
}

Point intersect(const Point &p1, const Point &p2, const Point &p3, const Point &p4){
	int a123 = (p2 - p1).cross(p3 - p1);
	int a124 = (p2 - p1).cross(p4 - p1);
	return (p4 * a123 - p3 * a124) / (a123 - a124);
}
```
    

## 逆序數對
```cpp=
int merge_sort(vector<int> &a,int l,int r){
    if(l==r) return 0;
    int mid = (l+r)/2;

    int lans = merge_sort(a,l,mid);
    int rans = merge_sort(a,mid+1,r);

    vector<int> b,c;

    for(int i = l;i<=mid;i++) b.pb(a[i]);
    for(int i = mid+1;i<=r;i++) c.pb(a[i]);

    int i=0,j=0,ans=0;
    for(int k=l;k<=r;k++){
        int Min = 0;
        if(i == b.size()) Min = c[j++];
        else if(j == c.size()) Min = b[i++];
        else if(b[i] <= c[j]) Min = b[i++];
        else {
            Min = c[j++];
            ans += mid - i - l + 1;
        }
        a[k] = Min;
    }
    return lans+rans+ans; 
}
```
    
## DP
### 01背包問題
```cpp=
int v[n];
int c[n];
int dp[101]={0};
int ans=0;
for(int i=0;i<n;i++){
    cin>> v[i] >> c[i];
    for(int j=100;j-v[i]>=0;j--){
        dp[j] = max(dp[j],dp[j-v[i]]+c[i]);
        ans = max(ans,dp[j]);
    }
}
cout << ans <<endl;
```

## GCD
>GCD求解
```cpp=
int gcd(int a, int b)
{ return a? gcd(b%a, a) : b;}
```

### 貝祖等式 (Bezout’s Thm)
>對整數 **a,b** , 存在整數 **x,y** 使得 **ax+by = gcd(a,b)**

### 擴展歐基里德
令 **g = gcd(0,g) = gcd(a,b)**
```cpp=
int gcd(int a, int b, int &x, int &y) {
  if(!a) {
    x = 0; // x 為任意整數 
    y = 1;
    return b;
  }

  int g = gcd(b%a, a, x, y);
  int xp = y - b/a * x, yp = x; // p := previous
  x = xp, y = yp;

  return g;
}
```
**exgcd** 可求模反元素 , 但僅於 b 跟 m 互質

## Karatsuba algorithm
```cpp=
// 不要爆掉，感覺都是對的
int getLength(long long value) {
	int counter = 0;
	while (value != 0) {
		counter++;
		value /= 10;
	}
	return counter;
}

long long karatsuba(long long x, long long y) {
	int xLength = getLength(x);
	int yLength = getLength(y);

	// the bigger of the two lengths
	int N = (int)(fmax(xLength, yLength));

	// if the max length is small it's faster to just flat out multiply the two nums
	if (N < 10)
		return x * y;

	//max length divided and rounded up
	N = (N/2) + (N%2);

	long long multiplier = pow(10, N);

	long long b = x/multiplier;
	long long a = x - (b * multiplier);
	long long d = y / multiplier;
	long long c = y - (d * N);

	long long z0 = karatsuba(a,c);
	long long z1 = karatsuba(a + b, c + d);
	long long z2 = karatsuba(b, d);


    return z0 + ((z1 - z0 - z2) * multiplier) + (z2 * (long long)(pow(10, 2 * N)));

}
               
//>
```
## 流
```cpp=
struct MaxFlow{
	struct edge{
		int to, cap, flow, rev;
	};
	vector<edge> v[MAXN];
	int s, t, dis[MAXN], cur[MAXN];
	int dfs(int u, int cap){
		if(u == t || !cap)
			return cap;
		for(int &i = cur[u]; i < v[u].size(); i++){
			edge &e = v[u][i];
			if(dis[e.to] == dis[u] + 1 && e.flow != e.cap){
				int df = dfs(e.to, min(e.cap - e.flow, cap));
				if(df){
					e.flow += df;
					v[e.to][e.rev].flow -= df;
					return df;
				}
			}
		}
		dis[u] = -1;
		return 0;
	}
	bool bfs(){
		memset(dis, -1, sizeof(dis));
		queue<int> q;
		q.push(s), dis[s] = 0;
		while(!q.empty()){
			int tmp = q.front();
			q.pop();
			for(auto &u : v[tmp]){
				if(!~dis[u.to] && u.flow != u.cap){
					q.push(u.to);
					dis[u.to] = dis[tmp] + 1;
				}
			}
		}
		return dis[t] != -1;
	}
	int maxflow(int _s, int _t){
		s = _s, t = _t;
		int flow = 0, df;
		while(bfs()){
			memset(cur, 0, sizeof(cur));
			while(df = dfs(s, INF)){
				flow += df;
			}
		}
		return flow;
	}
	void addedge(int st, int ed, int cap){
		v[st].push_back(edge{ed, cap, 0, v[ed].size()});
		v[ed].push_back(edge{st, 0, 0, v[st].size() - 1});
	}
};
```
## Matching
### 二分圖最大匹配 
>也可處理二分圖最小點覆蓋問題，而最小點覆蓋的補集就是一個最大點獨立集。
>>點覆蓋是圖的一個點集合，圖上所有邊都至少有一個端點在集合之中。
點獨立集是圖的一個點集合，集合中的點兩兩不相鄰。

>最小邊覆蓋之邊數等於點數扣掉最大匹配之邊數，而最大邊獨立集就是匹配問題。
>>最小邊覆蓋即為所有邊覆蓋中邊數最少的。
邊獨立集是圖的一個邊集合，且集合中的邊兩兩不相鄰。

>又，一側接上源點，一側接上匯點，即可利用網路流來解決最大二分匹配問題、最大（小）權二分匹配問題。
```cpp=
class BipartiteMatching {
    vector<bool> vis;
    vector<int> mx{}, my{};
    bool dfs(const vector<vector<int>>& adj, int x) {
        for (auto& y : adj[x]) {
            if (vis[y]) continue;
            vis[y] = true;
            if (my[y] == -1 || dfs(adj, my[y])) {
                mx[x] = y, my[y] = x;
                return true;
            }
        }
        return false;
    }
public:
    pair<vector<int>, vector<int>> operator()(const vector<vector<int>>& adj, int ny) {
        vis.resize(ny), mx.assign(adj.size(), -1), my.assign(ny, -1);
        int c{0};
        for (int x{0}; x < adj.size(); ++x) {
            if (mx[x] != -1) continue;
            fill(vis.begin(), vis.end(), false);
            if (dfs(adj, x)) ++c;
        }
        return {move(mx), move(my)};
    }
};
```
```cpp=
// Hopcroft–Karp algorithm - O(E√V)
class BipartiteMatching {
    vector<vector<int>> adjy;
    vector<int> dx, dy, mx, my;
    bool dfs(const vector<vector<int>>& adjx, int y) {
        for (auto& x : adjy[y]) {
            if (dx[x] + 1 != dy[y]) continue;
            if (mx[x] == -1 || dfs(adjx, mx[x])) {
                my[y] = x, mx[x] = y;
                dy[y] = -1; // remove from bfs forest
                return true;
            }
        }
        dy[y] = -1; // remove from bfs forest
        return false;
    }
    bool bfs(const vector<vector<int>>& adjx) {
        fill(dx.begin(), dx.end(), -1);
        fill(dy.begin(), dy.end(), -1);
        queue<int> qu{}, qu2{};
        for (int x{0}; x < adjx.size(); ++x) // use all unmatched x's as roots
            if (mx[x] == -1) qu.push(x), dx[x] = 0;
 
        bool found{false};
        while (!qu.empty() && !found) { // stop at the level found
            while (!qu.empty()) {
                auto x{qu.front()}; qu.pop();
                for (auto& y : adjx[x])
                    if (dy[y] == -1 && my[y] != x) {
                        dy[y] = dx[x] + 1;
                        if (my[y] == -1) found = true;
                        else qu2.push(my[y]), dx[my[y]] = dy[y] + 1;
                    }
            }
            qu.swap(qu2);
        }
        return found;
    }
public:
    pair<vector<int>, vector<int>> operator()(const vector<vector<int>>& adjx, int ny) {
        adjy.assign(ny, {});
        for (int x{0}; x < adjx.size(); ++x)
            for (auto& y : adjx[x]) adjy[y].push_back(x);
        dx.resize(adjx.size()), dy.resize(ny), mx.assign(adjx.size(), -1), my.assign(ny, -1);

        int c{0};
        while (bfs(adjx)) {
            for (int y{0}; y < ny; ++y)
                if (my[y] == -1 && dy[y] != -1) c += dfs(adjx, y);
        }
        return {move(mx), move(my)};
    }
};
```
### 一般圖最大匹配
```cpp=
// Edmonds' Blossom Algorithm - O(VEa(V))
class Matching {
    int label{0};
    vector<int> p{}, d{}, a{}, c1{}, c2{}, m{};
    // (alternating tree), (unvisited: -1, even: 0, odd: 1), (auxiliary for lca), (cross edge)
    /* DSU */
    vector<int> dsu_p{}, dsu_sz{}, dsu_b{};
    void dsu_reset() {
        fill(dsu_p.begin(), dsu_p.end(), -1);
        fill(dsu_sz.begin(), dsu_sz.end(), 1);
        iota(dsu_b.begin(), dsu_b.end(), 0);
    }
    int find(int x) { // collapsing find
        return dsu_p[x] == -1 ? x : dsu_p[x] = find(dsu_p[x]);
    }
    int base(int x) { return dsu_b[find(x)]; }
    void contract(int x, int y) { // weighted union
        auto rx{find(x)}, ry{find(y)};
        if (rx == ry) return ;
        auto b{dsu_b[rx]}; // treat x's base as new base
        if (dsu_sz[rx] < dsu_sz[ry]) swap(rx, ry);
        dsu_p[ry] = rx, dsu_sz[rx] += dsu_sz[ry], dsu_b[rx] = b;
    }
    /*******/
    queue<int> qu{}; // only even vertices
    void handle_one_side(int x, int y, int b) {
        for (int v{base(x)}; v != b; v = base(p[v])) {
            c1[v] = x, c2[v] = y, contract(b, v);
            if (d[v] == 1) qu.push(v); // odd vertex -> even vertex
        }
    }
    int lca(int u, int v) {
        ++label; // use next number in order not to clear a
        while (true) {
            if (a[u] == label) return u;
            a[u] = label;
            if (p[u] != -1) u = base(p[u]);
            swap(u, v);
        }
    }
    void augment(int r, int y) {
        if (r == y) return ;
        if (d[y] == 0) { // even vertex -> odd vertex
            // check d[y] == 0 instead of d[p[y]] == 1, so (m[y], y) is in the blossom
            augment(m[y], m[c1[y]]);
            augment(r, m[c2[y]]);
            m[c1[y]] = c2[y], m[c2[y]] = c1[y];
        } else {
            int x{p[y]};
            augment(r, m[x]);
            m[x] = y, m[y] = x;
        }
    }
    bool bfs(const vector<vector<int>>& adj, int r) {
        dsu_reset();
        fill(d.begin(), d.end(), -1);

        qu = {}, qu.push(r), d[r] = 0;
        while (!qu.empty()) {
            int x{qu.front()}; qu.pop();
            for (auto& y : adj[x]) {
                if (base(x) == base(y)) continue;
                if (d[y] == -1) {
                    p[y] = x, d[y] = 1;
                    if (m[y] == -1) { // augmenting path
                        augment(-1, y);
                        return true;
                    } else {
                        p[m[y]] = y, d[m[y]] = 0;
                        qu.push(m[y]);
                    }
                } else if (d[base(y)] == 0) { // blossom
                    int b{lca(base(x), base(y))};
                    handle_one_side(x, y, b);
                    handle_one_side(y, x, b);
                }
            }
        }
        return false;
    }
public:
    vector<int> operator()(const vector<vector<int>>& adj) {
        label = 0, p.assign(adj.size(), -1), d.resize(adj.size()), a.resize(adj.size()), c1.resize(adj.size()), c2.resize(adj.size()), m.assign(adj.size(), -1);
        dsu_p.resize(adj.size()), dsu_sz.resize(adj.size()), dsu_b.resize(adj.size());

        int c{0};
        /*
        // find random greey matching first
        vector<int> X(n); iota(X.begin(), X.end(), 0);
        random_shuffle(X.begin(), X.end());
        for (auto& x : X) {
            if (m[x] != -1) continue;
            for (auto& y : adj[x])
                if (m[y] == -1) {
                    m[x] = y, m[y] = x, ++c;
                    break;
                }
        }
         */
        for (int v{0}; v < adj.size(); ++v)
            if (m[v] == -1 && bfs(adj, v)) ++c;
        return move(m);
    }
};
```
### 二分圖最大權匹配
```cpp=
// Hungarian Algorithm - O(V^3)
class WeightedBipartiteMatching {
    int n;
    vector<long long> lx{}, ly{};
    vector<bool> vx{}, vy{};
    queue<int> qu{}; // only X's vertices
    vector<int> py{};
    vector<long long> dy{}; // smallest (lx[x] + ly[y] - w[x][y])
    vector<int> pdy{}; // & which x
    vector<int> mx{}, my{};
    void relax(const vector<vector<long long>>& w, int x) {
        for (int y{0}; y < n; ++y)
            if (lx[x] + ly[y] - w[x][y] < dy[y])
                dy[y] = lx[x] + ly[y] - w[x][y], pdy[y] = x;
    }
    void augment(int y) {
        while (y != -1) {
            int x{py[y]}, yy{mx[x]};
            mx[x] = y, my[y] = x;
            y = yy;
        }
    }
    bool bfs(const vector<vector<long long>>& w) {
        while (!qu.empty()) {
            int x{qu.front()}; qu.pop();
            for (int y{0}; y < n; ++y) {
                if (!vy[y] && lx[x] + ly[y] == w[x][y]) {
                    vy[y] = true, py[y] = x;
                    if (my[y] == -1) return augment(y), true;
                    int xx{my[y]};
                    qu.push(x), vx[xx] = true, relax(w, xx);
                }
            }
        }
        return false;
    }
    void relabel() {
        long long d{numeric_limits<long long>::max()};
        for (int y{0}; y < n; ++y) if (!vy[y]) d = min(d, dy[y]);
        for (int x{0}; x < n; ++x) if (vx[x]) lx[x] -= d;
        for (int y{0}; y < n; ++y) if (vy[y]) ly[y] += d;
        for (int y{0}; y < n; ++y) if (!vy[y]) dy[y] -= d;
    }
    bool expand(const vector<vector<long long>>& w) { // expand one layer with new equality edges
        for (int y{0}; y < n; ++y) {
            if (!vy[y] && dy[y] == 0) {
                vy[y] = true, py[y] = pdy[y];
                if (my[y] == -1) return augment(y), true;
                int xx{my[y]};
                qu.push(xx), vx[xx] = true, relax(w, xx);
            }
        }
        return false;
    }
public:
    pair<vector<int>, vector<int>> operator()(const vector<vector<long long>>& w) {
        n = w.size(), lx.assign(n, 0), ly.assign(n, 0), vx.resize(n), vy.resize(n), py.resize(n), dy.resize(n), pdy.resize(n), mx.assign(n, -1), my.assign(n, -1);
        for (int i{0}; i < n; ++i)
            for (int j{0}; j < n; ++j)
                lx[i] = max(lx[i], w[i][j]);
    
        for (int x{0}; x < n; ++x) {
            fill(vx.begin(), vx.end(), false);
            fill(vy.begin(), vy.end(), false);
            fill(dy.begin(), dy.end(), numeric_limits<long long>::max());

            qu = {}, qu.push(x), vx[x] = true, relax(w, x);
            while (!bfs(w)) {
                relabel();
                if (expand(w)) break;
            }
        }
     
        return {move(mx), move(my)};
    }
};
```
## 其他
### 一些 g++ 內建的 bit function
```cpp=
__builtin_clz()        //count leading 0s
__builtin_ctz()        //count trailing 0s
__builtin_popcount()   //count bit 1
```
### 凸包、旋轉尺
```cpp=
vector<Point> convexHull(vector<Point> v){
	int sz = 0;
	vector<Point> ans;
	sort(v.begin(), v.end(), cmp);
	for(size_t i = 0; i < v.size(); i++){
		while(sz >= 2 && (ans[sz - 1] - ans[sz - 2]).cross(v[i] - ans[sz - 2]) <= 0)
			ans.pop_back(), sz--;
		ans.push_back(v[i]), sz++;
	}
	for(int i = int(v.size()) - 2, t = sz + 1; i >= 0; i--){
		while(sz >= t && (ans[sz - 1] - ans[sz - 2]).cross(v[i] - ans[sz - 2]) <= 0)
			ans.pop_back(), sz--;
		ans.push_back(v[i]), sz++;
	}
	if(v.size() > 1)
		ans.pop_back();
	return ans;
}
    
double rotatingCaliper(const vector<Point> &v){
    double ans = 0;
    for(int i = 0; j = 0; i < v.size(); i++){
        while(abs2(v[i],v[j]) <= abs2(v[i], v[(j+1)%v.size()]))
            j = (j+1) % v.size();
        ans = max(ans, abs2(v[i], v[j]));
    }
}
```