作者: KSkun

后缀数组(倍增)算法原理及应用

后缀数组(倍增)算法原理及应用

概述

后缀数组是指对某一字符串的所有后缀按照字典序排序后的结果,这一信息常用于处理各种字符串问题,维护该信息的算法之一:倍增法的复杂度为O(n \log n)。尽管有更优秀的DC3法,但是倍增法更加通用好写,因此这里只介绍倍增法求后缀数组。

原理与实现

最初的思路

我们知道,如果处理出所有后缀,再使用快速排序等排序方式,复杂度是O(n^2 \log n)的。
我们希望优化复杂度,可以从排序算法上下手。比快排更优秀的排序方法也是有的,比如计数排序。计数排序需要消耗值域规模的空间,如果用在此处空间应该是O(n)的,不成问题。并且,计数排序并非基于比较的排序算法,之前复杂度中计入的O(n)比较字符串字典序现在也不需要了。那么,我们单次排序的复杂度降为O(n)了。
可是,计数排序需要有一个数值来代表每个字符串,基于数值来比较大小。对于代表字符串的数值,我们首先可以想到Hash。不过字符串长了Hash要取模无法比较大小不说,光是计算一遍Hash都需要O(n^2)时间,白白浪费了精力对排序进行优化。

倍增?

接下来是一个很神奇的操作,即倍增法的精髓所在:倍增。
我们可以考虑倍增关键字的长度,采用双关键字排序的方法。比如,当长度为1的时候,两个关键字就是每个后缀的第一个和第二个字符。每次对长度倍增,最多当长度变为 \displaystyle \lfloor \frac{n}{2} \rfloor的时候,所有排序就完成了。这样来看,复杂度就是O(n \log n)了。
接下来,我们要考虑具体的实现问题,比如,如何完成双关键字的计数排序。这一段,不放我们来直接阅读分析实现代码。

代码分析

建议在阅读下面的内容的时候配合后缀数组算法总结 | 远航休息栈一文中的若干示例,动手模拟一下,加深理解。
首先我们给出整个程序的代码。这里用struct简单封装了一下。

struct SuffixArray {
    char s[MAXN];
    int n, m, sa[MAXN], t[MAXN], t1[MAXN], c[MAXN];
    inline void init() {
        memset(this, 0, sizeof(SuffixArray));
    }
    inline void calsa() {
        int *x = t, *y = t1;
        for(int i = 1; i <= m; i++) c[i] = 0;
        for(int i = 1; i <= n; i++) c[x[i] = s[i]]++;
        for(int i = 1; i <= m; i++) c[i] += c[i - 1];
        for(int i = n; i >= 1; i--) sa[c[x[i]]--] = i;
        for(int k = 1; k <= n; k <<= 1) {
            int p = 0;
            for(int i = n - k + 1; i <= n; i++) y[++p] = i;
            for(int i = 1; i <= n; i++) if(sa[i] > k) y[++p] = sa[i] - k;
            for(int i = 1; i <= m; i++) c[i] = 0;
            for(int i = 1; i <= n; i++) c[x[y[i]]]++;
            for(int i = 1; i <= m; i++) c[i] += c[i - 1];
            for(int i = n; i >= 1; i--) sa[c[x[y[i]]]--] = y[i];
            std::swap(x, y);
            p = x[sa[1]] = 1;
            for(int i = 2; i <= n; i++) {
                x[sa[i]] = (y[sa[i - 1]] == y[sa[i]] && y[sa[i - 1] + k] == y[sa[i] + k]) 
                    ? p : ++p;
            }
            if(p >= n) break;
            m = p;
        }
    }
};
// Code by KSkun, 2018/5

变量声明

n代表字符串长度,m代表字符集大小。s是字符串本体,内容从下标1开始。
sa在排序中是一个排名数组,下标为排名,值为该排名对应的后缀首字母位置(后称“后缀编号”)。严格意义上说,这个数组实际上保存了每一次排序的结果,但如遇本次排序仍然存在两个后缀排序关键字相同,则长度大的(编号小的)后缀在前。SA计算完成后,该数组就是SA本体。
tt1是临时数组,用于存放后面的xy的信息。
c是计数排序的计数数组,下标为计数的值,值为该值出现的次数。排序过程中,我们会对其求前缀和以便计算排名。
x是一个中间量数组,意义为得到第一关键字的大小,对于一次排序,下标为代表后缀编号,值为象征对应后缀编号第一关键字大小的值(事实上可以视作排名)。
y是一个中间量数组,意义为第二关键字的排名,对于一次排序,下标为排名,值为第二关键字对应的后缀的编号。

