快速傅立叶之二

 

题意：求 ck=∑i=kn−1ai∗bi−kc_k=\sum\limits_{i=k}^{n-1}a_i*b_{i-k}

FFT

把 aa 数组翻过来。

ck=∑i=kn−1an−1−i∗bi−k=∑i=0n−1−kan−1−k−ibic_k=\sum\limits_{i=k}^{n-1}a_{n-1-i}*b_{i-k}=\sum\limits_{i=0}^{n-1-k}a_{n-1-k-i}b_i。

则 ci=(a∗b)n−i−1c_i=(a*b)_{n-i-1}

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
#include <bits/stdc++.h>
#define debug(x) cout << #x << ":\t" << (x) << endl;
using namespace std;
#define ll long long
const int N = 2e6 + 10;
const int INF = 0x3f3f3f3f;
const ll inf = 0x3f3f3f3f3f3f3f3f;
const ll mod = 1e9 + 7;
const double PI = acos(-1.0);
struct Complex {
	double x, y;
	Complex(double _x = 0.0, double _y = 0.0) {
		x = _x;
		y = _y;
	}
	Complex operator-(const Complex &b) const {
		return Complex(x - b.x, y - b.y);
	}
	Complex operator+(const Complex &b) const {
		return Complex(x + b.x, y + b.y);
	}
	Complex operator*(const Complex &b) const {
		return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
	}
};
void change(Complex y[], int len) {
	int i, j, k;
	for (i = 1, j = len / 2; i < len - 1; i++) {
		if (i < j) swap(y[i], y[j]);
		k = len / 2;
		while (j >= k) {
			j = j - k;
			k = k / 2;
		}
		if (j < k) j += k;
	}
}
/*
 * 做 FFT
 *len 必须是 2^k 形式
 *on == 1 时是 DFT，on == -1 时是 IDFT
 */
void fft(Complex y[], int len, int on) {	//0-len-1
	change(y, len);
	for (int h = 2; h <= len; h <<= 1) {
		Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
		for (int j = 0; j < len; j += h) {
			Complex w(1, 0);
			for (int k = j; k < j + h / 2; k++) {
				Complex u = y[k];
				Complex t = w * y[k + h / 2];
				y[k] = u + t;
				y[k + h / 2] = u - t;
				w = w * wn;
			}
		}
	}
	if (on == -1) {
		for (int i = 0; i < len; i++) {
			y[i].x /= len;
		}
	}
}
int n;
int a[N], b[N], c[N];
Complex A[N], B[N], C[N];
int main() {
	scanf("%d", &n);
	for (int i = 1; i <= n; i++)scanf("%d%d", &a[i], &b[i]);
	reverse(a + 1, a + n + 1);
	int len = 1, mx = 2 * n;
	while (len <= mx)len <<= 1;
	for (int i = 0; i < len; i++)
		A[i] = Complex(a[i] * 1.0, 0), B[i] = Complex(b[i] * 1.0, 0);
	fft(A, len, 1);
	fft(B, len, 1);
	for (int i = 0; i < len; i++)C[i] = A[i] * B[i];
	fft(C, len, -1);
	for (int i = 0; i < n; i++) {
		printf("%lld\n", (ll)(C[n - i + 1].x + 0.5));
	}
	return 0;
}

 

[Zjoi2014]力

 

题意：Fi=∑j<iqiqj(i−j)2−∑j>iqiqj(i−j)2Fi=\sum\limits_{j<i}\frac{q_i q_j}{(i-j)^2}-\sum\limits_{j>i}\frac{q_i q_j}{(i-j)^2}，求 Ei=FiqiE_i=\frac{F_i}{q_i}。

FFT

把 qiq_i 除掉，得到

Ei=∑j<iqj(i−j)2−∑j>iqj(i−j)2=(qi−112+qi−222+qi−332+⋯ )−(qi+112+qi+222+qi+332)\begin{align} E_i &=\sum_{j<i}\frac{q_j}{(i-j)^2}-\sum_{j>i}\frac{q_j}{(i-j)^2}\\ &= (\frac{q_{i-1}}{1^2}+\frac{q_{i-2}}{2^2}+\frac{q_{i-3}}{3^2}+\cdots)-(\frac{q_{i+1}}{1^2}+\frac{q_{i+2}}{2^2}+\frac{q_{i+3}}{3^2})\\ \end{align}

