BZOJ 4259 残缺的字符串 ——FFT

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

#define maxn 1200005
#define F(i,j,k) for (int i=j;i<=k;++i)
#define D(i,j,k) for (int i=j;i>=k;--i)

const double pi=acos(-1.0);

struct Complex{
	double x,y;
	Complex operator + (Complex a){Complex b; return b.x=x+a.x,b.y=y+a.y,b;}
	Complex operator - (Complex a){Complex b; return b.x=x-a.x,b.y=y-a.y,b;}
	Complex operator * (Complex a){Complex b; return b.x=x*a.x-y*a.y,b.y=x*a.y+y*a.x,b;}
}a[maxn],b[maxn],c[maxn];

int A[maxn],B[maxn],la,lb,n,m=1,len,rev[maxn],ans[maxn],cnt=0;
char s1[maxn],s2[maxn];

void FFT(Complex * x,int n,int f)
{
	F(i,0,n-1) if (rev[i]>i) swap(x[rev[i]],x[i]);
	for (int m=2;m<=n;m<<=1)
	{
		int mid=m>>1;
		Complex wn; wn.x=cos(2.0*pi/m*f); wn.y=sin(2.0*pi/m*f);
		for (int i=0;i<n;i+=m)
		{
			Complex w; w.x=1.0; w.y=0;
			F(j,0,mid-1)
			{
				Complex a=x[i+j],b=x[i+j+mid]*w;
				x[i+j]=a+b; x[i+j+mid]=a-b;
				w=w*wn;
			}
		}
	}
}

int main()
{
//	freopen("in.txt","r",stdin);
	scanf("%d%d",&lb,&la);
	scanf("%s",s2);
	scanf("%s",s1);
	n=la+lb+1;
	while (m<=n) m<<=1,len++; n=m;
	F(i,0,n-1)
	{
		int t=i,ret=0;
		F(j,1,len) ret<<=1,ret|=t&1,t>>=1;
		rev[i]=ret;
	}
	F(i,0,la-1)	{A[i]=s1[i]-‘a‘+1; if (s1[i]==‘*‘) A[i]=0;}
	F(i,0,lb-1)	{B[i]=s2[lb-1-i]-‘a‘+1;	if (s2[lb-1-i]==‘*‘) B[i]=0;}

	F(i,0,n-1) a[i].x=A[i]*A[i]*A[i],a[i].y=0;
	F(i,0,n-1) b[i].x=B[i],b[i].y=0;
	FFT(a,n,1); FFT(b,n,1);
	F(i,0,n-1) c[i]=a[i]*b[i];
	
	F(i,0,n-1) a[i].x=A[i],a[i].y=0;
	F(i,0,n-1) b[i].x=B[i]*B[i]*B[i],b[i].y=0;
	FFT(a,n,1); FFT(b,n,1);
	F(i,0,n-1) c[i]=c[i]+a[i]*b[i];
	
	F(i,0,n-1) a[i].x=2*A[i]*A[i],a[i].y=0;
	F(i,0,n-1) b[i].x=B[i]*B[i],b[i].y=0;
	FFT(a,n,1); FFT(b,n,1);
	F(i,0,n-1) c[i]=c[i]-a[i]*b[i];
	
	FFT(c,n,-1);
	F(i,0,n-1) c[i].x=c[i].x/n;
//	F(i,0,n-1) printf("%.3f ",c[i].x); printf("\n");
	F(i,lb-1,la-1) if (fabs(c[i].x)<0.4) ans[++cnt]=i-lb+2;
	printf("%d\n",cnt);
	F(i,1,cnt) printf("%d%c",ans[i],i==cnt?‘\n‘:‘ ‘);
}