FFT详解

FFT(快速傅里叶变换)

单位根

$$w_n = \cos(\frac{2\Pi}{n}) + \sin(\frac{2\Pi}{n})i$$

可理解为把复坐标轴上的单位圆分成了n份

性质1

$$w_n^k = w_{2n}^{2k}$$

证明

$$w_{2n}^{2k} = \cos(\frac{2\Pi*2k}{2n}) + \sin(\frac{2\Pi*2k}{2n})i = \cos(\frac{2\Pi*k}{n}) + \sin(\frac{2\Pi*k}{n})i = w_n^k$$

这里用到了复数相乘的性质，即模长相乘，夹角相加，因为模长都为1，所以不变，夹角相加。

性质2

$$w_n^{k+\frac{n}{2}} = -\ w_n^k$$

证明:

$$w_n^{\frac{n}{2}}=\cos(\pi) + \sin(\pi)i = -1$$ $$w_n^{k+\frac{n}{2}} =w_n^{\frac{n}{2}}*w_n^k = -\ w_n^k$$

性质3

$${(w_n^i)}^2 ={(w_n^{i+\frac{n}{2}})}^2 = w_n^{2i} = w_{\frac{n}{2}}^i$$

IFFT(系数转点值)

普通的一个多项式方程可以这样表示

$$\begin{pmatrix} x_0&x_1&x_2&\dots&x_{n-1} \end{pmatrix} \begin{pmatrix} a_1\\ a_2\\ a_3\\ \dots\\ a_{n-1} \end{pmatrix} =0$$

我们希望转换成点值，就需要在平面上寻找n个x，然后得到n个y。
就得到这个矩阵

$$\begin{pmatrix} x^1_0&x^1_1&x^1_2&\dots&x^1_{n-1}\\ x^2_0&x^2_1&x^2_2&\dots&x^2_{n-1}\\ x^3_0&x^3_1&x^3_2&\dots&x^3_{n-1}\\ \dots&\dots&\dots&\dots&\dots\\ x^{n-1}_0&x^{n-1}_1&x^{n-1}_2&\dots&x^n_{n-1}\\ \end{pmatrix} \begin{pmatrix} a_1\\ a_2\\ a_3\\ \dots\\ a_{n-1} \end{pmatrix} = \begin{pmatrix} y_1\\ y_2\\ y_3\\ \dots\\ y_{n-1} \end{pmatrix}$$

如图

现在我们需要通过n个蓝色框的x来求出红色框中的y，然后得到n个$(x_i,y_i)$的点，然后进行O(n)点乘即可。
现在就要用到刚刚所说的单位根，把$w_0 - w_{n-1}$代入函数中作为参数$x_0 - x_{n-1}$得到如下的方程组

$$\left\{ \begin{matrix} \begin{aligned} &f(w_n^0) =a_0+a_1*w_n^0+a_2*w_n^0+\dots+a_{n-1}*w_n^0 \\ &f(w_n^1) =a_0+a_1*w_n^1+a_2*w_n^1+\dots+a_{n-1}*w_n^1 \\ &f(w_n^2) =a_0+a_1*w_n^2+a_2*w_n^2+\dots+a_{n-1}*w_n^2 \\ &\vdots\\ &f(w_n^{n-1}) =a_0+a_1*w_n^{n-1}+a_2*w_n^{n-1}+\dots+a_{n-1}*w_n^{n-1} \\ \end{aligned} \end{matrix} \right.$$

然后我们对其中一个进行求解得到

$$\left\{ \begin{matrix} \begin{aligned} f(w_n^i)&=a_0+a_1w_n^i+a_2{(w_n^i)}^2+a_3{(w_n^i)}^3+\cdots++a_{n-1}{(w_n^i)}^{n-1}\\ &=a_0+a_2{(w_n^i)}^2+a_4{(w_n^i)}^4+\cdots+a_{n-2}{(w_n^i)}^{n-2}\\ &\qquad+w_n^i(a_1+a_3{(w_n^i)}^2 + a_5{(w_n^i)}^4 + a_{n-1}{(w_n^i)}^{n-2})\\ &=a_0+a_2(w_{\frac{n}{2}}^i)+a_4{(w_{\frac{n}{2}}^i)}^2+\cdots+a_{n-2}{(w_{\frac{n}{2}}^i)}^{\frac{n-2}{2}}\\ &\qquad+w_n^i(a_1+a_3{(w_{\frac{n}{2}}^i)} + a_5{(w_{\frac{n}{2}}^i)}^2 + a_{n-1}{(w_{\frac{n}{2}}^i)}^{\frac{n-2}{2}}) \end{aligned} \end{matrix} \right.$$

得到这个式子之后，我们就可以利用递归，一层层的求下去，然后log的复杂度求出我们想要的点值的y。
但如果是常规的递归的话我们需要每次传下去一个数组，类似归并的思想，空间上有所限制，然后就有了蝴蝶变换这个东西，观察下面这个图

这个图就是在递归的过程中的位置变换，我们发现，后序列相当于原序列的二进制反转后排序，然后根据这个规律就可以提前得知它应该在的位置，这样就不需要递归了。

rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));

