BZOJ4259 残缺的字符串 FFT

首先是csy博客的题解,但csy的并不十分优秀。

事实上,我们可以求出f[i] = sigma(A[j] * inv[B[i - j]]),当f[i] <= n时,匹配,否则不匹配。这样就只需要做一次FFT或者NTT即可。

 1 #include<cstdio>
 2 #include<cmath>
 3 #include<algorithm>
 4 #define N 1048576
 5 using namespace std;
 6 char sa[N],sb[N];int n,m,i,j,k,a[N],b[N],ans,q[N];
 7 struct comp{
 8   double r,i;comp(double _r=0,double _i=0){r=_r;i=_i;}
 9   comp operator+(const comp&x){return comp(r+x.r,i+x.i);}
10   comp operator-(const comp&x){return comp(r-x.r,i-x.i);}
11   comp operator*(const comp&x){return comp(r*x.r-i*x.i,r*x.i+i*x.r);}
12 }A[N],B[N],C[N];
13 const double pi=acos(-1.0);
14 void FFT(comp*a,int n,int t){
15   for(int i=1,j=0;i<n-1;i++){
16     for(int s=n;j^=s>>=1,~j&s;);
17     if(i<j)swap(a[i],a[j]);
18   }
19   for(int d=0;(1<<d)<n;d++){
20     int m=1<<d,m2=m<<1;
21     double o=pi/m*t;comp _w(cos(o),sin(o));
22     for(int i=0;i<n;i+=m2){
23       comp w(1,0);
24       for(int j=0;j<m;j++){
25         comp &A=a[i+j+m],&B=a[i+j],t=w*A;
26         A=B-t;B=B+t;w=w*_w;
27       }
28     }
29   }
30   if(t==-1)for(int i=0;i<n;i++)a[i].r/=n;
31 }
32 int main(){
33   scanf("%d%d%s%s",&m,&n,sa,sb);
34   for(i=0,j=m-1;i<j;i++,j--)swap(sa[i],sa[j]);
35   for(i=0;i<m;i++)if(sa[i]!='*')a[i]=sa[i]-'a'+1;
36   for(i=0;i<n;i++)if(sb[i]!='*')b[i]=sb[i]-'a'+1;
37   for(k=1;k<n+m;k<<=1);
38   for(i=0;i<k;i++)A[i]=comp(a[i]*a[i]*a[i],0),B[i]=comp(b[i],0);
39   for(FFT(A,k,1),FFT(B,k,1),i=0;i<k;i++)C[i]=C[i]+A[i]*B[i];
40   for(i=0;i<k;i++)A[i]=comp(a[i],0),B[i]=comp(b[i]*b[i]*b[i],0);
41   for(FFT(A,k,1),FFT(B,k,1),i=0;i<k;i++)C[i]=C[i]+A[i]*B[i];
42   for(i=0;i<k;i++)A[i]=comp(a[i]*a[i],0),B[i]=comp(b[i]*b[i],0);
43   for(FFT(A,k,1),FFT(B,k,1),i=0;i<k;i++)C[i]=C[i]-A[i]*B[i]*comp(2,0);
44   FFT(C,k,-1);
45   for(i=m-1;i<n;i++)if(C[i].r<0.5)q[++ans]=i-m+2;
46   for(printf("%d\n",ans),i=1;i<ans;i++)printf("%d ",q[i]);
47   if(ans)printf("%d",q[ans]);
48   return 0;
49 }

 NTT版本代码

  1 #include<bits/stdc++.h>
  2 using namespace std;
  3 
  4 const int MOD = 998244353;
  5 const int N = 1048576;
  6 const int M = 998244353;
  7 
  8 int WM[N + 2], IWM[N + 2];
  9 vector<int> IW[N + 2], W[N + 2];
 10 
 11 int Pw(int a, int b, int p)
 12 {
 13     int v = 1;
 14     for(; b; b >>= 1, a = 1ll * a * a % p)
 15         if (b & 1)
 16             v = 1ll * v * a % p;
 17     return v;
 18 }
 19 
 20 void ycl()
 21 {
 22     for (int m = 2; m <= N; m <<= 1)
 23     {
 24         WM[m] = Pw(3, (M - 1) / m, M), IWM[m] = Pw(3, (M - 1) / m * (m - 1), M);
 25         int o = 1;
 26         W[m].push_back(o);
 27         for (int i = 1; i < m; ++i)
 28             o = 1ll * o * WM[m] % M, W[m].push_back(o);
 29         o = 1;
 30         IW[m].push_back(o);
 31         for (int i = 1; i < m; ++i)
 32             o = 1ll * o * IWM[m] % M, IW[m].push_back(o);
 33     }
 34 }
 35 
 36 void NTT(int *a, int n, int f = 1)
 37 {
 38     int i, j, k, m, w, u, v;
 39     for (i = n >> 1, j = 1; j < n - 1; ++j)
 40     {
 41         if (i > j)
 42             swap(a[i], a[j]);
 43         for (k = n >> 1; k <= i; k >>= 1)
 44             i ^= k;
 45         i ^= k;
 46     }
 47     for (m = 2; m <= n; m <<= 1)
 48         for (i = 0; i < n; i += m)
 49             for (j = i; j < i + (m >> 1); ++j)
 50                 if (a[j] || a[j + (m >> 1)])
 51                 {
 52                     u = a[j];
 53                     v = 1ll * (f == 1 ? W[m][j - i] : IW[m][j - i]) * a[j + (m >> 1)] % M;
 54                     if ((a[j] = u + v) >= M)
 55                         a[j] -= M;
 56                     if ((a[j + (m >> 1)] = u - v) < 0)
 57                         a[j + (m >> 1)] += M;
 58                 }
 59     if (f == -1)
 60         for (w = Pw(n, M - 2, M), i = 0; i < n; ++i)
 61             a[i] = 1ll * a[i] * w % M;
 62 }
 63 
 64 char sa[N], sb[N], st[N];
 65 int a[5][N], b[5][N], c[N], q[N];
 66 
 67 int main()
 68 {
 69     ycl();
 70     int m, n;
 71     scanf("%d%d", &m, &n);
 72     scanf("%s%s", sa, sb);
 73     for (int i = 0, j = m - 1; i < j; ++i, --j)
 74         swap(sa[i],sa[j]);
 75     for (int i = 0; i < m; ++i)
 76         a[1][i] = sa[i] - 'a' + 1;
 77     for (int i = 0; i < n; ++i)
 78         b[1][i] = sb[i] - 'a' + 1;
 79     for (int i = 0; i < m; ++i)
 80         if(sa[i] == '*')
 81             a[1][i] = 0;
 82     for (int i = 0; i < n; ++i)
 83         if(sb[i] == '*')
 84             b[1][i] = 0;
 85     int K;
 86     for (K = 1; K < n + m; K <<= 1);
 87     for (int i = 2; i <= 3; ++i)
 88     {
 89         for (int j = 0; j < m; ++j)
 90             a[i][j] = a[i - 1][j] * a[1][j];
 91         for (int j = 0; j < n; ++j)
 92             b[i][j] = b[i - 1][j] * b[1][j];
 93     }
 94     for (int i = 1; i < 4; ++ i)
 95     {
 96         NTT(a[i], K);
 97         NTT(b[i], K);
 98     }
 99     for (int i = 0; i < K; ++ i)
100     {
101         c[i] = (c[i] + 1LL * a[3][i] * b[1][i] % MOD) % MOD;
102         c[i] = (c[i] + MOD - 2LL * a[2][i] * b[2][i] % MOD) % MOD;
103         c[i] = (c[i] + 1LL * a[1][i] * b[3][i] % MOD) % MOD;
104     }
105     NTT(c, K, -1);
106     
107     int ans(0);
108     for (int i = m - 1; i < n; ++i)
109         if (!c[i])
110             q[++ans] = i - m + 2;
111     printf("%d\n", ans);
112     for (int i = 1; i < ans; ++i)
113         printf("%d ", q[i]);
114     if (ans)
115         printf("%d\n", q[ans]);
116 
117     return 0;
118 }

 

posted @ 2018-04-28 16:34  Aseer  阅读(194)  评论(0编辑  收藏  举报