FFT/NTT總結+洛谷P3803 【模板】多項式乘法（FFT）（FFT/NTT）

前言

衆所周知，這兩個東西都是用來算多項式乘法的。

對於這種常人思惟難以理解的東西，就少些理解，多背板子吧！

所以只總結一下思路和代碼，什麼概念和推式子就靠巨佬們吧

推薦自爲風月馬前卒巨佬的概念和定理都很是到位的總結

推薦ppl巨佬的簡明易懂的總結

FFT

多項式乘法的蹊徑——點值表示法

通常咱們把兩個長度爲\(n\)的多項式乘起來，就相似於作豎式乘法，一位一位地乘再加到對應位上，是\(O(n^2)\)的

如何優化？直接看是沒有思路的，只好另闢蹊徑了。

多項式除了咱們經常使用的係數表示法\(y=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\)之外，還能夠用點值表示法。

所謂點值表示法，就是至關於用函數圖像上\(n\)個點的座標\((x_i,y_i)\)表示一個\(n\)次多項式，顯然這個表示是惟一的。

咱們能夠把係數表示轉化成點值表示。點值表示下的多項式怎麼相乘呢？就是選相同的\(x_i\)，把對應的\(y_i\)相乘。

固然，兩個長度爲\(n\)的多項式相乘獲得的是長度爲\(2n-1\)的多項式，須要\(2n-1\)個點值才能惟一表示，所以一開始兩個多項式也要選\(2n-1\)個點表示。

舉一個例子

\((x+1)(x+2)\)→ → → → \(x^2+3x+2\)
↓\(\qquad\qquad\qquad\qquad\qquad\qquad\)↑
↓（點值）\(\qquad\qquad\qquad\) （係數）↑
↓\(\qquad\qquad\qquad\qquad\qquad\qquad\)↑
\((0,1)(1,2)(2,3)\) （相乘）\(\qquad\quad\) ↑
\((0,2)(1,3)(2,4)\)→ → → →\((0,2)(1,6)(2,12)\)
但是，隨意選\(O(n)\)個點，計算\(y\)是\(O(n)\)的，總時間仍是\(O(n^2)\)的，把它還原成係數表達式更不輕鬆。

因此，若是選的點很普通，這只是一條蹊徑，並非一條捷徑。

神奇的單位根

介紹一個神奇的東西——\(n\)次單位根（記爲\(\omega\)）。

它是\(n\)個複數的集合，每個的\(n\)次方都等於\(1\)。其中的一個是\(e^{{2\pi i}\over n}\)，記爲\(\omega_n\)。

普及一下歐拉公式，\(e^{\theta i}=\cos\theta+(\sin\theta)i\)，\(\theta\)就是這個複數向量的旋轉角。顯然知足\((\omega_n)^n=e^{2\pi i}=1\)，那麼它的\(k\)次方\((\omega_n^k)^n=e^{2k\pi i}\)也等於\(1\)。

因而能夠看出，知足\(n\)次方等於\(1\)的複數的取值只會有\(n\)個，爲\(\omega_n^k(k\in[0,n-1])\)，由於會有\(\omega_n^{n+k}=\omega_n^k\)。

這\(n\)個複數向量在單位圓上呈放射狀。下面是算導的圖片：

介紹消去引理\(ω^{dk}_{dn}=ω^k_n\)，證實很容易的。

DFT&IDFT

單位根有什麼用呢？

看看咱們把\(\omega_n^k(k\in[0,n-1])\)分別帶入多項式求點值會發生什麼就知道了。這個過程叫DFT。

假設\(n\)爲偶數，那麼咱們能夠把它的奇偶項分開，用兩個多項式表示它

\(A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n−2}x^{{\frac n 2}−1}\)

\(A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n−1}x^{{\frac n 2}−1}\)

那麼\(A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)\)。

注意，下面的變化用到了\(ω^{2k}_n=ω^k_{\frac n 2}\)，\(ω^n_n=1\)和\(ω^{\frac n 2}_n=-1\)。

首先帶入單位根

