FFT

FFT

1.前言

FFT 要涉及很多前置基本概念，例如向量，复数等，在这里向量等简单概念不提。

2.1 复数

设 $a,b$ 为实数，$i^2=-1$ ，如果一个数 $z$ ，满足 $z=a+bi$ 的数叫复数，其中 $a$ 为实部，$b$ 为虚部，$i$ 为虚数单位，当 $b=0$ 时，称 $z$ 为实数；当 $b≠0$，且 $a=0$ 时，称 $z$ 为纯虚数。

2.2 复数域

复数域是实数域的代数闭包，即任何复系数多项式在复数域中总有根。

2.3 复平面

复平面是的每个点都表示一个复数，复数 $z=a+bi$ 在复平面中对应的坐标为 $(a,b)$

在复平面中， $x$ 轴又称为实轴；$y$ 轴又称为虚轴

$y$ 轴上有且仅有一个实点即为原点 $(0,0)$

复数集 $C$ 和复平面内所有点所构成的集合是一一对应的。

2.4 复数模

复数的实部与虚部的平方和的算术平方根称为该复数的模，记作 $|z|$

即对于复数 $z=a+bi$，它的模为 $|z|=\sqrt{a^2+b^2}$

复数 $z=a+bi$ 在复平面中的坐标为 $(a,b)$，到原点的距离正好是 $|z|=\sqrt{a^2+b^2}$

2.5 幅角

任意一个复数 $z$，都可以写成 $z=r×(cosθ+i×sinθ)$

其中 $r=|z|$，$θ$ 为 $z$ 的幅角，记做 $Arg(z)$

任意一个非 $0$ 的复数的幅角有无限多个取值，且这些值相差 $2π$ 的整数倍。

我们把满足 $−π≤θ<π$ 的辐角 $θ$ 的值，叫做辐角的主值，记作 $arg(z)$，$arg(z)$ 是唯一的。

如果放到复平面中，复数的幅角可以看作从 $x$ 轴正半轴到复数的转角

2.6 复数运算

2.6.1 平行四边形定则

复数的加法和向量的加法都满足平行四边形定则

2.6.2 加法

设复数 $z_1=a_1+b_1i,z_2=a_2+b_2i$

则 $z_1+z_2= a_1+b_1i+a_2+b_2i=(a_1+a_2)+(b_1+b_2)i$

2.6.1.1 乘法---代数定义

设复数 $z_1=a_1+b_1i,z_2=a_2+b_2i$

则 $z_1z_2=(a_1+b_1i)(a_2+b_2i)=(a_1a_2−b_1b_2)+(a_1b_2+a_2b_1)i$

2.6.1.2 乘法---几何意义

复数相乘，模长相乘，幅角相加

设复数 $z_1=r_1(cosθ_1+isinθ_1),z_2=r_2(cosθ_2+i×sinθ_2)$

$r_1$ 和 $r_2$ 分别为 $z_1$、$z_2$ 的模

则 $z_1z_2=r_1r_2(cos(θ_1+θ_2)+i×sin(θ_1+θ_2))$

3.1 单位根

在复平面中，以原点 $(0,0)$ 为圆心，$1$ 为半径作圆，将所作的圆称为单位圆

将圆 $n$ 等分，并将每个等分点标号

将 $n$ 等分后的 $1$ 号点记做 $ω_n$，称为 $n$ 次单位根

3.2 性质

性质1 ： 设 $n$ 等分后的 $j$ 号点的值为 $z_j$，则有 $z_j=z_1^j=ω_j^n$

性质2 ： $ω_{2n}^{2j}=ω_{n}^{j}$

性质3 ： $ω_{2n}^{n}=-1$

性质4 ： $ω_{2n}^{k+n}=-ω_{2n}^{k}$

性质5 ： $ω_{n}^{k+n}=ω_{n}^{k}$

4. 快速傅里叶 FFT

先再来了解了解多项式。

我们设定 $n$ 为 $2^k$，$k$ 为整数，不足则在多项式后面补 $0$

设一个 $n$ 次多项式 $A(x)=a_0+a_1x+a_2x^2⋯+a_{n−2}x^{n−2}+a_{n−1}x^{n−1}$

易知 $A(x)=(a_0+a_2x^2⋯+a_{n−2}x^{n−2})+(a_1x+a_3x^3⋯+a_{n−1}x^{n−1})$

设 :

$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)$

假设 $0≤k<n^2$，将 $x=ω_n^k$ 代入，得

$A(ω_{n}^{k})=A_1(ω_{n}^{2k})+ω_{n}^{k}A_2(ω_{n}^{2k})$

将 $x=ω_n^{k+\frac{n}{2}}$ 代入，得

$A(ω_{n}^{k+\frac{n}{2}})$

$=A_1(ω_{n}^{2k+n})+ω_{n}^{k+\frac{n}{2}}A(ω_{n}^{2k+n})$

$=A_1(ω_{n}^{2k})-ω_{n}^{k+\frac{n}{2}}A_2(ω_{n}^{2k})$

