BZOJ4503：两个串 （构造+FFT/随机化+NTT）

权限题传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=4503

---

题目大意：给你两个串S和T，询问T在S中出现了几次，以及在哪些位置出现。S中只包含小写字母，T中包含小写字母和’?’，’?’可以匹配任何字符。|T|<=|S|<=100000。

---

题目分析：26倍常数必定TLE。我们不妨先假设T中没有’?’，那么如何判断T和S中i开头的一段是否匹配？
for j=0->|T|-1 S[i+j]==T[j]
换句话说就是要求所有S[i+j]-T[j]==0。所以我们构造数组A，使得：

$$A[i]=\sum_{j=0}^{|T|-1}(S[i+j]-T[j])^2$$

如果A[i]为0，就意味着S中以i开头的串和T匹配。那么T中有’?’怎么办呢？我们不妨将’?’的权值设为0，然后重新构造：

$$A[i]=\sum_{j=0}^{|T|-1}(S[i+j]-T[j])^2T[j]$$

把这个式子裂项，发现有两项下标之差相等，另外一项是个定值。于是用两次FFT即可求出。时间复杂度$O(n\log(n))$。

---

上面的做法是网上众多dalao们用的方法，一眼看上去就感觉十分套路。然而我一开始接这题的时候并没有想到这种方法，我想了一个随机化算法。令M=1004535809（或者其它NTT模数），然后给’a’~’z’分别随机一个[1,M)的权值，再令’?’的权值为0（为什么是0，因为0乘以任何数都为0，它等于0的平方）。然后直接令$A[i]=\sum_{j=0}^{|T|-1}S[i+j]*T[j]$，这个用NTT算。如果最后A[i]等于$\sum_{j=0}^{|T|-1}T[j]^2$，那么就认为匹配成功。这样随机k次，只有每一次都匹配成功才能算进最终的答案（BZOJ上k=3的时候可以刚好卡进时限）。然后这个算法我自己用随机数据对拍了10min都没错，并且成功在B站上AC。至于其正确率我也不会算，我只是单纯地感觉从[1,M)中随机26个权值，冲突的概率应该是很小的，如果有路过的dalao知道怎么算的话欢迎提出QAQ。

---

CODE（构造+FFT）：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
using namespace std; 

const double pi=acos(-1.0);
const int maxn=1000000;

struct Complex
{
    double X,Y;
    Complex (double a=0.0,double b=0.0) : X(a),Y(b) {}
} ;

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);}

Complex A[maxn];
Complex B[maxn];

int Rev[maxn];
int N,Lg;

char S[maxn];
int ans[maxn];
int num=0;

int s[maxn];
int t[maxn];
int ls,lt;

void DFT(Complex *a,double f)
{
    for (int i=0; i<N; i++)
        if (i<Rev[i]) swap(a[i],a[ Rev[i] ]);

    for (int len=2; len<=N; len<<=1)
    {
        int mid=(len>>1);
        double ang=2.0*pi/((double)len);
        Complex e( cos(ang) , f*sin(ang) );

        for (Complex *p=a; p!=a+N; p+=len)
        {
            Complex wn(1.0,0.0);
            for (int i=0; i<mid; i++)
            {
                Complex temp=wn*p[mid+i];
                p[mid+i]=p[i]-temp;
                p[i]=p[i]+temp;
                wn=wn*e;
            }
        }
    }
}

void FFT()
{
    DFT(A,1.0);
    DFT(B,1.0);
    for (int i=0; i<N; i++) A[i]=A[i]*B[i];
    DFT(A,-1.0);
    for (int i=0; i<N; i++) A[i].X/=((double)N); 
}

