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

· · 题解

快速傅里叶变换

序言 & FFT 简介

FFT，即快速傅里叶变换，是离散傅里叶变换的快速算法，可以用来显著的加快多项式乘法（将原来朴素算法的 \Theta(n^2) 优化为 \Theta(n\log(n))），它的用处还有音频处理，图像滤波，调制解调等。

FFT 计算乘法的主要流程是先将系数表示法转化为点值表示法，然后再通过点值表示法计算完之后转化回系数表示法。

FFT 作为数论的一个重要板块，有着很多的前置知识，也有一定的思考难度，接下来就让我为大家分部分讲解 FFT 的前置知识与他本身。

前置知识

多项式的表示法

  1. 系数表示法 这种方法算得上是最常见的多项式表示法，一般的时候我们也用得最多，形如 f(x)=\sum^{n}_{k=0}{a_kx^k}。 以这种方式计算多项式乘法时间复杂度为 \Theta(n^2)。

  2. 点值表示法 这种方法类似于初中的给几个点求函数解析式的情景，一个 n 次多项式可以使用 n+1 个点来表示出来，即这个多项式被点 (x_1,y_1),(x_2,y_2),...,(x_{n+1},y_{n+1}) 唯一表示。 但是这样计算复杂度仍然为 \Theta(n^2)。

单位根

在讲单位根之前，我们先讲一下单位圆，单位圆就是以坐标轴为原点，半径为 1 的圆。 而单位根就是把一个圆 n 等分，取这 n 个点所表示的复数，就是这个圆的 n 次单位根。

根据定义，找单位根的方法就呼之欲出啦： 先将一个圆 n 等分，然后从 (1,0)（即 \omega^1_n）这个点开始逆时针编号，第 k 个点所对应的复数为 \omega^k_n。

注：为了便于实现，下文中的 n 将默认为 2

而根据复数的乘法法则可得：

(\omega^m_n)^k=\omega^{mk}_n

根据每个复数的幅角，可以计算出这个点的坐标，例如 \omega^k_n 所对应的点为 (\cos\frac{2k\pi}{n},\sin\frac{2k\pi}{n})，也是复数 \cos\frac{2k\pi}{n}+sin\frac{2k\pi}{n}\times i。

接下来是单位根的性质：

  1. \omega^{2k}_{2n}=\omega^{k}_{n}

证明：

\begin{aligned}\omega_{2n}^{2k} &= \cos\frac{4k\pi}{2n}+\sin\frac{4k\pi}{2n}\times i \\&= \cos\frac{2k\pi}{n}+\sin\frac{2k\pi}{n}\times i \\&= \omega_{n}^{k}\end{aligned}
  1. \omega^{k+n}_{n}=\omega^{k}_{n}

证明：

\begin{aligned}\omega_{n}^{k + n} &= \cos\frac{2(k + n)\pi}{n} + \sin\frac{2(k + n)\pi}{n} \times i \\&= \cos(\frac{2k\pi}{n} + 2\pi) + (\sin\frac{2k\pi}{n} + 2\pi) \times i \\&= \cos\frac{2k\pi}{n} + \sin\frac{2k\pi}{n} \times i \\&= \omega_{n}^{k}\end{aligned}
  1. \omega^{k+\frac{n}{2}}_{n}=-\omega^{k}_{n}

证明：

\begin{aligned}\omega^{k+\frac{n}{2}}_{n} &= \cos\frac{2(k+\frac{n}{2})\pi}{n}+\sin\frac{2(k+\frac{n}{2})\pi}{n}\times i \\&= \cos(\frac{2k\pi}{n}+\pi)+(\sin\frac{2k\pi}{n}+\pi)\times i \\&= -\cos\frac{2k\pi}{n}-\sin\frac{2k\pi}{n}\times i \\&= -\omega^{k}_{n}\end{aligned}

FFT 推导（正确性证明）

假设一个 n 次多项式：

