[TJOI2017]DNA，洛谷P3763，暗藏玄机的字符串FFT

最新推荐文章于 2023-01-11 00:02:18 发布

原创最新推荐文章于 2023-01-11 00:02:18 发布 · 535 阅读

1 ·

CC 4.0 BY-SA版权

数论专栏收录该内容

49 篇文章

订阅专栏

本文探讨了如何使用快速傅立叶变换(FFT)算法来解决字符串匹配问题，尤其是在字符集不限于{A,T,C,G}

正题

这题正解SAM？蒟蒻表示不懂。

拿起了FFT乱搞。

首先，题目没有说字符集只有{A,T,C,G};四种。

然后，我们就可以考虑FFT。

用 $F[i]$ 表示偏转i位之后，文本串和模式串的差异。

那么就有 $F[i]=\sum_{j=0}^{m-1}[a[i+j]!=b[j]]$

然后很明显，我们可以先令 $B[n-1-j]=b[j]$ ， $f[n+i-1]=F[i]$ 。

就可以写成 $f[n+i-1]=\sum_{j=0}^{m-1}[a[i+j]!=B[n-1-j]]$ 。

是不是很像卷积？

这是我们将分别令 ${A,T,C,G}$ 为1。

假设赋值后的数组是 $x,y$ 。

$\\f[n+i-1]=\sum_{j=0}^{m-1}(x[i+j]-y[n-1-j])^2 \\=\sum_{j=0}^{m-1}x[i+j]^2+y[n-1-j]^2-2*x[i+j]*y[n-1-j]$

前面两个值都是固定的，我们只需要用FFT求出最后的一项就可以了。

加上常数就是 $O(4*3n\log n+qn)$ ，q是一个大常数。

所以没办法，只能开 $O_2$ 来过

最后求出来的f实际上是差异的两倍，因为不同的两个字符各会被互相算一次。

所以最后和6比就好了。

#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
using namespace std;

const int maxn=1e5;
struct complex{
	double x,y;
	complex operator+(complex a)const {return (complex){x+a.x,y+a.y};}
	complex operator-(complex a)const {return (complex){x-a.x,y-a.y};}
	complex operator*(complex a)const {return (complex){x*a.x-y*a.y,x*a.y+y*a.x};}
}f[3*maxn+10],g[3*maxn+10];
int where[3*maxn+10];
int T;
char s[maxn+10];
int a[maxn+10],b[maxn+10],sum[maxn+10];
int n,m,l,lim;
int h[maxn+10];
const double Pi=acos(-1.0)*2.0;

int check(char c){
	if(c=='A') return 1;
	if(c=='T') return 2;
	if(c=='C') return 3;
	if(c=='G') return 4;
}

void dft(complex *now,int idft){
	for(int i=0;i<lim;i++) if(i>where[i]) swap(now[i],now[where[i]]);
	complex wn,w,a,b;
	for(int l=2;l<=lim;l*=2){
		wn=(complex){cos(Pi/l),sin(Pi/l)*idft};
		for(int i=0;i<lim;i+=l){
			w=(complex){1,0};
			for(int x=i,y=i+l/2;y<i+l;x++,y++,w=w*wn){
				a=now[x],b=w*now[y];
				now[x]=a+b;
				now[y]=a-b;
			}
		}
	}
}

int main(){
	scanf("%d",&T);
	while(T--){
		scanf("%s",s);n=strlen(s);
		for(int i=0;i<n;i++) a[i]=check(s[i]);
		scanf("%s",s);m=strlen(s);
		for(int i=0;i<m;i++) b[m-1-i]=check(s[i]);
		memset(h,0,sizeof(h));l=0,lim=1;
		while(lim<n+m-1) lim*=2,l++;
		for(int i=0;i<lim;i++) where[i]=((where[i>>1]>>1)|((i&1)<<(l-1)));
		for(int i=1;i<=4;i++){
			for(int j=0;j<lim;j++) f[j].x=f[j].y=g[j].x=g[j].y=0;
			for(int j=0;j<n;j++) {
				if(j!=0) sum[j]=sum[j-1];
				else sum[j]=0;
				if(a[j]==i) f[j].x=1,sum[j]++;
			}
			sum[n]=0;
			for(int j=0;j<m;j++) if(b[j]==i) g[j].x=1,sum[n]++;
			dft(f,1);dft(g,1);
			for(int j=0;j<lim;j++) f[j]=f[j]*g[j];
			dft(f,-1);
			for(int j=0;j<=n-m;j++) h[j]+=sum[j+m-1]-(j>=1?sum[j-1]:0)+sum[n];
			for(int j=m-1;j<=n-1;j++) 
				h[j-m+1]-=2*(int)(f[j].x/lim+0.5);
		}
		int ans=0;
		for(int j=0;j<=n-m;j++) ans+=(h[j]<=6);
		printf("%d\n",ans);
	}
}