很久以前就聽說了這東西，一直沒空學。最近重學多項式，就重新搞了一下。

MTT 主要解決的是任意模數（或者說是沒有模數）的多項式乘法，可以用于應對專門惡心人的毒瘤題。

首先，假設多項式次數 $10^5$,值域 $10^9$，那么樸素 FFT 光答案就高達 $10^{23}$，就算你不損失精度也轉不成 long long，況且中間值比這還要高幾個量級。

MTT 采用拆系數的方法，即對于多項式 $A(x)$，拆成兩個多項式 $M{A}_0(x)+A_1(x)$ 來減小精度壓力。其中 $M$ 為一較大整數，一般取 $\sqrt{MOD}$。要注意模數 $2e9$ 級別的時候不要偷懶用位運算，$2^{15}$ 或 $2^{16}$ 都容易炸，并且復雜度瓶頸不在這，不用卡這點常。

然后你只需要求 $(M{A}_0(x)+A_1(x))(M{B}_0(x)+B_1(x))=M^2{A}_0(x)B_0(x)+M(A_0(x)B_1(x)+A_1(x)B_0(x))+A_1(x)B_1(x)$ 就可以了，這樣答案是 $10^{14}$ 級別，開 long double 完全不虛。

這樣一共需要 $7$ 次 DFT，一般問題不大。接下來是利用高超的玄學技巧壓到 $4$ 次 DFT 的方法。

前置知識：普通 FFT 三次變兩次

對于一般的實系數多項式乘法 $A(x)?B(x)$，有一種優化方法是構造一個多項式 $A(x)+i?B(x)$，因為是實系數，所以直接把 $B(x)$ 的系數掛到虛部上就可以了。

對它做一次 DFT 然后自己乘自己，得到的是 $A^2(x)?B^2(x)+2i?A(x)B(x)$

這樣算出來的虛部除以 $2$ 就是答案。

該算法將大量用到該思想。

DFT

考慮我們對兩個實系數多項式 $A(x),B(x)$ 做 DFT，因為是實系數，我們把虛部在那里空著顯得很浪費，嘗試把這個優化成一次 DFT。

構造多項式 $A(x)+iB(x)$，花費一次 DFT 得到 $\mathrm{DFT}(A(x)+i(B(x)))$。因為 DFT 的本質是點值，所以這個也等于 $\mathrm{DFT}(A(x))+i\mathrm{DFT}(B(x))$。

這個顯然不是實部和虛部的關系了。現在我們需要不用 DFT 把這兩個分離開。

考慮共軛虛多項式 $A(x)?iB(x)$，如果它的 DFT 可以快速得到，我們就可以簡單地求出 $A(x)$ 和 $B(x)$ 的 DFT。

然后又是愉快的推式子時間。

設兩個多項式的系數為 $a_i,b_i$，代入一個單位根 ${\omega }_n^k$

$$\mathrm{DFT}(A(x)+iB(x))=A({\omega }_n^k)+iB({\omega }_n^k)$$

$$=\sum _{j=0}^{n?1}{a}_j{\omega }_n^{kj}+i{b}_j{\omega }_n^{kj}$$

推不動了，考慮硬算。令 ${\omega }_n^{kj}=x_j+i{y}_j$

$$\mathrm{DFT}(A(x)+iB(x))=\sum _{j=0}^{n?1}{a}_j(x_j+i{y}_j)+i{b}_j(x_j+i{y}_j)$$

$$=\sum _{j=0}^{n?1}(a_j{x}_j?b_j{y}_j)+i(a_j{y}_j+b_j{x}_j)$$

然后這不是三角函數公式嗎？

所以 DFT 的某一點值的本質就是兩個多項式對應項組成的向量 $(a_j,b_j)$ 旋轉同一角度后的向量和。

對于另一邊 $\mathrm{DFT}(A(x)?iB(x))$，其對應的向量是 $(a_j,?b_j)$ ，也就是關于 $x$ 軸對稱。那么我們把上面那個角度倒著旋轉，得到的向量和也是對稱的。

翻譯成人話， DFT 后翻轉再共軛就得到了共軛虛多項式的 DFT。

注意是長度為 $n$ 的翻轉，所以 $0$ 的位置需要特判。

這樣，我們用兩次 DFT 得到了 $A_0(x),A_1(x),B_0(x),B_1(x)$ 的 DFT。

