快速傅里叶变换FFT学习笔记
点值表示法
我们正常表示一个多项式的方式，形如 \(A(x)=a_0+a_1x+a_2x^2+...+a_nx^n\)，这是正常人容易看懂的，但是，我们还有一种表示法。
我们知道，\(n+1\)个点可以表示出一个 \(n\) 次的多项式。
于是，我们任意地取 \(n+1\) 个不同的值，表示 \(x\) ，求出的值与 \(x\) 对应，形成 \(n+1\)个点，这就可以表示。
复数
一种表示坐标的方法，对于坐标 \((x,y)\),可写作 \(x+iy\)，其中\(x\)为实部，\(y\)为虚部。
C++中有复数的模板，complex，可以直接作为变量类型使用。
运算规则，自然不用多说了，也就是直接拿式子算即可。
如果你不会，可以看看百度百科。
我们将点至原点的距离称为模长，将其与原点相连之后与 \(x\) 轴形成的一个夹角称为辐角，不过呢，对于第三第四象限的点，自然要加上一个 \(180°\)了。
我的表述自然不够专业，希望可以表述出这个意思吧。
复数的乘法可以理解为，模长相乘，辐角相加。
Tips：
此部分的证明来自 cjx 犇犇。
为啥是这样呢？证明如下:
那么这个点就是 \((ac-bd,bc+ad)\)，其模长:
那么我们应该可以看出来这个模长相乘了。
接下来是辐角相加，我们设原来两个辐角为 \(\theta_1,\theta_2\)，而模长为 \(t_1,t_2\)。
我们知道分母是乘积的模长，分子是横坐标，所以这个式子恰好就是乘积辐角对应的 \(\cos\) 值。
那么显然，辐角的值就是 \(\theta_1+\theta_2\)了。
把圆均分
如图是坐标轴上一个以原点为圆心的半径为 \(1\) 的圆。
我们定义\(\omega_n^k\)表示从\((1,0)\)开始，把圆均分为\(n\)份的第\(k\)个点的复数，其中\((1,0)\)为\(\omega_n^0\)。
那么，不难发现，\(\omega_n^k\)表示的点为 \((\cos(2\pi\frac k n),\sin(2\pi\frac n k))\)。
这是为什么呢？我们考虑它的辐角，由于其平分了一整个圆，所以其辐角为 \(360\frac k n°\)，转换为弧度后则为 \(2\pi\frac k n\)，且模长为 \(1\)，利用三角函数易得其坐标了。
简单推推式子不难发现 \((\omega_n^1)^k=\omega_n^k\)，这是利用了模长相乘，辐角相加，因为模长是 \(1\) ,怎么乘都是 \(1\)，于是辐角不断叠加，从定义上看是这样的。
- 性质1：\(\omega _{dn}^{dk}=\omega_n^k\)，根据定义可证。
- 性质2：\(\omega _{n}^{k+\frac n 2}=-\omega_n^k\)，两点对称。
这个东西有什么用呢？
离散傅里叶变换
离散傅里叶变换(Discrete Fourier Transform，简称DFT)的思想是利用 \(\omega_n^k\)将一个多项式转为点值表示法。
对于一个多项式\(A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\),我们按照前文所云，将所有的 \(\omega_n^k\)作为 \(x\) 代入。
于是我们得到了 \(n-1\) 个点，使用复数形式表示，成为一个数组 \((y_0,y_1,y_2,...,y_{n-1})\)的。
这被称为 \(A(x)\) 的傅里叶变换。
傅里叶逆变换
我们再将其作为一个多项式 \(B(x)=y_0+y_1x+....+y_{n-1}x^{n-1}\)。
对于这个多项式\(B(x)\)，代入所有的 \(\omega_n^{-k}\)，也就是\(\omega_n^k\)的倒数，得到 \((z_0,z_1,z_2,...,z_{n-1})\)。
易得：
关于后面那个等比数列，若 \(j=k\)，可得 \(1\)，否则用等比数列式子可知为 \(0\)。
因此：
所以我们可以求出原来的多项式了。
快速傅里叶变换
快速傅里叶变换(Fast Fourier Transform，简称FFT)，是在 DFT的基础上我们发现时间复杂度依然需要 \(O(n^2)\)，没有含金量，所以我们要给他含金量！
我们可以使用分治的思想，使得时间复杂度降至\(O(n\log n)\)。
对于\(A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\)，我们可以：
同理:
不错！利用这两个式子，我们可以在 \(O(n \log n)\) 的时间复杂度求出 \(A(x)\)。
不过注意这里一直有一个除二操作，为了方便，我们需要把多项式补成一个次数为\(2^x-1\) 的多项式。
可以写一个递归来求解。
注意这个取倒可以利用共轭复数，对于 \(a+ib\)，其共轭复数为 \(a-ib\)。
由于此处 \(a^2+b^2=1\) ,因此其共轭复数为其倒数。
void FFT(cp *a,LL n,bool inv)//inv 表示omega是否取倒
{
if(n==1)return;
static cp buf[N];
LL m=n/2;
for(int i=0;i<=m-1;i++)buf[i]=a[2*i],buf[i+m]=a[2*i+1];//奇偶分开
for(int i=0;i<=n-1;i++)a[i]=buf[i];
FFT(a,m,inv),FFT(a+m,m,inv);
for(int i=0;i<=m-1;i++)
{
cp x=omega(n,i);
if(inv)x=conj(x);//conj(x)求解共轭复数
buf[i]=a[i]+x*a[i+m],buf[i+m]=a[i]-x*a[i+m];
}
for(int i=0;i<=n-1;i++)a[i]=buf[i];
}
利用FFT求解多项式的乘积
这个还是十分简单的，直接把两个多项式转化为两个长度相同的次数为\(2^x-1\)的多项式。
求解出其傅里叶变换形式之后，对于每个点对应的复数，相乘即可。
为什么呢？
首先我们知道当前的 \(a_i\)表示的是 \(A(\omega_n^i)\)，\(b_i\)表示的是 \(B(\omega_n^i)\)，那么我们将两个值直接相乘。
因此多项式相乘以后，我们希望 \(a_i=A(\omega_n^i)B(\omega_n^i)\)，那么就是直接相乘喽。
给一个简单的实现：
、#include<bits/stdc++.h>
#define LL int
#define cp complex<double>
using namespace std;
const double PI=acos(-1.0000);
const int N=5e6+5;
cp omega(LL n, LL k)
{
return cp(cos(2*PI*k/n),sin(2*PI*k/n));
}
LL n,x,len,ans[N];
cp a[N],b[N];
//上文的 FFT 实现省去
int main()
{
scanf("%d",&n);
len=1;
while(len<2*n)len*=2;
for(int i=n-1;i>=0;i--)scanf("%1d",&x),a[i].real(x);
for(int i=n-1;i>=0;i--)scanf("%1d",&x),b[i].real(x);
FFT(a,len,0),FFT(b,len,0);
for(int i=0;i<=len-1;i++)a[i]*=b[i];
FFT(a,len,1);
for(int i=0;i<=len-1;i++)//进位
{
ans[i]+=floor(a[i].real()/len+0.5);
ans[i+1]+=ans[i]/10;
ans[i]%=10;
}
int i=len;
for(i;i>=0&&ans[i]==0;i--);//前导零
if(i==-1)len=0;
for(;i>=0;i--)printf("%d",ans[i]);
}
非递归FFT
这里有一个优化，我们发现每次递归有一个把 \(a_i\) 奇偶分开的过程，本质来看，就是将二进制末尾为 \(0\) 的数字与二进制末尾为 \(1\) 的数字分开。
我们不妨想一下，对于一个数 \(x\)，其位置可以根据其二进制确定，就是其二进制倒过来的数字。
我们先将每个 \(a_i\) 放置在对应的位置，然后向上逐渐合并。
void FFT(cp *a,bool inv)
{
LL lim=0;
while((1<<lim)<len)lim++;
for(int i=0;i<=len-1;i++)
{
LL t=0;
for(int j=0;j<lim;j++)
if((i>>j)&1)t|=(1<<(lim-j-1));//处理其翻转后的值
if(i<t)swap(a[i],a[t]);
}
static cp buf[N];
for(int l=2;l<=len;l*=2)
{
LL m=l/2;
for(LL j=0;j<=len-1;j+=l)
{
for(LL i=0;i<=m-1;i++)
{
cp x=omega(l,i+j);
if(inv)x=conj(x);
buf[i+j]=a[i+j]+x*a[i+j+m];
buf[i+j+m]=a[i+j]-x*a[i+j+m];
}
}
for(int j=0;j<=len-1;j++)a[j]=buf[j];
}
}
蝴蝶操作
这个东西其实就是想了个办法使得把工具人数组 buf 除掉了。
调调顺序即可。
void FFT(cp *a,bool inv)
{
LL lim=0;
while((1<<lim)<len)lim++;
for(int i=0;i<=len-1;i++)
{
LL t=0;
for(int j=0;j<lim;j++)
if((i>>j)&1)t|=(1<<(lim-j-1));
if(i<t)swap(a[i],a[t]);
}
for(int l=2;l<=len;l*=2)
{
LL m=l/2;
for(LL j=0;j<=len-1;j+=l)
{
for(LL i=0;i<=m-1;i++)
{
cp x=omega(l,i+j);
if(inv)x=conj(x);
x*=a[i+j+m];
a[i+j+m]=a[i+j]-x;
a[i+j]=a[i+j]+x;
}
}
}
}
一些小优化
对于 \(i\) 的二进制翻转可以先预处理出来。
然后 \(\omega_n^k\) 可以利用性质累乘，最后代码就长这样了：
#include<bits/stdc++.h>
#define LL int
#define cp complex<double>
using namespace std;
const double PI=acos(-1.0000);
const int N=5e6+5;
cp omega(LL n, LL k)
{
return cp(cos(2*PI*k/n),sin(2*PI*k/n));
}
LL n,len,lim,x,ans[N],r[N];
cp a[N],b[N];
void FFT(cp *a,bool inv)
{
for(int i=0;i<=len-1;i++)
{
LL t=r[i];
if(i<t)swap(a[i],a[t]);
}
for(int l=2;l<=len;l*=2)
{
LL m=l/2;
cp omg=omega(l,1);
for(LL j=0;j<=len-1;j+=l)
{
cp x(1,0);
for(LL i=0;i<=m-1;i++)
{
cp t=x;
if(inv)t=conj(t);
t*=a[i+j+m];
a[i+j+m]=a[i+j]-t,a[i+j]=a[i+j]+t;
x*=omg;
}
}
}
}
int main()
{
scanf("%d",&n);
len=1;
while(len<2*n)len*=2,lim++;
for(int i=0;i<=len-1;i++)
{
LL t=0;
for(int j=0;j<lim;j++)if((i>>j)&1)t|=(1<<(lim-j-1));
r[i]=t;
}
for(int i=n-1;i>=0;i--)scanf("%1d",&x),a[i].real(x);
for(int i=n-1;i>=0;i--)scanf("%1d",&x),b[i].real(x);
FFT(a,0),FFT(b,0);
for(int i=0;i<=len-1;i++)a[i]*=b[i];
FFT(a,1);
for(int i=0;i<=len-1;i++)
{
ans[i]+=floor(a[i].real()/len+0.5);
ans[i+1]+=ans[i]/10;
ans[i]%=10;
}
int i=len;
for(i;i>=0&&ans[i]==0;i--);
if(i==-1)len=0;
for(;i>=0;i--)printf("%d",ans[i]);
}