1. 程式人生 > >FFT--bzoj4503: 兩個串

FFT--bzoj4503: 兩個串

傳送門

先不管 ? ? ,考慮兩個串的對應位置,如果 t t s s

i i 處出現則一定滿足
j = i
i + l e n t 1
( s [ j ] t [ j i ] ) 2 = 0 \sum_{j=i}^{i+lent-1}(s[j]-t[j-i])^2=0

這個拆開然後翻轉 t t 就是卷積的形式

再考慮 ? ? ,發現只要把 ? ? 的值賦為 0 0 ,然後再在這個式子外乘上 t [ j ] t[j] 就好了:
j = i i + l e n t 1 ( s [ j ] t [ j i ] ) 2 × t [ j ] = 0 \sum_{j=i}^{i+lent-1}(s[j]-t[j-i])^2\times t[j]=0
拆開以後用兩次 F F T FFT ,注意兩次的變數要不一樣(因為清空時候莫名其妙出了問題,改成不同的變數就對了)

程式碼如下:

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define maxn 400005
using namespace std;
char s[maxn],t[maxn];
int lens,lent,limit=1,l,rev[maxn],cnt;
double f[maxn],sum;
const double Pi=acos(-1.0),eps=1e-2;

struct complex{
	double x,y;
	complex(double xx=0,double yy=0){x=xx,y=yy;}
}a[maxn],b[maxn],c[maxn],d[maxn];
complex operator +(complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
complex operator -(complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
complex operator *(complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

inline void FFT(complex *F,int type){
	for(int i=0;i<limit;i++)
		if(i<rev[i]) swap(F[i],F[rev[i]]);
	for(int mid=1;mid<limit;mid<<=1){
		complex Wn(cos(Pi/mid),1.0*type*sin(Pi/mid));
		for(int r=mid<<1,j=0;j<limit;j+=r){
			complex w(1,0);
			for(int k=0;k<mid;k++,w=w*Wn){
				complex x=F[j+k],y=w*F[j+mid+k];
				F[j+k]=x+y,F[j+mid+k]=x-y;
			}
		}
	}
}

int main(){
	scanf("%s%s",s,t); lens=strlen(s); lent=strlen(t);
	for(int i=0;i<lens;i++) a[i].x=s[i]-'a'+1,a[i].x=a[i].x*a[i].x;
	for(int i=0;i<lent;i++)
		if(t[i]=='?') b[lent-i-1].x=0;
		else b[lent-i-1].x=t[i]-'a'+1;
	for(int i=0;i<lent;i++) sum+=b[i].x*b[i].x*b[i].x;
	while(limit<lens+lent) limit<<=1,++l;
	for(int i=0;i<limit;i++)
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
	FFT(a,1); FFT(b,1);
	for(int i=0;i<=limit;i++) a[i]=a[i]*b[i];
	FFT(a,-1);
	for(int i=lent-1;i<lens;i++)
		f[i-lent+1]+=((a[i].x/limit)+0.5);
	for(int i=0;i<lens;i++) c[i].x=s[i]-'a'+1,c[i].x=c[i].x*2.0;
	for(int i=0;i<lent;i++) 
		if(t[lent-i-1]=='?') d[i].x=0;
		else d[i].x=t[lent-i-1]-'a'+1,d[i].x=d[i].x*d[i].x;
	FFT(c,1); FFT(d,1);
	for(int i=0;i<=limit;i++) c[i]=c[i]*d[i];
	FFT(c,-1);
	for(int i=lent-1;i<lens;i++) {
		f[i-lent+1]-=(c[i].x/limit+0.5);
		f[i-lent+1]+=sum;
		if(f[i-lent+1]<=eps && f[i-lent+1]>=-eps) cnt++;
	}
	printf("%d\n",cnt);
	for(int i=lent-1;i<lens;i++)
		if(f[i-lent+1]<=eps && f[i-lent+1]>=-eps) printf("%d\n",i-lent+1);
	return 0;
}