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 仍需多加练习，才能在考场上灵活应用