int main()
{
    freopen("4503.in","r",stdin);
    freopen("4503.out","w",stdout);

    scanf("%s",&S);
    ls=strlen(S);
    for (int i=0; i<ls; i++) s[i]=S[i]-'a'+1;
    scanf("%s",&S);
    lt=strlen(S);
    for (int i=0; i<lt; i++)
        if ( 'a'<=S[i] && S[i]<='z' ) t[i]=S[i]-'a'+1;
        else t[i]=0;

    N=1,Lg=0;
    while (N<ls+lt+2) N<<=1,Lg++;
    for (int i=0; i<N; i++)
        for (int j=0; j<Lg; j++)
            if (i&(1<<j)) Rev[i]|=(1<<(Lg-j-1));

    for (int i=0; i<ls; i++) A[i]=Complex((double)s[i]*s[i],0.0);
    for (int i=ls; i<N; i++) A[i]=Complex(0.0,0.0);
    for (int i=0; i<lt; i++) B[lt-i-1]=Complex((double)t[i],0.0);
    for (int i=lt; i<N; i++) B[i]=Complex(0.0,0.0);
    FFT();
    for (int i=0; i<=ls-lt; i++) ans[i]+=(int)floor( A[lt+i-1].X+0.5 );

    for (int i=0; i<ls; i++) A[i]=Complex((double)s[i],0.0);
    for (int i=ls; i<N; i++) A[i]=Complex(0.0,0.0);
    for (int i=0; i<lt; i++) B[lt-i-1]=Complex((double)t[i]*t[i],0.0);
    for (int i=lt; i<N; i++) B[i]=Complex(0.0,0.0);
    FFT();
    for (int i=0; i<=ls-lt; i++) ans[i]-=(2*(int)floor( A[lt+i-1].X+0.5 ));

    int temp=0;
    for (int i=0; i<lt; i++) temp+=(t[i]*t[i]*t[i]);
    for (int i=0; i<=ls-lt; i++) if (ans[i]+temp==0) num++;
    printf("%d\n",num);
    for (int i=0; i<=ls-lt; i++) if (ans[i]+temp==0) printf("%d\n",i);

    return 0;
}

CODE（随机化+NTT）：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
using namespace std; 

const long long M=1004535809;
const long long g=3;
const int maxn=1000000;
typedef long long LL;

const long long M1=998244353;
const long long M2=1000000007;
LL seed=0;

LL A[maxn];
LL B[maxn];

int Rev[maxn];
int N,Lg;

LL val[30];
bool vis[maxn];
int num=0;

char S[maxn];
int s[maxn];
int t[maxn];
int ls,lt;

LL Rand()
{
    seed=(seed*M1+M2)%M;
    return seed;
}

LL Pow(LL x,LL y)
{
    if (!y) return 1;
    LL temp=Pow(x,y>>1);
    temp=temp*temp%M;
    if (y&1) temp=temp*x%M;
    return temp;
}

void DFT(LL *a,int f)
{
    for (int i=0; i<N; i++)
        if (i<Rev[i]) swap(a[i],a[ Rev[i] ]);

    for (int len=2; len<=N; len<<=1)
    {
        int mid=(len>>1);
        LL e=Pow(g,(M-1)/len);
        if (f==-1) e=Pow(e,M-2);

        for (LL *p=a; p!=a+N; p+=len)
        {
            LL wn=1;
            for (int i=0; i<mid; i++)
            {
                LL temp=wn*p[mid+i]%M;
                p[mid+i]=(p[i]-temp+M)%M;
                p[i]=(p[i]+temp)%M;
                wn=wn*e%M;
            }
        }
    }
}

void NTT()
{
    DFT(A,1);
    DFT(B,1);
    for (int i=0; i<N; i++) A[i]=A[i]*B[i]%M;
    DFT(A,-1);
    LL inv=Pow(N,M-2);
    for (int i=0; i<N; i++) A[i]=A[i]*inv%M;
}

void Work()
{
    val[0]=0;
    for (int i=1; i<=26; i++) val[i]=Rand();

    for (int i=0; i<ls; i++) A[i]=val[ s[i] ];
    for (int i=ls; i<N; i++) A[i]=0;
    for (int i=0; i<lt; i++) B[lt-i-1]=val[ t[i] ];
    for (int i=lt; i<N; i++) B[i]=0;

    NTT();

    LL v=0;
    for (int i=0; i<lt; i++) v=(v+val[ t[i] ]*val[ t[i] ]%M)%M;
    for (int i=0; i<=ls-lt; i++)
        if (A[lt+i-1]!=v) vis[i]=false;
}

int main()
{
    //freopen("4503.in","r",stdin);
    //freopen("std.out","w",stdout);

    scanf("%s",&S);
    ls=strlen(S);
    for (int i=0; i<ls; i++) s[i]=S[i]-'a'+1;
    scanf("%s",&S);
    lt=strlen(S);
    for (int i=0; i<lt; i++)
        if ( 'a'<=S[i] && S[i]<='z' ) t[i]=S[i]-'a'+1;
        else t[i]=0;

    N=1,Lg=0;
    while (N<ls+lt+2) N<<=1,Lg++;
    for (int i=0; i<N; i++)
        for (int j=0; j<Lg; j++)
            if (i&(1<<j)) Rev[i]|=(1<<(Lg-j-1));

    for (int i=0; i<=ls-lt; i++) vis[i]=true;
    int k=3;
    for (int i=0; i<k; i++) Work();

    for (int i=0; i<=ls-lt; i++) if (vis[i]) num++;
    printf("%d\n",num);
    for (int i=0; i<=ls-lt; i++) if (vis[i]) printf("%d\n",i);

    return 0;
}