多项式乘法与快速傅里叶变换(FFT)

一、前置知识

本节介绍多项式前置知识、复数及单位根的相关内容。

1.1   多项式

设 \(A(x)\) 为一个 \(n\) 次多项式，则可以表示为 \[ A(x)=\sum\limits_{i=0}^n a_ix^i \] 其中，\(a_i\) 为多项式第 \(i\) 项的系数。

一个多项式在 \(x_0\) 处的取值 \(A(x_0)\) 为其在 \(x_0\) 上的一个点值。

一个 \(n\) 次多项式可以用 \(n+1\) 个点值表示出来。由 \(n+1\) 个对应位置上的点值能唯一表示一个多项式。

形式化地，一个多项式可以由 \(n+1\) 个点 \((x_i,y_i)\) 唯一确定，其中，\(y_i=\sum\limits_{j=0}^n a_jx_i^j\)。

1.2   复数

设 \(a,b\in \mathbf R\)，令 \(i^2=-1\)，称形如 \(a+bi\) 的数为复数。其中，称 \(a\) 为复数的实部，\(b\) 为复数的虚部。\(a=0\) 的数称为纯虚数。

在二维平面中，用横坐标表示实部，用纵坐标表示虚部，这样的平面为复平面。在复平面中，对于复数 \(a+bi\)，令其坐标为 \((a,b)\)，定义其模长为 \(r=\sqrt{a^2+b^2}\)，即该点到坐标原点的距离。称 \(\theta\) 表示该连线与 \(x\) 轴非负半轴的正夹角为辐角，不难发现，一个复数可以用有序数对 \((\theta,r)\) 唯一表示。

称圆心在原点，半径为 \(1\) 的圆为单位圆，则单位圆上的复数可以唯一表示为 \((\cos\theta,\sin\theta)\)，其中 \(\theta\) 为辐角。

复数的加法   对于两个复数 \(x=a+bi\) 与 \(y=c+di\)，定义 \(x+y=(a+c)+(b+d)i\)。

复数的乘法   对于两个复数 \(x=a+bi\) 与 \(y=c+di\)，定义 \(xy=(a+bi)(c+di)=(ac-bd)+(bc+ad)i\)。其几何意义为辐角相加，模长相乘。

在 C++ 中给出了复数的标准库 <complex>，通过声明 #include<complex> 引入头文件，通过 complex<_type> a; 声明一个实部虚部为 _type 类型的复数。该元素分别有两个值，即 a={_real,_imag}。其中，_real 值为该复数的实部，_imag 值为复数的虚部。通过调用 a.real() 与 a.imag() 分别访问该复数的实部与虚部。

1.3   单位根

对于是 \(2\) 的正整数幂次的整数 \(n\)，记单位辐角为 \(\frac{2\pi} n\)，定义 \(\omega_n^k\) 表示一个模长为 \(1\)，辐角为 \(\frac{2k\pi}n\) 的复数，其坐标为 \((\cos\frac{2k\pi}n,\sin\frac{2k\pi}n)\)​。

单位根有以下性质：

  • \(\omega_{n}^{k}=\cos\frac{2k\pi}n+i\sin\frac{2k\pi}n\)。根据定义可知。
  • \(\omega_{n}^{k+n}=\omega_{n}^{k}\)。相当于辐角多转 \(2\pi\)，根据诱导公式可知三角函数值相等。
  • \(\omega_{2n}^{2k}=\omega_{n}^{k}\)。根据定义可知。
  • \(\omega_{n}^{k+\frac n2}=-\omega_{n}^{k}\)，相当于辐角多转 \(\pi\)​，根据诱导公式可知三角函数值为相反数。
  • \((\omega_{n}^{k})^m=\omega_{n}^{mk}\)，可以由复数乘法的几何意义推得。

二、快速傅里叶变换(FFT)

接下来只讨论多项式项数 \(n\) 为 \(2\) 的正整数幂次的多项式函数。

对于 \(n\) 项多项式函数 \(A(x)=\sum\limits_{i=0}^{n-1} a_ix^i\)，定义向量 \((a_0,a_1,\cdots,a_{n-1})\) 为其系数。按照每项次数进行奇偶分类，有 \[ A(x)=(a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^3+\cdots+a_{n-1}x^{n-1}) \] 分别定义 \[ \begin{aligned} A_1(x)&=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac n 2 -1} \\ A_2(x)&=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac n 2 -1} \end{aligned} \] 不难有 \[ A(x)=A_1(x^2)+xA_2(x^2) \] 分别求其在 \(\omega_{n}^{k}(0\le k <\frac n2)\) 上的点值，有 \[ \begin{aligned} A(\omega_{n}^{k})&=A_1(\omega_{n}^{2k})+\omega_{n}^{k}A_2(\omega_{n}^{2k}) \\ &=A_1(\omega_{n/2}^{k})+\omega_{n}^{k}A_2(\omega_{n/2}^{k}) \end{aligned}\tag 1 \] 与 \[ \begin{aligned} A(\omega_{n}^{k+n/2})&=A_1(\omega_{n}^{2k+n})+\omega_{n}^{k+n/2}A_2(\omega_{n}^{2k+n}) \\ &=A_1(\omega_{n/2}^{k})-\omega_{n}^{k}A_2(\omega_{n/2}^{k}) \end{aligned} \tag 2 \] 发现上两式的形式类似，而我们将问题从 \(n\) 的范围递归地划分到 \(\frac n 2\) 的范围内。这帮助我们在 \(O(n\log n)\) 的时间复杂度内求出 \(n\) 个点的点值，这 \(n\) 个点为 \(\omega_{n-1}^{0}\sim \omega_{n-1}^{n-1}\)。

