学习笔记——多项式和 FFT

前言#

该来的还是来了，终于来啃这个东西了。以下是我对于多项式比较粗浅的理解。

多项式#

我们称 $F(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x$ 为一个 $n$ 次多项式。

多项式系数和系数表示法#

多项式的系数就是数列 $a_1,a_2,a_3,\cdots,a_n$，我们令 $F[i]$ 表示多项式 $F(x)$ 的 $i$ 次项的系数。

多项式的系数表示法，就是用多项式的系数来唯一表示一个多项式。

多项式单点求值和点值表示法#

其实多项式的形式非常像一个 $n$ 次函数对吧，那我们对某个值 $v$，可以求出多项式的具体的值。比如说有 $F(x)=2x^2+x-1$，那么对于 $v=2$，$F(2)=8+2-1=9$。

那点值表示法就非常简单了，就是考虑对于一个 $n$ 次多项式，当我们看成一个函数的时候，我们只要确定了 $n+1$ 个点对 $(x_i,y_i)$ 满足 $F(x_i)=y_i$，我们就可以确定这个多项式。具体可以考虑高斯消元的时候需要的方程组个数，容易证明。这样的话，我们用这 $n+1$ 个点对可以唯一表示一个多项式。

多项式乘法（卷积）#

多项式的乘法实际上是一种卷积。说白了也很简单，就是拆括号。

具体的，如果 $H=F*G$，那么 $H[k]=\sum\limits_{i=0}^kF[i]G[k-i]$。看起来非常高级，实际上就是小学就学的乘法分配律拆括号。

卷积的一般形式是 $H[k]=\sum\limits_{i\oplus j=k}F[i]G[j]$，其中 $\oplus$ 表示某种运算。那么多项式的乘法就是加法卷积。

FFT#

FFT，全名快速傅里叶变换，用于加速多项式乘法。

很容易发现，如果暴力卷积，求多项式乘法的复杂度显然是 $O(n^2)$ 的。这样的复杂度不让人满意。于是我们希望能够加速这个过程。FFT 能够把这个复杂度优化到 $O(n\log n)$。

那这是怎么做到的呢？

DFT&IDFT#

首先我们已经了解了多项式的两种表示法，现在我们来应用它们。

我们把一个多项式的系数表示法转化成点值表示法的过程，叫做 DFT，也叫多项式求值。

那么反过来，把点值表示法转化成系数表示法的过程叫做 IDFT，也叫多项式插值。

DFT 比较简单，就随便带几个值进去求值就行了，IDFT 就比较麻烦，似乎还需要高斯消元来解方程组，好像非常的困难。

而难以置信的是，FFT 的过程，是进行 DFT 然后再来一次 IDFT！

单位复根#

上面已经说了 FFT 的流程，其中需要进行 IDFT 仿佛比暴力卷积还慢。所以我们肯定不是高消暴力插值。也就是说，我们需要选取一些特殊的点来求值。

当然如果选一些看起来人畜无害的整数点，还是没办法解决最终 IDFT 的困难之处。于是我们选择单位复根，也就是 FFT 的精髓之一了。

复数#

大家都知道，$z=a+bi$，其中 $i^2=-1$，是复数的一般形式。大家都知道复平面上复数的加减乘完全切合于向量的加减乘。具体地，复数 $z=a+bi$ 用复平面上坐标 $(a,b)$ 表示。那类似与向量，我们把复数到原点的距离叫做模长（记为 $|z|$），把复数表示的坐标到原点的连线与 $x$ 轴形成的逆时针夹角叫做幅角。

尤其重要的，我们考虑两个复数相乘时，是模长相乘，幅角相加。证明略。

单位根#

$n$ 次单位根就是 $n$ 次幂为 $1$ 的复数。就是 $x^n=1$ 的复根。

我们把单位根当成这个很特殊的值来求多项式的值，为什么用单位根捏？它有一些很好的性质。

首先单位根的模长一定是 $1$。

证明

如果单位根 $z$，有 $|z|>1$，则 $|z^n|=|z|^n>1$，反之如果有 $|z|<1$，则 $|z^n|=|z|^n<1$，所以 $|z|=1$。

那么在复平面上，所有的单位根都在单位圆上。然后我们考虑找到这些单位根，很容易发现，$n$ 次单位根就是幅角为 $0,\dfrac{1}{n}\pi,\dfrac{2}{n}\pi,\cdots,\dfrac{n-1}{n}\pi$，然后由于 $n$ 次单位根是只有 $n$ 个的（算术基本定理），所以这就是所有的单位根。

接下来为了方便表述，我们给单位根一个符号就是 $\omega^m_n$，表示 $n$ 次单位根中的第 $m$ 个（从 $0$ 开始）。然后虽然只有 $n$ 个 $n$ 次单位根，后面可能会使 $m$ 大于 $n$ 或者小于 $0$，请自动 $\omega^m_n=\omega^{m\%n}_n$（这就是为什么下标从 $0$ 开始了）。

然后我们来搞点性质出来。

1. $\forall n,\omega^0_n=1$，显然。
2. $\omega^k_n=(\omega^1_n)^k$，显然。
3. $\omega^a_n\times\omega^b_n=\omega^{a+b}_n$，显然，是上一条的一般情况。
4. $\omega^{2m}_{2n}=\omega^m_n$，很好理解，考虑单位根是把单位圆等分为 $n$ 份，然后第 $i$ 个单位根就是从 $x$ 轴开始逆时针取 $i$ 份后的那个根。这样子分成 $2n$ 份取 $2m$ 份和分成 $n$ 份取 $m$ 份是一样的。同理可以拓展得到 $\omega^{km}_{kn}=\omega^m_n$。
5. 若 $n$ 是偶数，则 $\omega^{m+\frac{n}{2}}_n=-\omega^m_n$。证明：$\omega^{m+\frac{n}{2}}_n=\omega^m_n\times\omega^\frac{n}{2}_n$，然后考虑到 $\omega^\frac{n}{2}_n=-1$（很好想的，就是与负半轴重合的根），证毕。

分治加速 DFT#

接下来是用单位复根加速 DFT 的过程。

现在我们手上有一个多项式 $F(x)$。我们钦定它的项数是 $2^n$ 项。

然后我们把这个多项式的系数按奇偶分成两个 $\dfrac{n}{2}$ 项的多项式 $L(x)$ 和 $R(x)$，具体地有：

$$L(x)=F[0]+F[2]x+F[4]x^2+\cdots+F[n-2]x^{\frac{n}{2}-1}$$

$$R(x)=F[1]+F[3]x+F[5]x^2+\cdots+F[n-1]x^{\frac{n}{2}-1}$$

然后易得 $F(x)=L(x^2)+xR(x^2)$。

接下来掏出单位复根。代入 $\omega^m_n(m<\dfrac{n}{2})$。

此时，$F(\omega^m_n)=L((\omega^m_n)^2)+\omega^m_nR((\omega^m_n)^2)$。

又 $(\omega^m_n)^2=\omega^{2m}_n=\omega^m_\frac{n}{2}$。

所以 $F(\omega^m_n)=L(\omega^m_\frac{n}{2})+\omega^m_nR(\omega^m_\frac{n}{2})$。

然后我们再代一个 $\omega^{m+\frac{n}{2}}_n(m<\dfrac{n}{2})$。

此时：$$\begin{aligned}
F(\omega^{{m+\frac{n}{2}}_n)&=L((\omega}{2}}_n)^2)+\omega{2}}_nR((\omega^{{m+\frac{n}{2}}_n)}2)\&=L(\omega^{{2m+n}_n)+\omega}{2}}_nR(\omega^{{2m+n}_n)\&=L(\omega}_n)+\omega^{{m+\frac{n}{2}}_nR(\omega}n)\&=L(\omega^{m_\frac{n}{2})+\omega}{2}}nR(\omega^{m_\frac{n}{2})\&=L(\omega}m\frac{n}{2})-\omega^m_nR(\omegam\frac{n}{2})
\end{aligned}$$

