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

一、前置知识

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

1.1   多项式

设 $A(x)$ 为一个 $n$ 次多项式，则可以表示为

其中，$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})$ 为其系数。按照每项次数进行奇偶分类，有

分别定义

不难有

分别求其在 $\omega_{n}^{k}(0\le k <\frac n2)$ 上的点值，有

与

发现上两式的形式类似，而我们将问题从 $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$ 与 $a_k$ 之间的关系，则可反解出 $a_k$。对上面式子进行求解：

令 $S(x)=\sum\limits_{i=0}^{n-1}x^i$，有

两式相减得到

当 $\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}$ 处才有值，所以

即

这样，我们再次利用 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$ 为多项式 $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;
}