BZOJ.4892.[TJOI2017]DNA(后缀自动机/后缀数组)

题目链接

(Description)

给出两个串(S,T)，求(T)在(S)中出现了多少次。出现是指。可以有(3)次（(3)个字符）不匹配（修改使其匹配）。

(Solution)

一个套路的做法是构造多项式（CF528D），对每个字符c单独考虑，(f[i]=[S[i]可匹配c],g[i]=[T[i]==c])。
然后(F=f*g)，可以得到每个位置往后长(m)的串中有多少个位置(S,T)都匹配了(c)。如果某个位置匹配字符数(geq m-3)，则以它为左端点的串可行。
FFT/NTT实现，常数好也许能过。

SA做法：枚举(S)的每个位置(i)，设当前匹配(T)匹配到(j)，得到两个串的ht数组后我们可以(O(1))求出(LCP(suf[i],suf[j]))，直接(j+=LCP)就行了。
如果某个位置不匹配，可以至多用(3)次机会直接跳过去。所以实际枚举(j)的次数只有(5)。
复杂度(O(Tnlog n))。

SAM做法：得到parent树后，直接在上面DFS，如果能匹配则走，不能匹配则用一次次数。走了(m)步则加上该点的贡献(出现过多少次)。
复杂度(O(Tn))。

还有某种神奇的Hash做法。。好像复杂度比较优。

---

SAM：

//9224kb	1624ms
#include <cstdio>
#include <cstring>
#include <algorithm>
const int N=2e5+5;

struct Suffix_Automaton
{
	int n,Ans,tot,las,son[N][4],fa[N],len[N],cnt[N],tm[N],A[N],ref[233];
	char s[N];

	Suffix_Automaton() {tot=las=1;}
	void Insert(int c)
	{
		int np=++tot,p=las;
		len[las=np]=len[p]+1, cnt[np]=1;
		for(; p&&!son[p][c]; p=fa[p]) son[p][c]=np;
		if(!p) fa[np]=1;
		else
		{
			int q=son[p][c];
			if(len[q]==len[p]+1) fa[np]=q;
			else
			{
				int nq=++tot; len[nq]=len[p]+1;
				memcpy(son[nq],son[q],sizeof son[q]);
				fa[nq]=fa[q], fa[q]=fa[np]=nq;
				for(; son[p][c]==q; p=fa[p]) son[p][c]=nq;
			}
		}
	}
	void Build()
	{
		tot=las=1;
		ref['A']=0, ref['T']=1, ref['G']=2, ref['C']=3;
		memset(tm,0,sizeof tm);//! 你前缀和了→_→
		memset(cnt,0,sizeof cnt), memset(son,0,sizeof son);

		scanf("%s",s+1); int l=strlen(s+1);
		for(int i=1; i<=l; ++i) Insert(ref[s[i]]);
		for(int i=1; i<=tot; ++i) ++tm[len[i]];
		for(int i=1; i<=l; ++i) tm[i]+=tm[i-1];
		for(int i=1; i<=tot; ++i) A[tm[len[i]]--]=i;
		for(int i=tot,x=A[i]; i; x=A[--i]) cnt[fa[x]]+=cnt[x];
	}
	void DFS(int x,int use,int l)
	{
		if(l==n) return (void)(Ans+=cnt[x]);
		for(int i=0; i<4; ++i)
			if(son[x][i])
				if(ref[s[l]]==i) DFS(son[x][i],use,l+1);
				else if(use<3) DFS(son[x][i],use+1,l+1);
	}
	void Query()
	{
		scanf("%s",s), n=strlen(s);
		Ans=0, DFS(1,0,0), printf("%d
",Ans);
	}
}sam;

int main()
{
	int T; scanf("%d",&T);
	while(T--) sam.Build(), sam.Query();
	return 0;
}

---

SA：

//19768kb	5976ms(好慢...)
#include <cstdio>
#include <cstring>
#include <algorithm>
const int N=2e5+7;

int MAP[233],sa[N],sa2[N],rk[N],tm[N],ht[N],lg2[N],mn[18][N];
char s[N];

void Get_SA(int n)
{
	int *x=rk,*y=sa2,m=5;
	for(int i=0; i<=m; ++i) tm[i]=0;
	for(int i=1; i<=n; ++i) ++tm[x[i]=MAP[s[i]]];
	for(int i=1; i<=m; ++i) tm[i]+=tm[i-1];
	for(int i=n; i; --i) sa[tm[x[i]]--]=i;
	for(int k=1,p=0; k<n; k<<=1,m=p,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=0; i<=m; ++i) tm[i]=0;
		for(int i=1; i<=n; ++i) ++tm[x[i]];
		for(int i=1; i<=m; ++i) tm[i]+=tm[i-1];
		for(int i=n; i; --i) sa[tm[x[y[i]]]--]=y[i];

		std::swap(x,y), x[sa[1]]=p=1;
		for(int i=2; i<=n; ++i)
			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;
	}
	for(int i=1; i<=n; ++i) rk[sa[i]]=i;
	ht[1]=0;
	for(int i=1,k=0,p; i<=n; ++i)
	{
		if(rk[i]==1) continue;
		if(k) --k;
		p=sa[rk[i]-1];
		while(i+k<=n && p+k<=n && s[i+k]==s[p+k]) ++k;
		ht[rk[i]]=k;
	}
}
void Init_ST(int n)
{
	for(int i=1; i<=n; ++i) mn[0][i]=ht[i];
	for(int j=1; j<=lg2[n]; ++j)
		for(int i=1; i<=n; ++i)
			mn[j][i]=std::min(mn[j-1][i],mn[j-1][i+(1<<j-1)]);
}
inline int LCP(int l,int r)
{
	l=rk[l], r=rk[r]; if(l>r) std::swap(l,r);
	++l;
	int k=lg2[r-l+1];
	return std::min(mn[k][l],mn[k][r-(1<<k)+1]);
}

int main()
{
	MAP['A']=1, MAP['T']=2, MAP['C']=3, MAP['G']=4, MAP['Z']=5;
	lg2[1]=0;
	for(int i=2; i<=200005; ++i) lg2[i]=lg2[i>>1]+1;

	int T; scanf("%d",&T);
	while(T--)
	{
		int l,n;
		scanf("%s",s+1), s[l=strlen(s+1)+1]='Z';
		scanf("%s",s+l+1), n=strlen(s+1);
		Get_SA(n), Init_ST(n);
		int ans=0;
		for(int i=1,m=n-l,lim=l-m; i<=lim; ++i)
		{
			for(int j=1,t=0; t<=3; )
			{
				if(j>m) {++ans; break;}
				else if(s[i+j-1]!=s[l+j]) ++j, ++t;
				else j+=LCP(i+j-1,l+j);
			}
		}
		printf("%d
",ans);
	}
	return 0;
}