发现这个式子和上面那个就差了后面一个负号。

好这个东西有什么用呢？我们考虑我们的目标是求出 $F(x)$ 在所有单位根处的取值，那假如说我们已经知道了 $L(x)$ 和 $R(x)$ 的点值表示（用单位根），套用上面两个式子，我们可以在 $O(n)$ 内求出 $F(x)$ 的点值表示。

为了求 $L(x)$ 和 $R(x)$，递归下去就行了。

IDFT 和 FFT#

关于 IDFT，其求法就是把 DFT 中的 $\omega^1_n$ 换成 $\omega^{-1}_n$，然后最后每项除掉一个 $n$ 就可以了。

证明略，具体可以参考 command_block 的博客 中有详细的证明及其本质，这里就不展开了（反正记不住）。

那我们的 FFT 就是对于两个多项式，分别进行 DFT，然后把点值相乘，再 IDFT 就做完了。

对了，忘记提一嘴，我们为了分治，把整个补成了长度为 $2$ 的整次幂。这个需要预处理以下。

然后你可以写出暴力代码了，接下来讲一些优化。

蝴蝶变换优化#

我们考虑从奇偶分开的时候，暴力做需要很复杂的拷贝以及递归，很不爽，于是我们考虑最后分治是什么样子的。

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 4 2 6 1 5 3 7 的顺序进行处理，那么每次从中间分开就可以用循环处理这个东西了。那怎么求这个东西呢？其实容易发现，每个位置所对应的就是其二进制位翻转的结果，比如说 $1\to 4$ 就是 $001\to 100$。