然后我們就可以求出 $A_0(x)B_0(x),A_0(x)B_1(x),A_1(x)B_0(x),A_1(x)B_1(x)$ 的點值了。

IDFT

還是用一樣的方法抱團變換。

考慮構造多項式 $P(x)=A_0(x)B_0(x)+i{A}_0(x)B_1(x)$，注意是 IDFT 之后的，也就是這里設的是我們要求的答案。

因為我們知道 $A_0(x)B_0(x),A_0(x)B_1(x)$ 在所有單位根上的點值，所以可以直接求出 $P(x)$ 的點值，然后做 IDFT，實部虛部就是 $A_0(x)B_0(x)$ 和 $A_0(x)B_1(x)$。

另外兩個同理，需要 $2$ 次 IDFT。

總共做了 $4$ 次 DFT，比三模數快了 $8$ 倍。

注意輸出的時候答案已經很大了，要先模了再乘 $M$。

#include <iostream> #include <cstring> #include <cctype> #include <cstdio> #include <cmath> #define MAXN 262144+5 #define double long double using namespace std; inline int read() {int ans=0;char c=getchar();while (!isdigit(c)) c=getchar();while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();return ans; } typedef long long ll; const double Pi=acos(-1.0); int p; struct complex {double x,y;inline complex(const double& x=0,const double& y=0):x(x),y(y){} }; inline complex operator +(const complex& a,const complex& b){return complex(a.x+b.x,a.y+b.y);} inline complex operator -(const complex& a,const complex& b){return complex(a.x-b.x,a.y-b.y);} inline complex operator *(const complex& a,const complex& b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);} inline complex rev(const complex& a){return complex(a.x,-a.y);} complex rt[2][20]; int l,lim,r[MAXN]; inline void init(){lim=1<<l;for (int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));} inline void fft(complex* a,int type) {for (int i=0;i<lim;i++) if (i<r[i]) swap(a[i],a[r[i]]);for (int L=0;L<l;L++){int mid=1<<L,len=mid<<1;complex Wn=rt[type][L+1];for (int s=0;s<lim;s+=len){complex w(1,0);for (int k=0;k<mid;k++,w=w*Wn){complex x=a[s+k],y=w*a[s+mid+k];a[s+k]=x+y,a[s+mid+k]=x-y;}}}if (type) for (int i=0;i<lim;i++) a[i].x/=lim,a[i].y/=lim; } complex A[MAXN],B[MAXN],C[MAXN],D[MAXN]; int main() {freopen("test.in","r",stdin);freopen("test.out","w",stdout);for (int i=0;i<20;i++){double a=2*Pi/(1<<i);rt[0][i]=rev(rt[1][i]=complex(cos(a),sin(a)));}int n,m;n=read(),m=read(),p=read();for (l=0;(1<<l)<=n+m;l++);init();for (int i=0;i<=n;i++){int x=read();A[i]=complex(x>>15,x&((1<<15)-1));}for (int i=0;i<=m;i++){int x=read();C[i]=complex(x>>15,x&((1<<15)-1));}fft(A,0),fft(C,0);B[0]=rev(A[0]),D[0]=rev(C[0]);for (int i=1;i<lim;i++) B[i]=rev(A[lim-i]),D[i]=rev(C[lim-i]);for (int i=0;i<lim;i++) {complex a=A[i],b=B[i],c=C[i],d=D[i];A[i]=(a+b)*complex(0.5,0),B[i]=(a-b)*complex(0,-0.5);C[i]=(c+d)*complex(0.5,0),D[i]=(c-d)*complex(0,-0.5);}for (int i=0;i<lim;i++){complex a=A[i],b=B[i],c=C[i],d=D[i];A[i]=a*c+a*d*complex(0,1);B[i]=b*c+b*d*complex(0,1);}fft(A,1),fft(B,1);for (int i=0;i<=n+m;i++){ll ans=0;ll a=A[i].x+0.5,b=A[i].y+B[i].x+0.5,c=B[i].y+0.5;ans=(ans+((a%p)<<30))%p;ans=(ans+((b%p)<<15))%p;ans=(ans+c)%p;printf("%lld ",ans);}return 0; }

總結

以上是生活随笔為你收集整理的MTT 学习笔记的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。