前言

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

多项式

我们称 \(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\)。

那么在复平面上,所有的单位根都在单位圆上。然后我们考虑找到这些单位根,很容易发现,\(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)\),具体地有:

然后易得 \(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})\)。

此时,\(F(\omega^{m+\frac{n}{2}}_n)=L((\omega^{m+\frac{n}{2}}_n)^2)+\omega^{m+\frac{n}{2}}_nR((\omega^{m+\frac{n}{2}}_n)^2)\)。

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

好这个东西有什么用呢?我们考虑我们的目标是求出 \(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 70 2 4 6|1 3 5 70 4|2 6|1 5|3 70|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)<<' ';}