快速傅立叶变换（FFT）_idftkw

FFT

作用：快速求两个多项式的乘积/卷积

前置知识

多项式乘法（点表示法）：用任意n + 1个不同点，均可确定为一个n次多项式。

好处：如果有两个多项式A,B:
$$\begin{gathered} A(x)=a_0+a_1{x}^1+a_2{x}^2+a_3{x}^3+.... \\ B(x)=b_0+b_1{x}^1+b_2{x}^2+b_3{x}^3+.... \\ 多项式A点表示法为(x_1,A(x_1)),(x_2,A(x_2),(x_3,A(x_3)))... \\ 多项式B点表示法为(x_1,B(x_1)),(x_2,B(x_2),(x_3,B(x_3)))... \\ 那么C可以表示为A(x_i)\ast B(x_i); \\ 多项式C点表示法为(x_1,C(x_1)),(x_2,C(x_2),(x_3,C(x_3)))...（o(n)） \end{gathered}$$

复数（Complex）

定义：1.形如$a+bi$​​,其中a为实部，b为虚部。

​ 2.模长为：$\sqrt{a^2+b^2}$

​ 3.幅角：设逆时针方向为正方向，从x正半轴到向量的转角

运算法则：（满足结合律，交换律，分配律）

​ 1.加法满足平行四边形法则

​ 2.乘法：几何意义：幅角相加，模长相乘

​ 代数意义：$(ac-bd)+(ad+bc)i$

单位根

定义：$w_n^k$​​为n次单位根

性质：

• $w_n^j\ne w_n^i({\forall }_{i\ne j})$
• $w_n^k=cos\frac{2k\pi }{n}+sin\frac{2k\pi }{n}$
• $w_n^0=w_n^n=1$
• $w_{2n}^{2k}=w_n^k$​​（折半引理）
• $w_n^{k+\frac{n}{2}}=-w_n^k$​（消去引理）

离散傅立叶变换（Discrete Fourier Transform , DFT）

将系数矩阵转换为点表达式

设一个长度为n的数列$d[i]$,这个数列的第k项为$A(x)$在n次单位根第k次幂处的点值。

得：
$$d_k=\sum _{i=0}^{n-1}{a}_i\times w_n^i$$
但是这样仍然是$o(n^2)$​​，由单位根的性质，我们可得快速傅立叶变换(Fast Fourier transform , FFT)。

快速傅立叶变换(Fast Fourier transform , FFT)

对于$A(x)=a_0+a_1{x}^1+a_2{x}^2+a_3{x}^3+....+a_{n-1}{x}^{n-1}$​

按照奇偶分组：
$$\begin{gathered} A(x)=a_0+a_1{x}^1+a_2{x}^2+a_3{x}^3+....+a_{n-1}{x}^{n-1} \\ A(x)=(a_0+a_2{x}^2+...+a_{n-2}{x}^{n-2})+x(a_1+a_3{x}^2+...+a_{n-1}{x}^{n-2}) \\ 以x^2代x,不妨令： \\ {A}_1(x)=(a_0+a_{2}x+a_4{x}^2+...+a_{n-2}{x}^{\frac{n-2}{2}}) \\ {A}_2(x)=(a_1+a_{3}x+a_5{x}^2+...+a_{n-1}{x}^{\frac{n-2}{2}}) \\ 显然:A(x)=A_1(x^2)+x{A}_2(x^2) \end{gathered}$$
根据K的大小分类讨论： n

1. $k\in [0,\frac{n}{2}-1]$

$$\begin{gathered} A(w_n^k)=A_1(w_n^{2k})+w_n^k{A}_2(w_n^{2k}) \\ 由折半引理得：A(w_n^k)=A_1(w_{\frac{n}{2}}^k)+w_n^k{A}_2(w_{\frac{n}{2}}^k) \end{gathered}$$

