University of Engineering and Technology

1 Flow and Matching

University of Engineering and Technology - TurboDB
(22-23) Notebook 1.1 Maximum Flow (Dinic)
// In case we need to find Maximum flow of network with both minimum capacity and maximum capacity,
let s* and t* be virtual source and virtual sink.
Mục lục /*
Then, each edge (u->v) with lower cap l and upper cap r will be changed in to 3 edge:

1 Flow and Matching 1 - u->v whit capacity r-l

1.1 Maximum Flow (Dinic) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 - u->t* with capacity l
1.2 Dinic (Template by Tran Khoi Nguyen) . . . . . . . . . . . . . . . . . . . . . . 2
- s*->v with capacity l
*/
1.3 Min Cut (Template by Tran Khoi Nguyen) . . . . . . . . . . . . . . . . . . . . . 2
1.4 Maximum Matching (HopCroft - Karp) . . . . . . . . . . . . . . . . . . . . . . . 3 // We need add one other edge t->s with capacity Inf
Maximum Matching (Template by Tran Khoi Nguyen) . . . . . . . . . . . . . . . . . // Maximum Flow on original graph is the Maximum Flow on new graph: s*->t*
1.5 3
1.6 Maximum Matching (O(n2 )) . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 struct Edge
1.7 Min Cost Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 {
Min Cost Flow (Template by Tran Khoi Nguyen) . . . . . . . . . . . . . . . . . . . int u, v;
1.8 5
ll c;
Edge() {}
Edge(int u, int v, ll c)
2 Geometry 6 {
2.1 Pick’s Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 this->u = u;
2.2 Smallest Circle - Emo Welzl (Contain all points) . . . . . . . . . . . . . . . . . . . 6 this->v = v;
Closest pair of points in set . . . . . . . . . . . . . . . . . . . . . . . . . . . this->c = c;
2.3 7
}
2.4 Manhattan MST . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 };
struct Dinic
{
3 Numerical algorithms 8 const ll Inf = 1e17;
3.1 SQRT For Loop . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . 8 vector<vector<int>> adj;
vector<vector<int>::iterator> cur;
3.2 Rabin Miller - Prime Checker . . . . . . .. . . . . . . . . . . . . . . . . . . 9 vector<Edge> s;
3.3 Chinese Remain Theorem . . . . . . . . .. . . . . . . . . . . . . . . . . . . 9 vector<int> h;
3.4 Pollard Rho - Factorialize . . . . . . . . .. . . . . . . . . . . . . . . . . . . 9 int sink, t;
int n;
3.5 FFT . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . 10 Dinic(int n)
3.6 FFT (Mod 998244353) . . . . . . . . . .. . . . . . . . . . . . . . . . . . . 10 {
3.7 FFT Mod (Template by Tran Khoi Nguyen) . .. . . . . . . . . . . . . . . . . . . 11 this->n = n;
3.8 Count Primes . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . 11
adj.resize(n + 1);
h.resize(n + 1);
3.9 Interpolation (Mod a prime) . . . . . . . .. . . . . . . . . . . . . . . . . . . 12 cur.resize(n + 1);
3.10 Bignum . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . 12 s.reserve(n);
Bignum with FFT multiplication . . . . . .. . . . . . . . . . . . . . . . . . . }
3.11 13
void AddEdge(int u, int v, ll c)
3.12 Tonelli Shanks (Find square root modulo prime) . . . . . . . . . . . . . . . . . . . 14 {
3.13 x
Discrete Logarithm (Find x that a ≡ b (mod m)) . . . . . . . . . . . . . . . . . . 14 s.emplace_back(u, v, c);
3.14 k
Primitive Root (Exist k that g ≡ a (mod n) for all a) . . . . . . . . . . . . . . . . 14
adj[u].push_back(s.size() - 1);
s.emplace_back(v, u, 0);
3.15 Discrete Root (Find x that x ≡ a (mod n), n is a prime) .
k . . . . . . . . . . . . . . 15 adj[v].push_back(s.size() - 1);
3.16 Super Sieve of Primes (Code by RR) . . . . . . . . . . . . . . . . . . . . . . . 15 }
bool BFS()
{
4 Graph algorithms 16 fill(h.begin(), h.end(), n + 2);
4.1 Twosat (2-SAT) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
queue<int> pq;
h[t] = 0;
4.2 Eulerian Path . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 pq.emplace(t);
4.3 Biconnected Component Tree . . . . . . . . . . . . . . . . . . . . . . . . . . 17 while (pq.size())
Heavy Light Decomposition (Template by Tran Khoi Nguyen) . . . . . . . . . . . . . . {
4.4 17
int c = pq.front();
4.5 Check Odd Circle With DSU (Template by Tran Khoi Nguyen) . . . . . . . . . . . . . 17 pq.pop();
for (auto i : adj[c])
if (h[s[i ^ 1].u] == n + 2 && s[i ^ 1].c != 0)
5 String 18 {
5.1 Palindrome Tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 h[s[i ^ 1].u] = h[c] + 1;
Suffix Array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . if (s[i ^ 1].u == sink)
5.2 18
return true;
5.3 Suffix Array (Template by Tran Khoi Nguyen) . . . . . . . . . . . . . . . . . . . . 19 pq.emplace(s[i ^ 1].u);
5.4 Aho Corasick - Extended KMP . . . . . . . . . . . . . . . . . . . . . . . . . . 20 }
5.5 Aho Corasick (Template by Tran Khoi Nguyen) . . . . . . . . . . . . . . . . . . . 20
}
return false;
5.6 Suffix Tree (Template by Tran Khoi Nguyen) . . . . . . . . . . . . . . . . . . . . 21 }
5.7 Z Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 ll DFS(int v, ll flowin)
{
if (v == t)
6 Data structures 21 return flowin;
ll flowout = 0;
6.1 Ordered Set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 for (; cur[v] != adj[v].end(); ++cur[v])
6.2 Fenwick Tree (With Walk on tree) . . . . . . . . . . . . . . . . . . . . . . . . 22 {
6.3 Convex Hull Trick (Min) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 int i = *cur[v];
6.4 Dynamic Convex Hull Trick (Min) . . . . . . . . . . . . . . . . . . . . . . . . 22
if (h[s[i].v] + 1 != h[v] || s[i].c == 0)
continue;
6.5 SPlay Tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 ll q = DFS(s[i].v, min(flowin, s[i].c));
6.6 Hashing (Template by Tran Khoi Nguyen) . . . . . . . . . . . . . . . . . . . . . 23 flowout += q;

1
6.7 BIT 2D (Template by Tran Khoi Nguyen) . . . . . . . . . . . . . . . . . . . . . 24
if (flowin != Inf)
 University of Engineering and Technology
flowin -= q; if (d[e.v] == d[e.u] + 1)
s[i].c -= q; {
s[i ^ 1].c += q; ll amt = e.cap - e.flow;
if (flowin == 0) if (flow != -1 && amt > flow)
break; amt = flow;
} if (ll pushed = DFS(e.v, T, amt))
return flowout; {
} e.flow += pushed;
void BlockFlow(ll &flow) oe.flow -= pushed;
{ return pushed;
for (int i = 1; i <= n; ++i) }
cur[i] = adj[i].begin(); }
flow += DFS(sink, Inf); }
} return 0;
ll MaxFlow(int s, int t) }
{
this->sink = s; ll MaxFlow(int S, int T)
this->t = t; {
ll flow = 0; ll total = 0;
while (BFS()) while (BFS(S, T))
BlockFlow(flow); {
return flow; fill(pt.begin(), pt.end(), 0);
} while (ll flow = DFS(S, T))
}; total += flow;
}
return total;
}
};
1.2 Dinic (Template by Tran Khoi Nguyen)
struct Edge
{
1.3 Min Cut (Template by Tran Khoi Nguyen)
int u, v;
ll cap, flow;
Edge() {} typedef long long LL;
Edge(int u, int v, ll cap) : u(u), v(v), cap(cap), flow(0) {}
}; struct Edge
{
struct Dinic int u, v;
{ LL cap, flow;
int N; Edge() {}
vector<Edge> E; Edge(int u, int v, LL cap) : u(u), v(v), cap(cap), flow(0) {}
vector<vector<int>> g; };
vector<int> d, pt;
map<ll, map<ll, ll>> mp; bool chk[maxn];
Dinic(int N) : N(N), E(0), g(N), d(N), pt(N) {} vector<pll> mincut;
struct Dinic
void AddEdge(int u, int v, ll cap) {
{ int N;
if (u != v) vector<Edge> E;
{ vector<vector<int>> g;
E.emplace_back(u, v, cap); vector<int> d, pt;
g[u].emplace_back(E.size() - 1);
E.emplace_back(v, u, 0); Dinic(int N) : N(N), E(0), g(N), d(N), pt(N) {}
g[v].emplace_back(E.size() - 1);
} void AddEdge(int u, int v, LL cap)
} {
// cout <<u<<" "<<v<<" "<<cap<<endl;
bool BFS(int S, int T) if (u != v)
{ {
queue<int> q({S}); E.emplace_back(u, v, cap);
fill(d.begin(), d.end(), N + 1); g[u].emplace_back(E.size() - 1);
d[S] = 0; E.emplace_back(v, u, 0);
while (!q.empty()) g[v].emplace_back(E.size() - 1);
{ }
int u = q.front(); }
q.pop();
if (u == T) bool BFS(int S, int T)
break; {
for (int k : g[u]) queue<int> q({S});
{ fill(d.begin(), d.end(), N + 1);
Edge &e = E[k]; d[S] = 0;
if (e.flow < e.cap && d[e.v] > d[e.u] + 1) while (!q.empty())
{ {
d[e.v] = d[e.u] + 1; int u = q.front();
q.emplace(e.v); q.pop();
} if (u == T)
} break;
} for (int k : g[u])
return d[T] != N + 1; {
} Edge &e = E[k];
if (e.flow < e.cap && d[e.v] > d[e.u] + 1)
ll DFS(int u, int T, ll flow = -1) {
{ d[e.v] = d[e.u] + 1;
if (u == T || flow == 0) q.emplace(e.v);
return flow; }
for (int &i = pt[u]; i < g[u].size(); ++i) }
{ }
Edge &e = E[g[u][i]]; return d[T] != N + 1;

2
Edge &oe = E[g[u][i] ^ 1]; }
 University of Engineering and Technology
S.reserve(nx);
LL DFS(int u, int T, LL flow = -1) h.resize(ny + 5);
{ adj.resize(nx + 5);
if (u == T || flow == 0) match.resize(ny + 5, NoMatch);
return flow; }
for (int &i = pt[u]; i < g[u].size(); ++i) void Clear()
{ {
Edge &e = E[g[u][i]]; for (int i = 1; i <= nx; ++i)
Edge &oe = E[g[u][i] ^ 1]; adj[i].clear();
if (d[e.v] == d[e.u] + 1) S.clear();
{ fill(match.begin(), match.end(), NoMatch);
LL amt = e.cap - e.flow; }
if (flow != -1 && amt > flow) void AddEdge(int x, int y)
amt = flow; {
if (LL pushed = DFS(e.v, T, amt)) adj[x].emplace_back(y);
{ }
e.flow += pushed; bool BFS()
oe.flow -= pushed; {
return pushed; fill(h.begin(), h.end(), 0);
} queue<int> q;
} for (auto x : S)
} for (auto i : adj[x])
return 0; if (h[i] == 0)
} {
void dfs1(ll u) q.emplace(i);
{ h[i] = 1;
// cout <<u<<endl; }
chk[u] = 1; while (q.size())
for (int k : g[u]) {
{ int x, ypop = q.front();
Edge e = E[k]; q.pop();
if (e.cap - e.flow > 0) if ((x = match[ypop]) == NoMatch)
{ return true;
ll nxt = e.v; for (auto i : adj[x])
if (!chk[nxt]) if (h[i] == 0)
dfs1(nxt); {
} h[i] = h[ypop] + 1;
} q.emplace(i);
} }
void find_mincut(ll S) }
{ return false;
dfs1(S); }
for (int i = 0; i < E.size(); i += 2) void dfs(int v, int lv)
{ {
auto p = E[i]; for (auto i : adj[v])
if (chk[p.u] && !chk[p.v]) if (h[i] == lv + 1)
{ {
mincut.pb(make_pair(p.u, p.v)); if (match[i] == NoMatch)
} found = 1;
} else
} dfs(match[i], lv + 1);
if (found)
LL MaxFlow(int S, int T) {
{ match[i] = v;
LL total = 0; return;
while (BFS(S, T)) }
{ }
fill(pt.begin(), pt.end(), 0); }
while (LL flow = DFS(S, T)) int MaxMatch()
total += flow; {
} int ans(0);
return total; for (int i = 1; i <= nx; ++i)
} S.emplace_back(i);
}; while (BFS())
{
for (int i = S.size() - 1; ~i; --i)
{
found = 0;
1.4 Maximum Matching (HopCroft - Karp) dfs(S[i], 0);
if (found)
{
++ans;
// Trace to find vertex cover and independence set S[i] = S.back();
/* S.pop_back();
Y* = Set of vertices y such that exist an argument path from y to a vertex x which isn’t }
matched }
X* = Set of matched vertices x that x isn’t matched with a vertex in Y* }
return ans;
(X* v Y*) is vertex cover }
*/
};
struct HopCroft_Karp
{
const int NoMatch = -1;
vector<int> h, S, match;
vector<vector<int>> adj;
int nx, ny;
1.5 Maximum Matching (Template by Tran Khoi Nguyen)
bool found;
HopCroft_Karp(int nx = 0, int ny = 0)
{ struct bipartite_graph
this->nx = nx; {

3
this->ny = ny; int nx, ny;
 University of Engineering and Technology
vector<vector<int>> gr;
vector<int> match, x_not_matched;

vector<int> level;
1.6 Maximum Matching (O(n2 ))
bool found;
// start from 1
void init(int _nx, int _ny) // O(n^2)
{ struct Maximum_matching
nx = _nx; {
ny = _ny; int nx, ny, t;
gr.assign(nx, vector<int>()); vector<int> Visited, match;
x_not_matched.resize(nx); vector<vector<int>> a;
for (int i = 0; i < nx; ++i)
x_not_matched[i] = i; Maximum_matching(int nx = 0, int ny = 0)
match.assign(ny, -1); {
level.resize(ny); Assign(nx, ny);
} }

void add_edge(int x, int y) void Assign(int nx, int ny)

{ {
gr[x].push_back(y); this->nx = nx;
} this->ny = ny;
t = 0;
bool bfs() Visited.assign(nx + 5, 0);
{ match.assign(ny + 5, 0);
fill(level.begin(), level.end(), 0); a.resize(nx + 5, {});
queue<int> q; }
for (auto &x : x_not_matched)
for (auto &y : gr[x]) void AddEdge(int x, int y)
if (level[y] == 0) {
{ a[x].emplace_back(y);
level[y] = 1; }
q.push(y);
} bool visit(int u)
while (!q.empty()) {
{ if (Visited[u] != t)
int ypop = q.front(); Visited[u] = t;
q.pop(); else
int x = match[ypop]; return false;
if (x == -1)
return true; for (int i = 0; i < a[u].size(); i++)
for (auto &y : gr[x]) {
if (level[y] == 0) int v = a[u][i];
{ if (!match[v] || visit(match[v]))
level[y] = level[ypop] + 1; {
q.push(y); match[v] = u;
} return true;
} }
return false; }
} return false;
}
void dfs(int x, int lv)
{ int MaxMatch()
for (auto &y : gr[x]) {
if (level[y] == lv + 1) int ans(0);
{
level[y] = 0; for (int i = 1; i <= nx; i++)
if (match[y] == -1) {
found = true; t++;
else ans += visit(i);
dfs(match[y], lv + 1); }
if (found)
{ return ans;
match[y] = x; }
return; };
}
}
}

int max_matching()
{
1.7 Min Cost Flow
int res = 0;
while (bfs()) struct Edge
{ {
for (int i = sz(x_not_matched) - 1; i >= 0; --i) int u, v;
{ ll c, w;
found = false; Edge(const int &u, const int &v, const ll &c, const ll &w) : u(u), v(v), c(c), w(w) {}
dfs(x_not_matched[i], 0); };
if (found) struct MaxFlowMinCost
{ {
++res; const ll Inf = 1e17;
x_not_matched[i] = x_not_matched.back(); int n, source, sink;
x_not_matched.pop_back(); vector<ll> d;
} vector<int> par;
} vector<bool> inqueue;
} vector<Edge> s;
return res; vector<vector<int>> adj;
} MaxFlowMinCost(int n)
} man; {
this->n = n;

4
s.reserve(n * 2);
 University of Engineering and Technology
d.resize(n + 5); inline void amax(T &x, const T &y)
inqueue.resize(n + 5); {
par.resize(n + 5); if (x < y)
adj.resize(n + 5); x = y;
} }
void AddEdge(int u, int v, ll c, ll w) typedef vector<int> VI;
{ typedef ll Flow;
s.emplace_back(u, v, c, w); typedef ll Cost;
adj[u].emplace_back(s.size() - 1); const Flow FLOW_INF = 1LL << 60;
s.emplace_back(v, u, 0, -w); const Cost COST_INF = 1LL << 60;
adj[v].emplace_back(s.size() - 1);
} const int SIZE = 65;
bool SPFA() vector<pair<Cost, int>> B[SIZE];
{ Cost last;
fill(d.begin(), d.end(), Inf); int sz;
fill(par.begin(), par.end(), s.size());
fill(inqueue.begin(), inqueue.end(), 0); int bsr(Cost c)
d[sink] = 0; {
queue<int> q; if (c == 0)
q.emplace(sink); return 0;
inqueue[sink] = 1; return __lg(c) + 1;
while (q.size()) }
{
int c = q.front(); void init()
inqueue[c] = 0; {
q.pop(); last = sz = 0;
for (auto i : adj[c]) REP(i, SIZE)
if (s[i ^ 1].c > 0 && d[s[i].v] > d[c] + s[i ^ 1].w) B[i].clear();
{ }
par[s[i].v] = i ^ 1;
d[s[i].v] = d[c] + s[i ^ 1].w; void push(Cost cst, int v)
if (!inqueue[s[i].v]) {
{ assert(cst >= last);
q.emplace(s[i].v); sz++;
inqueue[s[i].v] = 1; B[bsr(cst ^ last)].emplace_back(cst, v);
} }
}
} pair<Cost, int> pop_min()
return (d[source] < Inf); {
} assert(sz);
pair<ll, ll> MaxFlow(int so, int t, ll k) if (B[0].empty())
{ {
source = so; int k = 1;
sink = t; while (k < SIZE && B[k].empty())
ll Flow(0), cost(0); k++;
while (k && SPFA()) assert(k < SIZE);
{ last = B[k][0].first;
ll q(Inf); EACH(e, B[k])
int v = source; amin(last, e->first);
while (v != sink) EACH(e, B[k])
{ B[bsr(e->first ^ last)].push_back(*e);
q = min(q, s[par[v]].c); B[k].clear();
v = s[par[v]].v; }
} assert(B[0].size());
pair<Cost, int> ret = B[0].back();
q = min(q, k); B[0].pop_back();
sz--;
cost += d[source] * q; return ret;
Flow += q; }
k -= q;
struct MinCostMaxFlow
v = source; {
while (v != sink) struct Edge
{ {
s[par[v]].c -= q; int dst;
s[par[v] ^ 1].c += q; Cost cst;
v = s[par[v]].v; Flow cap;
} int rev;
} };
return {Flow, cost}; typedef vector<vector<Edge>> Graph;
} Graph G;
}; bool negative_edge;
MinCostMaxFlow(int N) : G(N)
{
negative_edge = false;
}
1.8 Min Cost Flow (Template by Tran Khoi Nguyen) void add_edge(int u, int v, Cost x, Flow f)
{
if (u == v)
namespace MIN_COST_MAX_FLOW return;
{ if (x < 0)
negative_edge = true;
#define REP(i, n) for (int i = 0, i##_len = (n); i < i##_len; ++i) G[u].push_back((Edge){
#define EACH(i, c) for (__typeof((c).begin()) i = (c).begin(), i##_end = (c).end(); i != i##_end; ++i) v, x, f, (int)G[v].size()});
template <class T> G[v].push_back((Edge){
inline void amin(T &x, const T &y) u, -x, 0, (int)G[u].size() - 1});
{ }
if (y < x)
x = y; void bellman_ford(int s, vector<Cost> &h)
} {

5
template <class T> fill(h.begin(), h.end(), COST_INF);
 University of Engineering and Technology
vector<bool> in(G.size()); }
h[s] = 0; };
in[s] = true; };
VI front, back; using MinCostMaxFlow = MIN_COST_MAX_FLOW::MinCostMaxFlow;
front.push_back(s);
while (1)
{
if (front.empty())
{
if (back.empty())
return;
2 Geometry
swap(front, back);
reverse(front.begin(), front.end());
}
int v = front.back();
2.1 Pick’s Theorem
front.pop_back();
in[v] = false; Given a certain lattice polygon (all its vertices have integer coordinates in some 2D grid) with non-
EACH(e, G[v]) zero area.
if (e->cap)
{ We denote its area by $S$, the number of points with integer coordinates lying strictly inside the
int w = e->dst; polygon by $I$ and the number of points lying on polygon sides by $B$.
if (h[w] > h[v] + e->cst)
{ Then, the Pick’s formula states:
h[w] = h[v] + e->cst;
if (!in[w]) \begin{center}
{ $S=I+ \frac{B}{2} - 1$
back.push_back(w); \end{center}
in[w] = true;
}
}
}

}
}
2.2 Smallest Circle - Emo Welzl (Contain all points)
Flow flow;
Cost cost; #include <bits/stdc++.h>
Flow solve(int s, int t, Flow limit = FLOW_INF) using namespace std;
{ using ld = double;
flow = 0;
cost = 0; typedef pair<ld, ld> point;
vector<Cost> len(G.size()), h(G.size()); typedef pair<point, ld> circle;
if (negative_edge) #define X first
bellman_ford(s, h); #define Y second

vector<int> prev(G.size()), prev_num(G.size()); // O(n)

while (limit > 0) // Remember to change size of set points
{ // All point must be save in array a[] below
init();
push(0, s); namespace emowelzl
fill(len.begin(), len.end(), COST_INF); {
fill(prev.begin(), prev.end(), -1); const int N = 100005; // Size of set points
len[s] = 0; int n;
while (sz) point a[N];
{
pair<Cost, int> p = pop_min(); point operator+(point a, point b)
Cost cst = p.first; {
int v = p.second; return point(a.X + b.X, a.Y + b.Y);
if (cst > len[v]) }
continue; point operator-(point a, point b) { return point(a.X - b.X, a.Y - b.Y); }
for (int i = 0; i < (int)G[v].size(); i++) point operator/(point a, ld x) { return point(a.X / x, a.Y / x); }
{ ld abs(point a) { return sqrt(a.X * a.X + a.Y * a.Y); }
const Edge &f = G[v][i];
Cost tmp = len[v] + f.cst + h[v] - h[f.dst]; point center_from(ld bx, ld by, ld cx, ld cy)
if (f.cap > 0 && len[f.dst] > tmp) {
{ ld B = bx * bx + by * by, C = cx * cx + cy * cy, D = bx * cy - by * cx;
len[f.dst] = tmp; return point((cy * B - by * C) / (2 * D), (bx * C - cx * B) / (2 * D));
push(tmp, f.dst); }
prev[f.dst] = v;
prev_num[f.dst] = i; circle circle_from(point A, point B, point C)
} {
} point I = center_from(B.X - A.X, B.Y - A.Y, C.X - A.X, C.Y - A.Y);
} return circle(I + A, abs(I));
}
if (prev[t] == -1)
return flow; circle f(int n, vector<point> T)
for (int i = 0; i < (int)G.size(); i++) {
h[i] += len[i]; if (T.size() == 3 || n == 0)
{
Flow f = limit; if (T.size() == 0)
for (int v = t; v != s; v = prev[v]) return circle(point(0, 0), -1);
f = min(f, G[prev[v]][prev_num[v]].cap); if (T.size() == 1)
for (int v = t; v != s; v = prev[v]) return circle(T[0], 0);
{ if (T.size() == 2)
Edge &e = G[prev[v]][prev_num[v]]; return circle((T[0] + T[1]) / 2, abs(T[0] - T[1]) / 2);
e.cap -= f; return circle_from(T[0], T[1], T[2]);
G[e.dst][e.rev].cap += f; }
} random_shuffle(a + 1, a + n + 1);
limit -= f; circle Result = f(0, T);
flow += f; for (int i = 1; i <= n; i++)
cost += f * h[t]; if (abs(Result.X - a[i]) > Result.Y + 1e-9)
} {

6
return flow; T.push_back(a[i]);
 University of Engineering and Technology
Result = f(i - 1, T); void DAC(int l, int r)
T.pop_back(); {
} if (r - l <= 3)
return Result; {
} Bruteforce(l, r);
}; return;
}
int main() int mid = (l + r) / 2;
{
cin >> emowelzl::n; DAC(l, mid);
DAC(mid + 1, r);
for (int i = 1; i <= emowelzl::n; ++i)
cin >> emowelzl::a[i].X >> emowelzl::a[i].Y; vector<int> s;
for (int i = l; i <= r; ++i)
cout << emowelzl::f(emowelzl::n, {}).Y * 2; if (sq(a[i].x - a[mid].x) <= ans)
} s.push_back(i);

Brute(s);
}

2.3 Closest pair of points in set void calc()

{
sort(a + 1, a + n + 1, [&](const Point &a, const Point &b)
{ return a.x < b.x || (a.x == b.x && a.y < b.y); });
// Find pair of points that have closest distance ans = Inf;

#include <iostream> DAC(1, n);

#include <cstdio>
#include <algorithm> if (xa > ya)
#include <iomanip> swap(xa, ya);
#include <cmath>
#include <vector> cout << xa << " " << ya << "\n";
cout << fixed << setprecision(6) << sqrt((ld)ans);
using namespace std; }
};
using ll = long long;
using ld = long double; Point a[N];
int n;
const int N = 5e4 + 2;
const ll Inf = 1e17; void Read()
#define sq(x) ((x) * (x)) {
cin >> n;
struct Point for (int i = 1; i <= n; ++i)
{ {
ll x, y; cin >> a[i].x >> a[i].y;
int id; a[i].id = i;
}
Point(const ll &x = 0, const ll &y = 0) : x(x), y(y) {} }

Point operator-(const Point &a) const void Solve()

{ {
return Point(x - a.x, y - a.y); ClosestPoint::n = n;
} for (int i = 1; i <= n; ++i)
ll len() ClosestPoint::a[i] = a[i];
{
return x * x + y * y; ClosestPoint::calc();
} }
};
int32_t main()
namespace ClosestPoint {
{ ios_base::sync_with_stdio(0);
int n, xa, ya; cin.tie(0);
ll ans; cout.tie(0);
Point a[N]; Read();
Solve();
ll Bruteforce(int l, int r) }
{
for (; l < r; ++l)
for (int i = l + 1; i <= r; ++i)
if (ans > (a[l] - a[i]).len())
{
ans = (a[l] - a[i]).len();
2.4 Manhattan MST
xa = a[l].id;
ya = a[i].id; // Idea is to reduce number of edges which are candidates to be in the MST
} // Then apply Kruskal algorithm to find MST
return ans;
} #include <bits/stdc++.h>
using namespace std;
void Brute(vector<int> &s) using ll = long long;
{
sort(s.begin(), s.end(), [&](const int &x, const int &y) constexpr int N = 2e5 + 2;
{ return a[x].y < a[y].y; }); constexpr ll Inf = 1e17;
for (int i = 0; i < s.size(); ++i)
for (int j = i + 1; j < s.size() && sq(abs(a[s[i]].y - a[s[j]].y)) <= ans; ++j) namespace manhattanMST
if (ans > (a[s[i]] - a[s[j]]).len()) {
{ /// disjoint set union
ans = (a[s[i]] - a[s[j]]).len(); struct dsu
xa = a[s[i]].id; {
ya = a[s[j]].id; int par[N];
} dsu()
} {

7
memset(par, -1, sizeof par);
 University of Engineering and Technology
} sort(s.begin(), s.end());
s.resize(unique(s.begin(), s.end()) - s.begin());
int findpar(int v) sort(v.begin(), v.end());
{
return par[v] < 0 ? v : par[v] = findpar(par[v]); FenwickTreeMin f(n);
}
for (auto [num1, num2, cost, idx] : v)
bool Merge(int u, int v) {
{ num2 = Find(s, num2);
u = findpar(u);
v = findpar(v); int res = f.Get(num2).second;
if (u == v) if (res != -1)
return false; edges.emplace_back(res, idx, dist(res, idx));

if (par[u] < par[v]) f.Update(num2, make_pair(cost, idx));

swap(u, v); }
}
par[v] += par[u];
par[u] = v; void calc()
{
return true; edges.clear();
}
}; createEdge(1, -1, -1, 0, 1, 1); // R1
createEdge(-1, 1, 0, -1, 1, 1); // R2
// Fenwick Tree Min createEdge(-1, -1, 0, 1, 1, -1); // R3
createEdge(1, 1, -1, 0, 1, -1); // R4
struct FenwickTreeMin createEdge(-1, 1, 1, 0, -1, -1); // R5
{ createEdge(1, -1, 0, 1, -1, -1); // R6
pair<ll, int> a[N]; createEdge(1, 1, 0, -1, -1, 1); // R7
int n; createEdge(-1, -1, 1, 0, -1, 1); // R8

FenwickTreeMin(int n = 0) dsu f;
{
Assign(n); sort(edges.begin(), edges.end());
} vector<pair<int, int>> res;
ll ans(0);
void Assign(int n)
{ for (auto i : edges)
this->n = n; if (f.Merge(i.u, i.v))
fill(a, a + n + 1, make_pair(Inf, -1)); {
} ans += i.w;
res.emplace_back(i.u, i.v);
void Update(int p, pair<ll, int> v) }
{
for (; p <= n; p += p & -p) cout << ans << "\n";
a[p] = min(a[p], v);
} for (auto i : res)
cout << i.first << " " << i.second << "\n";
pair<ll, int> Get(int p) }
{ };
pair<ll, int> ans({Inf, -1});
int32_t main()
for (; p; p -= p & -p) {
ans = min(ans, a[p]); ios_base::sync_with_stdio(0);
return ans; cin.tie(0);
} cout.tie(0);
};
cin >> manhattanMST::n;
struct Edge
{ for (int i = 1; i <= manhattanMST::n; ++i)
int u, v; cin >> manhattanMST::x[i] >> manhattanMST::y[i];
ll w;
Edge(const int &u = 0, const int &v = 0, const ll &w = 0) : u(u), v(v), w(w) {} manhattanMST::calc();
bool operator<(const Edge &a) const }
{
return w < a.w;
}
};

int n;
ll x[N], y[N];
3 Numerical algorithms
vector<Edge> edges;

ll dist(int i, int j) 3.1 SQRT For Loop

{
return abs(x[i] - x[j]) + abs(y[i] - y[j]);
} // Calculate n/1 + n/2 + ... + n/n

#define Find(x, v) (lower_bound(x.begin(), x.end(), v) - x.begin() + 1) #define Cal(a, b) ((b) - (a) + 1)

void createEdge(int a1, int a2, int b1, int b2, int c1, int c2)
{ ll calc(ll n)
vector<array<ll, 4>> v; {
vector<ll> s; ll ans = 0;

for (int i = 1; i <= n; i++) for (ll i = 1; i <= n / i; ++i)

{ ans += Cal(n / (i + 1) + 1, n / i) * i;
v.push_back({a1 * x[i] + a2 * y[i],
b1 * x[i] + b2 * y[i], for (ll i = 1; i < n / i; ++i)
c1 * x[i] + c2 * y[i], ans += n / i;
i});
s.emplace_back(b1 * x[i] + b2 * y[i]); return ans;

8
} }
 University of Engineering and Technology
ll Pow(ll a, ll b, const ll &mod)
3.2 Rabin Miller - Prime Checker {
ll ans(1);
for (; b; b >>= 1)
{
// There is another version of Rabin Miller using random in the implementation of Pollard Rho
if (b & 1)
ans = Mul(ans, a, mod);
ll mul(ll a, ll b, ll mod)
a = Mul(a, a, mod);
{
}
a %= mod;
b %= mod;
return ans;
ll q = (ld)a * b / mod;
}
ll r = a * b - q * mod;
return (r % mod + mod) % mod;
ll calPhi(ll n)
}
{
ll pow(ll a, ll n, ll m)
ll ans = 1;
{
ll result = 1;
for (ll i = 2; i * i <= n; ++i)
a %= m;
if (n % i == 0)
while (n > 0)
{
{
while (n % i == 0)
if (n & 1)
{
result = mul(result, a, m);
n /= i;
n >>= 1;
ans *= i;
a = mul(a, a, m);
}
}
return result;
ans = ans / i * (i - 1);
}
}
pair<ll, ll> factor(ll n)
{
if (n != 1)
ll s = 0;
ans *= n - 1;
while ((n & 1) == 0)
{
return ans;
s++;
}
n >>= 1;
}
pair<ll, ll> solve(const vector<ll> &a, const vector<ll> &b, vector<ll> phi = {})
return {s, n};
{
}
assert(a.size() == b.size()); // Assume a and b have the same size
bool test(ll s, ll d, ll n, ll witness)
ll m = 1;
{
if (n == witness)
{
return true;
m = 1;
ll p = pow(witness, d, n);
for (auto i : b)
if (p == 1)
m *= i;
return true;
}
for (; s > 0; s--)
{
if (phi.empty())
if (p == n - 1)
{
return true;
for (auto i : b)
p = mul(p, p, n);
phi.emplace_back(calPhi(i));
}
}
return false;
}
ll r = 0;
bool miller(ll n)
{
for (int i = 0; i < (int)b.size(); ++i)
if (n < 2)
r = (r + Mul(Mul(a[i], m / b[i], m), Pow(m / b[i], phi[i] - 1, m), m)) % m;
return false;
if ((n & 1) == 0)
return make_pair(r, m);
return n == 2;
}
ll s, d;
};
tie(s, d) = factor(n - 1);
return test(s, d, n, 2) && test(s, d, n, 3) && test(s, d, n, 5) &&
test(s, d, n, 7) && test(s, d, n, 11) && test(s, d, n, 13) &&
test(s, d, n, 17) && test(s, d, n, 19) && test(s, d, n, 23);
}
3.4 Pollard Rho - Factorialize
// You can change code Rabin-Miller (preposition)
struct PollardRho
3.3 Chinese Remain Theorem {
ll n;
map<ll, int> ans;
using ll = long long; PollardRho(ll n) : n(n) {}
using ld = long double; ll random(ll u)
namespace CRT {
{ return abs(rand()) % u;
// b must contain distinct element }
// x % b_i = a_i ll mul(ll a, ll b, ll mod)
// x % m == (a1 * m2 * m3 * ... * m_k) * [(m2 * m3 * ... * mk) ^ -1 mod m_1] + (a2 * m1 * m3 * ... {
* m_k) / [(m1 * m3 * ... * mk) ^ -1 mod m2] + ... a %= mod;
// Call CRT(a, b, phi) [Default phi is empty] b %= mod;
// return {r, m} that x % m == r ll q = (ld)a * b / mod;
ll r = a * b - q * mod;
// In case of overflow, use this function return (r % mod + mod) % mod;
ll Mul(ll a, ll b, const ll &mod) }
{
ll q = (ld)a * b / mod; ll pow(ll a, ll b, ll m)
ll r = a * b - q * mod; {
ll ans = 1;
return (r % mod + mod) % mod; a %= m;

9
} for (; b; b >>= 1)
 University of Engineering and Technology
{ }
if (b & 1) ll p = 0;
ans = mul(ans, a, m);
a = mul(a, a, m); while (p == 0 || p == n)
} p = findfactor(n);
return ans;
} pollard_rho(n / p);
pollard_rho(p);
pair<ll, ll> factor(ll n) }
{ };
ll s = 0;
while ((n & 1) == 0)
{
s++;

}
n >>= 1; 3.5 FFT
return {s, n};
} using cd = complex<double>;
// Rabin - Miller const double PI = acos(-1);
bool miller(ll n) // invert == true means Interpolation
{ // invert == false means dft
if (n < 2) void fft(vector<cd> &a, bool invert)
return 0; {
if (n == 2) int n = a.size();
return 1; for (int i = 1, j = 0; i < n; i++)
ll s = 0, m = n - 1; {
while (m % 2 == 0) int bit = n >> 1;
{ for (; j & bit; bit >>= 1)
s++; j ^= bit;
m >>= 1; j ^= bit;
} if (i < j)
// 1 - 0.9 ^ 40 swap(a[i], a[j]);
for (int it = 1; it <= 40; it++) }
{
ll u = random(n - 2) + 2; for (int len = 2; len <= n; len <<= 1)
ll f = pow(u, m, n); {
if (f == 1 || f == n - 1) double ang = 2 * PI / len * (invert ? -1 : 1);
continue; cd wlen(cos(ang), sin(ang));
for (int i = 1; i < s; i++) for (int i = 0; i < n; i += len)
{ {
f = mul(f, f, n); cd w(1);
if (f == 1) for (int j = 0; j < len / 2; j++)
return 0; {
if (f == n - 1) cd u = a[i + j],
break; v = a[i + j + len / 2] * w;
} a[i + j] = u + v;
if (f != n - 1) a[i + j + len / 2] = u - v;
return 0; w *= wlen;
} }
return 1; }
} }
ll f(ll x, ll n) if (invert)
{ {
return (mul(x, x, n) + 1) % n; for (cd &x : a)
} x /= n;
// Find a factor }
ll findfactor(ll n) }
{
ll x = random(n - 1) + 2;
ll y = x;
ll p = 1;
while (p == 1)
{ 3.6 FFT (Mod 998244353)
x = f(x, n);
y = f(f(y, n), n);
p = __gcd(abs(x - y), n); constexpr int N = 1e5 + 5; // keep N double of n+m
}
return p; // Call init() before call mul()
}
// prime factorization constexpr ll mod = 998244353;
void pollard_rho(ll n)
{ ll Pow(ll a, ll b, ll mod)
if (n <= 1000000) {
{ ll ans(1);
for (int i = 2; i * i <= n; i++) for (; b; b >>= 1)
{ {
while (n % i == 0) if (b & 1)
{ ans = ans * a % mod;
ans[i]++; a = a * a % mod;
n /= i; }
} return ans;
} }
if (n > 1)
ans[n]++; namespace ntt
return; {
} const int N = ::N;
const long long mod = ::mod, rt = 3;
if (miller(n)) ll G[55], iG[55], itwo[55];
{ void add(int &a, int b)

10
ans[n]++; {
return; a += b;
 University of Engineering and Technology
if (a >= mod) for (static int k = 2; k < n; k *= 2)
a -= mod; {
} R.resize(n);
void init() rt.resize(n);
{ auto x = polar(1.0L, acos(-1.0L) / k);
int now = (mod - 1) / 2, len = 1, irt = Pow(rt, mod - 2, mod); for (int i = k; i < 2 * k; i++)
while (now % 2 == 0) rt[i] = R[i] = i & 1 ? R[i / 2] * x : R[i / 2];
{ }
G[len] = Pow(rt, now, mod); vector<int> rev(n);
iG[len] = Pow(irt, now, mod); for (int i = 0; i < n; i++)
itwo[len] = Pow(1 << len, mod - 2, mod); rev[i] = (rev[i / 2] | (i & 1) << L) / 2;
now >>= 1;
len++; for (int i = 0; i < n; i++)
} if (i < rev[i])
} swap(a[i], a[rev[i]]);
void dft(ll *x, int n, int fg = 1) // fg=1 for dft, fg=-1 for inverse dft for (int k = 1; k < n; k *= 2)
{ for (int i = 0; i < n; i += 2 * k)
for (int i = (n >> 1), j = 1, k; j < n; ++j) for (int j = 0; j < k; j++)
{ {
if (i < j) // C z = rt[j+k] * a[i+j+k]; // (25% faster if hand-rolled) /// include-line
swap(x[i], x[j]); auto x = (double *)&rt[j + k], y = (double *)&a[i + j + k]; /// exclude-line
for (k = (n >> 1); k & i; i ^= k, k >>= 1) C z(x[0] * y[0] - x[1] * y[1], x[0] * y[1] + x[1] * y[0]); /// exclude-line
; a[i + j + k] = a[i + j] - z;
i ^= k; a[i + j] += z;
} }
for (int m = 2, now = 1; m <= n; m <<= 1, now++) }
{
ll r = fg > 0 ? G[now] : iG[now]; typedef vector<ll> vl;
for (int i = 0, j; i < n; i += m) vl convMod(const vl &a, const vl &b)
{ {
ll tr = 1, u, v;
for (j = i; j < i + (m >> 1); ++j) vl res((int)a.size() + b.size() - 1);
{ int B = 32 - __builtin_clz(res.size()), n = 1 << B, cut = int(sqrt(M));
u = x[j];
v = x[j + (m >> 1)] * tr % mod; vector<C> L(n), R(n), outs(n), outl(n);
x[j] = (u + v) % mod;
x[j + (m >> 1)] = (u + mod - v) % for (int i = 0; i < (int)a.size(); i++)
mod; L[i] = C((int)a[i] / cut, (int)a[i] % cut);
tr = tr * r % mod;
} for (int i = 0; i < (int)b.size(); i++)
} R[i] = C((int)b[i] / cut, (int)b[i] % cut);
}
} fft(L), fft(R);

// Take two sequence a, b; for (int i = 0; i < n; i++)

// return answer in sequence a {
int j = -i & (n - 1);
void mul(ll *a, ll *b, int n, int m) outl[j] = (L[i] + conj(L[j])) * R[i] / (2.0 * n);
{ outs[j] = (L[i] - conj(L[j])) * R[i] / (2.0 * n) / 1i;
// a: 0,1,2,...,n-1; b: 0,1,2,...,m-1 }

int nn = n + m - 1; fft(outl), fft(outs);

if (n == 0 || m == 0)
{ for (int i = 0; i < (int)res.size(); i++)
memset(a, 0, nn * sizeof(a[0])); {
return; ll av = (ll)(real(outl[i]) + .5);
} ll cv = (ll)(imag(outs[i]) + .5);
int L, len; ll bv = (ll)(imag(outl[i]) + .5) + (ll)(real(outs[i]) + .5);
for (L = 1, len = 0; L < nn; ++len, L <<= 1) res[i] = ((av % M * cut + bv) % M * cut + cv) % M;
; }
if (n < L) return res;
memset(a + n, 0, (L - n) * sizeof(a[0])); }
if (m < L) };
memset(b + m, 0, (L - m) * sizeof(b[0]));
dft(a, L, 1); // dft(a)
dft(b, L, 1); // dft(b)
// Merge
for (int i = 0; i < L; ++i)
a[i] = a[i] * b[i] % mod;
3.8 Count Primes
// Interpolation
dft(a, L, -1); // To initialize, call init_count_primes() first.
for (int i = 0; i < L; ++i) // Function count_primes(n) will compute the number of
a[i] = a[i] * itwo[len] % mod; // prime numbers lower than or equal to n.
} //
}; // Time complexity: Around O(n ^ 0.75)

constexpr int N = 1e5 + 5; // keep N larger than max(sqrt(n) + 2)

bool prime[N];
int prec[N];
3.7 FFT Mod (Template by Tran Khoi Nguyen) vector<int> P;

ll rec(ll n, int k)
{
struct FFTmod if (n <= 1 || k < 0)
{ return 0;
typedef complex<double> C; if (n <= P[k])
const ll M = mod; return n - 1;
void fft(vector<C> &a)
{ if (n < N && ll(P[k]) * P[k] > n)
int n = a.size(), L = 31 - __builtin_clz(n); return n - 1 - prec[n] + prec[P[k]];

11
static vector<complex<long double>> R(2, 1);
static vector<C> rt(2, 1); // (^ 10% faster if double) const int LIM = 250;
 University of Engineering and Technology
static int memo[LIM * LIM][LIM]; prf[0] = suf[k + 1] = 1;
bool ok = n < LIM * LIM;
if (ok && memo[n][k]) for (int i = 1; i <= k; i++)
return memo[n][k]; prf[i] = prf[i - 1] * (x - i + mod) % mod;

ll ret = n / P[k] - rec(n / P[k], k - 1) + rec(n, k - 1); for (int i = k; i >= 1; i--)
suf[i] = suf[i + 1] * (x - i + mod) % mod;
if (ok)
memo[n][k] = ret; ll res = 0;
return ret;
} for (int i = 1; i <= k; i++)
void init_count_primes() {
{ if (!((k - i) & 1))
prime[2] = true; res = (res + prf[i - 1] * suf[i + 1] % mod * ifac[i - 1] % mod * ifac[k - i] % mod * a
for (int i = 3; i < N; i += 2) [i]) % mod;
prime[i] = true; else
for (int i = 3, j; i * i < N; i += 2) res = (res - prf[i - 1] * suf[i + 1] % mod * ifac[i - 1] % mod * ifac[k - i] % mod * a
if (prime[i]) [i] % mod + mod) % mod;
for (j = i * i; j < N; j += i + i) }
prime[j] = false;
return res;
for (int i = 1; i < N; ++i) }
if (prime[i]) };
P.push_back(i);

for (int i = 1; i < N; ++i)

prec[i] = prec[i - 1] + prime[i];
}
3.10 Bignum
ll count_primes(ll n)
{
if (n < N) /// M is the number of digits in the answer
return prec[n]; /// In case that we don’t use multiplication, let BASE be 1e17 or 1e18
int k = prec[(int)sqrt(n) + 1]; /// a = Bignum("5")
return n - 1 - rec(n, k) + prec[P[k]]; /// The operator / is only for integer, the result is integer too
}
using cd = complex<long double>;
const long double PI = acos(-1);
const int M = 2000;
const ll BASE = 1e8;
3.9 Interpolation (Mod a prime) const int gd = log10(BASE);
const int maxn = M / gd + 1;
struct Bignum
{
// You can change mod into other prime number int n;
// update k to the degree of polynomial ll a[maxn];
// Just work when we know a[1] = P(1), a[2] = P(2),..., a[k] = P(k) [The degree of P(x) is k-1] Bignum(ll x = 0)
// update() then build() then cal() {
memset(a, 0, sizeof a);
/* n = 0;
* Complexity: O(Nlog(mod), N) do
*/ {
a[n++] = x % BASE;
constexpr ll mod = 1e9 + 7; // Change mod here x /= BASE;
constexpr ll N = 1e5 + 5; // Change size here } while (x);
}
struct Interpolation Bignum(const string &s)
{ {
ll a[N], fac[N], ifac[N], prf[N], suf[N]; Convert(s);
int k; }
ll stoll(const string &s)
ll Pow(ll a, ll b) {
{ ll ans(0);
ll ans(1); for (auto i : s)
for (; b; b >>= 1) ans = ans * 10 + i - ’0’;
{ return ans;
if (b & 1) }
ans = ans * a % mod; void Convert(const string &s)
a = a * a % mod; {
} memset(a, 0, sizeof a);
n = 0;
return ans; for (int i = s.size() - 1; ~i; --i)
} {
int j = max(0, i - gd + 1);
void upd(int u, ll v) a[n++] = stoll(s.substr(j, i - j + 1));
{ i = j;
a[u] = v; }
} fix();
}
void build() void fix()
{ {
fac[0] = ifac[0] = 1; ++n;
for (int i = 1; i < N; i++) for (int i = 0; i < n - 1; ++i)
{ {
fac[i] = (long long)fac[i - 1] * i % mod; a[i + 1] += a[i] / BASE;
ifac[i] = Pow(fac[i], mod - 2); a[i] %= BASE;
} if (a[i] < 0)
} {
a[i] += BASE;
// Calculate P(x) --a[i + 1];

12
ll calc(int x) }
{ }
 University of Engineering and Technology
while (n > 1 && a[n - 1] == 0) return 3;
--n; }
} bool operator<(const Bignum &s) const
Bignum &operator+=(const Bignum &x) {
{ return com(s) == 1;
n = max(n, x.n); }
for (int i = 0; i < n; ++i) bool operator>(const Bignum &s) const
a[i] += x.a[i]; {
fix(); return com(s) == 2;
return *this; }
} bool operator==(const Bignum &s) const
Bignum &operator-=(const Bignum &x) {
{ return com(s) == 3;
for (int i = 0; i < x.n; ++i) }
a[i] -= x.a[i]; bool operator<=(const Bignum &s) const
fix(); {
return *this; return com(s) != 2;
} }
Bignum &operator*=(const Bignum &x) bool operator>=(const Bignum &s) const
{ {
vector<ll> c(x.n + n, 0); return com(s) != 1;
for (int i = 0; i < n; ++i) }
for (int j = 0; j < x.n; ++j) void read()
c[i + j] += a[i] * x.a[j]; {
n += x.n; string s;
for (int i = 0; i < n; ++i) cin >> s;
a[i] = c[i]; Convert(s);
fix(); }
return *this; void print()
} {
Bignum &operator/=(const ll &x) int i = n;
{ while (i > 0 && a[i] == 0)
ll r = 0ll; --i;
for (int i = n - 1; i > -1; --i) cout << a[i];
{ for (--i; ~i; --i)
r = r * BASE + a[i]; cout << setw(gd) << setfill(’0’)
a[i] = r / x; << a[i];
r %= x; }
} };
fix();
return *this;
}
Bignum operator+(const Bignum &s)
{
Bignum c; 3.11 Bignum with FFT multiplication
copy(a, a + n, c.a);
c.n = n;
c += s; // Replace function *= in Bignum implementation with below code:
return c; void fft(vector<cd> &a, bool invert)
} {
Bignum operator-(const Bignum &s) int n = a.size();
{ for (int i = 1, j = 0; i < n; i++)
Bignum c; {
copy(a, a + n, c.a); int bit = n >> 1;
c.n = n; for (; j & bit; bit >>= 1)
c -= s; j ^= bit;
return c; j ^= bit;
} if (i < j)
Bignum operator*(const Bignum &s) swap(a[i], a[j]);
{ }
Bignum c; for (int len = 2; len <= n; len <<= 1)
copy(a, a + n, c.a); {
c.n = n; double ang = 2 * PI / len * (invert ? -1 : 1);
c *= s; cd wlen(cos(ang), sin(ang));
return c; for (int i = 0; i < n; i += len)
} {
Bignum operator/(const ll &x) cd w(1);
{ for (int j = 0; j < len / 2; j++)
Bignum c; {
copy(a, a + n, c.a); cd u = a[i + j], v = a[i + j + len / 2] * w;
c.n = n; a[i + j] = u + v;
c /= x; a[i + j + len / 2] = u - v;
return c; w *= wlen;
} }
ll operator%(const ll &x) }
{ }
ll ans(0); if (invert)
for (int i = n - 1; ~i; --i) {
ans = (ans * BASE + a[i]) % x; for (cd &x : a)
return ans; x /= n;
} }
int com(const Bignum &s) const }
{ Bignum &operator*=(const Bignum &x)
if (n < s.n) {
return 1; int m = 1;
if (n > s.n) while (m < n + x.n)
return 2; m <<= 1;
for (int i = n - 1; i > -1; --i) vector<cd> fa(m), fb(m);
if (a[i] > s.a[i]) for (int i = 0; i < m; ++i)
return 2; {

13
else if (a[i] < s.a[i]) fa[i] = a[i];
return 1; fb[i] = x.a[i];
 University of Engineering and Technology
} 3.13 Discrete Logarithm (Find x that ax ≡ b (mod m))
fft(fa, false); /// dft
fft(fb, false); /// dft
for (int i = 0; i < m; i++)
fa[i] *= fb[i]; // Returns minimum x for which a ^ x % m = b % m.
fft(fa, true); /// Interpolation // Returns -1 if x isn’t exist
n = m;
for (int i = 0; i < n; ++i) using ll = int;
a[i] = round(fa[i].real());
fix(); ll DiscreteLogarithm(ll a, ll b, ll m)
return *this; {
} a %= m, b %= m;
ll k = 1, add = 0, g;
while ((g = __gcd(a, m)) > 1)
{
if (b == k)
return add;
3.12 Tonelli Shanks (Find square root modulo prime) if (b % g)
return -1;
b /= g, m /= g, ++add;
k = (k * 1ll * a / g) % m;
/* }
Takes as input an odd prime p and n < p and returns r such that r * r = n [mod p].
There’s exist r if and only if n ^ [(p-1) / 2] = 1 (mod p) ll n = sqrt((ld)m) + 1;
*/ ll an = 1;
for (ll i = 0; i < n; ++i)
using ll = int; // Change type of data here an = (an * 1ll * a) % m;

ll Pow(ll a, ll b, ll mod) unordered_map<ll, ll> vals;

{ for (ll q = 0, cur = b; q <= n; ++q)
ll ans(1); {
vals[cur] = q;
for (; b; b >>= 1) cur = (cur * 1ll * a) % m;
{ }
if (b & 1)
ans = ans * a % mod; for (ll p = 1, cur = k; p <= n; ++p)
a = a * a % mod; {
} cur = (cur * 1ll * an) % m;
if (vals.count(cur))
return ans; {
} ll ans = n * p - vals[cur] + add;
return ans;
ll tonelli_shanks(ll n, ll p) }
{ }
ll s = 0; return -1;
ll q = p - 1; }
while ((q & 1) == 0)
{
q /= 2;
++s;
}
if (s == 1)
3.14 Primitive Root (Exist k that g k ≡ a (mod n) for all a)
{
ll r = Pow(n, (p + 1) / 4, p);
if ((r * r) % p == n) /*
return r; g is a primitive root modulo n if and only if for any integer a such that gcd(a, n) = 1, there exists
return 0; an integer k such that:
} g^k = a (mod n).
// Find the first quadratic non-residue z by brute-force search
ll z = 1; Primitive root modulo n exists if and only if:
while (Pow(++z, (p - 1) / 2, p) != p - 1) - n is 1, 2, 4
; - n is power of an odd prime number (n = p ^ k)
ll c = Pow(z, q, p); - n is twice power of an odd prime number (n = 2 * p ^ k)
ll r = Pow(n, (q + 1) / 2, p);
ll t = Pow(n, q, p); This theorem was proved by Gauss in 1801.
ll m = s; */
while (t != 1)
{ using ll = int; // Change type of data here
ll tt = t;
ll i = 0; ll Pow(ll a, ll b, ll mod)
while (tt != 1) {
{ ll ans(1);
tt = (tt * tt) % p;
++i; for (; b; b >>= 1)
if (i == m) {
return 0; if (b & 1)
} ans = ans * a % mod;
ll b = Pow(c, Pow(2, m - i - 1, p - 1), p); a = a * a % mod;
ll b2 = (b * b) % p; }
r = (r * b) % p;
t = (t * b2) % p; return ans;
c = b2; }
m = i;
} ll GetPhi(ll n)
if ((r * r) % p == n) {
return r; ll ans(1);
return -1; // Can’t find
} for (ll i = 2; i * i <= n; ++i)
if (n % i == 0)
{

14
while (n % i == 0)
{
 University of Engineering and Technology
n /= i; // Sieve up to 10^9 by RR
ans *= i; namespace Sieve
} {
const int MAX = 1000000000LL;
ans = ans / i * (i - 1); const int WHEEL = 3 * 5 * 7 * 11 * 13;
} const int N_SMALL_PRIMES = 6536; // cnt primes less than 2^16
const int SIEVE_SPAN = WHEEL * 64; // one iteration of segmented sieve
if (n != 1) const int SIEVE_SIZE = SIEVE_SPAN / 128 + 1;
ans *= n - 1;
uint64_t ONES[64]; // ONES[i] = 1<<i
return ans; int small_primes[N_SMALL_PRIMES]; // primes less than 2^16
}
// each element of sieve is a 64-bit bitmask.
ll PrimitiveRoot(ll p) // Each bit (0/1) stores whether the corresponding element is a prime number.
{ // We only need to store odd numbers
vector<ll> fact; // -> 1st bitmask stores 3, 5, 7, 9, ...
ll phi = GetPhi(p); uint64_t si[SIEVE_SIZE];
ll n = phi; // for each ’wheel’, we store the sieve pattern (i.e. what numbers cannot be primes)
uint64_t pattern[WHEEL];
for (int i = 2; i * i <= n; ++i)
if (n % i == 0) inline void mark(uint64_t *s, int o) { s[o >> 6] |= ONES[o & 63]; }
{ inline int test(uint64_t *s, int o) { return (s[o >> 6] & ONES[o & 63]) == 0; }
fact.push_back(i);
while (n % i == 0) // update_sieve {{{
n /= i; void update_sieve(int offset)
} {
// copy each wheel pattern to sieve
if (n > 1) for (int i = 0, k; i < SIEVE_SIZE; i += k)
fact.push_back(n); {
k = std::min(WHEEL, SIEVE_SIZE - i);
for (ll res = 2; res <= p; ++res) memcpy(si + i, pattern, sizeof(*pattern) * k);
{ }
bool ok = true;
// Correctly mark 1, 3, 5, 7, 11, 13 as not prime / primes
for (int i = 0; i < fact.size() && ok; ++i) if (offset == 0)
ok &= Pow(res, phi / fact[i], p) != 1; {
si[0] |= ONES[0];
if (ok) si[0] &= ~(ONES[1] | ONES[2] | ONES[3] | ONES[5] | ONES[6]);
return res; }
}
// sieve for primes >= 17 (stored in ‘small_primes‘)
return -1; // can’t find for (int i = 0; i < N_SMALL_PRIMES; ++i)
} {
int j = small_primes[i] * small_primes[i];
if (j > offset + SIEVE_SPAN - 1)
break;
if (j > offset)
3.15 Discrete Root (Find x that xk ≡ a (mod n), n is a prime) else
j = (j - offset) >> 1;

{
j = small_primes[i] - offset % small_primes[i];
/* if ((j & 1) == 0)
Given a prime n and two integers a and k, find all x for which:\ j += small_primes[i];
- x ^ k = a (mod n) j >>= 1;
}
Notice: while (j < SIEVE_SPAN / 2)
- In case k = 2, let’s use Tonelli - Shanks {
- You must insert my implementation of Discrete Logarithm and Primitive Root to run algorithm mark(si, j);
*/ j += small_primes[i];
}
ll Pow(ll a, ll b, ll mod) }
{ }
ll ans(1); // }}}

for (; b; b >>= 1) void sieve()

{ {
if (b & 1) // init small primes {{{
ans = ans * a % mod; for (int i = 0; i < 64; ++i)
a = a * a % mod; ONES[i] = 1ULL << i;
}
// sieve to find small primes
return ans; for (int i = 3; i < 256; i += 2)
} {
if (test(si, i >> 1))
ll DiscreteRoot(ll a, ll k, ll n) {
{ for (int j = i * i / 2; j < 32768; j += i)
ll g = PrimitiveRoot(n); mark(si, j);
ll v = Pow(g, k, n); }
ll ans = DiscreteLogarithm(v, a, n); }
// store primes >= 17 in ‘small_primes‘ (we will sieve differently
if (ans == -1) // for primes 2, 3, 5, 7, 11, 13)
return -1; // Can’t find {
int m = 0;
return Pow(g, ans, n); for (int i = 8; i < 32768; ++i)
} {
if (test(si, i))
small_primes[m++] = i * 2 + 1;
}
}
3.16 Super Sieve of Primes (Code by RR) // }}}

15
// For primes 3, 5, 7, 11, 13: we initialize wheel pattern..
 University of Engineering and Technology
for (int i = 1; i < WHEEL * 64; i += 3) }
mark(pattern, i); if (low[u] == num[u])
for (int i = 2; i < WHEEL * 64; i += 5) {
mark(pattern, i); while (1)
for (int i = 3; i < WHEEL * 64; i += 7) {
mark(pattern, i); int v = stTarjan.back();
for (int i = 5; i < WHEEL * 64; i += 11) stTarjan.pop_back();
mark(pattern, i); root[v] = numComp;
for (int i = 6; i < WHEEL * 64; i += 13) if (u == v)
mark(pattern, i); break;
}
// Segmented sieve numComp++;
long long sum_primes = 2; }
}
for (int offset = 0; offset < MAX; offset += SIEVE_SPAN) bool solve()
{ {
update_sieve(offset); for (int i = 0; i < n; i++)
if (root[i] == -1)
for (uint32_t j = 0; j < SIEVE_SIZE; j++) tarjan(i);
{ for (int i = 0; i < n; i += 2)
uint64_t x = ~si[j]; {
while (x) if (root[i] == root[i ^ 1])
{ return 0;
uint32_t p = offset + (j << 7) + (__builtin_ctzll(x) << 1) + 1; color[i] = (root[i] < root[i ^ 1]);
if (p > offset + SIEVE_SPAN - 1) }
break; return 1;
if (p <= MAX) }
{ };
// p is a prime
}
x ^= (-x & x);
}

}
} 4.2 Eulerian Path
}
}; // Path that goes all edges
// Start from 1
struct EulerianGraph
{
vector<vector<pair<int, int>>> a;
int num_edges;
4 Graph algorithms EulerianGraph(int n)
{
a.resize(n + 1);
4.1 Twosat (2-SAT) }
num_edges = 0;

void add_edge(int u, int v, bool undirected = true)

// start from 0 {
// pos(V) is the vertex that represent V in graph a[u].push_back(make_pair(v, num_edges));
// neg(V) is the vertex that represent !V if (undirected)
// pos(V) ^ neg(V) = 1, use two functions below a[v].push_back(make_pair(u, num_edges));
// (U v V) <=> (!U -> V) <=> (!V -> U) num_edges++;
// You need do addEge(represent(U), represent(V)) }
// solve() == false mean no answer
// Want to get the answer ? vector<int> get_eulerian_path()
// color[pos(U)] = 1 means we choose U {
// otherwise, we don’t vector<int> path, s;
constexpr int N = 1e5 + 5; // Keep N double of n vector<bool> was(num_edges);
inline int pos(int u) { return u << 1; }
inline int neg(int u) { return u << 1 | 1; } s.push_back(1);
struct TwoSAT // start of eulerian path
{ // directed graph: deg_out - deg_in == 1
int n, numComp, cntTarjan; // undirected graph: odd degree
vector<int> adj[N], stTarjan; // for eulerian cycle: any vertex is OK
int low[N], num[N], root[N], color[N];
TwoSAT(int n) : n(n * 2) while (!s.empty())
{ {
memset(root, -1, sizeof root); int u = s.back();
memset(low, -1, sizeof low); bool found = false;
memset(num, -1, sizeof num); while (!a[u].empty())
memset(color, -1, sizeof color); {
cntTarjan = 0; int v = a[u].back().first;
stTarjan.clear(); int e = a[u].back().second;
} a[u].pop_back();
void addEdge(int u, int v) if (was[e])
{ continue;
adj[u ^ 1].push_back(v); was[e] = true;
adj[v ^ 1].push_back(u); s.push_back(v);
} found = true;
void tarjan(int u) break;
{ }
stTarjan.push_back(u); if (!found)
num[u] = low[u] = cntTarjan++; {
for (int v : adj[u]) path.push_back(u);
{ s.pop_back();
if (root[v] != -1) }
continue; }
if (low[v] == -1) reverse(path.begin(), path.end());

16
tarjan(v); return path;
low[u] = min(low[u], low[v]); }
 University of Engineering and Technology
}; return base;
if (left == right)
return st[id];
dosth(id, left, right);
ll mid = (left + right) / 2;
4.3 Biconnected Component Tree }
return min(get(id * 2, left, mid, x), get(id * 2 + 1, mid + 1, right, x));

ll cntnw = 0;
ll nchain = 1;
ll chainhead[maxn];
// Biconnected Component Tree ll chainid[maxn];
// 1 is the root of Tree ll id[maxn];
// n + i is the node that represent i-th bcc, its depth is even vector<pll> adj[maxn];
ll par1[maxn];
const int N = 3e5 + 5; // Change size to n + number of bcc (For safety, set N >= 2 * n) ll siz[maxn];
int n, nBicon, nTime; void hld(ll u, ll par)
int low[N], num[N]; {
vector<int> adj[N], nadj[N]; if (!chainhead[nchain])
vector<int> s; chainhead[nchain] = u;
cntnw++;
void dfs(int v, int p = -1) chainid[u] = nchain;
{ id[u] = cntnw;
low[v] = num[v] = ++nTime; ll nxt = -1;
s.emplace_back(v); for (auto p : adj[u])
{
for (auto i : adj[v]) ll to = p.ff;
if (i != p) if (to == par)
{ continue;
if (!num[i]) if (nxt == -1 || siz[nxt] < siz[to])
{ {
dfs(i, v); nxt = to;
low[v] = min(low[v], low[i]); }
}
if (low[i] >= num[v]) if (nxt != -1)
{ {
++nBicon; hld(nxt, u);
nadj[v].emplace_back(n + nBicon); }
for (auto p : adj[u])
int vertex; {
do ll to = p.ff;
{ if (to == par || to == nxt)
vertex = s.back(); continue;
s.pop_back(); nchain++;
hld(to, u);
nadj[n + nBicon].emplace_back(vertex); }
} while (vertex != i); }
} void update1(ll u, ll a, ll w)
} {
else ll p = chainid[u];
low[v] = min(low[v], num[i]); ll chk = chainid[a];
} while (1)
} {
if (p == chk)
{
update(1, 1, cntnw, id[a], id[u], w);
break;
4.4 Heavy Light Decomposition (Template by Tran Khoi }
update(1, 1, cntnw, id[chainhead[p]], id[u], w);
Nguyen) u = par1[chainhead[p]];
p = chainid[u];
}
}
ll st[4 * maxn];
ll la[4 * maxn];
void dosth(ll id, ll left, ll right)
{
if (left == right)
return;
4.5 Check Odd Circle With DSU (Template by Tran Khoi
st[id * 2] = min(st[id * 2], la[id]);
st[id * 2 + 1] = min(st[id * 2 + 1], la[id]); Nguyen)
la[id * 2] = min(la[id * 2], la[id]);
la[id * 2 + 1] = min(la[id * 2 + 1], la[id]);
la[id] = base; namespace ufs
} {
void update(ll id, ll left, ll right, ll x, ll y, ll w) struct node
{ {
if (x > right || y < left) int fa, val, size;
return; } t[30];
if (x <= left && y >= right) struct info
{ {
st[id] = min(st[id], w); int x, y;
la[id] = min(la[id], w); node a, b;
return; } st[30];
} inline void pre()
dosth(id, left, right); {
ll mid = (left + right) / 2; for (int i = 1; i <= n; i++)
update(id * 2, left, mid, x, y, w); t[i] = (node){i, 0, 1};
update(id * 2 + 1, mid + 1, right, x, y, w); }
} inline int find(int x)
ll get(ll id, ll left, ll right, ll x) {

17
{ while (t[x].fa != x)
if (x > right || x < left) x = t[x].fa;
 University of Engineering and Technology
return x; 5.2 Suffix Array
}
inline int dis(int x)
{ // string and array pos start from 0
int ans = 0; // but array sa and lcp start from 1
while (t[x].fa != x)
ans ^= t[x].val, x = t[x].fa; constexpr int N = 3e5 + 5; // change size to size of string;
return ans;
} struct SuffixArray
inline void link(int x, int y) {
{ string s;
int val = dis(x) ^ dis(y) ^ 1; int n, c[N], p[N], rp[N], lcp[N];
x = find(x);
y = find(y); //p[] : suffix array
if (t[x].size > t[y].size) // lcp[]: lcp array
swap(x, y);
t[x].fa = y; void Assign(const string &x)
t[x].val = val; {
t[y].size += t[x].size; s = x;
} s.push_back(’$’); // Change character here due to range of charater in string
n = s.size();
} Build();
using namespace ufs; s.pop_back();
n = s.size();
}

void Build()
{
5 String vector<int> pn(N), cn(N), cnt(N);

for (int i = 0; i < n; ++i)

++cnt[s[i]];
5.1 Palindrome Tree for (int i = 1; i <= 256; ++i)
cnt[i] += cnt[i - 1];
for (int i = 0; i < n; ++i)
p[--cnt[s[i]]] = i;
// base on idea odd palindrome, even palindrome
// 0-odd is the root of tree for (int i = 1; i < n; ++i)
struct node c[p[i]] = c[p[i - 1]] + (s[p[i]] != s[p[i - 1]]);
{
int len; int maxn = c[p[n - 1]];
node *child[26], *sufflink;
node() for (int i = 0; (1 << i) < n; ++i)
{ {
len = 0; for (int j = 0; j < n; ++j)
for (int i = 0; i < 26; ++i) p[j] = ((p[j] - (1 << i)) % n + n) % n;
child[i] = NULL; for (int j = 0; j <= maxn; ++j)
sufflink = NULL; cnt[j] = 0;
}
}; for (int j = 0; j < n; ++j)
struct PalindromeTree ++cnt[c[p[j]]];
{
node odd, even; for (int j = 1; j <= maxn; ++j)
PalindromeTree() cnt[j] += cnt[j - 1];
{
odd.len = -1; for (int j = n - 1; ~j; --j)
odd.sufflink = &odd; pn[--cnt[c[p[j]]]] = p[j];
even.len = 0;
even.sufflink = &odd; for (int j = 1; j < n; ++j)
} cn[pn[j]] = cn[pn[j - 1]] + (c[pn[j]] != c[pn[j - 1]] || c[(pn[j] + (1 << i)) % n] !=
void Assign(string &s) c[(pn[j - 1] + (1 << i)) % n]);
{
node *last = &even; maxn = cn[pn[n - 1]];
for (int i = 0; i < (int)s.size(); ++i)
{ for (int j = 0; j < n; ++j)
node *tmp = last; {
while (s[i - tmp->len - 1] != s[i]) p[j] = pn[j];
tmp = tmp->sufflink; c[j] = cn[j];
if (tmp->child[s[i] - ’a’]) }
{ }
last = tmp->child[s[i] - ’a’]; }
continue;
} void BuildLCP()
tmp->child[s[i] - ’a’] = new node; {
last = tmp->child[s[i] - ’a’]; for (int i = 1; i <= n; ++i)
last->len = tmp->len + 2; rp[p[i]] = i;
if (last->len == 1) for (int i = 0; i < n; ++i)
{ {
last->sufflink = &even; if (i)
continue; lcp[i] = max(lcp[i - 1] - 1, 0);
} if (rp[i] == n)
tmp = tmp->sufflink; continue;
while (s[i - tmp->len - 1] != s[i])
tmp = tmp->sufflink; while (lcp[i] < n - i && lcp[i] < n - p[rp[i] + 1] && s[i + lcp[i]] == s[p[rp[i] + 1] +
last->sufflink = tmp->child[s[i] - ’a’]; lcp[i]])
} ++lcp[i];
} }
}; }
} g;

18
 University of Engineering and Technology
5.3 Suffix Array (Template by Tran Khoi Nguyen) {
int v = sa[i];
if (v >= 1 && !ls[v - 1])
sa[buf[s[v - 1]]++] = v - 1;
struct suffix_array }
{ copy(sum_l.begin(), sum_l.end(), buf.begin());
vector<int> sa_naive(const vector<int> &s) for (auto i = n - 1; i >= 0; --i)
{ {
int n = (int)s.size(); int v = sa[i];
vector<int> sa(n); if (v >= 1 && ls[v - 1])
iota(sa.begin(), sa.end(), 0); sa[--buf[s[v - 1] + 1]] = v - 1;
sort(sa.begin(), sa.end(), [&](int l, int r) }
{ };
if(l == r) return false; vector<int> lms_map(n + 1, -1);
for(; l < n && r < n; ++ l, ++ r) if(s[l] != s[r]) return s[l] < s[r]; int m = 0;
return l == n; }); for (auto i = 1; i < n; ++i)
return sa; if (!ls[i - 1] && ls[i])
} lms_map[i] = m++;
vector<int> sa_doubling(const vector<int> &s) vector<int> lms;
{ lms.reserve(m);
int n = (int)s.size(); for (auto i = 1; i < n; ++i)
vector<int> sa(n), rank = s, tmp(n); if (!ls[i - 1] && ls[i])
iota(sa.begin(), sa.end(), 0); lms.push_back(i);
for (auto k = 1; k < n; k <<= 1) induce(lms);
{ if (m)
auto cmp = [&](int x, int y) {
{ vector<int> sorted_lms;
if (rank[x] != rank[y]) sorted_lms.reserve(m);
return rank[x] < rank[y]; for (auto v : sa)
int rx = x + k < n ? rank[x + k] : -1; if (lms_map[v] != -1)
int ry = y + k < n ? rank[y + k] : -1; sorted_lms.push_back(v);
return rx < ry; vector<int> rec_s(m);
}; int rec_sigma = 0;
sort(sa.begin(), sa.end(), cmp); rec_s[lms_map[sorted_lms[0]]] = 0;
tmp[sa[0]] = 0; for (auto i = 1; i < m; ++i)
for (auto i = 1; i < n; ++i) {
tmp[sa[i]] = tmp[sa[i - 1]] + (cmp(sa[i - 1], sa[i]) ? 1 : 0); int l = sorted_lms[i - 1], r = sorted_lms[i];
swap(tmp, rank); int end_l = (lms_map[l] + 1 < m) ? lms[lms_map[l] + 1] : n;
} int end_r = (lms_map[r] + 1 < m) ? lms[lms_map[r] + 1] : n;
return sa; bool same = true;
} if (end_l - l != end_r - r)
template <int THRESHOLD_NAIVE = 10, int THRESHOLD_DOUBLING = 40> same = false;
vector<int> sa_is(const vector<int> &s, int sigma) else
{ {
int n = (int)s.size(); for (; l < end_l; ++l, ++r)
if (n == 0) if (s[l] != s[r])
return {}; break;
if (n == 1) if (l == n || s[l] != s[r])
return {0}; same = false;
if (n == 2) }
{ if (!same)
if (s[0] < s[1]) ++rec_sigma;
return {0, 1}; rec_s[lms_map[sorted_lms[i]]] = rec_sigma;
else }
return {1, 0}; auto rec_sa = sa_is<THRESHOLD_NAIVE, THRESHOLD_DOUBLING>(rec_s, rec_sigma + 1);
} for (auto i = 0; i < m; ++i)
if (n < THRESHOLD_NAIVE) sorted_lms[i] = lms[rec_sa[i]];
return sa_naive(s); induce(sorted_lms);
if (n < THRESHOLD_DOUBLING) }
return sa_doubling(s); return sa;
vector<int> sa(n); }
vector<bool> ls(n); int n;
for (auto i = n - 2; i >= 0; --i) // data: sorted sequence of suffices including the empty suffix
ls[i] = (s[i] == s[i + 1]) ? ls[i + 1] : (s[i] < s[i + 1]); // rank[i]: position of the suffix i in the suffix array
vector<int> sum_l(sigma), sum_s(sigma); // lcp[i]: longest common prefix of data[i] and data[i + 1]
for (auto i = 0; i < n; ++i) // index start from 1
{ vector<int> data, rank, lcp;
if (!ls[i]) // O(n + sigma)
++sum_s[s[i]]; suffix_array(const vector<int> &s, int sigma) : n((int)s.size()), rank(n + 1), lcp(n)
else {
++sum_l[s[i] + 1]; assert(0 <= sigma);
} for (auto d : s)
for (auto i = 0; i < sigma; ++i) assert(0 <= d && d < sigma);
{ data = sa_is(s, sigma);
sum_s[i] += sum_l[i]; data.insert(data.begin(), n);
if (i + 1 < sigma) for (auto i = 0; i <= n; ++i)
sum_l[i + 1] += sum_s[i]; rank[data[i]] = i;
} for (auto i = 0, h = 0; i <= n; ++i)
auto induce = [&](const vector<int> &lms) {
{ if (h > 0)
fill(sa.begin(), sa.end(), -1); --h;
vector<int> buf(sigma); if (rank[i] == 0)
copy(sum_s.begin(), sum_s.end(), buf.begin()); continue;
for (auto d : lms) int j = data[rank[i] - 1];
{ for (; j + h <= n && i + h <= n; ++h)
if (d == n) if ((j + h == n) != (i + h == n) || j + h < n && s[j + h] != s[i + h])
continue; break;
sa[buf[s[d]]++] = d; lcp[rank[i] - 1] = h;
} }
copy(sum_l.begin(), sum_l.end(), buf.begin()); }

19
sa[buf[s[n - 1]]++] = n - 1; // O(n log n) time, O(n) space
for (auto i = 0; i < n; ++i) template <class T>
 University of Engineering and Technology
suffix_array(const vector<T> &s, bool prepare_lcp) : n((int)s.size()), rank(n + 1), lcp(n) // for convinient a->to[v] = the node x->to[v] that a match x and x->to[v] != NULL
{ // root -> suflink = root
vector<int> idx(n);
iota(idx.begin(), idx.end(), 0); void build()
sort(idx.begin(), idx.end(), [&](int l, int r) {
{ return s[l] < s[r]; }); queue<Node *> Q;
vector<int> s2(n); root->suflink = root;
int now = 0; Q.push(root);
for (auto i = 0; i < n; ++i)
{ while (!Q.empty())
if (i && s[idx[i - 1]] != s[idx[i]]) {
++now; Node *par = Q.front();
s2[idx[i]] = now; Q.pop();
} for (int c = 0; c < ALPHABET_SIZE; ++c)
data = sa_is(s2, now + 1); {
data.insert(data.begin(), n); if (par->to[c])
for (auto i = 0; i <= n; ++i) {
rank[data[i]] = i; par->to[c]->suflink = par == root ? root : par->suflink->to[c];
for (auto i = 0, h = 0; i <= n; ++i) Q.push(par->to[c]);
{ }
if (h > 0) else
--h; {
if (rank[i] == 0) par->to[c] = par == root ? root : par->suflink->to[c];
continue; }
int j = data[rank[i] - 1]; }
for (; j + h <= n && i + h <= n; ++h) }
if ((j + h == n) != (i + h == n) || j + h < n && s[j + h] != s[i + h]) }
break; };
lcp[rank[i] - 1] = h;
}
}
// RMQ must be built over lcp
// O(1)
template <class RMQ>
5.5 Aho Corasick (Template by Tran Khoi Nguyen)
int longest_common_prefix(int i, int j, const RMQ &rmq) const
{
assert(0 <= i && i <= n && 0 <= j && j <= n); struct aho_corasick
return i == j ? n - i : rmq.query(min(rank[i], rank[j]), max(rank[i], rank[j])); {
} struct tk
}; {
ll link;
ll nxt[27];
ll par;
char ch;
5.4 Aho Corasick - Extended KMP ll go[27];
ll val;
ll leaf;
constexpr int ALPHABET_SIZE = 26; tk(ll par = -1, char ch = ’a’) : par(par), ch(ch)
constexpr int firstCharacter = ’a’; {
memset(nxt, -1, sizeof(nxt));
struct Node memset(go, -1, sizeof(go));
{ val = -1;
Node *to[ALPHABET_SIZE]; link = -1;
Node *suflink; leaf = 0;
int ending_length; // 0 if is not ending }
};
Node() vector<tk> vt;
{ void init()
for (int i = 0; i < ALPHABET_SIZE; ++i) {
to[i] = NULL; vt.clear();
suflink = NULL; vt.pb({-1, ’a’});
ending_length = false; }
} ll add(string s, ll val)
}; {
ll nw = 0;
struct AhoCorasick for (auto to : s)
{ {
Node *root; if (vt[nw].nxt[to - ’a’ + 1] == -1)
{
AhoCorasick() vt[nw].nxt[to - ’a’ + 1] = vt.size();
{ vt.pb({nw, to});
root = new Node(); }
} nw = vt[nw].nxt[to - ’a’ + 1];
}
void add(const string &s) vt[nw].leaf = val;
{ return nw;
Node *cur_node = root; }
ll get_val(ll u)
for (char c : s) {
{ if (u == 0)
int v = c - firstCharacter; return 0;
if (vt[u].val == -1)
if (!cur_node->to[v]) {
cur_node->to[v] = new Node(); vt[u].val = vt[u].leaf + get_val(get_link(u));
}
cur_node = cur_node->to[v]; return vt[u].val;
} }
ll go(ll v, ll t)
cur_node->ending_length = s.size(); {
} if (vt[v].go[t] == -1)

20
{
// if a->to[v] == NULL if (vt[v].nxt[t] != -1)
 University of Engineering and Technology
{ ll &v = st[node].nxt[edge];
vt[v].go[t] = vt[v].nxt[t]; ll t = s[st[v].f + dis - 1];
} if (v == 0)
else {
{ v = st.size();
if (v == 0) st.emplace_back(node, n - dis, base);
vt[v].go[t] = 0; st[last].link = node;
else last = 0;
vt[v].go[t] = go(get_link(v), t); }
} else if (c == t)
} {
return vt[v].go[t]; st[last].link = node;
} return;
ll get_link(ll v) }
{ else
if (vt[v].link == -1) {
{ ll u = st.size();
if (vt[v].par == 0 || v == 0) st.emplace_back(node, st[v].f, dis - 1);
vt[v].link = 0; st[u].nxt[c] = st.size();
else st.emplace_back(u, n - 1, base);
vt[v].link = go(get_link(vt[v].par), vt[v].ch - ’a’ + 1); st[u].nxt[t] = v;
} st[v].f += (dis - 1);
return vt[v].link; st[v].len -= (dis - 1);
} v = u;
ll get(string s) st[last].link = u;
{ last = u;
ll nw = 0; }
ll ans = 0; if (node == 0)
for (auto to : s) dis--;
{ else
nw = go(nw, to - ’a’ + 1); node = st[node].link;
ans += get_val(nw); }
} }
return ans; };
}
};

5.7 Z Function
5.6 Suffix Tree (Template by Tran Khoi Nguyen)
// string start from 1
// f[i] = longest prefix match with s[i...i + f[i] - 1]
struct tk
{ constexpr int N = 2e5 + 5;
map<ll, ll> nxt;
ll par, f, len; void Build(string &s, int n, int f[N]) // n = size of string, f = z array
ll link; {
tk(ll par = -1, ll f = 0, ll len = 0) : par(par), f(f), len(len) int l(1), r(1);
{
nxt.clear(); f[1] = n;
link = -1;
} for (int i = 2; i <= n; ++i)
}; if (r < i)
{
struct Suffix_Tree l = r = i;
{ while (r <= n && s[r - i + 1] == s[r])
vector<tk> st; ++r;
ll node; f[i] = r - i;
ll dis; --r;
s }
ll n; else if (f[i - l + 1] < r - i + 1)
vector<ll> s; f[i] = f[i - l + 1];
void init() else
{ {
st.clear(); l = i;
node = 0; while (r <= n && s[r - i + 1] == s[r])
dis = 0; ++r;
st.emplace_back(-1, 0, base); f[i] = r - i;
n = 0; --r;
} }
}
void go_edge()
{
while (dis > st[st[node].nxt[s[n - dis]]].len)
{
node = st[node].nxt[s[n - dis]];
dis -= st[node].len; 6 Data structures
}
}

void add_char(ll c) 6.1 Ordered Set

{
ll last = 0;
s.pb(c); #include <ext/pb_ds/assoc_container.hpp>
n = s.size(); #include <ext/pb_ds/tree_policy.hpp>
dis++; using namespace __gnu_pbds;
while (dis > 0)
{ #define ordered_set tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update>

21
go_edge();
ll edge = s[n - dis]; /*
 University of Engineering and Technology
order_of_key (k) : Number of items strictly smaller than k. line.pop_back();
find_by_order(k) : K-th element in a set (counting from zero). if (!line.empty())
*/ point.pop_back();
}
else
break;
}
6.2 Fenwick Tree (With Walk on tree) else
{
if (ff(i, line.back()) <= ff(i, line[line.size() - 2]))
{
// This is equivalent to calculating lower_bound on prefix sums array line.pop_back();
// LOGN = log2(N) if (!line.empty())
point.pop_back();
struct FenwickTree }
{ else
int n, LOGN; break;
ll a[N]; // BIT array }
}
FenwickTree()
{ if (line.empty() || A[line.back()] != A[i])
memset(a, 0, sizeof a); {
} if (!line.empty())
point.emplace_back(ff(line.back(), i));
void Update(int p, ll v) line.emplace_back(i);
{ }
for (; p <= n; p += p & -p) }
a[p] += v; ll Get(int x)
} {
int j = lower_bound(point.begin(), point.end(), x) - point.begin();
ll Get(int p) return A[line[j - 1]] * x + B[line[j - 1]];
{ }
ll ans(0); };

for (; p; p -= p & -p)

ans += a[p];

return ans;
} 6.4 Dynamic Convex Hull Trick (Min)
int search(ll v)
{
ll sum = 0; struct Line
int pos = 0; {
for (int i = LOGN; i >= 0; i--) mutable ll k, m, p;
{ bool operator<(const Line& o) const
if (pos + (1 << i) <= n && sum + a[pos + (1 << i)] < v) {
{ if (k==o.k) return m>o.m;
sum += a[pos + (1 << i)]; return k > o.k;
pos += (1 << i); }
} bool operator<(ll x) const
} {
return pos + 1; return p < x;
//+1 because pos will be position of largest value less than v }
} };
}; struct LineContainer : multiset<Line, less<>>
{
static const ll inf = LLONG_MAX;
ll div(ll a, ll b)
{
return a / b - ((a ^ b) < 0 && a % b);
6.3 Convex Hull Trick (Min) }
bool isect(iterator x, iterator y)
{
// If you want to get maximum, sort coef (A) not decreasing and change B[line.back()] > B[i] into B[ if (y == end())
line.back()] < B[i] return x->p = inf, 0;
if (x->k == y->k)
struct ConvexHullTrick x->p = x->m < y->m ? inf : -inf;
{ else
vector<ll> A, B; x->p = div(y->m - x->m, x->k - y->k);
vector<int> line; return x->p >= y->p;
vector<ld> point; }
ConvexHullTrick(int n = 0) void add(ll k, ll m)
{ {
A.resize(n + 2, 0); auto z = insert({k, m, 0}), y = z++, x = y;
B.resize(n + 2, 0); while (isect(y, z))
point.emplace_back(-Inf); z = erase(z);
} if (x != begin() && isect(--x, y))
isect(x, y = erase(y));
ld ff(int x, int y) while ((y = x) != begin() && (--x)->p >= y->p)
{ isect(x, erase(y));
return (ld)1.0 * (B[y] - B[x]) / (A[x] - A[y]); }
} ll query(ll x)
{
void Add(int i) assert(!empty());
{ auto l = *lower_bound(x);
while ((int)line.size() > 1 || ((int)line.size() == 1 && A[line.back()] == A[i])) return l.k * x + l.m;
{ }
if (A[line.back()] == A[i]) };
{

22
if (B[line.back()] > B[i])
{
 University of Engineering and Technology
6.5 SPlay Tree x = x->R;
}
}
struct KNode return nil;
{ }
int Value;
int Size; void Split(QNode x, int k, QNode &a, QNode &b)
KNode *P, *L, *R; {
}; if (k == 0)
using QNode = KNode *; {
KNode No_thing_here; a = nil;
QNode nil = &No_thing_here, root; b = x;
return;
void Link(QNode par, QNode child, bool Right) }
{ QNode cur = The_kth(x, k);
child->P = par; Up_to_Root(cur);
if (Right) a = cur;
par->R = child; b = a->R;
else a->R = nil;
par->L = child; b->P = nil;
} Update(a);
}
void Update(QNode &a)
{ QNode Join(QNode a, QNode b)
a->Size = a->L->Size + a->R->Size + 1; {
} if (a == nil)
return b;
void Init() while (a->R != nil)
{ a = a->R;
nil->Size = 0; Up_to_Root(a);
nil->P = nil->L = nil->R = nil; Link(a, b, true);
root = nil; Update(a);
for (int i = 1; i <= n; ++i) return a;
{ }
QNode cur = new KNode;
cur->P = cur->L = cur->R = nil; void Print(QNode &a)
cur->Value = i; {
Link(cur, root, false); if (a->L != nil)
root = cur; Print(a->L);
Update(root); cout << (a->Value) << " ";
} if (a->R != nil)
} Print(a->R);
}
void Rotate(QNode x)
{
QNode y = x->P;
QNode z = y->P;
if (x == y->L)
{
6.6 Hashing (Template by Tran Khoi Nguyen)
Link(y, x->R, false);
Link(x, y, true);
} struct Hashing
else {
{ vector<vector<vector<ll>>> f;
Link(y, x->L, true); vector<ll> mod;
Link(x, y, false); vector<vector<ll>> mu;
} vector<ll> chr;
Update(y); ll num;
Update(x); ll base;
x->P = nil; void init()
if (z != nil) {
Link(z, x, z->R == y); num = 2;
} f.clear();
mod.clear();
void Up_to_Root(QNode x) mu.clear();
{ chr.clear();
while (1){ vector<ll> vt = {999244369, 999254351, 999154309, 989154311, 989254411, 997254397, 991294387,
QNode y = x->P; 991814399, 994114351, 994914359, 994024333};
QNode z = y->P; random_shuffle(vt.begin(), vt.end());
if(y == nil) base = 317;
break; for (int i = 1; i <= 26; i++)
if(z != nil){ chr.pb(abs((ll)(rnd())));
if((x == y->L) == (y == z->L)) for (int i = 0; i < num; i++)
Rotate(y); {
else f.emplace_back();
Rotate(x); mod.pb(vt[i]);
} vector<ll> pt;
Rotate(x); pt.pb(1);
} for (int j = 1; j < maxn; j++)
} {
pt.pb((pt.back() * base) % mod[i]);
QNode The_kth(QNode x, int k) }
{ mu.pb(pt);
while (true) }
{ }
if (x->L->Size == k - 1) ll add(string s)
return x; {
if (x->L->Size >= k) ll n = s.length();
x = x->L; ll id = f[0].size();
else for (int j = 0; j < num; j++)

23
{ {
k -= x->L->Size + 1; vector<ll> vt1;
 University of Engineering and Technology
vt1.pb(0); void fake_update(ll x, ll y)
for (int i = 1; i <= n; i++) {
{ for (int i = lower_bound(vt.begin(), vt.end(), x) - vt.begin() + 1; i <= vt.size(); i += i &
vt1.pb((vt1.back() * base + chr[s[i - 1] - ’a’]) % mod[j]); (-i))
} {
f[j].pb(vt1); node[i].pb(y);
} }
return id; }
} void fake_get(ll x, ll y)
pll get_hash(ll id, ll l, ll r) {
{ for (int i = lower_bound(vt.begin(), vt.end(), x) - vt.begin() + 1; i; i -= i & (-i))
return make_pair(((f[0][id][r] - f[0][id][l - 1] * mu[0][r - l + 1]) % mod[0] + mod[0]) % mod {
[0], ((f[1][id][r] - f[1][id][l - 1] * mu[1][r - l + 1]) % mod[1] + mod[1]) % mod[1]); node[i].pb(y);
} }
}; }
void update(ll x, ll y, ll val)
{
for (int i = lower_bound(vt.begin(), vt.end(), x) - vt.begin() + 1; i <= vt.size(); i += i &
(-i))
6.7 BIT 2D (Template by Tran Khoi Nguyen) {
for (int j = lower_bound(node[i].begin(), node[i].end(), y) - node[i].begin() + 1; j <=
node[i].size(); j += j & (-j))
{
struct BIT_2D f[i][j] = f[i][j] + val;
{ }
vector<ll> vt; }
vector<vector<ll>> f; }
vector<vector<ll>> node; ll get(ll x, ll y)
void init(vector<ll> vt1) {
{ ll ans = 0;
vt = vt1; for (int i = lower_bound(vt.begin(), vt.end(), x) - vt.begin() + 1; i; i -= i & (-i))
sort(vt.begin(), vt.end()); {
vt.resize(unique(vt.begin(), vt.end()) - vt.begin()); for (int j = lower_bound(node[i].begin(), node[i].end(), y) - node[i].begin() + 1; j; j -=
f = vector<vector<ll>>(vt.size() + 2, vector<ll>(0)); j & (-j))
node = vector<vector<ll>>(vt.size() + 2, vector<ll>(0)); {
} ans += f[i][j];
void rearrange() }
{ }
for (int i = 1; i <= vt.size(); i++) return ans;
{ }
sort(node[i].begin(), node[i].end()); };
node[i].resize(unique(node[i].begin(), node[i].end()) - node[i].begin());
f[i] = vector<ll>(node[i].size() + 2, 0);
}
}

24