POJ 3261 Milk Patterns 后缀数组

Openjudge原题链接

• 题意
  输入一个长度N(<=20000)的数组，求出其重复K次的最长可重叠子串
• 思路
  由于N可以达到20000，故考虑O(NlogN)的算法，于是想到后缀数组。
  假设取出了数组的全部后缀，那么重复K次的最长可重叠子串就是在全部后缀中出现K次的最长前缀。所以需要将后缀用倍增方法排序，并计算出相邻串的公共前缀长度（h数组），用RMQ分别查找[1...K],[2...K+1]...[N-K...N-1]的区间最小值，最大的一个就是重复了K次的可重叠子串长度。
• 难点在于求后缀数组和使用RMQ。
• 后缀数组的求法：sa[i]表示排第i的字符串是哪个，用rk[i] 第i个字符串排第几（相同的字符串有一样的排名）。根据所有后缀的前n个字符计算出rk数组，根据rk数组更新sa数组。然后根据后缀的前2n个字符计算出新的rk数组.
  为什么可以这样做呢？因为suffix(i,2n)=suffix(i,n)+suffix(i+n,n)。
  故前2n个字符的排名可以看成2元组(rk[i],rk[i+n])。用这个二元组进行排序即可。
  这样，每次计算rk和sa是O(n)的时间，而一共进行了O(logn)趟，故计算rk和sa耗时O(nlogn)。
• 然后，求h数组 h[k]=LCP(suffix[sa[k]],suffix[sa[k+1]])，这时候再一次利用后缀字符串的性质。
  假设我们已经求得h[k]=LCP(suffix[sa[k]],suffix[sa[k+1]]),假设sa[k]=i,sa[k+1]=j。
  那么LCP(suffix(i+1),suffix(j+1))=max(h[k]-1,0).找到i+1在sa中的位置t(t=rk(i+1)),sa[t+1]=p.
  即h[t]=LCP(suffix(i+1),suffix(p))
  而suffix(i)<suffix(j) => suffix(i+1)<suffix(j+1)
  而suffix(p)紧邻suffix(i+1)，故有
  suffix(i+1)<suffix(p)<=suffix(j+1)
  所以有h[t]>=max(h[k]-1,0),继续匹配则得到h[t]