$=A_1(ω_{n}^{2k})-ω_{n}^{k}A_2(ω_{n}^{2k})$

我们进行对比：

$$A(ω_{n}^{k})=A_1(ω_{n}^{2k})+ω_{n}^{k}A_2(ω_{n}^{2k})$$

$$A(ω_{n}^{k+\frac{n}{2}})=A_1(ω_{n}^{2k})-ω_{n}^{k}A_2(ω_{n}^{2k})$$

两者之差一个符号。

所以去在求 $A(ω_{n}^{k})$ 的时候可以同时 $O(1)$ 算出 $A(ω_{n}^{k+\frac{n}{2}})$

第一个式子，在 $k$ 取遍 $[0∼\frac{n}{2}−1]$ 时，第二个式子正好取遍 $[\frac{n}{2}∼n−1]$

所以原问题将缩小一半,子问题还是满足原问题的性质。时空复杂度 $O(nlog$ $n)$

代码实现最后说

5.快速傅里叶逆变换 IFFT

我们现在得到的是得到的是点值表示法下的乘积，但我们平时用系数表示的多项式，现在考虑快速地将点值表示法转换成系数表示法。

这个过程叫做傅里叶逆变换

$IFFT$ 是有一个结论的，由于这个结论过于复杂，我没听懂，不会证明。所以就背了算了。。。。

设一个 $A(x)=a_0+a_1x+a_2x^2⋯+a_{n−2}x^{n−2}+a_{n−1}x^{n−1}$

那么

$$a_i=\frac{\sum_{j=0}^{n-1}A(ω_n^j)ω_n^{-ij}}{n}$$

6.IFFT 代码实现部分

快速傅里叶逆变换

typedef struct complex{
    double a,b;
    complex(double x=0,double y=0){a=x,b=y;}
    complex operator+(const complex&t)const{return complex(a+t.a,b+t.b);}
    complex operator-(const complex&t)const{return complex(a-t.a,b-t.b);}
    complex operator*(const complex&t)const{return complex(a*t.a-b*t.b,a*t.b+b*t.a);}
}fastcomplex;

void FFT(fastcomplex a[],int n){
    if(n==1) return ;
    fastcomplex t[N];
    int m=n>>1;
    for(rint i=0;i<m;i++){
	t[i]=a[i<<1];
	t[i+m]=a[i<<1|1];
    }
    FFT(t,m),FFT(t+m,m);
    for(rint i=0,j=m;i<m;i++,j++){
        fastcomplex x(cos(2*pi*i/n),-sin(2*pi*i/n));
        a[i]=t[i]+x*t[j];
        a[j]=t[i]-x*t[j];
    }
    return ;
}

7.FFT 迭代实现——蝴蝶变化

观察下我们 FFT 中分治得到的结果，下面以 $n=8$ 为例

对比我们的原序列和操作之后得到的序列下标，并用二进制表示：

原序列：000,001,010,011,100,101,110,111

后序列：000,100,010,110,001,101,011,111

每个位置分治后的最终位置，为其二进制位翻转后得到的位置

同 IFFT ，这个东西我不会证明，，，，，背了算了

通过这个规律，可以 $O(n)$ 处理出每个数的二进制表示翻转之后的结果

对于一个数 $i$，我们将其二进制表示记做 $f[i]$，考虑如何求 $f[i]$

假设我们已经算出来了 $f[0∼i−1]$

若 $i$ 是偶数，那么 f[i]=f[i/2]+'0'

若 $i$ 是奇数，那么 f[i]=f[i/2]+'1'

二进制翻转代码:

int bit=log_2(n)-1;
for(rint i=0;i<n;i++){
    f[i]=f[i>>1]>>1|(i&1)<<bit;
}

这样，我们就可以把原序列中的每个数，先放在最终的位置上，再一步一步向上合并。就可以解决多项式乘法了P3803

#include<bits/stdc++.h>

#define rint register int
#define endl '\n'
#define complex New_complex

using namespace std;

const int N = 4e6 + 5;
const double pi = 2 * acos(-1);

struct New_complex{
    double a,b;
    complex(double x=0,double y=0){a=x,b=y;}
    complex operator+(const complex&t)const{return complex(a+t.a,b+t.b);}
    complex operator-(const complex&t)const{return complex(a-t.a,b-t.b);}
    complex operator*(const complex&t)const{return complex(a*t.a-b*t.b,a*t.b+b*t.a);}
};

int n,m,k,r[N]; // 存变换之后的下标，即上面的 f[i]
complex a[N],b[N];

int log_2(int x){
    int res=0;
    
    if(x&0xffff0000){
    	res+=16; 
	x>>=16;
    }
    
    if(x&0xff00){
	res+=8;
	x>>=8;
    }
    
    if(x&0xf0){
	res+=4;
	x>>=4;
    }
    
    if(x&0xc){
	res+=2;
	x>>=2;
    }
    
    if(x&2){
	res++;
    }
    
    return res;
}

void swap(complex &a,complex &b){
    static complex t;
    t=a,a=b,b=t;
    return ;
}

