快速傅里叶变换

【注】参考自算法导论。

1. 简介

快速傅里叶变换（FFT）是实现离散傅里叶变换（DFT）和离散傅里叶逆变换（IDFT）的快速算法，其时间复杂度为 O(nlog⁡n)O(n \log n)。DFT 在实际生活中有很多应用，比如通过离散傅里叶变换，可以将系数表示的多项式转为点值表示的多项式，从而使得多项式的乘法的复杂度由 O(n2)O(n^2) 降为 O(n)O(n) 。

2. 实现

快速傅里叶变换涉及到离散傅里叶变换的知识，FFT 能够实现的基础在于折半定理，折半定理提供了求大整数的 DFT 的一种思路：即通过将大整数不断划分为两部分，然后分别求出两部分的 DFT，最后在合并，本质是一种分治的思想。

3. 代码

3.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
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
#include <bits/stdc++.h>
using namespace std;

// 复数
struct Complex{
    double re, im;
    Complex(): re(0), im(0) {}
    Complex(double _re, double _im): re(_re), im(_im) {}
    Complex operator + (const Complex cp) const {
        return Complex(re+cp.re, im+cp.im);
    }
    Complex operator - (const Complex cp) const {
        return Complex(re-cp.re, im-cp.im);
    }
    Complex operator * (const Complex cp) const {
        return Complex(re*cp.re - im*cp.im, re*cp.im + im*cp.re);
    }
    Complex operator / (const double div) const {
        return Complex(re/div, im/div);
    }
    void operator = (const Complex cp) {
        re = cp.re;
        im = cp.im;
    }
};
// 快速傅里叶变换
struct FFT{
    #ifndef _FFT_
    #define _FFT_
    #define ll int
    #define MAXN 3000005
    #endif
    // 待计算
    ll cnt;         // 位数
    ll len;         // 数组大小（len = 2^cnt）
    ll r[MAXN];     // 位逆序置换数组
    const double Pi = acos(-1.0);
    FFT(): cnt(0), len(0) {}
    // 构建 2 的幂次的系数数组
    void build(Complex a[], ll curlen) {
        // 计算 len 所占位数
        cnt = 1;
        len = curlen;
        while((len >>= 1))   ++cnt;
        // 找到最小的大于等于 curlen 的 2 的幂值
        len = (ll)1 << cnt;    
        if(len-curlen == curlen) {
            --cnt;
            len = curlen;
        }
        // 超出 curlen 部分的系数补 0 
        for(ll i = curlen; i < len; ++i)    a[i] = Complex(0,0);
    }
    // 计算位逆序置换
    void bitReverse() {
        for(ll i = 0; i < len; ++i) r[i] = i;
        // 本质是动态规划
        // 利用 i>>1 的位逆序推导出 i 的位逆序
        for(ll i = 0; i < len; ++i) r[i] = (r[i>>1]>>1) | ((i&1)<<(cnt-1));
    }
    // 快速傅里叶变换
    void DFT(Complex y[]) {
        // 位逆序置换
        for(ll i = 0; i < len; ++i) {
            if(i < r[i])    swap(y[i], y[r[i]]);
        }
        // 合并 DFT
        // 枚举合并长度
        for(ll i = 1; i < len; i<<=1) {
            // 主单位复根
            Complex wn(cos(Pi/i), sin(Pi/i));
            // 枚举合并组首下标
            for(ll j = 0; j < len; j += (i<<1)) {
                Complex w(1,0);
                // 合并组 DFT
                for(ll k = j; k < j+i; ++k, w=w*wn) {
                    Complex u = y[k];
                    Complex v = w*y[k+i];
                    y[k] = u + v;
                    y[k+i] = u - v;
                }
            }
        }
    }
    // 快速傅里叶逆变换
    void IDFT(Complex a[]) {
        // 位逆序置换
        for(ll i = 0; i < len; ++i) {
            if(i < r[i])    swap(a[i], a[r[i]]);
        }
        // 合并 DFT
        // 枚举合并长度
        for(ll i = 1; i < len; i<<=1) {
            // 主单位复根
            Complex wn(cos(Pi/i), -sin(Pi/i));
            // 枚举合并组首下标
            for(ll j = 0; j < len; j += (i<<1)) {
                Complex w(1,0);
                // 合并组 DFT
                for(ll k = j; k < j+i; ++k, w=w*wn) {
                    Complex u = a[k];
                    Complex v = w*a[k+i];
                    a[k] = u + v;
                    a[k+i] = u - v;
                }
            }
        }
        for(ll i = 0; i < len; ++i) {
            a[i] = a[i]/len;
        }
    }
};