初始化

int *x = t, *y = t1;
for(int i = 1; i <= m; i++) c[i] = 0;
for(int i = 1; i <= n; i++) c[x[i] = s[i]]++;
for(int i = 1; i <= m; i++) c[i] += c[i - 1];
for(int i = n; i >= 1; i--) sa[c[x[i]]--] = i;

第1/8行分配tt1两个内存空间给x数组和y数组。
第2/9行清空用于计数排序的计数数组c
第3/10行,统计该次第一关键字(首字母)的出现次数。
第4/11行,计算前缀和。
第5/12行,计算第一次排序的sa。此次排序,仅区分首字母和后缀编号,首字母小的在前,首字母相同,后缀编号小(后缀长度大)的在前。

倍增:计算y数组

int p = 0;
for(int i = n - k + 1; i <= n; i++) y[++p] = i;
for(int i = 1; i <= n; i++) if(sa[i] > k) y[++p] = sa[i] - k;

第2/15行,有一些后缀拆分后的第二关键字为空,这个循环将这些后缀的排名提到最前,以表示空字符的最高优先级。
第3/16行,对于满的第二关键字排序,利用某一后缀的后缀也是原串的后缀这一特点,由上一次的sa计算而来。考虑那些排名靠前的后缀,如果往前补倍增长度位,那么这个后缀对应的长度前缀就变成了第二关键字。如我们有后缀aaabaa,长度为2,计算abaa的第二关键字的时候,实际上查到的是aa的排名,并且往前补了2位。

倍增:计数排序计算x数组

for(int i = 1; i <= m; i++) c[i] = 0;
for(int i = 1; i <= n; i++) c[x[y[i]]]++;
for(int i = 1; i <= m; i++) c[i] += c[i - 1];
for(int i = n; i >= 1; i--) sa[c[x[y[i]]]--] = y[i];

第2/18行,依然是在统计各个量的出现次数。y[i]表示的是第二关键字排名为i的后缀编号,而x[i]拿到的则是第一关键字排序时,后缀编号为i的后缀对应的象征第一关键字大小的值(在初始化过程中,我们将该值设置为首字母的值了,因此,对初始化还有一种理解,即y[i] = i,并没对第二关键字排序)。因此x[y[i]]的值反映的是代表第一关键字大小的值。
第4/20行,对i逆序统计,其用意是:当第一关键字相同的时候,第二关键字排名大的,计算后sa的值也会更大。

倍增:计算下一次的x数组

std::swap(x, y);
p = x[sa[1]] = 1;
for(int i = 2; i <= n; i++) {
    x[sa[i]] = (y[sa[i - 1]] == y[sa[i]] && y[sa[i - 1] + k] == y[sa[i] + k]) ? p : ++p;
}

第1/21行,把旧的x数组值换到y数组上。
第2/22行,规定当前排名第一的后缀x值为1,以后的后缀x值不会小于1。
第4/24行,考虑什么时候一个后缀和前面的后缀下一次的第一关键字(=这一次的第一+第二关键字,因为是倍增)相同,只有这一次的第一第二关键字都相等才可以判定下一次的第一关键字相等,否则应该比上一次的大。

倍增:边界条件

if(p >= n) break;
m = p;

第1/27行,当本次排序中已经产生了不小于n种不同的第一关键字的值,排序便没有必要继续下去了,因为已经区分开了这些后缀。
第2/28行,本次排序更新了关键字的种类数,因此要更新m值。

例题:【P3809】【模板】后缀排序 – 洛谷

把上面的板子打上去就可以了。

// Code by KSkun, 2018/5
#include <cstdio>
#include <cctype>
#include <cstring>

#include <algorithm>

const int MAXN = 1000005;

