学习笔记——多项式和 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)<<' ';
}