快速傅里叶变换学习笔记
一.单位根
满足\(x^n-1=0\)的\(x\)称为\(n\)次单位根
1.本原单位根
若 \(\omega\)是\(n\)次单位根,且\(\omega^0\),\(\omega^1\),\(\omega^2\),...,\(\omega^{n-1}\)恰好是所有的\(n\)次单位根，则称\(\omega\)为\(n\)次本原单位根,记作\(\omega_n\)
2.单位根的计算
在复数域\(\mathbb{C}\)上，\(e^{ix}=cosx+isinx\)（欧拉公式，证明的话左边将 \(\exp\) (ix) 展开成 \(\texttt{EGF}\) 的形式，右边同理即证）
重要的性质:\(n\)次单位根在复平面上等分单位圆
由"重要的性质"可以得到,\(\omega_n\)=\(\exp(\frac{2 \pi i}{n})\)=\(cos \frac{2\pi }{n}\)+\(isin \frac{2\pi }{n}\)
以此我们可以表示出所有的\(n\)次单位根，即:\(\omega_n^k\)=\(cos \frac{2\pi k}{n}\)+\(isin \frac{2\pi k}{n}\)
二.离散傅里叶变换DFT
前置芝士1:多项式的点值表达
对于一个多项式\(f(x)\)，我们代入每个\(x_i\),会得到对应的\(f(x_i)\),这些\((x_i,f(x_i))\)构成了这个多项式的点值表达，且唯一确定了这个多项式
前置芝士2：多项式相乘
两个多项式的乘积被定义为:
\(A(x)B(x)\)=\(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}{a_ib_jx^{i+j}}\)=\(\sum\limits_{i=1}^{2n}{c_kx^k}\)
其中\(c\)是\(a\)和\(b\)的卷积,计算两个多项式相乘的朴素算法是\(O(n^2)\)的
如果将两个多项式先转换为点值表达,再逐项相乘，时间复杂度只有\(O(n)\)
即 \((x,A(x))*(x,B(x))=(x,A(x)*B(x)）\)
那么，复杂度优化的瓶颈在于将多项式转化为点值表达（以及转回去），FFT解决的正是这个问题
1.未经优化的DFT
设\(a\)是长度为\(n\)的数列，对于\(0\leq k <n\),定义
\(b_k\)=\(\sum\limits_{i=0}^{n-1}{(a_i\omega^{ki})}\)
其中\(b\)即为\(a\)的\(DFT\)
这是，我们令多项式\(f(x)\)=\(\sum a_ix^i\),则\(b_k\)就是\(f(x)\)在\(\omega^k\)处的点值
当我们计算完两个多项式乘积的点值,如何将其转回多项式?
2.DFT的逆变换IDFT
\(a_k\)=\(\frac{1}{n}\sum\limits_{i=0}^{n-1}{b_i\omega^{-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)\)=\(\sum\limits_{i=0}^{n-1}a_ix^i\)
=\(\sum\limits_{i=0}^{m-1}a_{2i}x^{2i}\)+\(\sum\limits_{i=0}^{m-1}a_{2i+1}x^{2i+1}\)
=\(\sum\limits_{i=0}^{m-1}a_{2i}x^{2i}\)+\(x\sum\limits_{i=0}^{m-1}a_{2i+1}x^{2i}\)
我们令\(A_0(x)=\sum\limits_{i=0}^{m-1}a_{2i}x^{i}\),\(A_1(x)=\sum\limits_{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('\n')
#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('\n')
#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('\n')
#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;
}
完结撒花！！！