F(x)=a_0+a_1\times x+a_2\times x^2\cdots a_{n-1}\times x^{n-1}

将它的下标按奇偶性分类，并设：

F_1=a_0+a_2\times x+a_4\times x^2\cdots+a_{n-2}\times x^{\frac{n}{2}-1} F_2=a_1+a_3\times x+a_5\times x^2\cdots+a_{n-1}\times x^{\frac{n}{2}-1}

则有

F(x)=F_1(x^2)+xF_2(x^2)

分别将 \omega^k_n 与 \omega^{k+\frac{n}{2}}_{n} 代入得

\begin{aligned}F(\omega^k_n) &= F_1(\omega^{2k}_n)+\omega^k_nF_2(\omega^{2k}_n)\\&= F_1(\omega^{k}_{\frac{n}{2}})+\omega^k_nF_2(\omega^{k}_{\frac{n}{2}})\end{aligned}

和

F(\omega^{k+\frac{n}{2}}_{n})&=F_1(\omega^{2k+n}_{n})+\omega^{k+\frac{n}{2}}_{n}F_2(\omega^{2k+n}_{n})\\&=F_1(\omega^{2k}_{n})-\omega^k_nF_2(\omega^{2k+n}_{n})\\&=F_1(\omega^{k}_{\frac{n}{2}})-\omega^{k}_{n}F_2(\omega^{k}_{\frac{n}{2}})\end{aligned}

这两个式子长得极为相像，仅有一个常数的差距，意味着我们只需要求一半便可得知另一半的结果，如果我们将 F_1 与 F_2 也这么递归地计算下去，就可以得到递归版的代码，由于每层递归是 \Theta(n) 的并且由于每递归一次会少一半运算量，所以递归层数为 \Theta(\log(n))，总复杂度为 \Theta(n\log(n))。 但是，这样写常数巨大，十分的慢，这时候就要请出一个天才般的想法了，接下来，请观察分组前后的数组（以 n=8 为例）： 前：\{0,1,2,3,4,5,6,7\} 后：\{0,4,2,6,1,5,3,7\} 好像没什么规律？ 转化成二进制下试试！

\{(000)_2,(001)_2,(010)_2,(011)_2,(100)_2,(101)_2,(110)_2,(111)_2\} \{(000)_2,(100)_2,(010)_2,(110)_2,(001)_2,(101)_2,(011)_2,(111)_2\}

这不就是将分组前的数组每个数给翻转了吗，这样，我们就可以将翻转后的数组计算出来啦。 但，为什么可以这样呢？ 证明： 设 F(m,k) 表示一个序列中从 0 到 2^m-1 中第 k 个位置被分组后的位置。 现在将 k 按奇偶性分类讨论：

  1. 若 k 为奇数，则将其右移一位，并将它的二进制表示下的第 n-1 位修改成 1。

  2. 若 k 为偶数，则将其右移一位，并将它的二进制表示下的第 n-1 位修改成 0。

于是可得：

\begin{aligned}F(m,k) &=(k\&1)2^{n - 1}+((k\gg1)\&1)2^{n - 2}+\cdots+((k\gg n - 1)\&1)2^{0}\\&=\sum_{i = 0}^{n - 1}((k\gg i)\&1)2^{n - i - 1}\end{aligned}

即将 k 的二进制表示翻转。 于是就有了迭代版的代码：

for(int mid=1;mid<lim;mid<<=1){
        complex wn(cos(pie/mid),oper*sin(pie/mid));
        for(int R=mid<<1,j=0;j<lim;j+=R){
            complex w(1,0);
            for(int k=0;k<mid;k++,w=w*wn){
                complex y=x[j+k],z=w*x[j+mid+k];
                x[j+k]=y+z;
                x[j+mid+k]=y-z;
            }
        }
    }