• RMQ是一种用O(nlogn)时间预处理，O(1)时间求区间最小值的数据结构，用dp实现，其空间复杂度为O(nlogn)。
  用dp[i][j]表示[i,i+2^j-1]的区间最小值，那么转移方程就是
  dp[i][j]=min(dp[i][j-1],dp[i+(1<<(j-1))][j-1])
  那么为了查询区间[l,r]的最小值，只需要找到最大的k,s.t. (1<<k)<=(r-l+1)，则
  ans=min(dp[l][k],dp[r-(1<<k)+1][k].要注意这里的边界条件.
• AC代码

/*
    输入一个长度N(<=20000)的数组，求出其重复K次的最长可重叠子串
    由于N可以达到20000，故考虑O(NlogN)的算法，于是想到后缀数组。
    假设取出了数组的全部后缀，那么重复K次的最长可重叠子串就是在全
    部后缀中出现K次的最长前缀。所以需要将后缀用倍增方法排序，并计
    算出相邻串的公共前缀长度（h数组），用RMQ分别查找[1...K],
    [2...K+1]...[N-K...N-1]的区间最小值，最大的一个就是重复了K
    次的可重叠子串长度。
*/
#include<iostream>
#include<algorithm>
#include<cstring>
#include<string>
using namespace std;
const int N = 1e5+5;
const int LOGN = 20;
const int M = 1e6 + 5;
int a[N], sa[N], rk[N], h[N];// sa[i] 排第i的字符串是哪个; rk[i] 第i个字符串排第几
int cnt[M];
int t1[N], t2[N];//tmp array
int dp[N][LOGN];
int n, K;
void build_sa(int* s) {
    memset(cnt, 0, sizeof(cnt));
    int m = 0;
    for (int i = 1; i <= n; i++) {
        m = max(m, a[i]);//m for max
        cnt[s[i]]++;
        t1[i] = s[i];
    }
    for (int i = 1; i <= m; i++) {
        cnt[i] += cnt[i - 1];
    }
    int t = 1;
    for (int i = 1; i <= n; i++) {
        rk[i] = cnt[t1[i]];
    }
    for (int i = 1; i <= n; i++) {
        sa[cnt[t1[i]]--] = i;//初始化sa
    }
    /*a[0] = -1;
    for (int i = 1; i <= n; i++) {//初始化rk
        if (a[sa[i]] != a[sa[i - 1]])
            rk[sa[i]] = rk[sa[i - 1]] + 1;
        else rk[sa[i]] = rk[sa[i - 1]];
    }*/

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

        for (int i = 0; i <= m; i++) cnt[i] = 0;
        for (int i = 1; i <= n; i++) cnt[t1[i]]++;
        for (int i = 1; i <= m; i++) cnt[i] += cnt[i - 1];
        /* 双关键字排序 */
        for (int i = n; i >= 1; i--) sa[cnt[t1[t2[i]]]--] = t2[i];//wtf,根据rk(t1)更新sa,但这到底是在做什么！！

        for (int i = 1; i <= n; i++) swap(t1[i], t2[i]);//t2记录原先的t1
        p = 0; t1[sa[1]] = ++p;
        /* 计算t1为新的rk，即以前2^k个字符排序得到的rk */
        for (int i = 2; i <= n; i++)
            t1[sa[i]] = (t2[sa[i]] == t2[sa[i - 1]] &&
                t2[sa[i] + k] == t2[sa[i - 1] + k]) ? p : ++p;
        m = p; if (m >= n) break;//优化策略
    }

    for (int i = 1; i <= n; i++) rk[sa[i]] = i;
    /* 根据 rk,sa 计算h数组 */
    /* h[k]=LCP(suffix[sa[k-1]],suffix[sa[k]]) */
    int p = 0;
    for (int i = 1; i <= n; i++)
    {
        if (rk[i] == 1) p = 0;
        else
        {
            if (p) p--;
            while (i + p <= n && sa[rk[i] - 1] + p <= n
                && s[i + p] == s[sa[rk[i] - 1] + p]) p++;
        }
        h[rk[i]] = p;
    }
    /*求h数组 h[k]=LCP(suffix[sa[k]],suffix[sa[k+1]]) */
    /*int o, b;
    for (int i = 0; i < n; i++) {
        o = rk[i];//在sa数组中的位置
        if (o == n) continue;
        b = sa[o + 1];//b为要和i匹配的后缀的位置
        while (i+h[o]<=n&&b+h[o]<=n
            &&s[i + h[o]] == s[b + h[o]]) ++h[o];
        h[rk[i + 1]] = max(0, h[o] - 1);
    }*/
    /*
    for (int i = 1; i <= n; i++) {
        cout << rk[i] << " ";
    }
    cout << endl;
    for (int i = 1; i <= n; i++) {
        cout << h[i] << " ";
    }
    cout << endl;
    for (int i = 1; i <= n; i++) {
        for (int j = sa[i]; j <= n; j++) {
            cout << a[j];
        }
        cout << endl;
    }
    cout << endl;*/
}
/* dp[i][j]表示[i,i+2^j-1]的区间最小值              */
/* dp[i][j]=min(dp[i][j-1],dp[i+(1<<(j-1))][j-1])   */
void rmq_init() {
    int len = n;
    for (int i = 1; i <= len; i++)
        dp[i][0] = h[i];
    for (int j = 1; (1 << j) <= len; j++) {
        for (int i = 1; i + (1 << j) - 1 <= len; i++) {
            dp[i][j] = min(dp[i][j - 1], dp[i + (1 << (j - 1))][j - 1]);
        }
    }
}
/* 查询h[l...r]上的区间最小值    */      
int rmq(int l, int r) {
    int k = 0;
    while ((1 << (k + 1)) <= r - l + 1)
        k++;
    int ans = min(dp[l][k], dp[r - (1 << k) + 1][k]);
    return ans;
}
void solve() {
    rmq_init();
    int ans = -1;
    for (int i = 1; i <= n - K+1; i++) {
        ans = max(ans, rmq(i+1, i + K - 1));//h[i+1]...h[i+K-1]
    }
    cout << ans << endl;
}
int main() {
    cin >> n >> K;
    for (int i = 1; i <= n; i++) {
        cin >> a[i];
    }
    build_sa(a);
    solve();
    return 0;
}