Luogu4173 残缺的字符串 FFT

传送门

---

考虑如何使用FFT计算两个子串是否匹配。如果字符集比较小可以把每个字符都拿出来暴力做一遍，但是字符集比较大的时候复杂度就会有问题。这个时候可以考虑匹配函数。

先考虑没有通配符的情况。将\(A\)串翻转，然后设匹配函数\(chk(i,j) = (A_i - B_j)^2\)。不难知道\(A_i = B_j \Leftrightarrow chk(i,j) = 0\)。

又设\(C(x) = \sum\limits_{i=1}^m chk(m + 1 - i , x + i - 1)\)，那么\(B\)的以\(x\)为左端点的长度为\(m\)的子串能够跟\(A\)串匹配的充要条件就是\(C(x) = 0\)。

而\(C(x) = \sum\limits_{i=1}^m (A_{m+1-i} - B_{x+i-1})^2 = \sum\limits_{i=1}^m (A_{m+1-i}^2 + B_{x+i-1}^2) - 2\sum\limits_{i=1}^m A_{m+1-i}B_{x+i-1}\)。可以发现最后的式子是一个卷积的形式，而前面两项都可以前缀和。于是直接FFT算一下最后的卷积就可以\(O(NlogN)\)判断。

现在考虑通配符的情况，通配符可以与任意字符匹配，所以不妨设\(chk(i,j) = (A_i - B_j)^2 A_i B_j\)，其中如果某一个字符为通配符就令它的值为\(0\)，这样一个通配符匹配的贡献就永远都是\(0\)了。而

\[C(x) = \sum\limits_{i=1}^m chk(m + 1 - i , x+i-1) = \sum\limits_{i=1}^m (A_{m+1-i} - B_{x+i-1})^2 A_{m+1-i} B_{x+i-1} = \sum\limits_{i=1}^m A_{m+1-i}^3 B_{x+i-1} + \sum\limits_{i=1}^m A_{m+1-i}B_{x+i-1}^3 - 2\sum\limits_{i=1}^m A_{m+1-i}^2B_{x+i-1}^2\]

三个都用FFT算一遍就行了

#include<iostream>
#include<cstdio>
#include<cctype>
#include<algorithm>
#include<cstring>
#include<cmath>
//This code is written by Itst
using namespace std;

#define ld double
const int MAXN = (1 << 20) + 7;
struct comp{
    ld x , y;
    comp(ld _x = 0 , ld _y = 0) : x(_x) , y(_y){}
    comp operator +(comp a){return comp(x + a.x , y + a.y);}
    comp operator -(comp a){return comp(x - a.x , y - a.y);}
    comp operator *(comp a){return comp(x * a.x - y * a.y , x * a.y + y * a.x);}
}A[MAXN] , B[MAXN] , C[MAXN] , D[MAXN];
ld ans[MAXN];
const ld pi = acos(-1);
int M , N , need , dir[MAXN];
char s[MAXN];

inline void init(int len){
    need = 1;
    while(need < len)
        need <<= 1;
    for(int i = 1 ; i < need ; ++i)
        dir[i] = (dir[i >> 1] >> 1) | (i & 1 ? need >> 1 : 0);
}

void FFT(comp *arr , int type){
    for(int i = 1 ; i < need ; ++i)
        if(i < dir[i])
            swap(arr[i] , arr[dir[i]]);
    for(int i = 1 ; i < need ; i <<= 1){
        comp wn(cos(pi / i) , type * sin(pi / i));
        for(int j = 0 ; j < need ; j += i << 1){
            comp w(1 , 0);
            for(int k = 0 ; k < i ; ++k , w = w * wn){
                comp x = arr[j + k] , y = arr[i + j + k] * w;
                arr[j + k] = x + y; arr[i + j + k] = x - y;
            }
        }
    }
}

void work(){
    for(int i = 0 ; i < need ; ++i)
        C[i] = comp(A[i].x * A[i].x * A[i].x , 0);
    for(int i = 0 ; i < need ; ++i) D[i] = B[i];
    FFT(C , 1); FFT(D , 1);
    for(int i = 0 ; i < need ; ++i)
        C[i] = C[i] * D[i];
    FFT(C , -1);
    for(int i = M + 1 ; i <= N + 1 ; ++i)
        ans[i] += C[i].x / need;

    for(int i = 0 ; i < need ; ++i) C[i] = A[i];
    for(int i = 0 ; i < need ; ++i)
        D[i] = comp(B[i].x * B[i].x * B[i].x , 0);
    FFT(C , 1); FFT(D , 1);
    for(int i = 0 ; i < need ; ++i)
        C[i] = C[i] * D[i];
    FFT(C , -1);
    for(int i = M + 1 ; i <= N + 1 ; ++i)
        ans[i] += C[i].x / need;
    
    for(int i = 0 ; i < need ; ++i)
        C[i] = comp(A[i].x * A[i].x , 0);
    for(int i = 0 ; i < need ; ++i)
        D[i] = comp(B[i].x * B[i].x , 0);
    FFT(C , 1); FFT(D , 1);
    for(int i = 0 ; i < need ; ++i)
        C[i] = C[i] * D[i];
    FFT(C , -1);
    for(int i = M + 1 ; i <= N + 1 ; ++i)
        ans[i] -= 2 * C[i].x / need;
}

int main(){
#ifndef ONLINE_JUDGE
    freopen("in","r",stdin);
    //freopen("out","w",stdout);
#endif
    scanf("%d %d" , &M , &N);
    scanf("%s" , s + 1);
    reverse(s + 1 , s + M + 1);
    for(int i = 1 ; i <= M ; ++i)
        A[i].x = s[i] == '*' ? 0 : s[i] - 'a' + 1;
    scanf("%s" , s + 1);
    for(int i = 1 ; i <= N ; ++i)
        B[i].x = s[i] == '*' ? 0 : s[i] - 'a' + 1;
    init(M + N + 1);
    work();
    int cnt = 0;
    for(int i = M + 1 ; i <= N + 1 ; ++i)
        cnt += ans[i] <= 0.5;
    cout << cnt << endl;
    for(int i = M + 1 ; i <= N + 1 ; ++i)
        if(ans[i] <= 0.5)
            cout << i - M << ' ';
    return 0;
}

转载于:https://www.cnblogs.com/Itst/p/10438862.html