struct SuffixArray {
    char s[MAXN];
    int n, m, sa[MAXN], t[MAXN], t1[MAXN], c[MAXN];
    inline void init() {
        memset(this, 0, sizeof(SuffixArray));
    }
    inline void calsa() {
        int *x = t, *y = t1;
        for(int i = 1; i <= m; i++) c[i] = 0;
        for(int i = 1; i <= n; i++) c[x[i] = s[i]]++;
        for(int i = 1; i <= m; i++) c[i] += c[i - 1];
        for(int i = n; i >= 1; i--) sa[c[x[i]]--] = i;
        for(int k = 1; k <= n; k <<= 1) {
            int p = 0;
            for(int i = n - k + 1; i <= n; i++) y[++p] = i;
            for(int i = 1; i <= n; i++) if(sa[i] > k) y[++p] = sa[i] - k;
            for(int i = 1; i <= m; i++) c[i] = 0;
            for(int i = 1; i <= n; i++) c[x[y[i]]]++;
            for(int i = 1; i <= m; i++) c[i] += c[i - 1];
            for(int i = n; i >= 1; i--) sa[c[x[y[i]]]--] = y[i];
            std::swap(x, y);
            p = x[sa[1]] = 1;
            for(int i = 2; i <= n; i++) {
                x[sa[i]] = (y[sa[i - 1]] == y[sa[i]] && y[sa[i - 1] + k] == y[sa[i] + k]) 
                    ? p : ++p;
            }
            if(p >= n) break;
            m = p;
        }
    }
} sa;

int main() {
    scanf("%s", sa.s + 1);
    sa.n = strlen(sa.s + 1); sa.m = 256;
    sa.calsa();
    for(int i = 1; i <= sa.n; i++) {
        printf("%d ", sa.sa[i]);
    }
    return 0;
}

应用

height

上文提到过,x数组的值实际上是后缀编号对应的排名,我们将其内容复制一份进入rnk数组待用。
下面给出height数组的定义:代表该后缀与前一排名后缀的最长公共前缀(LCP)长度,下标为后缀编号。以下将该数组简写为hei数组。
我们只需要对每个后缀找到紧接着后面排名的后缀,来统计相同的部分即可。需要注意的是,当我们向后枚举排名的时候,实际上可以看做比较的后缀各删去了第一位,可以用上次的k-1来作为基础,然后再比较前k-1位之后有没有更多的相同的部分,这样可以优化时间。
以下代码就是求hei数组的代码,紧接着求sa的代码。

memcpy(rnk, x, sizeof(rnk));
int k = 0;
for(int i = 1; i <= n; i++) {
    int j = sa[rnk[i] + 1];
    if(k) k--;
    if(!j) continue;
    while(s[i + k] == s[j + k]) k++;
    hei[rnk[i]] = k;
}

求前缀信息

把字符串逆序一下扔进去求个SA就好了。

求前缀的LCS或后缀的LCP

要求两个排名不相邻的后缀的LCP怎么办?我们可以把上面的hei数组建一个ST表,查询排名区间的最小值即可。hei数组本身的含义就是LCP长度,这个性质有一定的传递性。在这里不太能说的明白,倒是可以采用观察一个sa的例子验证一下的做法。
至于前缀的LCS,其实你发现逆序了以后就变成一个后缀LCP问题了。
这个问题的一个经典应用题目是:[NOI2016]优秀的拆分

参考资料

[WC2011]最大XOR和路径 题解

[WC2011]最大XOR和路径 题解

题目地址:洛谷:【P4151】[WC2011]最大XOR和路径 – 洛谷、BZOJ:Problem 2115. — [Wc2011] Xor

题目描述

2606 1 - [WC2011]最大XOR和路径 题解

题意简述

给你一个n点m边的无向图,求一条1到n的路径使边权异或和最大。

输入输出格式

输入格式:
第一行包含两个整数 N 和 M, 表示该无向图中点的数目与边的数目。
接下来 M 行描述 M 条边,每行三个整数Si,Ti,Di,表示 Si 与 Ti 之间存在一条权值为 Di 的无向边。
图中可能有重边或自环。

输出格式:
仅包含一个整数,表示最大的XOR和(十进制结果),注意输出后加换行回车。

输入输出样例

输入样例#1:

5 7
1 2 2
1 3 2
2 4 1
2 5 1
4 5 3
5 3 4
4 3 2 

输出样例#1:

6

说明

2606 3 - [WC2011]最大XOR和路径 题解

题解

首先,我们考虑异或和最大的路径应该是怎么样的。我们考虑图中的环,显然我们可以从任意一个节点走到环上,绕环一圈,再走回来,这样只有这个环的权值异或和就被异或了进去。那么我们可以DFS处理一遍环的权值异或和来得到一个集合。
具体来说,DFS如果遇到一条边两端的点都已经访问过,那么就可以将dis[u]^dis[v]^边权记作环的权值异或和,这是由于dis[u]^dis[v]会将不在环上的边通过一次异或消去。
我们发现,我们只需要用一条1到n的路径再组合上这个集合的某一些权值使得异或和最大就行了。对于这个集合的处理,我们可以求线性基,这样就可以代表所有可能的情况的异或值。我们只需要任选一条1到n的路径并且使用线性基计算最大值就可以了。
任选一条路径为什么是正确的?我们考虑1到n的路径假如有2条,那么这两条路径会构成环,任选其中一条,再异或上这个环,实际上得到了另一条路径的异或和,而环的情况刚刚已经处理过,我们可以保证答案最优。

