[日志]9-25
我不是柳橙汁
2019-09-25 11:33:04
## 回来第一次复习FFT
今天就写一下FFT的小结
在此安利我所复习的博客，stO [自为风月马前卒的FFT详解](https://www.cnblogs.com/zwfymqz/p/8244902.html)
### 前置知识
1.复数相关
复数的表示法
$a+bi$
复数加减乘
$$(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$$
复数可在直角坐标系中表示
2.点值表示法
如一个多项式
他的系数表示法为
$$A_n=a_0+a_1x+a_2x^2+...+a_nx^n$$
假定这么一个多项式
$$A_n=2+5x+3x^2$$
将若干个$x$代入多项式中，得到一系列点$(x_i,y_i)$
比如这个多项式，我们将$x=0,x=1,x=2$代入
得到$(0,2),(1,10),(2,24)$这一系列点
我们也可以用这一系列点来表示这个多项式
为什么要用到这个呢
在正常的多项式乘法中
若是两个多项式相乘
需要将两个多项式中每个项都乘一遍，复杂度$o(n^2)$
若是采用点值表示法，只需要将每两个对应的点相乘，复杂度$o(n)$
3.1.**单位根**
以下内容默认n为2的整数次幂
在复数平面上，单位圆就是以原点为圆心，半径为1的一个圆
将这个圆n等分，按照分割线作n个向量
设幅度为正最小的那个向量为$\omega_n$
我们称之为n次单位根
可以知道其他n-1个向量分别为$\omega_n^2,\omega_n^3,...,\omega_n^n$
其中$\omega_n^n=\omega_n^0=1$
这些单位根的值可以通过欧拉公式计算
$$\omega_n^k=\cos \frac{2k\pi}{n}+i\sin \frac{2k\pi}{n}$$
在代数中，若$z^k=1$，我们称$z$为$k$次单位根
3.2.单位根的性质
$$\omega_n^k=\cos \frac{2k\pi}{n}+i\sin \frac{2k\pi}{n}$$
$$\omega_{2n}^{2k}=\omega_n^k$$
$$\omega_n^{k+\frac{n}{2}}=-\omega_n^k$$
$$\omega_n^n=\omega_n^0=1$$
$k\not= 0$时
$$\sum_{i=1}^n\omega_n^{ki}=0$$
$k=0$时
$$\sum_{i=1}^n\omega_n^{ki}=n$$
### 快速傅里叶变换(FFT)(Fast Fourier Transfromation)
我们设一个多项式
$$A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}^{n-1}$$
将其按照奇偶拆开
$$=(a_0+a_2x^2+....+a_{n-2}x^{n-2})+(a_1x+a_3x^3+...+a_{n-1}x^{n-1})$$
奇数项提出一个x
$$=(a_0+a_2x^2+....+a_{n-2}x^{n-2})$$
$$+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})$$
令
$$A_1(x)=a_0+a_2x+....+a_{n-2}x^{\frac{n}{2}-1}$$
$$A_2(x)=a_1+a_3x+...+a_{n-1}x^{\frac{n}{2}-1}$$
可知
$$A(x)=A_1(x^2)+xA_2(x^2)$$
将$\omega_n^k$代入
$$A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})$$
再将$\omega_n^{k+\frac{n}{2}}$代入
$$A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n})$$
$$=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})$$
我们将两个式子比较
$$A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})$$
$$A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})$$
发现枚举其中一个时，可以在$o(1)$的时间将另一个枚举出来
每次我们可以将问题缩小一半
这样就类似线段树的时间复杂度$o(nlogn)$
本身需要$o(n^2)$的时间，我们拿$o(n+2nlogn)$就完成了，得到了优化
### 快速傅里叶逆变换(IFFT)(Fast Fourier Inverse Transformation)
但是，我们常用的表示法并不是点值表示法
我们还应该将点值表示法转化为系数表示法
这个过程叫做**傅里叶逆变换**
首先
我们设$(y_1,y_2,...,y_n)$为$(a_1,a_2,...,a_n)$的点值表示
设另一个向量$(c_1,c_2,...,c_n)$满足
$$c_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i$$
$$=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^{i})^j)(\omega_n^{-k})^i$$
$$=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^{i})^j)(\omega_n^{-k})^i$$
$$=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^i)^{(j-k)}$$
于是我们知道
$j=k$时
$$\sum_{i=0}^{n-1}(\omega_n^i)^{(j-k)}=n$$
而$j\not=k$时
$$\sum_{i=0}^{n-1}(\omega_n^i)^{(j-k)}=0$$
于是
$$c_k=na_k$$
即
$$a_k=\frac{c_k}{n}$$
所以我们对乘起来的点值表示再做一次FFT
然后除去n就是所得到的多项式了
所用时间复杂度为$o(2nlogn+n+nlogn)$，即$o(nlogn)$
### 迭代实现
我们也可以尝试用递归的方法来完成FFT，但这样并不是最优的办法
如果我们事先就将需要转换的位置弄明白，这样可以直接转换
可以发现
反转后的序列的二进制数也反转过来
于是我们可以预处理我们需要得到的数列
就可以在$o(n)$的时间复杂度内完成了
### 板子题
[【模板】多项式乘法](https://www.luogu.org/problem/P3803)
```cpp
#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<iostream>
using namespace std;
#define fr(i,a,b) for(register long long i=a;i<=b;i++)
const double Pi=acos(-1.0);
struct complex{
double x,y;
complex(double xx=0,double yy=0){x=xx,y=yy;}
}a[5000010],b[5000010];
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);}
long long n,m,limit,l;
long long r[5000010];
inline long long v_in(){
char ch=getchar();long long sum=0,f=1;
while(!isdigit(ch)){if(ch=='-') f=-1;ch=getchar();}
while(isdigit(ch)) sum=(sum<<3)+(sum<<1)+(ch^48),ch=getchar();
return sum*f;
}
void FFT(complex *A,long long type){
fr(i,0,limit-1) if(i<r[i]) swap(A[i],A[r[i]]);
for(register long long mid=1;mid<limit;mid<<=1){
complex Wn(cos(Pi/mid),type*sin(Pi/mid));
for(register long long j=0;j<limit;j+=(mid<<1)){
complex W(1,0);
for(register long long k=0;k<mid;k++,W=W*Wn){
complex x=A[j+k],y=W*A[j+mid+k];
A[j+k]=x+y;
A[j+mid+k]=x-y;
}
}
}
}
void inpt(){
n=v_in(),m=v_in();
fr(i,0,n) a[i].x=v_in();
fr(i,0,m) b[i].x=v_in();
limit=1,l=0;
while(limit<n+m+1) limit<<=1,l++;
fr(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
void work(){
FFT(a,1);
FFT(b,1);
fr(i,0,limit) a[i]=a[i]*b[i];
FFT(a,-1);
}
int main(){
inpt();
work();
fr(i,0,n+m) printf("%lld ",(int)(a[i].x/limit+0.5));
printf("\n");
return 0;
}
```