快速傅里叶变换 学习笔记
快速傅里叶变换 学习笔记
- 前言：这篇学习笔记以\(\text{Dew}\)的意会
yy为主，有些地方会比较简略，不过该有的证明应该还是会有的。
多项式的表示法
-
系数表示法
\(f(x)=a_0+a_1x+a_2x^2+\dots+a_nx^n\)
-
点值表示法
\(n+1\)个不同的点可以确定一个\(n\)次的多项式。
即一个\(n\)次多项式可以被\((x_1,y_1),(x_2,y_2),\dots,(x_n,y_n)\)全部表示出来
-
考虑如何求出点值表示法的所有点
复数与单位根
-
定义复数运算
-
\(a+bi+c+di=a+c+(b+d)i\)
\(a+bi-(c+di)=a-c+(b-d)i\)
几何定义同向量
-
\((a+bi)(c+di)=ac-bd+(ad+bc)i\)
几何定义：模长相乘，幅角相加
struct complex { double x,y; complex(){} complex(double x,double y){this->x=x,this->y=y;} complex friend operator +(complex n1,complex n2){return complex(n1.x+n2.x,n1.y+n2.y);} complex friend operator -(complex n1,complex n2){return complex(n1.x-n2.x,n1.y-n2.y);} complex friend operator *(complex n1,complex n2){return complex(n1.x*n2.x-n1.y*n2.y,n1.x*n2.y+n1.y*n2.x);} }; -
-
单位根
-
按照乘法的几何定义，我们发现在模长为\(1\)的单位圆上，两个负数相乘的模长不变，于是考虑单位圆上的表示法。
-
在半径为\(1\)的单位圆上，向量的坐标为\((\cos \theta,i\sin \theta)\)，我们可以稍稍证明一下幅角相加了。
设有负数\((\cos \theta_1,i\sin \theta_1)\),\((\cos \theta_2,i\sin \theta_2)\)，那么它们相乘有
\(\cos \theta_1\cos \theta_2-\sin\theta_1\sin\theta_2+i(\cos\theta_1\sin\theta_2+\sin\theta_1\cos\theta_2)\)
\(=\cos(\theta_1+\theta_2)+i\sin(\theta_1+\theta_2)\)
-
将复平面的单位圆划分成\(n\)份，定义第\(k\)份的向量对应的复数为\(w_n^k\)。特殊的\(w_n^0=w_n^n=1\)。
为了方便和必要，以下的\(n\)均认为是\(2\)的正整次数幂。
-
欧拉公式：
\(w_n^k\)的坐标表示为\((\cos\frac{2k\pi}{n},i\sin\frac{2k\pi}{n})\)
-
性质
- \((w_n^k)^2=w_n^{2k}=w_{\frac{n}{2}}^k\)
- \(w_n^{k+\frac{n}{2}}=-w_n^k\)
-
快速傅里叶变换
对多项式\(F(x)=a_0+a_1x+a_2x^2+\dots+a_{n-1}x^{n-1}\)
\(F(x)=a_0+a_2x^2+\dots+a_{n-2}x^{n-2}+x(a_1+a_3x^2+\dots a_{n-1}x^{n-2})\)
设\(FL(x)=a_0+a_2x+\dots+a_{n-2}x^{\frac{n}{2}-1}\),\(FR(x)=a_1+a_3x+\dots+a_{n-1}x^{\frac{n}{2}-1}\)
则\(F(x)=FL(x^2)+xFR(x^2)\)
取\(k<\frac{n}{2}\)，并代入\(x=w_n^k\)
\(F(w_n^k)=FL(w^k_{\frac{n}{2}})+w_n^kFR(w^k_{\frac{n}{2}})\)
\(F(w_n^{k+\frac{n}{2}})=FL(w^k_\frac{n}{2})-w_n^kFR(w_\frac{n}{2}^k)\)
显然可以子问题求解，于是我们得到了\(O(nlogn)\)求点值表示法的方法了。
void FFT(int len,int *a,int typ)
{
if(len==1) return;
complex a1[len>>1],a2[len>>1];
for(int i=0;i<n;i+=2)
a1[i>>1]=a[i],a2[i>>1]=a[i+1];
FFT(len>>1,a1,typ);
FFT(len>>1,a2,typ);
complex wn=complex(cos(2*pi/len),typ*sin(2*pi/len)),w=(1,0);
for(int i=0;i<len>>1;i++.w=w*wn)
{
a[i]=a1[i]+w*a2[i];
a[i+(len>>1)]=a1[i]-w*a2[i];
}
}
快速傅里叶逆变换
显然我们可以有
\(\begin{bmatrix}1&1&1&\cdots&1\\1 &w_n^1&w_n^2&\cdots&w_n^{n-1}\\\vdots& \vdots &\vdots&\ddots&\vdots\\1& w_n^{n-1}&w_n^{2(n-1)}&\cdots&w_n^{(n-1)(n-1)}\\\end{bmatrix}\times \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\a_{n-1}\\ \end{bmatrix}=\begin{bmatrix} y_0 \\ y_1 \\ \vdots \\y_{n-1}\\\end{bmatrix}\)
现在我们有\(y\)和\(w\)那些东西，要求\(a\)，按道理我们要找到一个\(w\)那个的逆矩阵然后乘在右边，然后我们有结论
\(\begin{bmatrix}1&1&1&\cdots&1\\1 &w_n^1&w_n^2&\cdots&w_n^{n-1}\\\vdots& \vdots &\vdots&\ddots&\vdots\\1& w_n^{n-1}&w_n^{2(n-1)}&\cdots&w_n^{(n-1)(n-1)}\\\end{bmatrix}\times \frac{1}{n}\begin{bmatrix}1&1&1&\cdots&1\\1&w_n^{-1}&w_n^{-2}&\cdots&w_n^{-(n-1)}\\\vdots &\vdots&\vdots&\ddots&\vdots\\1& w_n^{-(n-1)}&w_n^{-2(n-1)}&\cdots&w_n^{-(n-1)(n-1)}\\\end{bmatrix}=I\)
这个结合前面单位根的性质很容易意会到的。
然后发现这个也可以看做一个系数表示法转点值表示法的过程，而且和前面的系数矩阵很像，所以只需要改个参数就搞到啦。
\(Code:\)(递归版)
#include <cstdio>
#include <cmath>
const int N=(1<<21)+10;
struct complex
{
double x,y;
complex(){}
complex(double x,double y){this->x=x,this->y=y;}
complex friend operator +(complex n1,complex n2){return complex(n1.x+n2.x,n1.y+n2.y);}
complex friend operator -(complex n1,complex n2){return complex(n1.x-n2.x,n1.y-n2.y);}
complex friend operator *(complex n1,complex n2){return complex(n1.x*n2.x-n1.y*n2.y,n1.x*n2.y+n1.y*n2.x);}
}a[N],b[N];
int n,m;
const double pi=3.1415926535897;
void FFT(int len,complex *a,int typ)
{
if(len==1) return;
complex a1[len>>1],a2[len>>1];
for(int i=0;i<len;i+=2)
a1[i>>1]=a[i],a2[i>>1]=a[i+1];
FFT(len>>1,a1,typ);
FFT(len>>1,a2,typ);
complex wn=complex(cos(2*pi/len),typ*sin(2*pi/len)),w=complex(1,0);
for(int i=0;i<len>>1;i++,w=w*wn)
{
a[i]=a1[i]+w*a2[i];
a[i+(len>>1)]=a1[i]-w*a2[i];
}
}
int main()
{
scanf("%d%d",&n,&m);
for(int i=0;i<=n;i++) scanf("%lf",&a[i].x);
for(int i=0;i<=m;i++) scanf("%lf",&b[i].x);
int len=1;while(len<=n+m) len<<=1;
FFT(len,a,1),FFT(len,b,1);
for(int i=0;i<=len;i++) a[i]=a[i]*b[i];
FFT(len,a,-1);
for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/len+0.5));
return 0;
}
\(Butterfly\)优化
- 递归版一看就好慢的。
如上图，
| 原数列数字 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
|---|---|---|---|---|---|---|---|---|
| 二进制 | 000 | 001 | 010 | 011 | 100 | 101 | 110 | 111 |
| 最底层数列数字 | 0 | 4 | 2 | 6 | 1 | 5 | 3 | 7 |
|---|---|---|---|---|---|---|---|---|
| 二进制 | 000 | 100 | 010 | 110 | 001 | 101 | 011 | 111 |
然后我们发现二进制表示反过来了。
于是我们可以得到一个转换方法。
for(int i=0;i<len;i++) turn[i]=turn[i>>1]>>1|((i&1)<<L);
然后直接从下往上做就不需要进行递归啦
Code:(非递归版)
#include <cstdio>
#include <algorithm>
#include <cmath>
const int N=(1<<21)+10;
struct complex
{
double x,y;
complex(){}
complex(double x,double y){this->x=x,this->y=y;}
complex friend operator +(complex n1,complex n2){return complex(n1.x+n2.x,n1.y+n2.y);}
complex friend operator -(complex n1,complex n2){return complex(n1.x-n2.x,n1.y-n2.y);}
complex friend operator *(complex n1,complex n2){return complex(n1.x*n2.x-n1.y*n2.y,n1.x*n2.y+n1.y*n2.x);}
}a[N],b[N],tmpx,tmpy,wn,w;
int n,m,len=1,L=-1,turn[N];
const double pi=3.1415926535897;
void FFT(complex *a,int typ)
{
for(int i=0;i<len;i++)
if(turn[i]>i) std::swap(a[i],a[turn[i]]);
for(int le=1;le<len;le<<=1)
{
wn=complex(cos(pi/le),typ*sin(pi/le));
for(int dl=le<<1,p=0;p<len;p+=dl)
{
w=complex(1,0);
for(int k=0;k<le;k++,w=w*wn)
{
tmpx=a[p+k],tmpy=w*a[p+k+le];
a[p+k]=tmpx+tmpy;
a[p+k+le]=tmpx-tmpy;
}
}
}
}
int main()
{
scanf("%d%d",&n,&m);
for(int i=0;i<=n;i++) scanf("%lf",&a[i].x);
for(int i=0;i<=m;i++) scanf("%lf",&b[i].x);
while(len<=n+m) len<<=1,++L;
for(int i=0;i<len;i++) turn[i]=turn[i>>1]>>1|((i&1)<<L);
FFT(a,1),FFT(b,1);
for(int i=0;i<=len;i++) a[i]=a[i]*b[i];
FFT(a,-1);
for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/len+0.5));
return 0;
}