代码

// Code by KSkun, 2018/5
#include <cstdio>
#include <cctype>

#include <algorithm>
#include <vector>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    register char c = fgc();
    while(!isdigit(c)) {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(isdigit(c)) {
        res = (res << 1) + (res << 3) + c - '0';
        c = fgc();
    }
    return res * neg;
}

const int MAXN = 50005;

int n, m;

struct Edge {
    int to;
    LL w;
};

std::vector<Edge> gra[MAXN];

bool vis[MAXN];
LL dis[MAXN];
std::vector<LL> cir;

inline void dfs(int u) {
    vis[u] = true;
    for(int i = 0; i < gra[u].size(); i++) {
        int v = gra[u][i].to;
        if(vis[v]) {
            cir.push_back(dis[v] ^ dis[u] ^ gra[u][i].w);
        } else {
            dis[v] = dis[u] ^ gra[u][i].w;
            dfs(v);
        }
    }
}

LL mat[65];

int main() {
    n = readint(); m = readint();
    int u, v; LL w;
    for(int i = 1; i <= m; i++) {
        u = readint(); v = readint(); w = readint();
        gra[u].push_back(Edge {v, w});
        gra[v].push_back(Edge {u, w});
    }
    dfs(1);
    for(int i = 0; i < cir.size(); i++) {
        for(int j = 63; j >= 0; j--) {
            if(cir[i] & (1ll << j)) { // 这里不是1ll会爆int
                if(!mat[j]) {
                    mat[j] = cir[i];
                    break;
                } else {
                    cir[i] ^= mat[j];
                }
            }
        }
    }
    LL ans = dis[n];
    for(int i = 63; i >= 0; i--) {
        ans = std::max(ans, ans ^ mat[i]);
    }
    printf("%lld", ans);
    return 0;
}
[NOI2010]航空管制 题解

[NOI2010]航空管制 题解

题目地址:洛谷:【P1954】[NOI2010]航空管制 – 洛谷、BZOJ:Problem 2535. — [Noi2010]Plane 航空管制2

题目描述

世博期间,上海的航空客运量大大超过了平时,随之而来的航空管制也频频发生。最近,小X就因为航空管制,连续两次在机场被延误超过了两小时。对此,小X表示很不满意。
在这次来烟台的路上,小X不幸又一次碰上了航空管制。于是小X开始思考关于航空管制的问题。
假设目前被延误航班共有n个,编号为1至n。机场只有一条起飞跑道,所有的航班需按某个顺序依次起飞(称这个顺序为起飞序列)。定义一个航班的起飞序号为该航班在起飞序列中的位置,即是第几个起飞的航班。
起飞序列还存在两类限制条件:

  • 第一类(最晚起飞时间限制):编号为i的航班起飞序号不得超过ki;
  • 第二类(相对起飞顺序限制):存在一些相对起飞顺序限制(a, b),表示航班a的起飞时间必须早于航班b,即航班a的起飞序号必须小于航班b的起飞序号。

小X思考的第一个问题是,若给定以上两类限制条件,是否可以计算出一个可行的起飞序列。第二个问题则是,在考虑两类限制条件的情况下,如何求出每个航班在所有可行的起飞序列中的最小起飞序号。

输入输出格式

输入格式:
输入文件plane.in第一行包含两个正整数n和m,n表示航班数目,m表示第二类限制条件(相对起飞顺序限制)的数目。
第二行包含n个正整数k1, k2, …, kn。
接下来m行,每行两个正整数a和b,表示一对相对起飞顺序限制(a, b),其中1≤a,b≤n, 表示航班a必须先于航班b起飞。

输出格式:
输出文件plane.out由两行组成。
第一行包含n个整数,表示一个可行的起飞序列,相邻两个整数用空格分隔。输入数据保证至少存在一个可行的起飞序列。如果存在多个可行的方案,输出任意一个即可。
第二行包含n个整数t1, t2, …, tn,其中ti表示航班i可能的最小起飞序号,相邻两个整数用空格分隔。