2.$k+\frac{n}{2}\in [\frac{n}{2},n-1]$​
$$\begin{gathered} A(w_n^{k+\frac{n}{2}})=A_1(w_n^{2k+n})+w_n^{k+\frac{n}{2}}{A}_2(w_n^{2k+n}) \\ 由消去引理得：A(w_n^{k+\frac{n}{2}})=A_1(w_{\frac{n}{2}}^k)-w_n^k{A}_2(w_{\frac{n}{2}}^k) \end{gathered}$$
综上得： 若已知 $w_{\frac{n}{2}}^i$ 那么可以快速求得$A(x)$点表达式。

而$A_1(x),A_2(x)$,为$A(x)$的一半规模，可以递归求解问题，那么时间复杂度为$o(nlogn)$

离散傅立叶逆变换（Inverse Discrete Fourier Transform , IDFT）

既然有了$o(nlogn)$的将系数表达式转化为点表达式的方法，那么我们自然也想用一个理想的方式将点表达式转化为系数表达式。

已知多项式$A(x)={\sum }_{i=0}^{n-1}{a}_i{x}^i$,和其点表达式$d_k={\sum }_{i=0}^{n-1}{a}_i\times w_n^{ik}$,求$a_i$​.

结论：
$$c_k=\sum _{i=0}^{n-1}{d}_i(w_n^{-k}{)}^i$$
证明：
$$\begin{gathered} 构造一多项式F(x)=d_0+d_{1}x+d_2{x}^2+...+d_{n-1}{x}^{n-1} \\ 设C_k为x=w_n^{-k}下的点表示,即c_k=\sum _{i=0}^{n-1}{d}_i(w_n^{-k}{)}^i \\ 即多项式B(x)=d_0,d_{1}x,d_2{x}_2,...,d_{n-1}{x}^{n-1}在w_n^0,w_n^{-1},....,w_{n-1}^{-(n-1)}的点表达式 \\ (c_0,c_1,c_2,c_3,...,c_{n-1})满足: \\ {c}_k=\sum _{i=0}^{n-1}{d}_i(w_n^{-k}{)}^i\rightleftharpoons d_i=\sum _{j=0}^{n-1}{a}_j\times (w_n^i{)}^j \\ {c}_k=\sum _{i=0}^{n-1}{d}_i(w_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j(w_n^i{)}^j)(w_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j(w_n^j{)}^i)(w_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j(w_n^j{)}^i(w_n^{-k}{)}^i) \\ =\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j(w_n^{j-k}{)}^i \\ =\sum _{j=0}^{n-1}{a}_j(\sum _{i=0}^{n-1}(w_n^{j-k}{)}^i) \\ 设S(x)=\sum _{i=0}^{n-1}{x}_i \\ 代入x=w_n^k，得：S(w_n^k)=1+(w_n^k)+(w_n^k{)}^2+...+(w_n^k{)}^{n-1} \\ 分情况讨论得： \\ 当k\ne 0时，两边同时\times w_n^k,得 \\ {w}_n^{k}S(w_n^k)=w_n^k+(w_n^k{)}^2+(w_n^k{)}^3+...+(w_n^k{)}^n \\ {w}_n^{k}s(w_n^k)-s(w_n^k)=(w_n^k{)}^n-1 \\ s(w_n^k)=\frac{(w_n^k{)}^n-1}{w_n^k-1} \\ s(w_n^k)=\frac{(w_n^n{)}^k-1}{w_n^k-1} \\ s(w_n^k)=\frac{0}{w_n^k-1} \\ s(w_n^k)=0 \\ 当k=0时，s(w_n^k)=1+1+...+1=n \\ 由于c_k=\sum _{j=0}^{n-1}{a}_j(\sum _{i=0}^{n-1}(w_n^{j-k}{)}^i) \\ 当且仅当j=k时,c_k=n{a}_k \\ {a}_k=\frac{c_k}{n} \\ {c}_k=\sum _{i=0}^{n-1}{d}_i{w}_n^{-ki} \\ {a}_k=\frac{1}{n}\sum _{i=0}^{n-1}{d}_i{w}_n^{-ki} \end{gathered}$$
这个是DFT的逆变换，称为IDFT,可以通过这个式子求出多项式的系数表示法.

