# # 快速傅里叶变换

## # 多项式的系数表示和点值表示

$f(x)=\sum_{i=0}^na_ix^i$

## # 点值表示与多项式乘法的关系

$\{(x_0,f(x_0)g(x_0)),(x_1,f(x_1)g(x_1)),\cdots,(x_n,f(x_n)g(x_n))\}$

## # FFT的实现

### # FFT

$\omega_n^n=1$

$\omega_n^i=\omega_{2n}^{2i}$

$\omega_{2n}^{n+i}=-\omega_{2n}^i$

$f(x)=a_0+a_1x^1+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7$

\begin{aligned} f(x) &=(a_0+a_2x^2+a_4x^4+a_6x^6)+x(a_1+a_3x^2+a_5x^4+a_7x^6) \\ &=G(x^2)+xH(x^2) \end{aligned}

$\text{DFT}(f(x))=\text{DFT}(G(x^2))+x\text{DFT}(H(x^2))$

\begin{aligned} \text{DFT}(f(\omega_n^k))&=\text{DFT}(G(\omega_n^{2k}))+\omega_n^k\text{DFT}(H(\omega_n^{2k})) \\ &=\text{DFT}(G(\omega_{n/2}^k))+\omega_n^k\text{DFT}(H(\omega_{n/2}^k)) \end{aligned}

\begin{aligned} \text{DFT}(f(\omega_n^{k+n/2}))&=\text{DFT}(G(\omega_n^{2k+n}))+\omega_n^k\text{DFT}(H(\omega_n^{2k+n})) \\ &=\text{DFT}(G(\omega_{n/2}^{k+n/2}))+\omega_n^{k+n/2}\text{DFT}(H(\omega_{n/2}^{k+n/2})) \\ &=\text{DFT}(G(\omega_{n/2}^k))-\omega_n^k\text{DFT}(H(\omega_{n/2}^k)) \end{aligned}

$T(n)=2T(n/2)$

### # 模板题：洛谷 P3803 - 多项式乘法（FFT） (opens new window)

#include <cmath>
#include <complex>
#include <iostream>
#define MAXN (1 << 22)

using namespace std;
typedef complex<double> cd;
const cd I{0, 1};
cd tmp[MAXN], a[MAXN], b[MAXN];
void fft(cd *f, int n, int rev) {
if (n == 1)
return;
for (int i = 0; i < n; ++i)
tmp[i] = f[i];
for (int i = 0; i < n; ++i) {
if (i & 1)
f[n / 2 + i / 2] = tmp[i];
else
f[i / 2] = tmp[i];
}
cd *g = f, *h = f + n / 2;
fft(g, n / 2, rev), fft(h, n / 2, rev);
cd omega = exp(I * (2 * M_PI / n * rev)), now = 1;
for (int k = 0; k < n / 2; ++k) {
tmp[k] = g[k] + now * h[k];
tmp[k + n / 2] = g[k] - now * h[k];
now *= omega;
}
for (int i = 0; i < n; ++i)
f[i] = tmp[i];
}
int main() {
int n, m;
cin >> n >> m;
int k = 1 << (32 - __builtin_clz(n + m + 1));
for (int i = 0; i <= n; ++i)
cin >> a[i];
for (int j = 0; j <= m; ++j)
cin >> b[j];
fft(a, k, 1);
fft(b, k, 1);
for (int i = 0; i < k; ++i)
a[i] *= b[i];
fft(a, k, -1);
for (int i = 0; i < k; ++i)
a[i] /= k;
for (int i = 0; i < n + m + 1; ++i)
cout << (int)round(a[i].real()) << " ";
}

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

### # 蝴蝶变换

$\{a_0,a_2,a_4,a_6\},\{a_1,a_3,a_5,a_7\}$

$\{a_0,a_4\},\{a_2,a_6\},\{a_1,a_5\},\{a_3,a_7\}$

$\{a_0\},\{a_4\},\{a_2\},\{a_6\},\{a_1\},\{a_5\},\{a_3\},\{a_7\}$

### # FFT的非递归实现

