學習筆記——多項式

前言

該來的還是來了,終于來啃這個東西了。以下是我對于多項式比較粗淺的理解。

多項式

我們稱 \(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)<<' ';}