快速傅立叶逆变换（Inverse Fast Fourier transform,IFFT）

我们考虑的问题是如何快速求解：$A(x)={\sum }_{i=0}^{n-1}{b}_i{x}^i$在$w_n^{-ki}于[0,n)$​处的点值。

从图上，我们得到：
$$w_n^{-k}=w_n^{k+n}$$

只要先算出$A(X)$在$w_n$​在各个幂次的弧度，然后反过来计算就可以
$$a_k=\frac{1}{n}\sum _{i=0}^{n}A(w_n^k)(w_n^{-k}{)}^{-1}$$
至此通过FFT,IFFT我们可以$o(nlogn)$的时间复杂度由系数表达式变成点表达式，由点表达式变成系数表达式。

蝴蝶变换

由于递归写法常数太大，我们采用神奇的蝴蝶变换，就可以让他十分高效。
$$\begin{gathered} 第一层：a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7 \\ 第二层：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 \end{gathered}$$
可以发现:
$$\begin{gathered} 原来：01234567 \\ 变换：04261537 \\ 实际上就是二进制位翻转 \end{gathered}$$

$$递推式：rev[i]=(rev[i>>1]>>1)|(i\&1)<<(bit-1)$$

OK，到这就完全结束了。

例题

例题：多项式乘法

y总代码（注释版)：

#include <bits/stdc++.h>
using namespace std;

const int N = 300010;
const double PI = acos(-1);

int n, m;
struct Complex
{
    double x, y;
  	//重载复数加法
    Complex operator+ (const Complex& t) const
    {
        return {x + t.x, y + t.y};
    }
  	//重载复数减法 
    Complex operator- (const Complex& t) const
    {
        return {x - t.x, y - t.y};
    }
  	//重载复数乘法
    Complex operator* (const Complex& t) const
    {
        return {x * t.x - y * t.y, x * t.y + y * t.x};
    }
}a[N], b[N];//A，B，两个多项式
//rev存的是逆转后的序列，也就是递归后的位置。
//bit为最高次幂。
//tot为总共多少项，补充到2的整数次项。
int rev[N], bit, tot;

void fft(Complex a[], int inv)
{
  	//先将序列按照rev排序
    for (int i = 0; i < tot; i ++ )
        if (i < rev[i])
            swap(a[i], a[rev[i]]);
  	//枚举每个分块的长度，开始迭代，一开始是1.
    for (int mid = 1; mid < tot; mid <<= 1)
    {
      	//w1是单位根，乘inv是为了方便倒着作fft
        auto w1 = Complex({cos(PI / mid), inv * sin(PI / mid)});
      	//从第i块开始枚举，每次加两块
        for (int i = 0; i < tot; i += mid * 2)
        {
          	//从w0开始枚举
            auto wk = Complex({1, 0});
          	//每次都++,算出x,y，存下来
            for (int j = 0; j < mid; j ++, wk = wk * w1)
            {
                auto x = a[i + j], y = wk * a[i + j + mid];
                a[i + j] = x + y, a[i + j + mid] = x - y;
            }
        }
    }
}

int main()
{
    scanf("%d%d", &n, &m);
    for (int i = 0; i <= n; i ++ ) scanf("%lf", &a[i].x);
    for (int i = 0; i <= m; i ++ ) scanf("%lf", &b[i].x);
    while ((1 << bit) < n + m + 1) bit ++;
    tot = 1 << bit;
    for (int i = 0; i < tot; i ++ )
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
    fft(a, 1), fft(b, 1);
    for (int i = 0; i < tot; i ++ ) a[i] = a[i] * b[i];
    fft(a, -1);
    for (int i = 0; i <= n + m; i ++ )
        printf("%d ", (int)(a[i].x / tot + 0.5));

    return 0;
}
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

参考：acwing，
【学习笔记】超简单的快速傅里叶变换（FFT）（含全套证明）