构造数组 a[i]=q[i]a[i]=q[i]，b[i]=1i2b[i]=\frac{1}{i^2}

对于前半部分，直接就是卷积形式。(a∗b)i−1(a*b)_{i-1}。

对于后半部分，两个下标都递增，和上题类似，把 aa 数组翻过来，(reva∗b)n−i−2(reva*b)_{n-i-2}。

下标由于应该从0开始，所以要注意一下。

darkbzoj上没有spj。0.o

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
#include <bits/stdc++.h>
#define debug(x) cout << #x << ":\t" << (x) << endl;
using namespace std;
#define ll long long
const int N = 2e6 + 10;
const int INF = 0x3f3f3f3f;
const ll inf = 0x3f3f3f3f3f3f3f3f;
const ll mod = 1e9 + 7;
const double PI = acos(-1.0);
struct Complex {
	double x, y;
	Complex(double _x = 0.0, double _y = 0.0) {
		x = _x;
		y = _y;
	}
	Complex operator-(const Complex &b) const {
		return Complex(x - b.x, y - b.y);
	}
	Complex operator+(const Complex &b) const {
		return Complex(x + b.x, y + b.y);
	}
	Complex operator*(const Complex &b) const {
		return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
	}
};
void change(Complex y[], int len) {
	int i, j, k;
	for (i = 1, j = len / 2; i < len - 1; i++) {
		if (i < j) swap(y[i], y[j]);
		k = len / 2;
		while (j >= k) {
			j = j - k;
			k = k / 2;
		}
		if (j < k) j += k;
	}
}
/*
 * 做 FFT
 *len 必须是 2^k 形式
 *on == 1 时是 DFT，on == -1 时是 IDFT
 */
void fft(Complex y[], int len, int on) {	//0-len-1
	change(y, len);
	for (int h = 2; h <= len; h <<= 1) {
		Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
		for (int j = 0; j < len; j += h) {
			Complex w(1, 0);
			for (int k = j; k < j + h / 2; k++) {
				Complex u = y[k];
				Complex t = w * y[k + h / 2];
				y[k] = u + t;
				y[k + h / 2] = u - t;
				w = w * wn;
			}
		}
	}
	if (on == -1) {
		for (int i = 0; i < len; i++) {
			y[i].x /= len;
		}
	}
}
int n;
Complex a[N], b[N], c[N];
double ans[N];
int main() {
	scanf("%d", &n);
	for (int i = 0; i < n; i++)scanf("%lf", &a[i].x);
	int mx = 2 * n;
	int len = 1;
	while (len <= mx)len <<= 1;
	for (int i = 0; i < n; i++)b[i].x = 1.0 / (i + 1) / (i + 1);
	fft(a, len, 1);
	fft(b, len, 1);
	for (int i = 0; i < len; i++)c[i] = a[i] * b[i];
	fft(c, len, -1);
	for (int i = 1; i < n; i++)ans[i] = c[i - 1].x;
	fft(a, len, -1);
	reverse(a, a + n);
	fft(a, len, 1);
	for (int i = 0; i < len; i++)c[i] = a[i] * b[i];
	fft(c, len, -1);
	for (int i = 0; i < n - 1; i++)ans[i] -= c[n - i - 2].x;
	for (int i = 0; i < n; i++)printf("%.3lf\n", ans[i]);
	return 0;
}

 

两个串

 

题意：给定两个字符串，问 T 在 S 中所有出现位置，T 中有通配符 ‘?’。

FFT

设 S 长度 n，T 长度 m。

如果没有通配符，可以用

Ak=∑i=0m−1(S[k+i]−T[i])2A_k=\sum_{i=0}^{m-1}(S[k+i]-T[i])^2\\

判断 k 开始的 S 串是否与 T 匹配。

加入通配符，则无论 S[k+i]S[k+i] 是什么，S[k+i]−T[i]S[k+i]-T[i] 都要为0，显然不行，但是可以使得每一项都乘以 T[i]T[i]，且令 T[i]T[i] 为通配符时 T[i]=0T[i]=0，否则为非零。