#include <cmath>
#include <complex>
#include <iostream>
#define MAXN (1 << 22)

using namespace std;
typedef complex<double> cd;
const cd I{0, 1};
cd a[MAXN], b[MAXN];
void change(cd *f, int n) {
int i, j, k;
for (int i = 1, j = n / 2; i < n - 1; i++) {
if (i < j)
swap(f[i], f[j]);
k = n / 2;
while (j >= k) {
j = j - k;
k = k / 2;
}
if (j < k)
j += k;
}
}
void fft(cd *f, int n, int rev) {
change(f, n);
for (int len = 2; len <= n; len <<= 1) {
cd omega = exp(I * (2 * M_PI / len * rev));
for (int j = 0; j < n; j += len) {
cd now = 1;
for (int k = j; k < j + len / 2; ++k) {
cd g = f[k], h = now * f[k + len / 2];
f[k] = g + h, f[k + len / 2] = g - h;
now *= omega;
}
}
}
if (rev == -1)
for (int i = 0; i < n; ++i)
f[i] /= n;
}
int main() {
int n, m;
cin >> n >> m;
int k = 1 << (32 - __builtin_clz(n + m + 1));
for (int i = 0; i <= n; ++i)
cin >> a[i];
for (int j = 0; j <= m; ++j)
cin >> b[j];
fft(a, k, 1);
fft(b, k, 1);
for (int i = 0; i < k; ++i)
a[i] *= b[i];
fft(a, k, -1);
for (int i = 0; i < n + m + 1; ++i)
cout << (int)round(a[i].real()) << " ";
}

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

import kotlin.math.*

}

fun nearest(n: Int): Int {
return 1.shl(ceil(log2(n.toDouble())).toInt() + 1)
}

data class Complex(val real: Double, val imag: Double) {
operator fun plus(other: Complex): Complex {
return Complex(real + other.real, imag + other.imag)
}

operator fun minus(other: Complex): Complex {
return Complex(real - other.real, imag - other.imag)
}

operator fun times(other: Double): Complex {
return Complex(real * other, imag * other)
}

operator fun times(other: Complex): Complex {
return Complex(real * other.real - imag * other.imag, real * other.imag + imag * other.real)
}
}

val rev = IntArray(nearest(1000002))

fun change(f: Array<Complex>, n: Int) {
if (rev[1] == 0)
for (i in 0 until n) {
rev[i] = rev[i shr 1] shr 1
if (i % 2 == 1)
rev[i] = rev[i] or (n shr 1)
}
for (i in 0 until n)
if (i < rev[i]) {
val tmp = f[i]
f[i] = f[rev[i]]
f[rev[i]] = tmp
}
}

fun fft(f: Array<Complex>, n: Int, rev: Int) {
change(f, n)
var len = 2
while (len <= n) {
val omega = Complex(cos(2.0 * PI / len * rev), sin(2.0 * PI / len * rev))
var l = 0
while (l < n) {
var now = Complex(1.0, 0.0)
for (i in l until l + len / 2) {
val g = f[i]
val h = now * f[i + len / 2]
f[i] = g + h
f[i + len / 2] = g - h
now *= omega
}
l += len
}
len *= 2
}
if (rev == -1)
for (i in 0 until n)
f[i] = f[i] * (1.0 / n)
}

fun main() {
val k = nearest(n + m + 1)
val a = Array(k) { Complex(0.0, 0.0) }
a[i] = Complex(v.toDouble(), 0.0)
val b = Array(k) { Complex(0.0, 0.0) }
b[i] = Complex(v.toDouble(), 0.0)
fft(a, k, 1)
fft(b, k, 1)
for (i in 0 until k)
a[i] = a[i] * b[i]
fft(a, k, -1)
val ans = a.slice(0 until n + m + 1).map { round(it.real).toInt() }
println(ans.joinToString(" "))
}

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

## # 学习资源

### #Matters Computational (opens new window)

• 第二十一章 快速傅里叶变换

## # 练习题

### #SPOJ - ADAMATCH (opens new window)

#include <cmath>
#include <complex>
#include <cstring>
#include <iostream>
#include <vector>
#define MAXN (1 << 22)

