详解 FFT

Createsj · 2020-04-07 19:07:56 · 个人记录

前置知识：多项式，分治。

应用场景：多项式乘法。

温馨提示：本文证明不必掌握，仅供想要了解的人阅读。

〇、导入

您一定算过多项式乘法吧！有的时候，这算起来比较麻烦，比如：

用 \LaTeX 表示就更麻烦了。

在实际应用上，有时面对的多项式甚至多达上万项！这个时候再人工手算效率过低，且容易出错。幸好，我们已经有了计算机，能够用非常快的速度算出结果！

暴力算法是很容易想到的：

for(int i=0;i<n;++i)
    for(int j=0;j<m;++j)
        c[i+j]+=a[i]*b[j];

但它的时间复杂度为 \Theta(n^2) 级，如果遇到上万的数据就容易被卡了。

有没有更快的方法呢？当然有！那就是我们现在要讲的 FFT ！

一、知识补充

1. 多项式

1-1 多项式的一般表达

我们通常用 F(x) 来表示一个多项式，定义一个多项式只需用 F(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_0 ，如 F(x)=x^2+3x-5 。可以把它理解为函数，比如 F(2) 就是将 x=2 代入多项式 F(x) 后的值。

1-2 多项式的点值表达

我们知道，在平面直角坐标系中，n+1 个不重合的点可以唯一确定一个一元 n 次多项式。

所以我们可以用 n+1 个点值来表示一个一元 n 次多项式！

那么如何通过点值表达来计算多项式乘法呢？

设已知两个一元 n 次多项式 F(x),G(x) 的点值表达，W(x)=F(x)\times G(x) 。很显然，多项式 W(x) 在 x=i 时的点值为 W(i)=F(i)\times G(i) 。

但是，就像前面所说的那样，n+1 个不重合的点可以唯一确定一个一元 n 次多项式，而 W(x) 显然是个一元 2n 次多项式，n+1 个点值是根本不够的。怎么办？

很简单，只需先分别算出 F(i),G(i) 的 2n+1 个点值就行了！

1-3 多项式乘法的本质——卷积

接下来我们来思考一下乘法的本质。

容易想到，积的每一项的系数都可以表示为：

对于

这样的式子，我们称之为卷积。（了解更多）

可以看出，多项式乘法其实就是加法卷积。

2. 复数

2-1 定义

学复数之前，我们所接触的数仅限于实数范围内。

从自然数到整数、有理数再到实数，数的天地越来越大。现在，我们将数系扩充到复数范围内：

复数

• 实数
  • 有理数
  • 整数
    • 正数
    • 0
    • 负数
  • 分数
  • 无理数
• 虚数

我们用 \mathrm i 来表示虚数单位 \sqrt{-1} ，复数用 a+b\mathrm i(a,b\in\mathbf R) 表示，其中 a 为实部，b 为虚部。

我们定义复数 a+b\mathrm i 与 a-b\mathrm i 互为共轭复数。（通常记复数 z 的共轭复数为 \overline z）

复数 z=a+b\mathrm i 可以分类如下：

2-2 几何意义

当我们还在学实数的时候，数轴是这样的：

现在又多了复数。于是，我们需要用一个平面直角坐标系来表示数轴了：

这个建立了直角坐标系来表示复数的平面叫做复平面，x 轴叫做实轴，y 轴叫做虚轴。显然，在实轴上的点都表示实数；除了原点外，虚轴上的点都表示纯虚数。复平面内的点 (a,b) 表示虚数 a+b\mathrm i 。

由上图可以看出，复数 z=a+bi 、点 Z(a,b) 和 平面向量 \overrightarrow{OZ} 一一对应。我们将向量 \overrightarrow{OZ} 的模 r 叫做复数 z 的模，记作 \lvert z\rvert 。由模的定义可知：

为方便起见，我们常把复数 z=a+b\mathrm i 说成点 Z 或说成向量 \overrightarrow{OZ} ，并且规定，相等的向量表示同一个复数。

辐角

当 z\not=0 时，向量 \overrightarrow{OZ} 与实轴正向的夹角称为复数 z 的辐角，记作 \operatorname{Arg}z 。辐角的符号规定为：由正实轴依反时针方向转到 \overrightarrow{OZ} 为正，依顺时针方向转到 \overrightarrow{OZ} 为负。

显然一个非零复数 z 的辐角有无穷多个值，它们相差 \pi 的整数倍，但 \operatorname{Arg}z 中只有一个值 \theta_0 满足条件 -\pi\le\theta_0\le\pi ，称为复数的主辐角，记为 \arg z ，于是

这里我们使用的是弧度制，单位为 \mathrm{rad}，可以省略不写。\pi~\mathrm{rad}=180\degree.

2-3 四则运算

2-3-1 公式

2-3-2 几何意义

复数的加减法与向量的加减法对应，很容易理解。

复数相乘时，模长相乘，辐角相加。

证明：

设有两个复数 z_1=a+b\mathrm i,z_2=c+d\mathrm i 分别对应点 A(a,b),B(c,d) ，它们的积为 z_1z_2=(ac-bd)+(ad+bc)\mathrm i ，对应点 C(ac-bd,ad+bc) 。

• 模长相乘： \begin{aligned} \because~& \left\{ \begin{aligned} \lvert z_1\rvert &=\sqrt{a^2+b^2},\\ \lvert z_2\rvert&=\sqrt{c^2+d^2},\\ \end{aligned} \right.\\ \therefore~& \begin{aligned} \lvert z_1\rvert\lvert z_2\rvert&=\sqrt{a^2+b^2}\times\sqrt{c^2+d^2}\\ &=\sqrt{(a^2+b^2)(c^2+d^2)}, \end{aligned}\\ \text{又}\because~& \begin{aligned} \lvert z_1z_2\rvert&=\sqrt{(ac-bd)^2+(ad+bc)^2}\\ &=\sqrt{a^2c^2+b^2d^2+a^2d^2+b^2c^2}\\ &=\sqrt{(a^2+b^2)(c^2+d^2)}, \end{aligned}\\ \therefore~& \lvert z_1z_2\rvert=\lvert z_1\rvert\lvert z_2\rvert. \end{aligned}
• 辐角相加：设有点 P(1,0) ，连接 OP,PB,AC 。 \begin{aligned} \because~ &\begin{aligned} AC^2&=[a-(ac-bd)]^2+[b-(ad+bc)]^2\\ &=a^2c^2+a^2d^2+b^2c^2+b^2d^2+a^2+b^2-2a^2c-2b^2c\\ &=a^2(c^2+d^2)+b^2(c^2+d^2)+a^2+b^2-2(a^2+b^2)c\\ &=(a^2+b^2)(c^2+d^2)+(a^2+b^2)(1-2c)\\ &=(a^2+b^2)(c^2+d^2-2c+1), \end{aligned}\\ &\begin{aligned} OA^2PB^2&=(a^2+b^2)[(1-c)^2+(-d)^2]\\ &=(a^2+b^2)(c^2+d^2-2c+1), \end{aligned}\\ \therefore~& AC:PB=OA,\\ \text{又}\because~& OC=OA\cdot OB,OP=1,\\ \therefore~& OC:OB=OA:OP=AC:PB=OA,\\ \therefore~&\Delta OAC\sim\Delta OBP,\\ \therefore~&\angle AOC=\angle BOP,\\ \therefore~&\arg z_1z_2=\angle AOC+\angle AOP=\arg z_1+\arg z_2. \end{aligned}

给出一个图方便大家验证：

3. 单位根

3-1 定义

我们称 n 次幂为 1 的复数为 n 次单位根，即方程

的复数解。

因为 n 次方程有且只有 n 个复数根，所以 n 次单位根一共有 n 个。

为了求出单位根，我们在复平面上放一个单位圆。单位圆即圆心在原点上且半径为 1 的圆。如图：

显然，单位圆是由所有模长为 1 的复数表示的点构成的。

关于一个复数 z 和 z^n 的关系：

• 当 \lvert z\rvert>1 时，\lvert z^n\rvert=\lvert z\rvert^n>1 ，
• 当 \lvert z\rvert<1 时，\lvert z^n\rvert=\lvert z\rvert^n<1 ，
• 当 \lvert z\rvert=1 时，\lvert z^n\rvert=\lvert z\rvert^n=1 。

所以，只有模长为 1 的复数可能为 n 次单位根，即所有单位根表示的点都在单位圆上。可以发现，n 次单位根是且只能是辐角为 \frac{2k\pi}n(0\le k<n,k\in\mathbf Z)（即 \frac 1 n 圆周）的复数。

所以这 n 个 n 次单位根 n 等分单位圆。

3-2 性质

0. 首先已知 \arg z\omega_n^1\pi\equiv(\arg z+\frac{2\pi}n)\ (\bmod\ 2\pi)\quad(z\not=0) ，即一个非零复数每乘一次 \omega_n^1 都相当于辐角增加 \frac 1 n 圆周。
1. 证明： $$ \begin{aligned} \because~&\arg \omega_n^k=\arg \omega_n^{k\bmod n}\\ \therefore~&\omega_n^k=\omega_n^{k\bmod n}. \end{aligned} $$
2. 证明：由性质 $0$ 可知 $\arg(\omega_n^1)^k=k\arg\omega_n^1=\frac{2k\pi}n=\arg\omega_n^k$ ，即 $(\omega_n^1)^k$ 相当于辐角为 $\frac k n$ 圆周的复数，即 $\omega_n^k$ 。
3. 证明：由性质 $0,1$ 可得。
4. 证明：由性质 $2$ 可得。
5. 证明：因为 $\arg\omega_n^k=\frac{2k\pi}{n},\arg\omega_{pn}^{pk}=\frac{2pk\pi}{pn}=\frac{2k\pi}{n}$ ，即 $\omega_{pn}^{pk}$ 的辐角为 $\frac{pk}{pn}$ 圆周，$\omega_n^k$ 的辐角为 $\frac k n$ 圆周，所以 $\arg\omega_n^k=\arg\omega_{pn}^{pk}$ 。
6. 证明：由性质 $3$ 可得 $\omega_n^{(k+n/2)}=\omega_n^k\times\omega_n^{(n/2)}=\omega_n^k\times-1=-\omega_n^k\quad(2|n)$ 。

3-3 单位根反演

在这里只是为了验证 FFT 的正确性，不想了解可以跳过。其中

是不是跟 C++ 的 bool 值一样？

• \frac 1 n\sum\limits_{i=0}^{n-1}\omega_n^{(a-b)i} =[a=b]\quad(a,b<n),

  证明：当 a=b 时，

  \frac 1 n\sum\limits_{i=0}^{n-1}\omega_n^{(a-b)i} =\frac 1 n\sum\limits_{i=0}^{n-1}\omega_n^0 =1,

  当 a\not=b 时，

  \frac 1 n\sum\limits_{i=0}^{n-1}\omega_n^{(a-b)i} =\frac 1 n\times\frac{\omega_n^{(a-b)n}-1}{\omega_n^{a-b}-1 } =\frac 1 n\times\frac{\omega_n^0-1}{\omega_n^{a-b}-1 } =0.
• \frac 1 a\sum\limits_{i=0}^{a-1}\omega_a^{bi} =[a|b],

  证明：当 a|b 时，

  \frac 1 a\sum\limits_{i=0}^{a-1}\omega_a^{bi} =\frac 1 a\sum\limits_{i=0}^{a-1}\omega_a^0 =1,

  当 a\not|b 时，

  \frac 1 a\sum\limits_{i=0}^{a-1}\omega_a^{bi} =\frac 1 a\sum\limits_{i=0}^{a-1}\omega_a^{bi\bmod a} =\frac 1 a\times\frac{\omega_a^{ab}-1}{\omega_a^b-1} =\frac 1 a\times\frac{\omega_a^0-1}{\omega_a^b-1} =0.

补充完了前置知识 ，感觉身体被掏空，让我们开始学习 FFT 的核心思想吧！

二、核心

利用 FFT 计算多项式乘法大致分为两步：DFT（离散傅里叶变换） 和 IDFT（离散傅里叶逆变换） 。其中 DFT 表示将多项式的系数表达转换为点值表达，IDFT 则是它的逆运算，即将多项式的点值表达转换为系数表达。FFT 通过分治将 DFT 加速到 \Theta(n\log n) ，能过 10^6 级的数据。

1. 思想

1-1 DFT

设有一一元 n-1(n=2^t,t\in\mathbf N) 次（即有 n 项）多项式

将它的每一项安装次数的奇偶分成两部分：

设有两个多项式

那么

将 \omega_n^k(k<n/2) 代入 F(x) 中，得

将 \omega_n^{k+n/2}(k<n/2) 代入 F(x) 中，得

如果我们知道 F_l(x) 和 F_r(x) 在 \omega_{n/2}^0,\omega_{n/2}^1,\cdots,\omega_{n/2}^{n/2-1} 上的点值表达，那么可以用

来求 F(x) 在 \omega_n^0,\omega_{n/2}^1,\cdots,\omega_n^{n/2-1} 上的点值表达，用

来求 F(x) 在 \omega_n^{n/2},\omega_{n/2}^1,\cdots,\omega_n{n-1} 上的点值表达。

也就是说，如果我们知道 F_l(x) 和 F_r(x) 在 \omega_{n/2}^0,\omega_{n/2}^1,\cdots,\omega_{n/2}^{n/2-1} 上的点值表达，能够 \Theta(n) 求出 F(x) 在 \omega_n^0,\omega_{n/2}^1,\cdots,\omega_n^{n-1} 上的点值表达。

但问题在于：我们不知道 F_l(x) 和 F_r(x) 在 \omega_{n/2}^0,\omega_{n/2}^1,\cdots,\omega_{n/2}^{n/2-1} 上的点值表达啊？

这时候就要用到分治了！

可以想到，用同样的方式处理 F_l(x) 和 F_r(x)，直到处理的多项式为常数（即只有一项）时就不需要再操作了。

这种通过分治加速 DFT 的方法即为 FFT 。

1-2 IDFT

设多项式 F(x)=\sum\limits_{i=0}^{n-1}a_ix^i 在 \omega_n^k 上的点值表达为 b_k ，则

只需记住

就可以了！

所以，IDFT 相当于 DFT 用 \omega_n^{-k} 代替 \omega_n^k 后再除以 n 。

证明：（需要用到单位根反演）

2. 初步实现

2-1 多项式

我们需要一个复数数组 a[n] 来储存一个 n-1 次多项式。用 a[i] 储存一个多项式的 i 次系数或在 \omega_n^i 上的点值表达。

前面说过，n=2^t,t\in\mathbf N ，即 n 为 2 的整数次幂。但在实际应用中，n 很可能不为 2 的整数次幂。怎么办？

很简单，如果它不为 2 的整数次幂，那就将它向上补到 2 的整数次幂，再将多出来的项的系数赋值为 0 ！

2-2 复数

STL 自带复数库 complex ，但为防止被卡，最好自己手写一个复数结构体。在 FFT 中我们只需要用到加减乘即可。

struct Complex
{
    double a,b;
    Complex(){}
    Complex(const double &ia,const double &ib){a=ia,b=ib;}
    // 重载运算符
    inline Complex operator+(const Complex &x){return Complex(a+x.a,b+x.b);}
    inline Complex operator-(const Complex &x){return Complex(a-x.a,b-x.b);}
    inline Complex operator*(const Complex &x){return Complex(a*x.a-b*x.b,a*x.b+b*x.a);}
};

2-3 单位根

这里需要一点三角函数的知识。

由单位根的性质可知

由此我们可以通过递推求出 \omega_n^k 。

另外，为了更高的精度，我们这样定义 \pi ：

const double pi=acos(-1.0);

2-4 递归实现 FFT

对比公式和代码有助于更好地理解。

void fft(Complex *&a,const int n,const bool &inv)
{
    if(n==1)    return;
    const int bn=n>>1;  // n/2
    Complex *la=a,*ra=a+bn;
    for(int k=0;k<n;++k)    tmp[k]=bg[k];
    // 根据次数的奇偶分成两部分
    for(int k=0;k<bn;++k)   la[k]=tmp[k<<1],ra[k]=tmp[(k<<1)|1];
    // 分治
    fft(bg,ra,inv),fft(ra,ed,inv);
    // w1 为第 1 个 n 次单位根，wk 为第 k 个 n 次单位根
    Complex w1(cos(2.0*pi/n),sin(2.0*pi/n)),wk(1.0,0.0),t;
    if(inv) w1.b=-w1.b;
    for(int k=0;k<bn;++k)
    {
        t=wk*ra[k];     // 由于复数乘法较慢，这里我用一个临时变量储存积
        tmp[k]=la[k]+t;
        tmp[k+bn]=la[k]-t;
        wk=wk*w1;
    }
    for(int k=0;k<n;++k)    a[k]=tmp[k];
}

2-5 蝴蝶变换

我们来看分治过程中各项的原下标：

0 1 2 3 4 5 6 7
0 2 4 6|1 3 5 7
0 4|2 6|1 5|3 7
0|4|2|6|1|5|3|7

下标变化如下：

0->0
1->4
2->2
3->6
4->1
5->5
6->3
7->7

可以发现是两两交换：

0<->0
1<->4
2<->2
3<->6
5<->5
7<->7

然而貌似并没有什么用……

转换成二进制试试：

000<->000
001<->100
010<->010
011<->110
101<->101
111<->111

发现玄机了吧？ ????

显然就是将下标的二进制翻转了！

用递推搞定：

for(int i=0;i<N;++i)
    rev[i]=(rev[i>>1]>>1)|((i&1)?(n>>1):0);

仔细看就能搞懂了。

2-6 合并优化

对于

for(int k=0;k<bn;++k)
{
    t=wk*ra[k];
    tmp[k]=la[k]+t;
    tmp[k+bn]=la[k]-t;
    wk=wk*w1;
}
for(int k=0;k<n;++k)    a[k]=tmp[k];

这部分，由于 la=a,ra=a+bn ，所以可以替换为

for(int k=0;k<bn;++k)
{
    t=wk*ra[k];
    // 注意：下面两条语句的顺序已交换，请按照此顺序编写
    a[k+bn]=la[k]-t;
    a[k]=la[k]+t;
    wk=wk*w1;
}

这样就可以避免大量的数组拷贝。

现在 FFT 实现如下：

void fft(Complex *&a,const int n,const bool &inv)
{
    if(n==1)    return;
    const int bn=n>>1;
    Complex *la=a,*ra=a+bn;
    fft(bg,ra,inv),fft(ra,ed,inv);
    Complex w1(cos(2.0*pi/n),sin(2.0*pi/n)),wk(1.0,0.0),t;
    if(inv) w1.b=-w1.b;
    for(int k=0;k<bn;++k)
    {
        t=wk*ra[k];
        a[k+bn]=la[k]-t;
        a[k]=la[k]+t;
        wk=wk*w1;
    }
}

请在 FFT 之前交换二进制翻转的下标：

for(int i=0;i<n;++i)
    if(i<rev[i])
        swap(a[i],a[rev[i]]);

2-7 迭代实现 FFT

从下往上合并，这样可以减少许多不必要的操作，使得时间得到进一步优化：

void fft(Complex *a,const bool &inv)
{
    // 蝴蝶变换
    for(int i=0;i<N;++i)
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    // 枚举 n
    for(int n=2;n<=N;n<<=1)
    {
        const int bn=n>>1;
        Complex w1(cos(2.0*pi/n),sin(2.0*pi/n));
        if(inv) w1.b=-w1.b;
        for(int l=0;l<N;l+=n)
        {
            Complex wk(1.0,0.0),t,*la=a+l,*ra=a+l+bn;
            for(int k=0;k<bn;++k)
            {
                t=wk*ra[k];
                // 这里用 la 替代 a
                la[k+bn]=la[k]-t;
                la[k]=la[k]+t;
                wk=wk*w1;
            }
        }
    }
}

三、应用

费了那么大的劲弄出来的 FFT，总不能仅仅是 DFT 后再 IDFT 吧？那不就什么也没干。让我们来应用吧！

1. 多项式乘法

模板题：P3803 【模板】多项式乘法（FFT）。

还记得之前说过的吗？

设已知两个一元 n 次多项式 F(x),G(x) 的点值表达，W(x)=F(x)\times G(x) 。很显然，多项式 W(x) 在 x=i 时的点值为 W(i)=F(i)\times G(i) 。

所以，我们分别对两个多项式 DFT 一遍，然后将两者的点值的积 IDFT 一下就行了！

我的 AC 代码：

#include <cmath>
#include <cstdio>
const double pi=acos(-1.0);
int N,M,rev[0x200005];
struct Complex
{
    double a,b;
    Complex(){}
    Complex(const double &ia,const double &ib){a=ia,b=ib;}
    inline Complex operator+(const Complex &x){return Complex(a+x.a,b+x.b);}
    inline Complex operator-(const Complex &x){return Complex(a-x.a,b-x.b);}
    inline Complex operator*(const Complex &x){return Complex(a*x.a-b*x.b,a*x.b+b*x.a);}
}fa[0x200005],ga[0x200005];
void swap(Complex &x,Complex &y){Complex t=x;x=y;y=t;}
int getint()
{
    char c;
    for(;c=getchar(),c<'0' || c>'9';);
    return c&0xf;
}
void putint(const int num)
{
    if(num>=10)
        putint(num/10);
    putchar(num%10|'0');
}
void fft(Complex *a,const bool &inv)
{
    for(register int i=0;i<N;++i)   if(i<rev[i])    swap(a[i],a[rev[i]]);
    for(register int n=2,l,k;n<=N;n<<=1)
    {
        const register int bn=n>>1;
        register Complex w1(cos(2.0*pi/n),sin(2.0*pi/n));
        if(inv) w1.b=-w1.b;
        for(l=0;l<N;l+=n)
        {
            register Complex wk(1.0,0.0),t,*la=a+l,*ra=a+l+bn;
            for(k=0;k<bn;++k)
            {
                t=wk*ra[k];
                ra[k]=la[k]-t;
                la[k]=la[k]+t;
                wk=wk*w1;
            }
        }
    }
}
inline void init()
{
    register int i;
    scanf("%d%d",&N,&M);
    for(i=0;i<=N;++i)   fa[i].a=getint();
    for(i=0;i<=M;++i)   ga[i].a=getint();
    for(M+=N,N=1;N<=M;N<<=1);
    for(i=0;i<N;++i)    rev[i]=(rev[i>>1]>>1)|((i&1)?(N>>1):0);
}
inline void solve()
{
    fft(fa,false),fft(ga,false);
    for(register int i=0;i<N;++i)   fa[i]=fa[i]*ga[i];
    fft(fa,true);
    for(register int i=0;i<=M;++i)  putint((int)(fa[i].a/N+0.49)),putchar(' ');
}
int main()
{
    init();
    solve();
    return 0;
}

2. 高精度乘法

模板题：P1919 【模板】A*B Problem升级版（FFT快速傅里叶）。

高精度乘法其实就是将多项式乘法稍稍改动了而已。

设有两个自然数 a,b ，两个数字分别有 n,m 位。设存在一元 n 次多项式 F(x) 和一元 m 次多项式 G(x) ，使 F(10)=a,G(10)=b ，则 a\times b=(F\times G)(10) 。

The End

参考：

• 傅里叶变换(FFT)学习笔记
• 【多项式】FFT
• 百度百科
• 高中《数学》