但这还没完，还有要将点值表示法转回点值表示法的步骤没有完成，这就需要 IFFT（快速傅里叶逆变换）登场了。 使用矩阵来理解 IFFT 是一种比较简单的理解方式，从矩阵的角度看我们刚才的计算，实际上就是在进行如下操作：

\begin{bmatrix}F_0\\F_1\\F_2\\\vdots\\F_{n-1}\end{bmatrix}=\begin{bmatrix}1&1&1&\cdots&1\\1&\omega^{1}_{n}&\omega^{2}_{n}&\cdots&\omega^{n-1}_{n}\\1&\omega^{2}_{n}&\omega^{4}_{n}&\cdots&\omega^{2(n-1)}_{n}\\\vdots&\vdots&\vdots&\ddots&\vdots\\1&\omega^{n-1}_{n}&\omega^{2(n-1)}_{n}&\cdots&\omega^{(n-1)^2}_{n}\end{bmatrix}\begin{bmatrix}a_0\\a_1\\a_2\\\vdots\\a_{n-1}\end{bmatrix}

现在，我们已经完成了求左边式子的值这个任务，而中间的值在多项式的点值表示法中也是一一对应的，接下来只需要在两边左乘中间矩阵的逆矩阵即可。（因为矩阵乘法不满足交换律，所以把一个矩阵乘在另一个矩阵的左边和右边结果不同，乘在左边就是左乘）

由于这个矩阵的元素非常特殊，它的逆矩阵也有特殊的性质，就是每一项取倒数，再除以变换的长度 n，就能得到它的逆矩阵。

——引用自 oi wiki。 在代码实现上就是将原来 FFT 的单位根替换为它的倒数，然后计算出结果再除以 n 就好了。

实现代码

#include<bits/stdc++.h>
using namespace std;
constexpr int M=1e7+5e3;
const double pie=acos(-1.0);
struct comp{// 手打复数，避免被卡（其实一般也没有人来卡这个吧）。
    double real,img;
    comp(double xx=0,double yy=0){
        real=xx;
        img=yy;
    }
    comp operator + (const comp b)const{ 
        return comp(real+b.real , img+b.img);
    }
    comp operator - (const comp b)const{ 
        return comp(real-b.real , img-b.img);
    }
    comp operator * (const comp b)const{ 
        return comp(real*b.real-img*b.img , real*b.img+img*b.real);
    }
}f[M],g[M];
int n,m,l,lim=1;
int r[M];
void fft(comp *F,int oper){// oper 是操作，传入 1 是FFT，-1 是 IFFT。
    for(int i=0;i<lim;i++){
        if(i>r[i]-1){
            continue;
        }
        swap(F[i],F[r[i]]);
    }
    for(int mid=1;mid<lim;mid<<=1){
        comp wn(cos(pie/mid),oper*sin(pie/mid));
        for(int R=mid<<1,j=0;j<lim;j+=R){
            comp w(1.0,0.0);
            for(int k=0;k<mid;k++,w=w*wn){
                comp y=F[j+k],z=w*F[j+mid+k];
                F[j+k]=y+z;
                F[j+mid+k]=y-z;
            }
        }
    }
}
int main(){
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++){
        scanf("%lf",&f[i].real);
    }
    for(int i=0;i<=m;i++){
        scanf("%lf",&g[i].real);
    }
    while(lim<n+m+1){
        lim<<=1;
        l++;
    }
    for(int i=0;i<lim;i++){
        r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
    }
    fft(f,1);
    fft(g,1);
    for(int i=0;i<lim;i++){
        f[i]=f[i]*g[i];
    }
    fft(f,-1);
    for(int i=0;i<n+m+1;i++){
        printf("%d ",(int)(f[i].real/lim+0.5));// 防止精度爆炸。
    }
    return 0;
}

参考资料

oi wiki （部分内容参考与代码参考）。

attcak大佬的文章 （代码参考）。

_不会dp不改名_大佬的文章 （位逆序置换证明参考）。