输入输出样例

输入样例#1:

5 5
4 5 2 5 4
1 2
3 2
5 1
3 4
3 1

输出样例#1:

3 5 1 4 2
3 4 1 2 1

输入样例#2:

5 0
3 3 3 5 5

输出样例#2:

3 2 1 5 4
1 1 1 4 4

说明

【样例说明】
在样例1 中:
起飞序列3 5 1 4 2满足了所有的限制条件,所有满足条件的起飞序列有:
3 4 5 1 2 3 5 1 2 4 3 5 1 4 2 3 5 4 1 2
5 3 1 2 4 5 3 1 4 2 5 3 4 1 2
由于存在(5, 1)和(3, 1)两个限制,航班1只能安排在航班5和3之后,故最早起飞时间为3,其他航班类似。
在样例2 中:
虽然航班4、5没有相对起飞顺序限制,但是由于航班1、2、3都必须安排在前3个起飞,所以4、5最早只能安排在第4个起飞。
【数据范围】
对于30%数据:n≤10;
对于60%数据:n≤500;
对于100%数据:n≤2,000,m≤10,000。

题解

第一问,考虑把限制条件建个图,一定是一个DAG,直接进行拓扑排序,只不过拓扑的队列要改成优先队列,优先级按最晚时间限制从小到大排序,这样出来的顺序就是正确的了。
但是考虑第二问的时候,发现就算建反图DFS还可能出现其他需要满足的条件,例如样例2中前三个飞机必须在前三飞。此时我们考虑倒着找,建反图按限制从大到小的顺序拓扑,我们枚举每个飞机u,每次都拓扑一遍,遇到u的时候直接把它设置为不可以被加入队列中,这样就可以无视掉u的前置,直到条件不能被满足位置。倒着拓扑的意义,实质上是在求满足了u以后剩下最多的点。用n减去这个数字就得到了我们想要的答案。
这个做法的复杂度是O(n^2 \log n)的,在洛谷得开O2才能跑过。

代码

// Code by KSkun, 2018/5
#include <cstdio>
#include <cctype>
#include <cstring>

#include <algorithm>
#include <queue>
#include <vector>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF 
        : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    register char c = fgc();
    while(!isdigit(c)) {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(isdigit(c)) {
        res = (res << 1) + (res << 3) + c - '0';
        c = fgc();
    }
    return res * neg;
}

const int MAXN = 20005;

std::vector<int> gra[MAXN];
int deg[MAXN];

inline void addedge(int u, int v) {
    gra[u].push_back(v); deg[v]++;
}

int n, m, lim[MAXN];

struct Node {
    int u, lim;
    inline bool operator<(const Node &rhs) const {
        return lim < rhs.lim;
    }
};

int now[MAXN];
int ans[MAXN];

inline void toposort() {
    std::priority_queue<Node> pq;
    memcpy(now, deg, sizeof(deg));
    for(int i = 1; i <= n; i++) {
        if(!now[i]) pq.push(Node {i, lim[i]});
    }
    int cnt = n;
    while(!pq.empty()) {
        int u = pq.top().u; pq.pop();
        ans[cnt--] = u;
        for(int i = 0; i < gra[u].size(); i++) {
            int v = gra[u][i];
            if(!--now[v]) {
                pq.push(Node {v, lim[v]});
            }
        }
    }
}

inline int toposort1(int uu) {
    std::priority_queue<Node> pq;
    memcpy(now, deg, sizeof(deg));
    now[uu] = n;
    for(int i = 1; i <= n; i++) {
        if(!now[i]) pq.push(Node {i, lim[i]});
    }
    for(int i = n; i; i--) {
        if(pq.empty() || pq.top().lim < i) return i;
        int u = pq.top().u; pq.pop();
        for(int i = 0; i < gra[u].size(); i++) {
            int v = gra[u][i];
            if(!--now[v]) pq.push(Node {v, lim[v]});
        }
    }
}

int main() {
    n = readint(); m = readint();
    for(int i = 1; i <= n; i++) {
        lim[i] = readint();
    }
    for(int i = 1, u, v; i <= m; i++) {
        u = readint(); v = readint();
        addedge(v, u);
    }
    toposort();
    for(int i = 1; i <= n; i++) {
        printf("%d ", ans[i]);
    }
    printf("\n");
    for(int i = 1; i <= n; i++) {
        printf("%d ", toposort1(i));
    }
    return 0;
}