Ak=∑i=0m−1(S[k+i]−T[i])2⋅T[i]A_k=\sum_{i=0}^{m-1}(S[k+i]-T[i])^2\cdot T[i]\\

平方展开得到

Ak=∑i=0m−1S2[k+i]⋅T[m−1−i]−2∑i=0m−1S[k+i]⋅T2[m−1−i]+∑i=0m−1T3[m−1−i]A_k=\sum_{i=0}^{m-1}S^2[k+i]\cdot T[m-1-i]-2\sum_{i=0}^{m-1}S[k+i]\cdot T^2[m-1-i]+\sum_{i=0}^{m-1}T^3[m-1-i]\\

最后一项是常数项，前两项是卷积。

这里 S2,T2S^2,T^2 是对应项的平方，不是多项式乘积。

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
#include <bits/stdc++.h>
#define debug(x) cout << #x << ":\t" << (x) << endl;
using namespace std;
#define ll long long
const int N = 2e6 + 10;
const int INF = 0x3f3f3f3f;
const ll inf = 0x3f3f3f3f3f3f3f3f;
const ll mod = 1e9 + 7;
const double PI = acos(-1.0);
struct Complex {
	double x, y;
	Complex(double _x = 0.0, double _y = 0.0) {
		x = _x;
		y = _y;
	}
	Complex operator-(const Complex &b) const {
		return Complex(x - b.x, y - b.y);
	}
	Complex operator+(const Complex &b) const {
		return Complex(x + b.x, y + b.y);
	}
	Complex operator*(const Complex &b) const {
		return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
	}
};
void change(Complex y[], int len) {
	int i, j, k;
	for (i = 1, j = len / 2; i < len - 1; i++) {
		if (i < j) swap(y[i], y[j]);
		k = len / 2;
		while (j >= k) {
			j = j - k;
			k = k / 2;
		}
		if (j < k) j += k;
	}
}
/*
 * 做 FFT
 *len 必须是 2^k 形式
 *on == 1 时是 DFT，on == -1 时是 IDFT
 */
void fft(Complex y[], int len, int on) {	//0-len-1
	change(y, len);
	for (int h = 2; h <= len; h <<= 1) {
		Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
		for (int j = 0; j < len; j += h) {
			Complex w(1, 0);
			for (int k = j; k < j + h / 2; k++) {
				Complex u = y[k];
				Complex t = w * y[k + h / 2];
				y[k] = u + t;
				y[k + h / 2] = u - t;
				w = w * wn;
			}
		}
	}
	if (on == -1) {
		for (int i = 0; i < len; i++) {
			y[i].x /= len;
		}
	}
}
int n, m;
char s[N], t[N];
int S[N], T[N];
Complex a1[N], a2[N], b1[N], b2[N], c1[N], c2[N];
ll C;
vector<int>vc;
int main() {
	scanf("%s%s", s, t);
	n = (int)strlen(s), m = (int)strlen(t);
	for (int i = 0; i < n; i++)S[i] = s[i] - 'a' + 1;
	for (int i = 0; i < m; i++) {
		if (t[i] == '?')T[i] = 0;
		else T[i] = t[i] - 'a' + 1;
		C += T[i] * T[i] * T[i];
	}
	reverse(T, T + m);
	int len = 1, mx = n + m;
	while (len <= mx)len <<= 1;
	for (int i = 0; i < n; i++) {
		a1[i] = Complex(S[i], 0); a2[i] = Complex(S[i] * S[i], 0);
		b1[i] = Complex(T[i], 0); b2[i] = Complex(T[i] * T[i], 0);
	}
	fft(a1, len, 1); fft(a2, len, 1);
	fft(b1, len, 1); fft(b2, len, 1);
	for (int i = 0; i < len; i++)c1[i] = a2[i] * b1[i];
	for (int i = 0; i < len; i++)c2[i] = a1[i] * b2[i];
	fft(c1, len, -1); fft(c2, len, -1);
	for (int i = 0; i <= n - m; i++) {
		ll tmp1 = (ll)(c1[m + i - 1].x + 0.5);
		ll tmp2 = (ll)(c2[m + i - 1].x + 0.5);
		if (tmp1 - 2 * tmp2 + C == 0)vc.push_back(i);
	}
	printf("%d\n", (int)vc.size());
	for (int u : vc)printf("%d\n", u);
	return 0;
}

 