using namespace std;
typedef complex<double> cd;
const cd I{0, 1};
cd a[MAXN], b[MAXN];
void change(cd *f, int n) {
for (int i = 1, j = n / 2; i < n - 1; i++) {
if (i < j)
swap(f[i], f[j]);
int k = n / 2;
while (j >= k) {
j = j - k;
k = k / 2;
}
if (j < k)
j += k;
}
}
void fft(cd *f, int n, int rev) {
change(f, n);
for (int len = 2; len <= n; len <<= 1) {
cd omega = exp(I * (2 * M_PI / len * rev));
for (int j = 0; j < n; j += len) {
cd now = 1;
for (int k = j; k < j + len / 2; ++k) {
cd g = f[k], h = now * f[k + len / 2];
f[k] = g + h, f[k + len / 2] = g - h;
now *= omega;
}
}
}
if (rev == -1)
for (int i = 0; i < n; ++i)
f[i] /= n;
}
int main() {
string s, r;
cin >> s >> r;
int n = s.size(), m = r.size();
int k = 1 << (32 - __builtin_clz(n + m + 1));
vector<int> cnt(k);
for (char c : "ACGT") {
memset(a, 0, sizeof(a));
memset(b, 0, sizeof(b));
for (int i = 0; i < n; ++i)
a[i] = s[i] == c;
for (int i = 0; i < m; ++i)
b[i] = r[m - i - 1] == c;
fft(a, k, 1);
fft(b, k, 1);
for (int i = 0; i < k; ++i)
a[i] *= b[i];
fft(a, k, -1);
for (int i = 0; i < k; ++i)
cnt[i] += (int)round(a[i].real());
}
int ans = m;
for (int i = m - 1; i < n; ++i)
ans = min(ans, m - cnt[i]);
cout << ans;
}

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

### #SPOJ - TSUM (opens new window)

#include <cmath>
#include <complex>
#include <iostream>
#include <vector>

#define MAXN 131072
#define OFFSET 20000

using namespace std;
typedef complex<double> cd;
const cd I{0, 1};

void change(vector<cd> &f, int n) {
for (int i = 1, j = n / 2; i < n - 1; i++) {
if (i < j)
swap(f[i], f[j]);
int k = n / 2;
while (j >= k) {
j = j - k;
k = k / 2;
}
if (j < k)
j += k;
}
}

void fft(vector<cd> &f, int n, int rev) {
change(f, n);
for (int len = 2; len <= n; len <<= 1) {
cd omega = exp(I * (2 * M_PI / len * rev));
for (int j = 0; j < n; j += len) {
cd now = 1;
for (int k = j; k < j + len / 2; ++k) {
cd g = f[k], h = now * f[k + len / 2];
f[k] = g + h, f[k + len / 2] = g - h;
now *= omega;
}
}
}
if (rev == -1)
for (int i = 0; i < n; ++i)
f[i] /= n;
}

int main() {
int n;
cin >> n;
vector<cd> a(MAXN), a2(MAXN);
vector<int> a3(MAXN);
for (int i = 0; i < n; ++i) {
int m;
cin >> m;
a[m + OFFSET] = cd{1, 0};
a2[(m + OFFSET) << 1] = cd{1, 0};
a3[(m + OFFSET) * 3] = 1;
}
vector<cd> tot(a), b(a);
fft(tot, MAXN, 1);
fft(b, MAXN, 1);
fft(a2, MAXN, 1);
for (int i = 0; i < MAXN; ++i)
tot[i] *= b[i] * b[i], a2[i] *= b[i];
fft(tot, MAXN, -1);
fft(a2, MAXN, -1);
for (int i = 0; i < MAXN; ++i) {
int cnt1 = round(tot[i].real()); // ABC, with permutation
int cnt2 = round(a2[i].real());  // AAB, no permutation
int cnt3 = a3[i];                // AAA
int cnt = (cnt1 - cnt2 * 3 + cnt3 * 2) / 6;
if (cnt > 0)
cout << i - OFFSET * 3 << " : " << cnt << endl;
}
}

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