这个东西怎么求快呢？类似于数位 dp，我们令 $rev[i]$ 表示 $i$ 翻转的结果。然后 $rev[i]=rev[i>>1]>>1$，就是把 $i$ 二进制除了最后一位翻转的结果作为它的最后几位，然后判断这个数的奇偶性在最高位加 $1$ 就行了。

Code#

代码实现的话，首先要实现一个复数。（不要用自带的 STL，太慢了）

struct CP{
    double a,b;
    CP(double _a=0,double _b=0){a=_a;b=_b;}
    CP friend operator+(CP x,CP y){return CP(x.a+y.a,x.b+y.b);}
    CP friend operator-(CP x,CP y){return CP(x.a-y.a,x.b-y.b);}
    CP friend operator*(CP x,CP y){return CP(x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a);}
};

然后剩下的就按照上面说的做就行了。

• P3803 【模板】多项式乘法（FFT）

#include<bits/stdc++.h>
#define ll long long
#define inf (1<<30)
#define INF (1ll<<60)
#define pii pair<int,int>
#define pll pair<ll,ll>
#define mkp make_pair
#define fi first
#define se second
#define rep(i,j,k) for(int i=(j);i<=(k);i++)
#define per(i,j,k) for(int i=(j);i>=(k);i--)
using namespace std;
const int MAXN=(1<<22)+10;
struct CP{
    double a,b;
    CP(double _a=0,double _b=0){a=_a;b=_b;}
    CP friend operator+(CP x,CP y){return CP(x.a+y.a,x.b+y.b);}
    CP friend operator-(CP x,CP y){return CP(x.a-y.a,x.b-y.b);}
    CP friend operator*(CP x,CP y){return CP(x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a);}
}F[MAXN],G[MAXN];
int rev[MAXN];
void makerev(int n){
    rep(i,1,n-1){
        rev[i]=rev[i>>1]>>1;
        if(i&1) rev[i]|=(n>>1);
    }
}
void FFT(CP f[],int n,int fl){
    makerev(n);
    rep(i,0,n-1) if(i<rev[i]) swap(f[i],f[rev[i]]);
    for(int h=2;h<=n;h<<=1){
        CP w=CP(cos(M_PI*2/h),sin(M_PI*2/h)*fl);
        for(int i=0;i<n;i+=h){
            CP nw=CP(1,0);
            rep(j,i,i+h/2-1){
                CP L=f[j],R=f[j+h/2]*nw;
                f[j]=L+R;f[j+h/2]=L-R;
                nw=nw*w;
            }
        }
    }if(fl==-1) rep(i,0,n-1) f[i].a/=n;
}
int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    int n,m,x;cin>>n>>m;
    rep(i,0,n) cin>>x,F[i]=CP(x,0);
    rep(i,0,m) cin>>x,G[i]=CP(x,0);
    int len=1;while(len<n+m) len<<=1;
    len<<=1;FFT(F,len,1);FFT(G,len,1);
    rep(i,0,len-1) F[i]=F[i]*G[i];
    FFT(F,len,-1);
    rep(i,0,n+m) cout<<int(F[i].a+0.5)<<' ';
}