常用算法模板
大量早期代码，码风很差 ......
有些东西放这儿只是为了 CF 复制方便。
倍增后缀数组，SA Doubling
这里贴了新的板子，包含了求 height(LCP) 和 ST 表维护任意 LCP 的过程，稍微写一下主要的问题。
\(sa_i\) 是第 \(i\) 名的后缀，\(rk_i\) 是 \({\rm suffix}(i)\) 的排名。
radixSort() 是基数排序过程，第一关键字是 \(rk\)，第二关键字是 \(tp\)，其中 \(tp_i\) 表示第二关键字第 \(i\) 大的后缀是谁，排序时只需要按第一关键字做前缀和，第二关键字从后往前加入即可。
先按字符排序，随后倍增。按照当前 SA 的顺序加入第二关键字，排序后计算新的 rank。注意第一和第二关键字都相同算作一样，如果两者中任一或全部都没有第二关键字，视为不同，因为长度不同。
计算 height 就十分简单了，按照原串中的顺序依次计算，性质显然。
为了使用 ST 表，需要全局预处理出 \(\log_2 x\)。
复杂度 \(\mathcal O(n \log n)\)
int lg[N];
void preLog2(){
for(int i=2; i<Len; i<<=1) lg[i] = 1;
for(int i=3; i<Len; ++i) lg[i] += lg[i - 1];
}
struct SuffixArray{
char s[Len];
int sa[Len], rk[Len], ht[Len], st[lgLen][Len], tp[Len];
int n, m;
void radixSort(){
static int pos[N];
fill_n(pos, m + 1, 0);
for(int i=0; i<n; i++) ++pos[rk[i]];
for(int i=1; i<=m; i++) pos[i] += pos[i - 1];
for(int i=n-1; i>=0; i--) sa[--pos[rk[tp[i]]]] = tp[i];
}
void buildSA(){
m = 26;
for(int i=0; i<n; i++) rk[i] = s[i] - 'a', tp[i] = i;
radixSort();
for(int w = 1, p = 0; p + 1 < n; m = p, w <<= 1){
p = 0;
for(int i=n-w; i<n; i++) tp[p++] = i;
for(int i=0; i<n; i++)
if(sa[i] >= w) tp[p++] = sa[i] - w;
radixSort();
copy_n(rk, n, tp);
rk[sa[0]] = 0; p = 0;
for(int i=1; i<n; i++)
if(tp[sa[i]] == tp[sa[i - 1]] &&
sa[i] + w < n && sa[i-1] + w < n &&
tp[sa[i] + w] == tp[sa[i - 1] + w])
rk[sa[i]] = p;
else rk[sa[i]] = ++p;
}
int p = 0;
for(int i=0; i<n; i++){
if(p) --p;
if(rk[i] == 0) continue;
int j = sa[rk[i] - 1];
while(s[i + p] == s[j + p]) ++p;
ht[rk[i]] = p;
}
}
void buildST(){
for(int i=0; i<n; i++) st[0][i] = ht[i];
for(int i=1; i<=lg[n]; i++)
for(int j=0; j+(1<<i)-1 < n; j++)
st[i][j] = min(st[i-1][j], st[i-1][j + (1 << (i-1))]);
}
int lcp(int x, int y){
if(x == y) return n - x;
x = rk[x], y = rk[y];
if(x > y) swap(x, y);
++x;
int k = lg[y - x + 1];
return min(st[k][x], st[k][y - (1 << k) + 1]);
}
}
exgcd
int exgcd(int a, int b, int &x, int &y){
if(b == 0){
x = 1; y = 0;
return a;
}
int result = exgcd(b, a%b, x, y);
int tp = x; x = y; y = tp - a/b*y;
return result;
}
tarjan 缩点
void tarjan(int x){
low[x] = dfn[x] = ++id;
stk[++top] = x; instk[x] = true;
for(int i=hd[x]; i; i=nt[i]){
if(dfn[to[i]] == 0) tarjan(to[i]), low[x] = min(low[x], low[to[i]]);
else if(instk[to[i]]) low[x] = min(low[x], dfn[to[i]]);
}
if(low[x] == dfn[x]){
++color; col[x] = color;
while(stk[top] != x) instk[stk[top]] = false, col[stk[top--]] = color;
top--; instk[x] = false;
}
}
tarjan 割点/点双连通分量
以下代码统计每个 BCC 中的割点数量，摘自题目Luogu P3225 [HNOI2012]矿场搭建
void tarjan(int x, int fa){
dfn[x] = low[x] = ++ID;
stk[++top] = x;
int ch = 0, last = 0;
for(int i=hd[x]; i; i=nt[i]){
int v = to[i];
if(!dfn[v]){
++ch;
tarjan(v, x);
low[x] = min(low[x], low[v]);
if(x == fa && ch >= 2){
cut[x] = true;
if(ch == 2 && last != 0) ++cutCnt[last]; //是根且有2个或以上孩子，是割点
}
if(low[v] >= dfn[x]){ //是根或没有返回边，得到一个点双
if(x != fa) cut[x] = true; //不是根就是割点
++num; last = num; cnt[num] = cutCnt[num] = 0;
while(stk[top] != v) cutCnt[num] += (int)(cut[stk[top--]]), ++cnt[num];
cutCnt[num] += (int)(cut[stk[top--]]); ++cnt[num]; //注意加入v
cutCnt[num] += (int)cut[x]; ++cnt[num]; //注意加入自己
}
}
else if(v != fa)
low[x] = min(low[x], dfn[v]);
}
}
ISAP 最大流
template<class cap_t, class val_t, size_t node_c, size_t edge_c>
struct maxflow_graph {
static_assert(sizeof(val_t) >= sizeof(cap_t), "result may overflow");
struct edge {
int v;
cap_t f;
} E[edge_c * 2];
vector<int> G[node_c];
size_t cur[node_c];
int ec = 0;
int s, t, nc;
int h[node_c], gap[node_c + 1];
void clear() {
ec = 0;
memset(cur, 0, sizeof(cur));
for (size_t i = 0; i < node_c; ++i) G[i].clear();
}
void link(int x, int y, cap_t flow) {
if (!flow) return;
G[x].push_back(ec); E[ec++] = {y, flow};
G[y].push_back(ec); E[ec++] = {x, 0};
}
bool bfs() {
memset(h, -1, sizeof(int) * nc);
h[t] = 0;
int *q = new int[node_c], *ql = q - 1, *qr = q;
*qr = t;
while (ql != qr) {
int x = *(++ql);
++gap[h[x]];
for (int i : G[x]) {
if (h[E[i].v] == -1 && !E[i].f)
h[E[i].v] = h[x] + 1, *(++qr) = E[i].v;
}
}
delete[] q;
return h[s] != -1;
}
val_t augment(int x, val_t flow) {
if (x == t) return flow;
val_t res = 0;
for (size_t &i = cur[x]; i != G[x].size(); ++i) {
edge &e = E[G[x][i]];
if (!e.f || h[e.v] + 1 != h[x]) continue;
cap_t aug = min<val_t>(e.f, flow);
aug = augment(e.v, aug);
res += aug; flow -= aug;
e.f -= aug; E[G[x][i] ^ 1].f += aug;
if (!flow) return res;
}
--gap[h[x]];
if (!gap[h[x]]) h[s] = nc;
++gap[++h[x]];
cur[x] = 0;
return res;
}
val_t maxflow() {
val_t flow = 0;
if (!bfs()) return 0;
while (h[s] < nc) flow += augment(s, numeric_limits<val_t>::max());
return flow;
}
};
最小费用最大流
又长又慢
建议写 Dijkstra 版，不要用这个！
struct MinCostMaxFlow{
int to[Edge], nt[Edge], hd[Node], f[Edge], w[Edge], h[Node], cur[Node], fromE[Node], ec = 0;
int S, T;
void clear(){memset(hd, -1, sizeof(hd)); ec = 0;}
MinCostMaxFlow(){clear();}
void link(int x, int y, int flow, int cost){
to[ec] = y; nt[ec] = hd[x]; hd[x] = ec; f[ec] = flow; w[ec] = cost; ++ec;
to[ec] = x; nt[ec] = hd[y]; hd[y] = ec; f[ec] = 0; w[ec] = -cost; ++ec;
}
bool spfa(){
deque<int> q;
memset(h, 0x3f, sizeof(h));
h[S] = 0; q.push_back(S);
while(!q.empty()){
int x = q.front(); q.pop_front();
for(int i=hd[x]; i!=-1; i=nt[i]){
if(f[i] && h[to[i]] > h[x] + w[i]){
h[to[i]] = h[x] + w[i];
fromE[to[i]] = i;
if(!q.empty() && h[to[i]] <= h[q.front()]) q.push_front(to[i]);
else q.push_back(to[i]);
}
}
}
return h[T] != 0x3f3f3f3f;
}
bool vis[Node];
int dfs(int x, int flow){
if(x == T || !flow) return flow;
int res = 0;
vis[x] = true;
for(int i=hd[x]; i!=-1; i=nt[i]){
if(vis[to[i]] || h[to[i]] != h[x] + w[i] || !f[i]) continue;
int aug = min(flow, f[i]);
if((aug = min(aug, dfs(to[i], aug)))){
f[i] -= aug; f[i^1] += aug;
flow -= aug; res += aug;
}
if(!flow) break;
}
vis[x] = false;
if(!res) h[x] = 0x3f3f3f3f;
return res;
}
pair<int, int> mcmf(){
int flow = 0, cost = 0;
while(spfa()){
int fl = dfs(S, 0x3f3f3f3f);
flow += fl; cost += fl * h[T];
}
return make_pair(flow, cost);
}
}G;
回文自动机 PAM
namespace pam{
char s[N];
int c[N][26], fa[N], len[N], cnt[N], last, pos = 0, nc = 0;
void init() {
len[0] = 0, len[1] = -1;
fa[0] = 1;
s[0] = -1;
nc = 1;
last = 0;
}
int getfail(int x) {
while (s[pos] != s[pos - len[x] - 1])
x = fa[x];
return x;
}
void extend(int x) {
s[++pos] = x;
int now = getfail(last);
if(!c[now][x]) {
++nc;
len[nc] = len[now] + 2;
fa[nc] = c[getfail(fa[now])][x];
c[now][x] = nc;
}
last = c[now][x];
++cnt[last];
}
}; // namespace pam
模素数组合数 / 二项式系数
int fac[N], ifac[N];
void prefac(int n) {
fac[0] = 1;
for (int i = 1; i <= n; ++i) fac[i] = mul(fac[i - 1], i);
ifac[n] = inv(fac[n]);
for (int i = n - 1; i >= 0; --i) ifac[i] = mul(ifac[i + 1], i + 1);
}
int binom(int n, int m) {
return mul(fac[n], ifac[m], ifac[n - m]);
}
SAIS 线性后缀数组
封装了求 LCP，但注意 sais() 函数需要保证字符串是 int 类型，且最后存在极小字符，最后极小字符被排在开头。
（即 0 下标变 1 下标）
struct suffix_array {
int sa[N], pos[N];
int st[logN][N];
const int *ht = st[0];
void sais(const int *a, int *sa, int n, int w) {
bool *type = new bool[n];
int *buc = new int[w], *pl = new int[w], *pr = new int[w];
int n1 = 0, *p1;
type[n - 1] = true;
for (int i = n - 1; i; --i) {
type[i - 1] = a[i - 1] == a[i] ? type[i] : a[i - 1] < a[i];
n1 += type[i] && !type[i - 1];
}
p1 = new int[n1];
for (int i = 1, *p = p1; i < n; ++i)
if (type[i] && !type[i - 1]) *(p++) = i;
memset(buc, 0, sizeof(int) * w);
for (int i = 0; i < n; ++i)
++buc[a[i]];
for (int i = 0, sum = 0; i < w; ++i) {
sum += buc[i];
pr[i] = sum; pl[i] = sum - buc[i];
}
auto induce = [&]() {
memcpy(buc, pl, sizeof(int) * w);
for (int i = 0; i < n; ++i)
if (int j = sa[i] - 1; j >= 0 && !type[j]) sa[buc[a[j]]++] = j;
memcpy(buc, pr, sizeof(int) * w);
for (int i = n - 1; i >= 0; --i)
if (int j = sa[i] - 1; j >= 0 && type[j]) sa[--buc[a[j]]] = j;
};
memset(sa, -1, sizeof(int) * n);
memcpy(buc, pr, sizeof(int) * w);
for (int i = 1; i < n; ++i)
if (type[i] && !type[i - 1]) sa[--buc[a[i]]] = i;
induce();
int id = 0;
int *s1 = new int[n1], *sa1 = new int[n1];
for (int i = 0, *p = sa1; i < n; ++i)
if (type[sa[i]] && sa[i] && !type[sa[i] - 1]) *(p++) = sa[i];
sa[sa1[0]] = 0;
for (int i = 0; i < n1 - 1; ++i) {
int p = sa1[i], q = sa1[i + 1];
if (a[p] != a[q]) {
sa[q] = ++id;
continue;
}
++p, ++q;
while (a[p] == a[q] && type[p] == type[q] && !(type[p] && !type[p - 1]) && !(type[q] && !type[q - 1]))
++p, ++q;
id += !((type[p] && !type[p - 1]) && (type[q] && !type[q - 1]));
sa[sa1[i + 1]] = id;
}
for (int i = 0; i < n1; ++i)
s1[i] = sa[p1[i]];
if (id + 1 == n1)
for (int i = 0; i < n1; ++i)
sa1[s1[i]] = i;
else
sais(s1, sa1, n1, id + 1);
memcpy(buc, pr, sizeof(int) * w);
memset(sa, -1, sizeof(int) * n);
for (int i = n1 - 1; i >= 0; --i)
sa[--buc[a[p1[sa1[i]]]]] = p1[sa1[i]];
induce();
delete[] type; delete[] buc; delete[] pl; delete[] pr; delete[] p1; delete[] s1; delete[] sa1;
}
void init(char *s) {
int n = strlen(s);
int *a = new int[n + 1];
copy(s, s + n, a);
a[n] = 0;
sais(a, sa, n + 1, 128);
delete[] a;
for (int i = 1; i <= n; ++i)
pos[sa[i]] = i;
int h = 0;
for (int i = 0; i < n; ++i) {
h = !h ? 0 : h - 1;
int j = sa[pos[i] - 1];
while (s[i + h] == s[j + h]) ++h;
st[0][pos[i]] = h;
}
for (int i = 1, li = __lg(n); i <= li; ++i)
for (int j = 1, li = n - (1 << i) + 1; j <= li; ++j)
st[i][j] = min(st[i - 1][j], st[i - 1][j + (1 << (i - 1))]);
}
int lcp(int i, int j) {
i = pos[i], j = pos[j];
if (i > j) swap(i, j);
int k = __lg(j - i);
return min(st[k][i + 1], st[k][j - (1 << k) + 1]);
}
}
复数 FFT
using cp = complex<double>;
const double pi = acos(-1);
cp root[N];
cp MUL(const cp &a, const cp &b) {
double x = a.real(), y = a.imag(), u = b.real(), v = b.imag();
return cp(x * u - y * v, x * v + y * u);
}
void init_roots() {
for (int i = N / 2; i; i >>= 1) {
double theta = pi / i;
for (int j = 0; j < i; ++j)
root[i + j] = polar(1.0, j * theta);
}
}
void DIF(cp *a, int n) {
for (int l = n, k = l / 2; k; l >>= 1, k >>= 1) {
cp *w = root + k;
for (cp *i = a, *t = a + n; i != t; i += l)
for (int j = 0; j < k; ++j) {
cp t = i[j];
i[j] += i[j + k];
i[j + k] = MUL((t - i[j + k]), w[j]);
}
}
}
void DIT(cp *a, int n) {
for (int l = 2, k = 1; k < n; l <<= 1, k <<= 1) {
cp *w = root + k;
for (cp *i = a, *t = a + n; i != t; i += l)
for (int j = 0; j < k; ++j) {
cp t = MUL(i[j + k], w[j]);
i[j + k] = i[j] - t;
i[j] += t;
}
}
double iv = 1.0 / n;
reverse(a + 1, a + n);
for (int i = 0; i < n; ++i)
a[i] *= iv;
}
待补充...
浙公网安备 33010602011771号