\(A(ω^k_n)=A^{[0]}(ω^{2k}_n)+ω^k_nA^{[1]}(ω^{2k}_n)\\\qquad=A^{[0]}(ω^k_{\frac n 2})+ω^k_nA^{[1]}(ω^k_{\frac n 2})(k<\frac n 2)\)

那\(k\geq\frac n 2\)時又會發生什麼呢？把它變成\(\frac n 2+k\)

\(A(ω^{\frac n 2+k}_n)=A^{[0]}(ω^{n+2k}_n)+ω^{\frac n 2+k}_nA^{[1]}(ω^{n+2k}_n)\\\qquad\quad=A^{[0]}(ω^{2k}_n)+ω^{\frac n 2}_nω^k_nA^{[1]}(ω^{2k}_n)\\\qquad\quad =A^{[0]}(ω^k_{\frac n 2})-ω^k_nA^{[1]}(ω^k_{\frac n 2})\)

也就是說，若是咱們要求一個長度爲\(n\)的多項式取\(\omega_n^k(k\in[0,n-1])\)的\(n\)個點值，咱們只須要求出兩個長度爲\(n/2\)的多項式取\(\omega_{\frac n 2}^k(k\in[0,\frac n 2-1])\)的\(\frac n 2\)個點值，再經過上述兩個式子合併結果。這徹底就是個遞歸過程，很容易寫一個函數來計算。

能夠由算法寫出DFT的複雜度\(T(n)=2T(\frac n 2)+O(n)=O(n\log n)\)

---

固然，別忘了還原成係數表示，這個過程叫作IDFT。

蒟蒻以爲這裏理解太麻煩了，所以再也不證實IDFT的過程，想了解的參見其它的總結。

只須要記住，IDFT的\(\omega_n\)是\(e^{-\frac{2\pi i}n}\)，最後的結果除以\(n\)，其它過程同DFT，能夠寫在一個函數裏。具體看下面的代碼：

void FFT(R complex<double>*a,R int n,R int op){//op=1爲DFT，-1爲IDFT
    if(!n)return;//爲了方便，n的意義與上面不同，這裏的n是a0、a1的長度
    complex<double>a0[n],a1[n];
    for(R int k=0;k<n;++k)
        a0[k]=a[k<<1],a1[k]=a[k<<1|1];//奇偶項分離
    FFT(a0,n>>1,op);FFT(a1,n>>1,op);//遞歸處理
    R complex<double>wn(cos(PI/n),sin(PI/n)*op),w(1,0);//單位根
    for(R int k=0;k<n;++k,w*=wn)//k從到n/2
        a[k]=a0[k]+w*a1[k],a[k+n]=a0[k]-w*a1[k];
}

遞歸版過程

引入例題：洛谷P3803 【模板】多項式乘法（FFT）

因爲係數很小，咱們沒必要擔憂精度的問題了（固然float是使不得的）

遞歸版代碼：

#include<cstdio>
#include<cmath>
#include<complex>
#include<iostream>
#define R register
#define G c=getchar()
using namespace std;
const int N=4.2e6;
const double PI=acos(-1);//自定義π
complex<double>f[N],g[N];
inline int in(){
    R char G;
    while(c<'-')G;
    return c&15;
}
void FFT(R complex<double>*a,R int n,R int op){//同上
    if(!n)return;
    complex<double>a0[n],a1[n];
    for(R int i=0;i<n;++i)
        a0[i]=a[i<<1],a1[i]=a[i<<1|1];
    FFT(a0,n>>1,op);FFT(a1,n>>1,op);
    R complex<double>W(cos(PI/n),sin(PI/n)*op),w(1,0);
    for(R int i=0;i<n;++i,w*=W)
        a[i]=a0[i]+w*a1[i],a[i+n]=a0[i]-w*a1[i];
}
int main(){
    R int n,m,i;
    scanf("%d%d",&n,&m);
    for(i=0;i<=n;++i)f[i]=in();
    for(i=0;i<=m;++i)g[i]=in();
    for(m+=n,n=1;n<=m;n<<=1);//把長度補到2的冪，沒必要擔憂高次項的係數，由於默認爲0
    FFT(f,n>>1,1);FFT(g,n>>1,1);//DFT
    for(i=0;i<n;++i)f[i]*=g[i];//點值直接乘
    FFT(f,n>>1,-1);//IDFT
    for(i=0;i<=m;++i)printf("%.0f ",fabs(f[i].real())/n);//注意結果除以n，當心「-0」
    puts("");
    return 0;
}