这个用了动态规划,i代表位置，然后对照上图即可理解
然后就可以使用点进行乘法了

IDFT

我们得到点值后需要转化为多项式的系数值，其实过程是一样的

$$\begin{pmatrix} a_1\\ a_2\\ a_3\\ \dots\\ a_{n-1} \end{pmatrix} = \begin{pmatrix} y_1\\ y_2\\ y_3\\ \dots\\ y_{n-1} \end{pmatrix} {\begin{pmatrix} x^1_0&x^1_1&x^1_2&\dots&x^1_{n-1}\\ x^2_0&x^2_1&x^2_2&\dots&x^2_{n-1}\\ x^3_0&x^3_1&x^3_2&\dots&x^3_{n-1}\\ \dots&\dots&\dots&\dots&\dots\\ x^{n-1}_0&x^{n-1}_1&x^{n-1}_2&\dots&x^n_{n-1}\\ \end{pmatrix}}^{-1}$$

也就是把矩阵变为逆矩阵，在图上的话可以想象一开始是顺时针转，现在逆时针转动,实际上得到的矩阵是

所以转回来的时候需要乘以1/n

FFT板子

struct cp
{
    double r,i;
    complex(double _r = 0,double _i = 0)
    {
        r = _r; i = _i;
    }
    complex operator +(const complex &b)
    {
        return complex(r+b.r,i+b.i);
    }
    complex operator -(const complex &b)
    {
        return complex(r-b.r,i-b.i);
    }
    complex operator *(const complex &b)
    {
        return complex(r*b.r-i*b.i,r*b.i+i*b.r);
    }
};
void FFT(cp *s,int inv){
	for(int i=0;i<n;i++){
		if(i<r[i])swap(s[i],s[r[i]]);
	}
	for(int len=2;len<=n;len*=2){
		cp wn = cp(cos(pi*2.0/len),inv*sin(pi*2.0/len));
		for(int st=0;st<n;st+=len){
			cp w = cp(1.0,0.0);
			for(int i=st;i<st+len/2;i++,w*=wn){
				cp x=s[i];//偶次幂 
				cp y=w*s[i+len/2];//奇次幂 
				s[i]=x+y;
				s[i+len/2]=x-y;
			} 
		}
	}
}

洛谷P3803

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define pi acos(-1)
const int N = 2e6+10;
int f[N],g[N],r[N];
struct cp
{
    double r,i;
    cp(double _r = 0,double _i = 0)
    {
        r = _r; i = _i;
    }
    cp operator +(const cp &b)
    {
        return cp(r+b.r,i+b.i);
    }
    cp operator -(const cp &b)
    {
        return cp(r-b.r,i-b.i);
    }
    cp operator *(const cp &b)
    {
        return cp(r*b.r-i*b.i,r*b.i+i*b.r);
    }
}F[N],G[N];
void FFT(cp *s,int n,int inv){
	int bit=0;
    while ((1<<bit)<n)bit++;
	for(int i=0;i<n;i++){
			if(i<r[i])swap(s[i],s[r[i]]);
	}
	for(int len=2;len<=n;len*=2){
		cp wn = cp(cos(pi*2.0/len),inv*sin(pi*2.0/len));
		for(int st=0;st<n;st+=len){
			cp w = cp(1.0,0.0);
			for(int i=st;i<st+len/2;i++,w=w*wn){
				cp x=s[i];//偶次幂 
				cp y=w*s[i+len/2];//奇次幂 
				s[i]=x+y;
				s[i+len/2]=x-y;
			} 
		}
	}
}
int main(){
	int n,m;
	scanf("%d %d",&n,&m);
	for(int i=0;i<=n;i++)scanf("%d",&f[i]),F[i].r=f[i];
	for(int i=0;i<=m;i++)scanf("%d",&g[i]),G[i].r=g[i];
	int len = 1;
	while(len<=n+m)len<<=1;
	int bit=0;
    while ((1<<bit)<len)bit++;
    for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(bit-1));
	FFT(F,len,1);
	FFT(G,len,1);
	for(int i=0;i<=len;i++){
		F[i]=F[i]*G[i];
	}
	FFT(F,len,-1);
	for(int i=0;i<=n+m;i++){
		printf("%d ",(int)(F[i].r/len+0.5));
	}
}