上述过程称为快速傅里叶变换。

这时候我们已经可以通过递归的办法实现，但由于递归实现常数较大，下面介绍迭代实现的办法。

因为每一次递归都要奇偶分组，下图展示的原始系数数列下标与分组后的序列下标的关系。

1101696-20180212074250859-1560811086.png

发现后序列的二进制表示是原序列的二进制表示的反串（高低位颠倒）。我们可以按照这个性质对系数数列重新排序后迭代地计算。

即 \(trans_i\) 表示后序列第 \(i\) 个位置上对应的原序列下标，可以通过 \(trans_i=trans_{\lfloor n/2\rfloor}\times 2+(i\bmod 2)^{n-1}\) 得到。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
// typedef complex<double> cd;
void FFT(cd *a){ // a为多项式的系数向量
    for(int i=0;i<lim;i++) if(i<trans[i]) swap(a[i],a[trans[i]]); //下标置换
    for(int k=1;k<lim;k<<=1){ //枚举段长，合并两个 n'=k 的区间
        cd e={cos(pi/k),sin(pi/k)}; //定义单位根
        for(int i=0;i<lim;i+=k*2){
            cd w={1,0}; //当前的复数幂
            for(int j=0;j<k;j++){
                cd x=a[i+j],y=w*a[i+j+k];
                a[i+j]=x+y;
                a[i+j+k]=x-y;
                w*=e;
            }
        }
    }
}

三、快速傅里叶逆变换(IFFT)

对于多项式 \(A(x)=\sum\limits_{i=0}^{n-1}a_ix^i\)，其中 \((a_0,a_1,\cdots,a_{n-1})\) 为其系数向量。对其进行快速傅里叶变换后得到 \((y_0,y_1,y_2,\cdots,y_{n-1})\)，称为其 FFT 向量。通过 FFT 向量反解出系数向量的过程是快速傅里叶逆变换(IFFT)。

考虑多项式 \(B(x)=\sum\limits_{i=0}^{n-1}y_ix^i\)，即系数为 FFT 向量的新多项式。对于其在 \(\omega_{n}^{-k}\) 处的点值 \[ c_k=\sum\limits_{i=0}^{n-1}y_i(\omega_{n}^{-k})^i \] 我们试图找到 \(c_k\) 与 \(a_k\) 之间的关系，则可反解出 \(a_k\)。对上面式子进行求解： \[ \begin{aligned} c_k&=\sum_{i=0}^{n-1}y_i\left(\omega_{n}^{-k} \right)^i \\ &=\sum_{i=0}^{n-1}\left(\sum_{j=0}^{n-1}a_j(\omega_{n}^{i})^j \right)\left(\omega_{n}^{-k} \right)^i \\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j\left(\omega_{n}^{j} \right)^i\left(\omega_{n}^{-k} \right)^i \\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j\left(\omega_{n}^{j-k} \right)^i \\ &=\sum_{j=0}^{n-1}a_j\left(\sum_{i=0}^{n-1}\left(\omega_{n}^{j-k} \right)^i \right) \end{aligned} \tag 3 \] 令 \(S(x)=\sum\limits_{i=0}^{n-1}x^i\)，有 \[ \begin{aligned} S(\omega_{n}^{k})&=1+\omega_{n}^{k}+(\omega_{n}^{k})^2+\cdots+(\omega_{n}^{k})^{n-1} \\ \omega_{n}^{k}S(\omega_{n}^{k})&=\omega_{n}^{k}+(\omega_{n}^{k})^2+(\omega_{n}^{k})^2+\cdots+(\omega_{n}^{k})^n \end{aligned} \] 两式相减得到 \[ S(\omega_{n}^{k})=\dfrac{(\omega_{n}^{k})^n-1}{\omega_{n}^{k}-1}=\dfrac{1-1}{\omega_{n}^{k}-1} \] 当 \(\omega_{n}^{k}\not=1\)，即 \(k\not=0\) 时，\(S(\omega_{n}^{k})=0\)；