好記又好寫的遞歸版結果怎樣呢？

Fast Fast TLE!只有77分。

最主要的緣由在於，空間調用太多了。

蝴蝶操做和Rader排序

針對效率過低的問題，咱們繼續觀察FFT實現過程當中的更多特色。

觀察這一句代碼a[k]=a0[k]+w*a1[k],a[k+n]=a0[k]-w*a1[k];，在操做過程當中，取出了a0[k]和a1[k]的值，經過計算又求出了a[k]和a[k+n]的值。咱們把這樣的一次運算叫作「蝴蝶操做」。

這樣的操做有什麼特色呢？咱們試着將a0和a1合併成一個數組a，每一次蝴蝶操做後，取出了a[k]和a[k+n]的值，又求出了a[k]和a[k+n]的值。最後，整個數組都完成了求值，而並無用到兩個數組！

以\(n=8\)爲例，看看遞歸過程的結構

其實，咱們徹底能夠遞推求解！假設\(a\)數組已經變成了第四層的樣子，先對a0和a四、a2和a六、a1和a五、a3和a7蝴蝶操做，變成第三層；再對a0和a二、a4和a六、a1和a三、a5和a7蝴蝶操做，變成第二層；最後對a0和a一、a2和a三、a4和a五、a6和a7蝴蝶操做，變成第一層，FFT就完成了，而空間複雜度僅爲\(O(n)\)。這個過程能夠用循環來控制。

剩下的問題就是把初始的數組變成最後一層的樣子了。先別急着寫一個遞歸函數暴力把位置換過去。來觀察一下最後序列的編號的二進制表示000,100,010,110,001,101,011,111，是否是與原來000,001,010,011,100,101,110,111相比，每一個位置上的二進制位都反過來了？這樣的變化叫作Rader排序。

咱們能夠\(O(n)\)將Rader排序的映射關係求出。設\(i\)Rader排序後的數爲\(r_i\)，咱們能夠經過\(r_{\frac i 2}\)遞推求出，具體方法能夠看看代碼&動動腦筋qwq

#include<cmath>
#include<cstdio>
#define R register
#define I inline
using namespace std;
const int N=4.2e6;
const double PI=acos(-1);
int n,r[N];
struct C{//手寫complex，比STL快一點點
    double r,i;
    I C(){r=i=0;}
    I C(R double x,R double y){r=x;i=y;}
    I C operator+(R C&x){return C(r+x.r,i+x.i);}
    I C operator-(R C&x){return C(r-x.r,i-x.i);}
    I C operator*(R C&x){return C(r*x.r-i*x.i,r*x.i+i*x.r);}
    I void operator+=(R C&x){r+=x.r;i+=x.i;}
    I void operator*=(R C&x){R double t=r;r=r*x.r-i*x.i;i=t*x.i+i*x.r;}
}f[N],g[N];
I int in(){
    R char c=getchar();
    while(c<'-')c=getchar();
    return c&15;
}
I void FFT(R C*a,R int op){
    R C W,w,t,*a0,*a1;
    R int i,j,k;
    for(i=0;i<n;++i)//根據映射關係交換元素
        if(i<r[i])t=a[i],a[i]=a[r[i]],a[r[i]]=t;
    for(i=1;i<n;i<<=1)//控制層數
        for(W=C(cos(PI/i),sin(PI/i)*op),j=0;j<n;j+=i<<1)//控制一層中的子問題個數
            for(w=C(1,0),a1=i+(a0=a+j),k=0;k<i;++k,++a0,++a1,w*=W)
                t=*a1*w,*a1=*a0-t,*a0+=t;//蝴蝶操做
}
int main(){
    R int m,i,l=0;
    scanf("%d%d",&n,&m);
    for(i=0;i<=n;++i)f[i].r=in();
    for(i=0;i<=m;++i)g[i].r=in();
    for(m+=n,n=1;n<=m;n<<=1,++l);
    for(i=0;i<n;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));//遞推求r
    FFT(f,1);FFT(g,1);
    for(i=0;i<n;++i)f[i]*=g[i];
    FFT(f,-1);
    for(i=0;i<=m;++i)printf("%.0f ",fabs(f[i].r)/n);
    puts("");
    return 0;
}

