多項式
系數表示法
設\(f(x)\)為一個\(n-1\)次多項式，則 \(f(x)=\sum\limits_{i=0}^{n-1}a_i*x^i\)
其中\(a_i\)為\(f(x)\)的系數，用這種方法計算兩個多項式相乘（逐位相乘）復雜度為\(O(n^2)\)
點值表示法
根據小學知識，一個\(n-1\)次多項式可以唯一地被\(n\)個點確定
即，如果我們知道了對於一個多項式的\(n\)個點\((x_1,y_1),(x_2,y_2)……(x_n,y_n)\)
那么這個多項式唯一滿足，對任意\(1\le i \le n\)，滿足\(y_i=\sum\limits_{j=0}^{n-1}a_j*x_i^j\)
那么用點值實現多項式相乘是什么復雜度呢？
首先我們需要選\(n\)個點，每個點需要求出其在多項式中的值，復雜度為\(O(n^2)\)
然后把兩個點值表示的多項式相乘，由於\(c(x_i)=a(x_i)*b(x_i)\)，復雜度為\(O(n)\)
最后插值法用點值求出系數，復雜度為\(O(n^2)\)（我還不會插值）
考慮如果可以快速實現點值轉系數和系數轉點值，豈不是可以快速計算多項式乘法（說的簡單，你倒是告訴我怎么快速轉化啊）
前置芝士
復數
定義虛數單位\(i=\sqrt{-1}\)，\(a,b\)為實數，則形如\(a+bi\)的數叫復數
其中\(a\)為復數的實部，\(b\)為復數的虛部
在復平面中，復數可以被表示為向量，所以和向量具有很多相似的性質，我們也可以用向量來理解復數，但是復數具有更多性質，比如作為一個數代入多項式
其中模長定義為\(\sqrt{a^2+b^2}\),幅度定義為\(x\)軸正半軸到向量轉角的有向角
復數運算法則：
加減法與向量相同，重點是乘法：
幾何定義為：模長相乘，幅角相加
代數定義為：
單位根
我們首先定義圓心為坐標原點，半徑為\(1\)的圓叫做單位圓
我們將這個圓\(n\)等分，得到\(n\)個圓上的點，以這\(n\)個圓上點的橫坐標作為實部，縱坐標作為虛部，就得到了\(n\)個復數
網上扒的圖片，侵刪\(QwQ\)
首先我們不自找麻煩，以\((1,0)\)作為這\(n\)個點的起點，記作\({\omega _n}^0\),逆時針方向第\(k\)個點記作\({\omega _n}^k\)
根據模長相乘幅角相加，我們可以看出\(\omega _n^k\)是\(\omega _n^0\)的\(k\)次方，其中\(\omega _n^1\)被稱為\(n\)次單位根
根據幅角，我們可以計算出\(\omega _n^k\)表示的復數為\(cos{\frac{k}{n}2\pi}+i*sin{\frac{k}{n}2\pi}\)
單位根具有一些性質：
\(1.\omega _n^k=cos{\frac{k}{n}2\pi}+i*sin{\frac{k}{n}2\pi}=e^{\frac{2\pi ki}{n}}\)
證明：這個第一步到第二步由定義得出，第二步到第三步由歐拉公式得出
\(2.\omega _{2n}^{2k}=\omega_{n}^{k}\)
證明：\(\omega _{2n}^{2k}=e^{\frac{2\pi 2ki}{2n}}=e^{\frac{2\pi ki}{n}}=\omega_{n}^{k}\)
\(3.\omega _{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}\)
證明：\(\omega _{n}^{k+\frac{n}{2}}=\omega _{2n}^{2k+n}=\omega _{2n}^{2k}*\omega _{2n}^{n}=\omega_{n}^{k}*(cos\pi+i*sin\pi)=-\omega_{n}^{k}\)
\(4.\omega _{n}^{0}=\omega _{n}^{n}=1\)
證明：不用證了吧……
正文之前
這段話有可能有助於您理解本算法：
傅立葉這個大神仙根本就沒見過計算機長什么樣，所以他提出的傅立葉變換和逆變換只是一種將系數轉點值和將點值轉系數的方法，沒有任何降低復雜度的功效，至於快速傅立葉變換是后人再研究傅立葉變換發現的一種加速方法，是對\(DFT\)和\(IDFT\)的優化
離散傅立葉變換（DFT）
假設\(f(x)=\sum\limits_{i=0}^{n-1}a_i*x_i\)
\(DFT(a)=(f(1),f(\omega _{n}),f(\omega _{n}^{2}),……,f(\omega _{n}^{n-1}))\)
通俗點說，就是對於一個系數表示法的多項式\(f(x)\)，將\((1,\omega _{n},\omega _{n}^{2},……,\omega _{n}^{n-1})\)帶入求出該多項式的點值表示法
離散傅立葉逆變換（IDFT）
將\(f(x)\)在\(n\)個\(n\)次單位根處的點值表示轉化為系數表示
這里就可以回答，為什么我們要讓\(n\)次單位根作為\(x\)代入多項式
假設\((y_0,y_1,y_2,……,y_{n-1})\)是多項式\(A(x)=\sum\limits_{i=0}^{n-1}a_i*x_i\)的離散傅立葉變換
我們另有一個多項式\(B(x)=\sum\limits_{i=0}^{n-1}=y_i*x_i\)
將上述\(n\)次單位根的倒數\((1,\omega _{n}^{-1},\omega _{n}^{-2},……,\omega _{n}^{-(n-1)})\)代入\(B(x)\)得到新的離散傅立葉變換\((z_0,z_1,z_2,……,z_{n-1})\)
則我們發現
對於\(\sum\limits_{i=0}^{n-1}(\omega _n^{j-k})^i\)我們單獨考慮：
當\(j-k=0\)時
答案為\(n\)
當\(j\ne k\)時
等比數列求和得到\(\frac{(\omega _n^{j-k})^n-1}{\omega _n^{j-k}-1}=\frac{(\omega _n^n)^{j-k}-1}{\omega _n^{j-k}-1}=\frac{1^{j-k}-1}{\omega _n^{j-k}-1}=0\)
所以
即
得出結論：對於以\(A(x)\)的離散傅立葉變換作為系數的多項式\(B(x)\)，取單位根的倒數\((1,\omega _{n}^{-1},\omega _{n}^{-2},……,\omega _{n}^{-(n-1)})\)作為\(x\)代入，再將結果除以\(n\)即為\(A(x)\)的系數
這個結論實現了將多項式點值轉化為系數
快速傅立葉變換
我們一頓分析最后發現復雜度……仍然是\(O(n^2)\)
我們學這破玩意的意義不就是降低復雜度嘛，所以我們接下來講怎么降復雜度
我們先設\(A(x)=\sum\limits_{i=0}^{n-1}a_i*x_i\)
我們將\(A(x)\)的下標按奇偶分類，得到
設兩個多項式
那么就可以發現\(A(x)=A_1(x^2)+xA_2(x^2)\)
將\(x=\omega _n^{k}(k<\frac{n}{2})\)代入
將\(x=\omega _n^{k+\frac{n}{2}}(k<\frac{n}{2})\)代入
我們一點也不驚喜地發現，只要求出\(A_1(x)\)和\(A_2(x)\)在\((\omega _{\frac{n}{2}}^0,\omega _{\frac{n}{2}}^1,……,\omega _{\frac{n}{2}}^{\frac{n}{2}-1})\)的點值表示，就可以\(O(n)\)地求出\(A(x)\)在\((1,\omega _{n},\omega _{n}^{2},……,\omega _{n}^{n-1})\)
所以我們可以遞歸實現\(O(nlogn)\)求解多項式乘法了
注意：我們假設\(f*g=h\)，那么對於\(f\)和\(g\)都要直接求出大於\(n+m+1\)個的\(2^k\)個點值（由於分治要求，點數一定是\(2\)的整次冪）
\(code\)
#include<bits/stdc++.h>
using namespace std;
namespace red{
#define mid ((l+r)>>1)
#define eps (1e-8)
inline int read()
{
int x=0;char ch,f=1;
for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
if(ch=='-') f=0,ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return f?x:-x;
}
const int N=5e6+10;
const double pi=acos(-1.0);
int n,m;
struct complex
{
double x,y;
complex(double tx=0,double ty=0){x=tx,y=ty;}
inline complex operator + (const complex t) const
{
return complex(x+t.x,y+t.y);
}
inline complex operator - (const complex t) const
{
return complex(x-t.x,y-t.y);
}
inline complex operator * (const complex t) const
{
return complex(x*t.x-y*t.y,x*t.y+y*t.x);
}
}a[N],b[N];
inline void fft(int limit,complex *a,int inv)
{
if(limit==1) return;
complex a1[limit>>1],a2[limit>>1];
for(int i=0;i<limit;i+=2)
{
a1[i>>1]=a[i],a2[i>>1]=a[i+1];
}
fft(limit>>1,a1,inv);
fft(limit>>1,a2,inv);
complex Wn=complex(cos(2.0*pi/limit),inv*sin(2.0*pi/limit)),w=complex(1,0);
for(int i=0;i<(limit>>1);++i,w=w*Wn)
{
a[i]=a1[i]+w*a2[i];
a[i+(limit>>1)]=a1[i]-w*a2[i];
}
}
inline void main()
{
n=read(),m=read();
for(int i=0;i<=n;++i) a[i].x=read();
for(int i=0;i<=m;++i) b[i].x=read();
int limit=1;
while(limit<=n+m) limit<<=1;
fft(limit,a,1);
fft(limit,b,1);
for(int i=0;i<=limit;++i)
{
a[i]=a[i]*b[i];
}
fft(limit,a,-1);
for(int i=0;i<=n+m;++i) printf("%d ",(int)(a[i].x/limit+0.5));
}
}
signed main()
{
red::main();
return 0;
}
然而我們發現好像有點慢
迭代優化
眾所周知遞歸比較慢，我們有沒有什么方法可以用迭代代替遞歸呢？
扒圖時間到
通過一頓找規律，我們根本不能發現，每個數字在分治后的位置就是它所在位置的二進制翻轉
這個規律也有一個好聽的名字，叫蝴蝶定理
那么我們只要預處理出每個數字在最后一次遞歸中的位置，然后自底向上合並，豈不是可以擺脫遞歸
#include<bits/stdc++.h>
using namespace std;
namespace red{
#define eps (1e-8)
inline int read()
{
int x=0;char ch,f=1;
for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
if(ch=='-') f=0,ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return f?x:-x;
}
const int N=5e6+10;
const double pi=acos(-1.0);
int n,m;
int limit=1,len;
int pos[N];
struct complex
{
double x,y;
complex(double tx=0,double ty=0){x=tx,y=ty;}
inline complex operator + (const complex t) const
{
return complex(x+t.x,y+t.y);
}
inline complex operator - (const complex t) const
{
return complex(x-t.x,y-t.y);
}
inline complex operator * (const complex t) const
{
return complex(x*t.x-y*t.y,x*t.y+y*t.x);
}
}a[N],b[N],buf[N];
inline void fft(complex *a,int inv)
{
for(int i=0;i<limit;++i)
if(i<pos[i]) swap(a[i],a[pos[i]]);
for(int mid=1;mid<limit;mid<<=1)
{
complex Wn(cos(pi/mid),inv*sin(pi/mid));
for(int r=mid<<1,j=0;j<limit;j+=r)
{
complex w(1,0);
for(int k=0;k<mid;++k,w=w*Wn)
{
buf[j+k]=a[j+k]+w*a[j+k+mid];
buf[j+k+mid]=a[j+k]-w*a[j+k+mid];
}
}
for(int i=0;i<limit;++i) a[i]=buf[i];
}
}
inline void main()
{
n=read(),m=read();
for(int i=0;i<=n;++i) a[i].x=read();
for(int i=0;i<=m;++i) b[i].x=read();
while(limit<=n+m) limit<<=1,++len;
for(int i=0;i<limit;++i)
pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
fft(a,1);
fft(b,1);
for(int i=0;i<=limit;++i) a[i]=a[i]*b[i];
fft(a,-1);
for(int i=0;i<=n+m;++i) printf("%d ",(int)(a[i].x/limit+0.5));
}
}
signed main()
{
red::main();
return 0;
}
蝴蝶操作
考慮這里
for(int r=mid<<1,j=0;j<limit;j+=r)
{
complex w(1,0);
for(int k=0;k<mid;++k,w=w*Wn)
{
buf[j+k]=a[j+k]+w*a[j+k+mid];
buf[j+k+mid]=a[j+k]-w*a[j+k+mid];
}
}
for(int i=0;i<limit;++i) a[i]=buf[i];
之所以加\(buf\)數組是因為兩次賦值\(a\)的值會變化，我們可以提前存儲\(a\)數組的值，然后優化掉\(buf\)數組
for(int k=0;k<mid;++k,w=w*Wn)
{
complex x=a[j+k],y=w*a[j+k+mid];
a[j+k]=x+y;
a[j+k+mid]=x-y;
}
三次變兩次優化
觀察到上面的代碼我們跑了三次肥肥兔，現在我們有一種方法可以少跑一次
假設我們求\(f(x)*g(x)\)
設復多項式\(h(x)=f(x)+i*g(x)\)，實部為\(f(x)\)，虛部為\(g(x)\)
那么\(h(x)^2=(f(x)+i*g(x))^2=f(x)^2-g(x)^2+i*2*f(x)*g(x)\)
我們只要把\(h(x)^2\)的虛部除以\(2\)就得到了結果
完全版：
#include<bits/stdc++.h>
using namespace std;
namespace red{
#define eps (1e-8)
inline int read()
{
int x=0;char ch,f=1;
for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
if(ch=='-') f=0,ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return f?x:-x;
}
const int N=5e6+10;
const double pi=acos(-1.0);
int n,m;
int limit=1,len;
int pos[N];
struct complex
{
double x,y;
complex(double tx=0,double ty=0){x=tx,y=ty;}
inline complex operator + (const complex t) const
{
return complex(x+t.x,y+t.y);
}
inline complex operator - (const complex t) const
{
return complex(x-t.x,y-t.y);
}
inline complex operator * (const complex t) const
{
return complex(x*t.x-y*t.y,x*t.y+y*t.x);
}
}a[N];
inline void fft(complex *a,int inv)
{
for(int i=0;i<limit;++i)
if(i<pos[i]) swap(a[i],a[pos[i]]);
for(int mid=1;mid<limit;mid<<=1)
{
complex Wn(cos(pi/mid),inv*sin(pi/mid));
for(int r=mid<<1,j=0;j<limit;j+=r)
{
complex w(1,0);
for(int k=0;k<mid;++k,w=w*Wn)
{
complex x=a[j+k],y=w*a[j+k+mid];
a[j+k]=x+y;
a[j+k+mid]=x-y;
}
}
}
}
inline void main()
{
n=read(),m=read();
for(int i=0;i<=n;++i) a[i].x=read();
for(int i=0;i<=m;++i) a[i].y=read();
while(limit<=n+m) limit<<=1,++len;
for(int i=0;i<limit;++i)
pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
fft(a,1);
for(int i=0;i<=limit;++i) a[i]=a[i]*a[i];
fft(a,-1);
for(int i=0;i<=n+m;++i) printf("%d ",(int)(a[i].y/limit/2+0.5));
}
}
signed main()
{
red::main();
return 0;
}
注意三次變兩次優化會令精度誤差平方，請根據題目值域考慮是否使用
參考博客
%%%attack
%%%rabbithu