当 \(\omega_{n}^{k}=1\)，即 \(k=0\) 时，根据定义有 \(S(\omega_{n}^{k})=n\)。

对于 \((3)\) 式，用 \(S(\omega_{n}^{j-k})\) 表示 \(\sum\limits_{i=0}^{n-1}(\omega_{n}^{j-k})^i\)，因为其在 \(\omega_{n}^{0}\) 处才有值，所以 \[ \begin{aligned} c_k&=\sum_{j=0}^{n-1}a_jS(\omega_{n}^{j-k}) \\ &=na_k \end{aligned} \] 即 \[ a_k=\dfrac {c_k} n\tag 4 \] 这样，我们再次利用 FFT 求出多项式 \(B(x)\) 在 \(\omega_{n}^{-k}\) 上的点值 \(c_k\)，可以反解出 \(a_k\)，时间复杂度 \(O(n\log n)\)。

不难发现 IFFT 与 FFT 有许多相同之处，并且需要利用 FFT 求出点值 \(c_k\)，于是求解 IFFT 的时候直接更改 \(\omega_{n}^{k}\) 为 \(\omega_{n}^{-k}\) 即可利用 FFT 求解出 \(c_k\)，再除以 \(n\) 得到 \(a_k\)。

四、多项式乘法

对于两个多项式 \(A(x)=\sum\limits_{i=0}^n a_ix^i\) 和 \(B(x)=\sum\limits_{i=0}^m b_ix^i\)，定义多项式乘法 \(C(x)=A(x)\times B(x)\) 为 \[ c_{i}=\sum_{j=0}^{i}a_jb_{i-j} \] 其中 \(c_i\) 为多项式 \(C(x)\) 的第 \(i\) 项系数。所得的多项式 \(C(x)\) 是一个 \(n+m\) 次多项式，共有 \(n+m+1\) 项。

朴素的多项式乘法时间复杂度为 \(O(n^2)\) 级别，利用 FFT 可以在 \(O(n\log n)\) 的复杂度内求解，大致思路如下图所示。

算法流程如下：

  1. 分别对 \(A(x),B(x)\) 多项式求出其在 \(\omega_{n}^{k}\) 处的点值，利用 FFT 做到时间复杂度 \(O(n\log n)\)；
  2. 每个位置上的点值对应相乘，所得点值就是多项式 \(C(x)\) 在对应位置上的点值；
  3. 根据点值利用 IFFT 反解出 \(c_k\)，时间复杂度 \(O(n\log n)\)。

关键点在于将系数表示转换为点值表示，这样在对应位置上的点值只与当前位置有关，做到了位置分离，进而便于求解。

需要注意的是，在求解 FFT 与 IFFT 中，多项式项数 \(n\) 需要是 \(2\) 的正整数次幂，这要求我们在计算的时候需要将下标上届扩展到大于 \(n+m\) 的 \(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
// Luogu P3803 【模板】多项式乘法（FFT）
#include<iostream>
#include<complex>
#include<cmath>
using namespace std;
#define MAXN 4000005
typedef complex<double> cd;
const double pi=acos(-1.);
int n,m,lim,trans[MAXN];
cd a[MAXN],b[MAXN],c[MAXN];
void FFT(cd *a,int op){ //通过添加系数 op 将 FFT 与 IFFT 的前半部分合二为一
    for(int i=0;i<lim;i++) if(i<trans[i]) swap(a[i],a[trans[i]]);
    for(int k=1;k<lim;k<<=1){
        cd e={cos(pi/k),sin(pi/k)*op};
        for(int i=0;i<lim;i+=k*2){
            cd w={1,0};
            for(int j=0;j<k;j++){
                cd x=a[i+j],y=w*a[i+j+k];
                a[i+j]=x+y;
                a[i+j+k]=x-y;
                w*=e;
            }
        }
    }
}
inline int read(){
    int w=0;
    char ch=getchar();
    while(!isdigit(ch)) ch=getchar();
    while(isdigit(ch)){
        w=w*10+ch-'0';
        ch=getchar();
    }
    return w;
}
int main(){
    n=read(),m=read();
    for(int i=0;i<=n;i++) cin>>a[i];
    for(int i=0;i<=m;i++) cin>>b[i];
    int len=0;
    lim=1;while(lim<=n+m) lim<<=1,len++; //扩展 lim 到 2 的幂次
    for(int i=0;i<lim;i++) trans[i]=((trans[i>>1]>>1)|((i&1)<<(len-1)));
    FFT(a,1);
    FFT(b,1);
    for(int i=0;i<lim;i++) c[i]=a[i]*b[i];
    FFT(c,-1);
    for(int i=0;i<lim;i++) c[i]/=lim; // IFFT 的最后一步
    for(int i=0;i<=n+m;i++) printf("%d ",int(c[i].real()+0.5)); // 避免精度误差
    printf("\n");
    return 0;
}