NTT

原根

FFT的全過程都是基於小數運算，一旦數據較大就會丟失精度。咱們如何在整數範圍內尋找和單位根\(\omega_n^n=1\)性質同樣的數呢？

顯然，若是不在模意義下這樣的式子是不成立的。

引入原根，設爲\(g\)，若是\(g^n=1(\mod p)\)中\(n\)的最小正整數解爲歐拉函數\(\phi(p)\)，則稱\(g\)爲\(p\)的一個原根。

爲了理解NTT，瞭解定義就夠了，固然更多內容能夠參考YL的總結

用\(\omega_n\)替換\(g\)，根據定義顯然也知足消去引理\(ω^{dk}_{dn}=ω^k_n\)。那麼放到多項式乘法裏沒問題。

一個問題在於，遞推過程當中如何根據求出不一樣子問題的\(ω\)？顯然由於\(g^{p-1}=\omega_n^n=1\)，有\(\omega_n=g^{\frac{p-1}n}\)

固然，多項式的長度都被補成了\(2\)的冪，這就要求\(\phi(p)\)的質因子中含有大量的\(2\)；另外，IDFT中的單位根要對原根求逆，最後除以\(n\)也須要求逆，爲了方便還要求\(p\)是質數。這樣，能被用做NTT的模數的數很是少，常見的就是\(998244353(998244352=2^{23}×7×17)\)，\(3\)是它的一個原根。

上面那題的NTT代碼

#include<cstdio>
#define I inline
#define RG register
#define RL RG LL
#define R RG int
typedef long long LL;
const int N=4.2e6,YL=998244353;//一塊兒來%YL
int n,f[N],g[N],r[N];
I int in(){
    RG char c=getchar();
    while(c<'-')c=getchar();
    return c&15;
}
I int qpow(RL b,R k){//快速冪
    RL a=1;
    for(;k;k>>=1,b=b*b%YL)
        if(k&1)a=a*b%YL;
    return a;
}
I void NTT(R*a,R op){
    R i,j,k,d,wn,w,t,*a0,*a1;
    for(i=0;i<n;++i)
        if(i<r[i])t=a[i],a[i]=a[r[i]],a[r[i]]=t;
    for(i=1;i<n;i<<=1){
        wn=qpow(3,(YL-1)/(d=i<<1));//原根變換
        if(op&2)wn=qpow(wn,YL-2);//注意要求逆
        for(j=0;j<n;j+=d)
            for(w=1,a1=(a0=a+j)+i,k=0;k<i;++k,++a0,++a1,w=(LL)w*wn%YL)
                t=(LL)*a1*w%YL,*a1=(*a0-t+YL)%YL,*a0=(*a0+t)%YL;
    }
}
int main(){
    R m,i,l=0;
    scanf("%d%d",&n,&m);
    for(i=0;i<=n;++i)f[i]=in();
    for(i=0;i<=m;++i)g[i]=in();
    for(m+=n,n=1;n<=m;n<<=1,++l);
    for(i=0;i<n;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    NTT(f,1);NTT(g,1);
    for(i=0;i<n;++i)f[i]=(LL)f[i]*g[i]%YL;
    NTT(f,-1);
    R inv=qpow(n,YL-2);//一樣注意要求逆
    for(i=0;i<=m;++i)printf("%lld ",(LL)f[i]*inv%YL);
    puts("");
    return 0;
}