一.单位根
满足(x^n-1=0)的(x)称为(n)次单位根
1.本原单位根
若 (omega)是(n)次单位根,且(omega^0),(omega^1),(omega^2),...,(omega^{n-1})恰好是所有的(n)次单位根，则称(omega)为(n)次本原单位根,记作(omega_n)
重要的性质:(n)次单位根在复平面上等分单位圆
2.单位根的计算
在复数域(mathbb{C})上，(e^{ix}=cosx+isinx)
由前面"重要的性质"可以得到,(omega_n)=(exp(2pi in))=(cos frac{2pi }{n})+(isin frac{2pi }{n})
以此我们可以表示出所有的(n)次单位根，即:(omega_n^k)=(cos frac{2pi k}{n})+(isin frac{2pi k}{n})
二.离散傅里叶变换DFT
前置芝士1:多项式的点值表达
对于一个多项式(f(x))，我们代入每个(x_i),会得到对应的(f(x_i)),这些((x_i,f(x_i)))构成了这个多项式的点值表达，且唯一确定了这个多项式
前置芝士2：多项式相乘
两个多项式的乘积被定义为:
(A(x)B(x))=(sumlimits_{i=1}^{n}sumlimits_{j=1}^{n}{a_ib_jx^{i+j}})=(sumlimits_{i=1}^{2n}{c_kx^k})
其中(c)是(a)和(b)的卷积,计算两个多项式相乘的朴素算法是(O(n^2))的
如果将两个多项式先转换为点值表达,再相乘，时间复杂度只有(O(n))
那么，复杂度优化的瓶颈在于将多项式转化为点值表达，FFT解决的正是这个问题
1.未经优化的DFT
设(a)是长度为(n)的数列，对于(0leq k <n),定义
(b_k)=(sumlimits_{i=0}^{n-1}{(a_iomega^{ki})})
其中(b)即为(a)的(DFT)
这是，我们令多项式(f(x))=(sum a_ix^i),则(b_k)就是(f(x))在(omega^k)处的点值
当我们计算完两个多项式乘积的点值,如何将其转回原多项式?
2.DFT的逆变换IDFT
(a_k)=(frac{1}{n}sumlimits_{i=0}^{n-1}{b_iomega^{-ki}})
这个柿子与DFT的相似度极高的特性使得我们不需要重新写一个函数来处理IDFT，只需要提前预处理单位根的倒数（即共轭复数），利用FFT转化，最后全部除以 (n) 即可
四.快速傅里叶变换FFT
前置芝士:单位根的性质
((1)omega_{2n}^{2k}=omega_n^k)
((2)omega_{2n}^{n+k}=- omega_{2n}^k)
1.FFT
啥是FFT？就是利用DFT的奇偶性质快速计算DFT的一个分治算法
设(n=2m),将(A(x))按次数奇偶分类，有:
(A(x))=(sumlimits_{i=0}^{n-1}a_ix^i)
=(sumlimits_{i=0}^{m-1}a_{2i}x^{2i})+(sumlimits_{i=0}^{m-1}a_{2i+1}x^{2i+1})
=(sumlimits_{i=0}^{m-1}a_{2i}x^{2i})+(xsumlimits_{i=0}^{m-1}a_{2i+1}x^{2i})
我们令(A_0(x)=sumlimits_{i=0}^{m-1}a_{2i}x^{i}),(A_1(x)=sumlimits_{i=0}^{m-1}a_{2i+1}x^{i})
则(A(x)=A_0(x^2)+xA_1(x^2))
通过上式我们可以得出，如果我们得到了(A_0(x))(A_1(X))在(omega_m^0),(omega_m^1),(omega_m^2),...,(omega_m^{n-1})处的点值，则可在(O(n))内得到(A(x))在(omega_m^0),(omega_m^1),(omega_m^2),...,(omega_m^{n-1})处的点值
而(A_0(x),A_1(x))是可以分治计算的,递归深度不会超过(log n)
综上，我们可以在(O(n log n))内完成对(A(x))的点值求值
五.一些优化
1.位逆序置换
我们发现，令(rev(x))表示(x)经过二进制反转后的数，且令(b_i=a_{rev(i)}),则每次对(a_i)进行的操作在(b_i)中变为了对相邻两个序列进行的操作，那么我们只需要预先处理出(b_i)，直接递归向上不断还原即可
它在代码中一般这样子体现：
fr(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
2.蝴蝶操作
因为double的乘法比较慢，我们在代码中有这样一段：
Complex x=arr[j+k];
arr[j+k]=x+w*arr[j+mid+k],arr[j+mid+k]=x-w*arr[j+mid+k];
我们发现w*arr[j+mid+k]这东西被算了两次，所以应该这么写
Complex x=arr[j+k],y=w*arr[j+mid+k];
arr[j+k]=x+y,arr[j+mid+k]=x-y;
然后就优化完了/cy
所以这难道不是常数优化？
六.代码实现
P3803 【模板】多项式乘法（FFT）
#include<bits/stdc++.h>
using namespace std;
namespace my_std
{
typedef long long ll;
typedef double db;
#define pf printf
#define pc putchar
#define fr(i,x,y) for(register int i=(x);i<=(y);++i)
#define pfr(i,x,y) for(register int i=(x);i>=(y);--i)
#define go(u) for(int i=head[u];i;i=e[i].nxt)
#define enter pc('
')
#define space pc(' ')
#define fir first
#define sec second
#define MP make_pair
const int inf=0x3f3f3f3f;
const ll inff=1e15;
inline int read()
{
int sum=0,f=1;
char ch=0;
while(!isdigit(ch))
{
if(ch=='-') f=-1;
ch=getchar();
}
while(isdigit(ch))
{
sum=sum*10+(ch^48);
ch=getchar();
}
return sum*f;
}
inline void write(int x)
{
if(x<0)
{
x=-x;
pc('-');
}
if(x>9) write(x/10);
pc(x%10+'0');
}
inline void writeln(int x)
{
write(x);
enter;
}
inline void writesp(int x)
{
write(x);
space;
}
}
using namespace my_std;
const int maxn=1e7+50;
struct Complex
{
double x,y;
Complex(double xx=0,double yy=0){x=xx,y=yy;}
}a[maxn],b[maxn];
double PI=acos(-1.0);
Complex operator + (Complex a,Complex b){return Complex(a.x+b.x,a.y+b.y);}
Complex operator - (Complex a,Complex b){return Complex(a.x-b.x,a.y-b.y);}
Complex operator * (Complex a,Complex b){return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int N,M,limit=1,l,r[maxn];
inline void fft(Complex *arr,int pd)
{
fr(i,0,limit-1) if(i<r[i]) swap(arr[i],arr[r[i]]);
for(int mid=1;mid<limit;mid<<=1)
{
Complex Wn(cos(PI/mid),pd*sin(PI/mid));
for(int j=0,R=mid<<1;j<limit;j+=R)
{
Complex w(1,0);
for(int k=0;k<mid;++k,w=w*Wn)
{
Complex x=arr[j+k],y=w*arr[j+mid+k];
arr[j+k]=x+y,arr[j+mid+k]=x-y;
}
}
}
}
int main(void)
{
N=read(),M=read();
fr(i,0,N) a[i].x=read();
fr(i,0,M) b[i].x=read();
//system("pause");
while(limit<=M+N) limit<<=1,++l;
fr(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
fft(a,1),fft(b,1);
fr(i,0,limit) a[i]=a[i]*b[i];
fft(a,-1);
fr(i,0,M+N) writesp((int)(a[i].x/limit+0.5));
return 0;
}
快速数论变换NTT
#include<bits/stdc++.h>
using namespace std;
namespace my_std
{
typedef long long ll;
typedef double db;
#define pf printf
#define pc putchar
#define fr(i,x,y) for(register ll i=(x);i<=(y);++i)
#define pfr(i,x,y) for(register ll i=(x);i>=(y);--i)
#define go(u) for(ll i=head[u];i;i=e[i].nxt)
#define enter pc('
')
#define space pc(' ')
#define fir first
#define sec second
#define MP make_pair
const ll inf=0x3f3f3f3f;
const ll inff=1e15;
inline ll read()
{
ll sum=0,f=1;
char ch=0;
while(!isdigit(ch))
{
if(ch=='-') f=-1;
ch=getchar();
}
while(isdigit(ch))
{
sum=sum*10+(ch^48);
ch=getchar();
}
return sum*f;
}
inline void write(ll x)
{
if(x<0)
{
x=-x;
pc('-');
}
if(x>9) write(x/10);
pc(x%10+'0');
}
inline void writeln(ll x)
{
write(x);
enter;
}
inline void writesp(ll x)
{
write(x);
space;
}
}
using namespace my_std;
const ll maxn=1e7+50,G=3,mod=998244353;
ll N,M,limit=1,a[maxn],b[maxn],l,r[maxn];
inline ll ksmod(ll a,ll b)
{
ll ans=1;
while(b)
{
if(b&1) ans=(ans*a)%mod;
a=(a*a)%mod;
b>>=1;
}
return ans%mod;
}
inline void NTT(ll *a,ll pd)
{
fr(i,0,limit-1) if(i<r[i]) swap(a[i],a[r[i]]);
for(ll i=1;i<limit;i<<=1)
{
ll gn=ksmod(G,(mod-1)/(i<<1));
for(ll j=0;j<limit;j+=(i<<1))
{
ll g=1;
for(ll k=0;k<i;++k,g=(g*gn)%mod)
{
ll x=a[j+k],y=g*a[j+k+i]%mod;
a[j+k]=(x+y)%mod,a[j+k+i]=(x-y+mod)%mod;
}
}
}
if(pd==1) return ;
ll inv=ksmod(limit,mod-2);
reverse(a+1,a+limit);
fr(i,0,limit-1) a[i]=a[i]*inv%mod;
}
int main(void)
{
N=read(),M=read();
fr(i,0,N) a[i]=read();
fr(i,0,M) b[i]=read();
while(limit<=M+N) limit<<=1,++l;
fr(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
NTT(a,1),NTT(b,1);
fr(i,0,limit-1) a[i]=(a[i]*b[i])%mod;
NTT(a,-1);
fr(i,0,M+N) writesp(a[i]%mod);
return 0;
}
洛谷1919【模板】A*B Problem升级版（FFT快速傅里叶）
#include<bits/stdc++.h>
using namespace std;
namespace my_std
{
typedef long long ll;
typedef double db;
#define pf printf
#define pc putchar
#define fr(i,x,y) for(register int i=(x);i<=(y);++i)
#define pfr(i,x,y) for(register int i=(x);i>=(y);--i)
#define go(u) for(int i=head[u];i;i=e[i].nxt)
#define enter pc('
')
#define space pc(' ')
#define fir first
#define sec second
#define MP make_pair
const int inf=0x3f3f3f3f;
const ll inff=1e15;
inline int read()
{
int sum=0,f=1;
char ch=0;
while(!isdigit(ch))
{
if(ch=='-') f=-1;
ch=getchar();
}
while(isdigit(ch))
{
sum=sum*10+(ch^48);
ch=getchar();
}
return sum*f;
}
inline void write(int x)
{
if(x<0)
{
x=-x;
pc('-');
}
if(x>9) write(x/10);
pc(x%10+'0');
}
inline void writeln(int x)
{
write(x);
enter;
}
inline void writesp(int x)
{
write(x);
space;
}
}
using namespace my_std;
const int maxn=2e6+50;
struct Complex
{
double x,y;
Complex(double xx=0,double yy=0){x=xx,y=yy;}
}a[maxn],b[maxn];
double PI=acos(-1.0);
Complex operator + (Complex a,Complex b){return Complex(a.x+b.x,a.y+b.y);}
Complex operator - (Complex a,Complex b){return Complex(a.x-b.x,a.y-b.y);}
Complex operator * (Complex a,Complex b){return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int N,M,limit=1,l,r[maxn],ans[maxn];
char sa[maxn],sb[maxn];
inline void fft(Complex *arr,int pd)
{
fr(i,0,limit-1) if(i<r[i]) swap(arr[i],arr[r[i]]);
for(int mid=1;mid<limit;mid<<=1)
{
Complex Wn(cos(PI/mid),pd*sin(PI/mid));
for(int j=0,R=mid<<1;j<limit;j+=R)
{
Complex w(1,0);
for(int k=0;k<mid;++k,w=w*Wn)
{
Complex x=arr[j+k],y=w*arr[j+mid+k];
arr[j+k]=x+y,arr[j+mid+k]=x-y;
}
}
}
}
int main(void)
{
scanf("%s%s",sa,sb);
int lena=0,lenb=0,hhh=strlen(sa),hhhh=strlen(sb);
pfr(i,hhh-1,0) a[lena++].x=sa[i]-48;
pfr(i,hhhh-1,0) b[lenb++].x=sb[i]-48;
while(limit<lena+lenb) limit<<=1,++l;
fr(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
fft(a,1),fft(b,1);
fr(i,0,limit) a[i]=a[i]*b[i];
fft(a,-1);
int tot=0;
fr(i,0,limit)
{
ans[i]+=(int)(a[i].x/limit+0.5);
if(ans[i]>=10) ans[i+1]+=ans[i]/10,ans[i]%=10,limit+=(i==limit);
}
while(!ans[limit]&&limit>=1) --limit;
++limit;
while(--limit>=0) write(ans[limit]);
return 0;
}
完结撒花！！！