void FFT(complex a[],int f){
    for(rint i=0;i<k;i++)
	if(i<r[i])
	    swap(a[i],a[r[i]]);// 将 a 变换成分治后的序列
	
    complex t,w,x,y;// 要用到的四个变量，t 存单位根，w 存单位根的 i 次方，x 和 y 做临时变量
	
    for(rint len=2;len<=k;len<<=1){ // 枚举区间长度
	t=complex(cos(pi/len),f*sin(pi/len));  // 找到单位根
	for(rint l=0;l<k;l+=len){  // 枚举所有区间左端点
	    w=complex(1,0);  // 单位根的 0 次方为 1
            for(rint i=0,j=len>>1;j<len;i++,j++,w=w*t){  // 分治处理
		x=a[l+i];
		y=w*a[l+j];
		a[l+i]=x+y;
		a[l+j]=x-y; // 先用 x 和 y 表示 a[l+i] 和 w*a[l+j],然后再计算,不然会覆盖掉前面的
	    }
	}
    }
}

signed main(){
    cin>>n>>m;
    for(rint i=0;i<=n;i++)
        cin>>a[i].a;
    for(rint i=0;i<=m;i++)
        cin>>b[i].a;
		
    int bit=log_2(n+m);
    k=1<<bit+1; //这里用 k 存补后的多项式的次数
	
    for(rint i=0;i<k;i++)
        r[i]=r[i>>1]>>1|(i&1)<<bit; // 先将 r 数组处理好
	
    FFT(a,1),FFT(b,1);
	
    for(rint i=0;i<k;i++)
        a[i]=a[i]*b[i];
		
    FFT(a,-1);
	
    for(rint i=0;i<=n+m;i++)
        printf("%d ",(int)(a[i].a/k+0.5));
	   
    return 0;
}

8.三两优化

PS:这个优化名字是瞎叫的，它没有学名。。。。

设 $c_k=a_k+ib_k$

那么 $c_k^2=a_k^2−b_k^2+2ia_kb_k$

那么我们就可以把两个多项式存到一个数组里面，然后做两次 FFT

另外，注意到上面重复调用了很多次三角函数，我们可以直接预处理出 $ω_n^i$
，可以减少很多调用三角函数的次数

#include<bits/stdc++.h>

#define rint register int
#define endl '\n'
#define complex New_complex

using namespace std;

const int N = 4e6 + 5;
const double pi = 2 * acos(-1);

struct New_complex{
    double a,b;
    complex(double x=0,double y=0){a=x,b=y;}
    complex operator+(const complex&t)const{return complex(a+t.a,b+t.b);}
    complex operator-(const complex&t)const{return complex(a-t.a,b-t.b);}
    complex operator*(const complex&t)const{return complex(a*t.a-b*t.b,a*t.b+b*t.a);}
};

int n,m,k,r[N]; 
complex a[N],w[N];// w 预处理 omega^k , a 存两个多项式

int log_2(int x){
    int res=0;
    
    if(x&0xffff0000){
    	res+=16; 
	x>>=16;
    }
    
    if(x&0xff00){
	res+=8;
	x>>=8;
    }
    
    if(x&0xf0){
	res+=4;
	x>>=4;
    }
    
    if(x&0xc){
	res+=2;
	x>>=2;
    }
    
    if(x&2){
	res++;
    }
    
    return res;
}

void swap(complex &a,complex &b){
    static complex t;
    t=a,a=b,b=t;
    return ;
}

void FFT(bool type){
    for(rint i=0;i<k;i++)
	if(i<r[i])
	     swap(a[i],a[r[i]]);
			
    static complex x,y;
	
    for(rint len=2;len<=k;len<<=1){
	for(rint l=0;l<k;l+=len){
	    for(rint i=0,j=len>>1;j<len;i++,j++){
		x=a[l+i];
	        y=w[k/len*i];
		if(type)
		    y=y*a[l+j];
		else 
		    y=complex(y.a,-y.b)*a[l+j];	
		a[l+i]=x+y;
		a[l+j]=x-y;
            }	    	
	}		
     }

     return ;
}

int main(){
    cin>>n>>m;
    for(rint i=0;i<=n;i++)
        cin>>a[i].a;
    for(rint i=0;i<=m;i++)
	cin>>a[i].b;
		
    int bit=log_2(n+m);
    k=1<<bit+1;
	
    for(rint i=0;i<k;i++)
	r[i]=r[i>>1]>>1|(i&1)<<bit;
	
    w[0].a=1;
	
    complex t(cos(pi/k),sin(pi/k));
	
    for(rint i=1;i<k;i++)
	w[i]=w[i-1]*t;
	   
    FFT(true);
	
    for(rint i=0;i<k;i++)
        a[i]=a[i]*a[i];
		
    FFT(false);
	
    for(rint i=0;i<=n+m;i++)
	printf("%d ",(int)(a[i].b/k/2.0+0.5));
		
    return 0;
}

好了，就这么多，FFT 仍需多加练习，才能在考场上灵活应用