3.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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
#include <bits/stdc++.h>
using namespace std;

// 复数
struct Complex{
    double re, im;
    Complex(): re(0), im(0) {}
    Complex(double _re, double _im): re(_re), im(_im) {}
    Complex operator + (const Complex cp) const {
        return Complex(re+cp.re, im+cp.im);
    }
    Complex operator - (const Complex cp) const {
        return Complex(re-cp.re, im-cp.im);
    }
    Complex operator * (const Complex cp) const {
        return Complex(re*cp.re - im*cp.im, re*cp.im + im*cp.re);
    }
    Complex operator / (const double div) const {
        return Complex(re/div, im/div);
    }
    void operator = (const Complex cp) {
        re = cp.re;
        im = cp.im;
    }
};
// 大整数乘法（快速傅里叶变换实现）
struct BigIntegersMultiply{
    #ifndef _BIGINTEGERSMULTIPLY_
    #define _BIGINTEGERSMULTIPLY_
    #define ll int
    #define MAXN 3000005
    #endif
    const double Pi = acos(-1.0);
    // 待计算
    ll cnt;         // 位数
    ll len;         // 最终乘积数组大小（len = 2^cnt）
    ll r[MAXN];     // 位逆序置换数组
    ll ans[MAXN];   // 最终答案（逆序存储）
    // 待输入
    ll len1, len2;  // 初始系数数组大小
    Complex x1[MAXN], x2[MAXN]; // 两个大整数乘数（逆序储存）

    BigIntegersMultiply(): cnt(0), len(0) {}
    // 构建 2 的幂次的系数数组
    void build() {
        // 计算两个乘数最小的公共 2 的幂次 len 及所占位数 cnt
        cnt = 0, len = 1;
        while(len < (len1<<1) || len < (len2<<1))   ++cnt, len <<= 1;
        // 超出 curlen 部分的系数补 0 
        for(ll i = len1; i < len; ++i)  x1[i] = Complex(0,0);
        for(ll i = len2; i < len; ++i)  x2[i] = Complex(0,0);
    }
    // 计算位逆序置换
    void bitReverse() {
        for(ll i = 0; i < len; ++i) r[i] = i;
        // 本质是动态规划
        // 利用 i>>1 的位逆序推导出 i 的位逆序
        for(ll i = 0; i < len; ++i) r[i] = (r[i>>1]>>1) | ((i&1)<<(cnt-1));
    }
    // 快速傅里叶变换
    void DFT(Complex x[]) {
        // 位逆序置换
        for(ll i = 0; i < len; ++i) {
            if(i < r[i])    swap(x[i], x[r[i]]);
        }
        // 合并 DFT
        // 枚举合并长度
        for(ll i = 1; i < len; i<<=1) {
            // 主单位复根
            Complex wn(cos(Pi/i), sin(Pi/i));
            // 枚举合并组首下标
            for(ll j = 0; j < len; j += (i<<1)) {
                Complex w(1,0);
                // 合并组 DFT
                for(ll k = j; k < j+i; ++k, w=w*wn) {
                    Complex u = x[k];
                    Complex v = w*x[k+i];
                    x[k] = u + v;
                    x[k+i] = u - v;
                }
            }
        }
    }
    // 快速傅里叶逆变换
    void IDFT(Complex x[]) {
        // 位逆序置换
        for(ll i = 0; i < len; ++i) {
            if(i < r[i])    swap(x[i], x[r[i]]);
        }
        // 合并 DFT
        // 枚举合并长度
        for(ll i = 1; i < len; i<<=1) {
            // 主单位复根
            Complex wn(cos(Pi/i), -sin(Pi/i));
            // 枚举合并组首下标
            for(ll j = 0; j < len; j += (i<<1)) {
                Complex w(1,0);
                // 合并组 DFT
                for(ll k = j; k < j+i; ++k, w=w*wn) {
                    Complex u = x[k];
                    Complex v = w*x[k+i];
                    x[k] = u + v;
                    x[k+i] = u - v;
                }
            }
        }
        for(ll i = 0; i < len; ++i) {
            x[i] = x[i]/len;
        }
    }
    // 大整数乘法
    void bigMultiply() {
        build();
        bitReverse();
        DFT(x1), DFT(x2);
        for(ll i = 0; i < len; ++i) x1[i] = x1[i]*x2[i];
        IDFT(x1);
        for(ll i = 0; i < len; ++i) ans[i] = ll(x1[i].re + 0.5);
        for(ll i = 0; i < len; ++i) ans[i+1] += ans[i]/10, ans[i] %= 10;
        len = len1 + len2 - 1;
        while(len && !ans[len]) --len;
        ++len;
    }
};