POJ 3261 Milk Patterns(後綴數組+二分答案)
題意
給定一個長度為 \(n?\) 的由正整數構成的串和一個整數 \(K?\),求至少出現 \(K?\) 次的串的最長長度(可以重疊)。
\(1 \leq n\leq 20000?\)
\(2\leq K\leq n\)
思路
變量及其含義
後綴數組(SA)是一個比較靈活的處理字符串有關問題的算法,它能對一個串的所有後綴進行排序。我們用串 abaaab 作例子,來分析這個算法。
sa(suffix array)
算法名由來的數組,可見該數組是核心。\(sa[i]\) 的定義是排第 \(i\) 的後綴開頭在哪個位置。
sa | suffix |
---|---|
3 | aaab |
4 | aab |
5 | ab |
1 | abaaab |
6 | b |
2 | baaab |
比如例串的 \(sa\) 數組如上表
rk(rank)
\(rk\) 數組是 \(sa\) 數組的逆運算,\(sa\) 是排名映射到位置,而 \(rk\) 則是位置映射到排名,即 \(rk[sa[i]]=i\) 。
H(height)
\(H\) 數組也就是常說的高度數組,它是在求解後綴數組問題時一個重要的數組,滿足 \(H[i]=\text{lcp}(sa[i],sa[i-1])\) ,其中 \(\text{lcp}\) 表示最長公共前綴。高度數組的含義也就是排名相鄰的兩個後綴的最長公共前綴長度。
這麽說可能不太形象,那麽我們把上面的表格再拓展一列。
sa | H | suffix |
---|---|---|
3 | / | aaab |
4 | 2 | aab |
5 | 1 | ab |
1 | 1 | abaaab |
6 | 0 | b |
2 | 1 | baaab |
具體流程
一般後綴數組使用 \(O(n\log n)\) 的倍增算法進行構造,當然有一個更加優秀的叫做 \(\text{DC3}?\) 的算法(我不會)是線性的,聽說較難實現。
基數排序
為了能實現較快的排序,有必要實現一個線性的排序算法,而基排就是一個相當靈活的線形排序算法。基數排序這裏不在詳細介紹,這裏僅僅把它挖的更深一點,方便理解後綴數組中的基排。我們以它的一次排序(其實就是計數排序)來演示。
比如直接把長度為 \(n?\) ,數值屬於 \([1,n]?\)
FOR(i,1,m)c[i]=0;
FOR(i,1,n)c[a[i]]++;
FOR(i,2,m)c[i]+=c[i-1];
DOR(i,n,1)b[c[a[i]]--]=a[i];
而直接按照下標放置(即把 \(a_i?\) 當作比較關鍵字,對一個\(1?\) 到 \(n?\) 的全排列進行排序),只用在第 \(4?\) 行把 \(a[]?\) 去掉即可。
FOR(i,1,m)c[i]=0;
FOR(i,1,n)c[a[i]]++;
FOR(i,2,m)c[i]+=c[i-1];
DOR(i,n,1)b[c[a[i]]--]=i;
再拓展一下,把一個長度為 \(n?\) 的 \(a?\) 數組,按比較關鍵字 \(f?\) 數組進行排序,放入 \(b?\) 數組中,其中 \(f_i\in[1,m]?\)。
FOR(i,1,m)c[i]=0;
FOR(i,1,n)c[f[a[i]]]++;
FOR(i,2,m)c[i]+=c[i-1];
DOR(i,n,1)b[c[f[a[i]]]--]=a[i];
利用基排屬於穩定排序的順序,可以在此基礎上繼續優先度更高的排序。
現在對基排的理解有沒有更深?九種排序什麽的要不僅僅會打板子,還要知道靈活運用才行。
倍增過程
回憶一下我們學過的倍增。如果可以依次求出 \(1,2,3,4\) 直到 \(k\) 的答案,那麽有時稍作修改,通過依次求 \(2^0,2^1,2^2,2^3\) 的答案,就可以一直到 \(2^k\) 的答案。
那麽構造後綴數組的倍增也是如此,我們想得到後綴的排序,把問題抽象一下,就是求從每個字符開始長度為 \(k?\) 的串的順序(\(k>=n?\) ,越界則該位置無窮小),我們把問題轉化成給 \(n?\) 個長為 \(k?\) 的串排序( \(k?\) 為 \(2?\) 的正整數次冪),通過長為 \(k?\) 的串的順序合並出長為 \(2k?\) 的串的順序。而長為 \(1?\) 則是直接對每個字符排序。
每次合並的過程就是給每第 \(i\) 個位置一個 \((rk[i],rk[i+k])\) 的二元組,按這個排序得到一個新的順序,就是長為 \(2k\) 的串的順序。
更加具體的流程只能通過代碼+註釋講解了。
memset(tmp,0,sizeof(tmp));
int *x=tmp[0],*y=tmp[1],*c=tmp[2]; //c是基數排序的計數數組
FOR(i,1,m)c[i]=0;
FOR(i,1,n)c[x[i]=s[i]]++;
FOR(i,2,m)c[i]+=c[i-1];
DOR(i,n,1)sa[c[x[i]]--]=i; //先解決2^0的排序
for(int k=1;k<=n;k<<=1)
{ //x在這裏表示長為k的串的rk,y在這裏表示二元組的後面那一維的順序
int p=0;
FOR(i,n-k+1,n)y[++p]=i; //i+k>n,故rk[i+k]為極小值
FOR(i,1,n)if(sa[i]>k)y[++p]=sa[i]-k;
FOR(i,1,m)c[i]=0; //對y數組以x為關鍵字排序,放在sa中
FOR(i,1,n)c[x[y[i]]]++;
FOR(i,2,m)c[i]+=c[i-1];
DOR(i,n,1)sa[c[x[y[i]]]--]=y[i];
std::swap(x,y); //其實這裏就是把x數組賦給y,並丟掉原來的x
p=x[sa[1]]=1;
FOR(i,2,n)x[sa[i]]=y[sa[i]]==y[sa[i-1]]&&y[sa[i]+k]==y[sa[i-1]+k]?p:++p; //x在這裏表示長為2k的串的rk,通過已經得到的順序sa算出
if(p==n)break; //如果rk已經兩兩不同
m=p;
}
FOR(i,1,n)rk[sa[i]]=i; //逆運算,得到rk
高度數組
- 定理:\(H[rk[i+1]]>=H[rk[i]]-1\)
用人類的語言闡述一下,就是在原字符串 \(i+1\) 位置在後綴數組對應位置中的 \(H\) 最少只可能比 \(i\) 的少 \(1\) 。
這個定理還是比較顯然的,設原字符串的 \(i\) 位置 \(\text{lcp}\) 為 \(l\) ,那麽 \(i+1\) 位置保底也有一個 \(\text{lcp}\) 長度為 \(l-1\) (就是與 \(i\) 位置形成 \(\text{lcp}\) 的那個後綴去掉開頭一位),至於有沒有更長就不知道了。
那我們只用按原串的順序,求出後綴數組對應位置的高度即可,通過 \(\text{k- -}\) 與 \(\text{k++}\) 的次數分析,可得復雜度 \(O(n)\) 。
int k=0;
FOR(i,1,n)
{
if(k)k--;
if(rk[i]==1)continue;
int j=sa[rk[i]-1];
while(i+k<=n&&j+k<=n&&s[i+k]==s[j+k])k++;
H[rk[i]]=k;
}
如何用SA解題
利用 \(sa,rk,H\) 三個關鍵的數組,將字符串問題轉化成這三個數組上的問題,是後綴數組解題的一般套路。
後綴的前綴就是子串,是後綴結構的核心原理。
我們此題題,求至少出現 \(K\) 次的串的最長長度為例進行分析。
首先求出後綴數組,至少出現 \(K\) 次,就意味著一段後綴數組上的一段長度大於 \(K\) 的區間,它們相鄰的 \(\text{lcp}\)( \(H\) 值)最大。那麽很自然的想到二分這個最長長度,檢查後綴數組上是否存在一個長度為 \(K\) 的區間滿足 \(H\) 的最小值大於等於這個長度即可。
代碼
#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#define FOR(i,x,y) for(int i=(x),i##END=(y);i<=i##END;++i)
#define DOR(i,x,y) for(int i=(x),i##END=(y);i>=i##END;--i)
template<typename T,typename _T>inline bool chk_min(T &x,const _T y){return y<x?x=y,1:0;}
template<typename T,typename _T>inline bool chk_max(T &x,const _T y){return x<y?x=y,1:0;}
typedef long long ll;
const int N=2e4+5;
int sa[N],rk[N],H[N],tmp[3][N];
int disc[N],D;
int s[N];
int n,K;
void get_SA(int m)
{
memset(tmp,0,sizeof(tmp));
int *x=tmp[0],*y=tmp[1],*c=tmp[2];
FOR(i,1,m)c[i]=0;
FOR(i,1,n)c[x[i]=s[i]]++;
FOR(i,2,m)c[i]+=c[i-1];
DOR(i,n,1)sa[c[x[i]]--]=i;
for(int k=1;k<=n;k<<=1)
{
int p=0;
FOR(i,n-k+1,n)y[++p]=i;
FOR(i,1,n)if(sa[i]>k)y[++p]=sa[i]-k;
FOR(i,1,m)c[i]=0;
FOR(i,1,n)c[x[y[i]]]++;
FOR(i,2,m)c[i]+=c[i-1];
DOR(i,n,1)sa[c[x[y[i]]]--]=y[i];
std::swap(x,y);
p=x[sa[1]]=1;
FOR(i,2,n)x[sa[i]]=y[sa[i]]==y[sa[i-1]]&&y[sa[i]+k]==y[sa[i-1]+k]?p:++p;
if(p==n)break;
m=p;
}
FOR(i,1,n)rk[sa[i]]=i;
int k=0;
FOR(i,1,n)
{
if(k)k--;
if(rk[i]==1)continue;
int j=sa[rk[i]-1];
while(i+k<=n&&j+k<=n&&s[i+k]==s[j+k])k++;
H[rk[i]]=k;
}
}
bool check(int len)
{
int cnt=0;
FOR(i,1,n)
{
cnt++;
if(cnt>=K)return true;
if(H[i+1]<len)cnt=0;
}
return false;
}
int main()
{
scanf("%d%d",&n,&K);
FOR(i,1,n)scanf("%d",&s[i]);
FOR(i,1,n)disc[++D]=s[i];
std::sort(disc+1,disc+D+1);
D=std::unique(disc+1,disc+D+1)-disc-1;
FOR(i,1,n)s[i]=std::lower_bound(disc+1,disc+D+1,s[i])-disc;
get_SA(n);
int l=0,r=n;
while(l<r)
{
int mid=(l+r+1)>>1;
if(check(mid))
l=mid;
else r=mid-1;
}
printf("%d\n",l);
return 0;
}
POJ 3261 Milk Patterns(後綴數組+二分答案)