bzoj4259: 残缺的字符串

Posted 2023-02-03 thy_asdf

传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=4259

思路：FFT是一种高效的字符串匹配算法

只要想到FFT就好办了

令为‘*’的位置值为0

那么两个字符串匹配当且仅当

Σ(a[i]-b[i])^2*a[i]*b[i]=0

式子的含义很好理解，前面就是字符相同，后面就是该位是否为通配符

为了FFT能做，我们要把A串反过来才有一个卷积的形式

设F[i]就是B串以i结尾的子串B[i-m+1,i]与A匹配得到的值

f[i]=Σ(a[j]-b[i-j])^2*a[j]*b[i-j]

    =Σ(a[j]^3*b[i-j]-2*(a[j])^2*(b[i-j])^2)+a[j]*b[i-j]^3)

分成三个部分，分别做FFT，加起来就是答案了

（公式写的比较挫，可以去看claris的博客http://www.cnblogs.com/clrs97/p/4814499.html）

#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
const int maxn=300010,maxt=(1<<19)+10;
using namespace std;
int n,m,ans[maxn],a[maxn],b[maxn];char s1[maxn],s2[maxn];
struct comdouble r,i;A[maxt],B[maxt],C[maxt];
com operator +(com a,com b)return (com)a.r+b.r,a.i+b.i;
com operator -(com a,com b)return (com)a.r-b.r,a.i-b.i;
com operator *(com a,com b)return (com)a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r;
com operator /(com a,int b)return (com)a.r/b,a.i/b;
com operator *(com a,int b)return (com)a.r*b,a.i*b;
double sqr(int x)return x*x;
double cube(int x)return x*x*x;

struct FFT
	int rev[maxt],n;
	void clear(com a[])memset(a,0,sizeof(*a)*(n+1));
	void init(int nn)
		for (n=1;n<nn;n<<=1);
		for (int i=1;i<n;i++) rev[i]=rev[i-(i&-i)]+(n>>1)/(i&-i);
	
	void fft(com a[],int op)
		for (int i=1;i<n;i++) if (rev[i]<i) swap(a[rev[i]],a[i]);
		for (int sz=2;sz<=n;sz<<=1)
			com w0=(com)cos(2*M_PI/sz),sin(2*M_PI*op/sz);
			for (int bg=0;bg<n;bg+=sz)
				com w=(com)1,0;
				for (int pos=bg;pos<bg+(sz>>1);pos++)
					com x=a[pos],y=w*a[pos+(sz>>1)];
					a[pos]=x+y,a[pos+(sz>>1)]=x-y;
					w=w*w0;
				
			
		
		if (op==-1) for (int i=0;i<n;i++) a[i]=a[i]/n;
	
T;

void init()
	scanf("%d%d",&m,&n);
	scanf("%s%s",s1,s2);
	for (int i=0;i<m;i++) a[i]=s1[i]=='*'?0:s1[i]-'a'+1;
	for (int i=0;i<n;i++) b[i]=s2[i]=='*'?0:s2[i]-'a'+1;
	for (int i=0;i<(m>>1);i++) swap(a[i],a[m-i-1]);
	//for (int i=0;i<m;i++) printf("%d ",(int)a[i]);puts("");
	//for (int i=0;i<n;i++) printf("%d ",(int)b[i]);puts("");


void work()
	T.init(n);
	
	T.clear(A),T.clear(B);
	for (int i=0;i<m;i++) A[i]=(com)a[i],0;
	for (int i=0;i<n;i++) B[i]=(com)cube(b[i]),0;
	T.fft(A,1),T.fft(B,1);
	for (int i=0;i<T.n;i++) C[i]=C[i]+A[i]*B[i];
	//for (int i=0;i<T.n;++i) printf("r=%.3lf i=%.3lf\\n",C[i].r,C[i].i);puts("");
	
	T.clear(A),T.clear(B);
	for (int i=0;i<m;i++) A[i]=(com)sqr(a[i]),0;
	for (int i=0;i<n;i++) B[i]=(com)sqr(b[i]),0;
	T.fft(A,1),T.fft(B,1);
	for (int i=0;i<T.n;i++) C[i]=C[i]-A[i]*B[i]*2;
	//for (int i=0;i<T.n;++i) printf("r=%.3lf i=%.3lf\\n",C[i].r,C[i].i);puts("");
	
	T.clear(A),T.clear(B);
	for (int i=0;i<m;i++) A[i]=(com)cube(a[i]),0;
	for (int i=0;i<n;i++) B[i]=(com)b[i],0;
	T.fft(A,1),T.fft(B,1);
	for (int i=0;i<T.n;i++) C[i]=C[i]+A[i]*B[i];
	//for (int i=0;i<T.n;++i) printf("r=%.3lf i=%.3lf\\n",C[i].r,C[i].i);puts("");
	
	T.fft(C,-1);
	int cnt=0;
	for (int i=0;i<n-m+1;i++)
		if (fabs(C[m-1+i].r)<1e-1)
			ans[++cnt]=i+1;
	printf("%d\\n",cnt);
	for (int i=1;i<cnt;i++) printf("%d ",ans[i]);
	printf("%d\\n",ans[cnt]);


int main()
	init(),work();
	return 0;