残缺的字符串

 

题意：同上题，两串都有通配符。

FFT

同上题，三项卷积。

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
#include <bits/stdc++.h>
#define debug(x) cout << #x << ":\t" << (x) << endl;
using namespace std;
#define ll long long
const int N = 8e5 + 10;
const int INF = 0x3f3f3f3f;
const ll inf = 0x3f3f3f3f3f3f3f3f;
const ll mod = 1e9 + 7;
const double PI = acos(-1.0);
struct Complex {
	double x, y;
	Complex(double _x = 0.0, double _y = 0.0) {
		x = _x;
		y = _y;
	}
	Complex operator-(const Complex &b) const {
		return Complex(x - b.x, y - b.y);
	}
	Complex operator+(const Complex &b) const {
		return Complex(x + b.x, y + b.y);
	}
	Complex operator*(const Complex &b) const {
		return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
	}
};
void change(Complex y[], int len) {
	int i, j, k;
	for (i = 1, j = len / 2; i < len - 1; i++) {
		if (i < j) swap(y[i], y[j]);
		k = len / 2;
		while (j >= k) {
			j = j - k;
			k = k / 2;
		}
		if (j < k) j += k;
	}
}
/*
 * 做 FFT
 *len 必须是 2^k 形式
 *on == 1 时是 DFT，on == -1 时是 IDFT
 */
void fft(Complex y[], int len, int on) {	//0-len-1
	change(y, len);
	for (int h = 2; h <= len; h <<= 1) {
		Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
		for (int j = 0; j < len; j += h) {
			Complex w(1, 0);
			for (int k = j; k < j + h / 2; k++) {
				Complex u = y[k];
				Complex t = w * y[k + h / 2];
				y[k] = u + t;
				y[k + h / 2] = u - t;
				w = w * wn;
			}
		}
	}
	if (on == -1) {
		for (int i = 0; i < len; i++) {
			y[i].x /= len;
		}
	}
}
int n, m;
char s[N], t[N];
int S[N], T[N];
Complex a1[N], a2[N], a3[N], b1[N], b2[N], b3[N], c1[N], c2[N], c3[N];
vector<int>vc;
int main() {
	scanf("%d%d", &m, &n);
	scanf("%s%s", t, s);
	for (int i = 0; i < n; i++) {
		if (s[i] == '*')S[i] = 0;
		else S[i] = s[i] - 'a' + 1;
	}
	for (int i = 0; i < m; i++) {
		if (t[i] == '*')T[i] = 0;
		else T[i] = t[i] - 'a' + 1;
	}
	reverse(T, T + m);
	int len = 1, mx = n + m;
	while (len <= mx)len <<= 1;
	for (int i = 0; i < n; i++) {
		a1[i] = Complex(S[i], 0); a2[i] = Complex(S[i] * S[i], 0); a3[i] = Complex(S[i] * S[i] * S[i], 0);
		b1[i] = Complex(T[i], 0); b2[i] = Complex(T[i] * T[i], 0); b3[i] = Complex(T[i] * T[i] * T[i], 0);
	}
	fft(a1, len, 1); fft(a2, len, 1); fft(a3, len, 1);
	fft(b1, len, 1); fft(b2, len, 1); fft(b3, len, 1);
	for (int i = 0; i < len; i++)c1[i] = a3[i] * b1[i];
	for (int i = 0; i < len; i++)c2[i] = a1[i] * b3[i];
	for (int i = 0; i < len; i++)c3[i] = a2[i] * b2[i];
	fft(c1, len, -1); fft(c2, len, -1); fft(c3, len, -1);
	for (int i = 0; i <= n - m; i++) {
		ll tmp1 = (ll)(c1[m + i - 1].x + 0.5);
		ll tmp2 = (ll)(c2[m + i - 1].x + 0.5);
		ll tmp3 = (ll)(c3[m + i - 1].x + 0.5);
		if (tmp1 + tmp2 - 2 * tmp3 == 0)vc.push_back(i + 1);
	}
	printf("%d\n", (int)vc.size());
	for (int i = 0; i < (int)vc.size(); i++)printf("%d%c", vc[i], " \n"[i == (int)vc.size() - 1]);
	return 0;
}