BZOJ - 4259 残缺的字符串(FFT)

Posted on | In ,

思路

参考了Claris老师的博客。该题最精髓的操作就是将形如A[i]B[i]的式子做一个反转，形成A[n-i]B[i]的卷积形式，这样就能够往FFT上靠。似乎BZOJ一道权限题“快速傅里叶变换”也是使用了这样一个思路，所以这应该算一个比较经典的技巧，在知道这一步后这道题就是FFT的模板题了。

一点小吐槽

  1. 为什么洛谷这题只有5组数据…最后我不得不拿了BZOJ的数据重测了下。
  2. O2优化这么牛批的嘛，不开跑了4.4s开了只跑了1s多啊？如果比5组的总时间，那么速度上快了将近10倍。。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
#include <bits/stdc++.h>
#define ll long long
#define ull unsigned long long
#define INTINF 0x7fffffff
#define LLINF 0x7fffffffffffffff
#define ms(a,val) memset((a),(val),(sizeof(a)))
#define sqr(x) ((x)*(x))
#define all(x) (x).begin(),(x).end()
using namespace std;
ll quickpow(ll a,ll b,ll MOD){ll ans=1;a%=MOD;while(b){if(b&1ll)ans=ans*a%MOD;a=a*a%MOD,b>>=1;}return ans;}
ll gcd(ll a,ll b){return a%b==0?b:gcd(b,a%b);}
//head

const int MAXN=1048576;//1<<n

const double PI=acos(-1);
struct Complex{
    double real,image;
    Complex(){real=0,image=0;}
    Complex(double real){this->real=real,image=0;}
    Complex(double real,double image){this->real=real,this->image=image;}
    Complex operator +(const Complex &c){
        return Complex(c.real+real,c.image+image);
    }
    Complex operator -(const Complex &c){
        return Complex(real-c.real,image-c.image);
    }
    Complex operator *(const Complex &c){
        return Complex(c.real*real-c.image*image,c.real*image+c.image*real);
    }
    Complex operator *(double x){
        return Complex(real*x,image*x);
    }
    Complex operator /(double x){
        return Complex(real/x,image/x);
    }
    Complex rev(){
        return Complex(real,-image);
    }
};
void fft(Complex *a, int n,Complex *com) {
	int k = 0;
	while ((1 << k)<n)k++;
	for (int i = 0;i<n;i++) {
		int t = 0;
		for (int j = 0;j<k;j++)if (i&(1 << j))t |= (1 << (k - j - 1));
		if (i<t)swap(a[i], a[t]);
	}
	for (int l = 2;l <= n;l *= 2) {
		int m = l / 2;
		for (Complex *p = a;p<a + n;p += l) {
			for (int i = 0;i<m;i++) {
				Complex t = com[n / l*i] * p[m + i];
				p[m + i] = p[i] - t;
				p[i]=p[i]+t;
			}
		}
	}
}
Complex omega[MAXN], iomega[MAXN];
void init(int n) {
	for (int i = 0;i<n;i++) {
		omega[i] = Complex(cos(2 * PI / n*i), sin(2 * PI / n*i));
		iomega[i] = omega[i].rev();
	}
}
void dft(Complex *a, int n) {
	fft(a, n,omega);
}
void idft(Complex *a, int n) {
	fft(a, n,iomega);
	for (int i = 0;i<n;i++)a[i] =a[i]/ n;
}
Complex a[MAXN], b[MAXN],c[MAXN];
int A[MAXN],B[MAXN];
int ans[MAXN];
char str[MAXN];

int main() {
    #ifndef ONLINE_JUDGE
    freopen("input.in","r",stdin);
    freopen("output.out","w",stdout);
    #endif // ONLINE_JUDGE
    int m,n;
    while(~scanf("%d%d",&m,&n)){
        //ms(ans,0),ms(A,0),ms(B,0);
        scanf("%s",str);
        for(int i=0;i<m;i++)if(isalpha(str[m-i-1]))A[i]=str[m-i-1]-'a'+1;
        scanf("%s",str);
        for(int i=0;i<n;i++)if(isalpha(str[i]))B[i]=str[i]-'a'+1;
        int len=1;
        while(len<(m+n))len<<=1;
        init(len);
        for(int i=0;i<len;i++)a[i]=Complex(A[i]*A[i]*A[i],0),b[i]=Complex(B[i],0);
        dft(a,len),dft(b,len);
        for(int i=0;i<len;i++)c[i]=a[i]*b[i];
        for(int i=0;i<len;i++)a[i]=Complex(A[i],0),b[i]=Complex(B[i]*B[i]*B[i],0);
        dft(a,len),dft(b,len);
        for(int i=0;i<len;i++)c[i]=c[i]+a[i]*b[i];
        for(int i=0;i<len;i++)a[i]=Complex(A[i]*A[i],0),b[i]=Complex(B[i]*B[i],0);
        dft(a,len),dft(b,len);
        for(int i=0;i<len;i++)c[i]=c[i]-a[i]*b[i]*2.0;
        idft(c,len);
        int cnt=0;
        for(int i=m-1;i<n;i++)if(fabs(c[i].real)<0.5)ans[cnt++]=i-m+2;
        printf("%d\n",cnt);
        for(int i=0;i<cnt;i++){
            if(i)printf(" ");
            printf("%d",ans[i]);
        }
        printf("\n");
    }
    return 0;
}