P4173 残缺的字符串

题意

先考虑没有残缺位置的情况:

先将两个字符匹配形式化:
定义(C(x,y)=(A_x-B_y)^2)
如果有匹配函数(P(x)=sumlimits_{i=0}^{m-1}C(i,x-m+i+1)=0),那么(B)(x)开始向后(m)个字符和(A)匹配。

之后将(A)翻转,得到(A'_i=A_{m-i-1})。于是匹配函数可以写成:(P(x)=sumlimits_{i=0}^{m-1}(A'_{m-i-1}-B_{x-m+i+1})^2)

展开这个式子:
(P(x)=sumlimits_{i=0}^{m-1}(A'_{m-i-1})^2+sumlimits_{i=0}^{m-1}(B_{x-m+i+1})^2-2sumlimits_{i=0}^{m-1}A'_{m-i-1}B_{x-m+i+1})
我们设(T=sumlimits_{i=0}^{m-1}(A'_i)^2,f(x)=sumlimits_{i=0}^{x}B_i^2,g(x)=sumlimits_{i+j=x}A'(i)B(j)),上式即为:
(=T+f(x)-f(x-m)-2g(x))

显然(T)(f(x))是可以预处理出来的,(g(x))是个卷积的形式,我们直接(FFT)就好了。

现在考虑有残缺位置的情况:

我们将残缺位置设为(0),之后重新定义两字符匹配:(C(x,y)=(A_x-B_y)^2A_xB_y)

与上面相同,我们开始推式子:
(A'_i=A_{m-i-1})
(P(x)=sumlimits_{i=0}^{m-1}(A'_{m-i-1}-B_{x-m+i+1})^2A'_{m-i-1}B_{x-m+i+1})
(=sumlimits_{i=0}^{m-1}(A'_{m-i-1})^3B_{x-m+i+1}+sumlimits_{i=0}^{m-1}A'_{m-i-1}(B_{x-m+i+1})^3-2sumlimits_{i=0}^{m-1}(A'_{m-i-1}B_{x-m+i+1})^2)
发现这三个式子全都是卷积的形式:
(=sumlimits_{i+j=x}(A'_i)^3B_j+sumlimits_{i+j=x}A'_i(B_j)^3-2sumlimits_{i+j=x}(A_i)^2(B_i)^2)

于是我们做(6)(NTT)最后再来一遍(INTT)即可。

code:

#pragma GCC optimize("Ofast")
#include<bits/stdc++.h>
using namespace std;
#define re register
typedef long long ll;
const int maxn=2400010;
const ll mod=1004535809;
const ll G=3;
const ll invG=334845270;
int n,m;
int a[maxn],b[maxn];
char s1[maxn],s2[maxn];
vector<int>ans;
inline int read()
{
	char c=getchar();re int res=0,f=1;
	while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
	while(c>='0'&&c<='9')res=res*10+c-'0',c=getchar();
	return res*f;
}
int lim=1,len;
int pos[maxn];
ll A[maxn],B[maxn],Ans[maxn];
inline ll power(ll x,ll k)
{
    ll res=1;
    while(k)
    {
        if(k&1)res=res*x%mod;
        x=x*x%mod;k>>=1;
    }
    return res;
}
inline void NTT(ll* a,int op)
{
    for(int i=0;i<lim;i++)if(i<pos[i])swap(a[i],a[pos[i]]);
    for(int l=1;l<lim;l<<=1)
    {
        ll wn=power(op==1?G:invG,(mod-1)/(l<<1));
        for(int i=0;i<lim;i+=l<<1)
        {
            ll w=1;
            for(int j=0;j<l;j++,w=w*wn%mod)
            {
                int x=a[i+j],y=w*a[i+l+j]%mod;
                a[i+j]=(x+y)%mod;a[i+l+j]=(x-y+mod)%mod;
            }
        }
    }
    if(op==1)return;
    ll inv=power(lim,mod-2);
    for(int i=0;i<lim;i++)a[i]=a[i]*inv%mod;
}
int main()
{
	m=read(),n=read();
	scanf("%s%s",s1,s2);
	reverse(s1,s1+m);
	for(re int i=0;i<m;i++)a[i]=(s1[i]=='*')?0:s1[i]-'a'+1;
	for(re int i=0;i<n;i++)b[i]=(s2[i]=='*')?0:s2[i]-'a'+1;
	while(lim<n*2)lim<<=1,len++;
	for(re int i=0;i<lim;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
	//first
	for(re int i=0;i<m;i++)A[i]=a[i]*a[i]*a[i];
	for(re int i=m;i<lim;i++)A[i]=0;
	for(re int i=0;i<n;i++)B[i]=b[i];
	for(re int i=n;i<lim;i++)B[i]=0;
	NTT(A,1);NTT(B,1);
	for(re int i=0;i<lim;i++)Ans[i]=(Ans[i]+A[i]*B[i]%mod)%mod;
	//second
	for(re int i=0;i<m;i++)A[i]=a[i];
	for(re int i=m;i<lim;i++)A[i]=0;
	for(re int i=0;i<n;i++)B[i]=b[i]*b[i]%mod*b[i]%mod;
	for(re int i=n;i<lim;i++)B[i]=0;
	NTT(A,1);NTT(B,1);
	for(re int i=0;i<lim;i++)Ans[i]=(Ans[i]+A[i]*B[i]%mod)%mod;
	//third
	for(re int i=0;i<m;i++)A[i]=a[i]*a[i]%mod;
	for(re int i=m;i<lim;i++)A[i]=0;
	for(re int i=0;i<n;i++)B[i]=b[i]*b[i]%mod;
	for(re int i=n;i<lim;i++)B[i]=0;
	NTT(A,1);NTT(B,1);
	for(re int i=0;i<lim;i++)Ans[i]=(Ans[i]-2*A[i]%mod*B[i]%mod+mod)%mod;
	NTT(Ans,-1);
	for(re int i=m-1;i<n;i++)if(!Ans[i])ans.push_back(i-m+2);
	printf("%d
",(int)ans.size());
	for(re unsigned int i=0;i<ans.size();i++)printf("%d ",ans[i]);
	return 0;
}
原文地址:https://www.cnblogs.com/nofind/p/12146564.html