[学习笔记] 快速傅里叶变换 (FFT)

0. 前置芝士

0.1. 复数

0.1.1. 定义

形如 (z=a+bi)(a,b) 均为实数)的数称为复数,其中 (a) 称为实部,(b) 称为虚部。其中 (i) 是虚数单位,满足 (i^2=-1)

另外,可以将其想象成平面直角坐标系上的点/向量 ((a,b))

定义 ( heta) 为复数 (z) 的辐角(由上可想象),(r=sqrt{a^2+b^2}) 是模长。

那么 (z) 可以表示为 (r imes (cos heta+isin heta))(显然 (r imes cos heta) 就是 (a))。

还可以用欧拉公式((e^{xi}=cos x+isin x))表示:(z=r imes e^{xi})

不过下文就会省略 (r)(因为是单位圆)。

0.1.2. 计算

相加/减就是实部与虚部相加/减。

相乘就是 (z_1z_2=r_1 imes r_2 imes e^{(x+y)i}),即 模长相乘,幅角相加

这里再提一嘴,如果想让复数转动一定角度,我们只用乘上此角度对应的 单位复数 以保障模长不变。

0.2. 多项式表示法

0.2.1. 系数多项式

[f(x)=sum_{i=0}^{n-1}a_ix^i ]

0.2.2. 点值多项式

选取 (n)(x_i) 代入 (f(x)) 得到 (y_i),每一对表示为 ((x_i,y_i))。因为 (n) 个未知系数 (a_i) 需要 (n) 个方程即可解。

0.3. 单位复根

0.3.1. 定义

由欧拉公式可知,(e^{xi}) 对应的复数的辐角为 (x)(弧度制),故 (e^{2pi i}) 对应 (2pi)

假设我们需要把 (2pi) 分成 (n) 等份(如下图,(n=8)),设一份对应着 (omega_n)

我们需要将 (omega_n) 的辐角加 (n-1) 个它本身才是 (2pi)

想到什么了吗?相乘即是模长相乘,幅角相加

所以有:

[omega_n^n=e^{2pi i} ]

[omega_n=e^{frac{2pi i}{n}} ]

0.3.2. 性质

  1. (omega_n^j=omega_{n imes k}^{j imes k})。这就是 (e) 的指数分子分母同时乘上 (k),所以相等。
  2. (omega_n^j=-omega_n^{j+frac{n}{2}})。容易发现就是辐角加上 (pi)
  3. (omega_n^n=1)。显然 (e^{2pi i}=1 imes (cos 2pi +isin 2pi)=1)

1. 正文

1.1. (mathtt{FFT}) 可以干什么

求两个多项式 (f,g) 相乘:((a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}) imes (b_0+b_1x+b_2x^2+...+b_{m-1}x^{m-1}))

显然我们需要计算 (x)(0)(n+m-2) 次方项的系数。

虽然我们得到和最终要求的都是系数多项式,可我们发现点值多项式合并只需要 (y) 相乘,整体复杂度是 (mathcal O(n)) 的。

如果我们能快速地转换系数多项式与点值多项式,就能快速求解原问题。

朴素是 (mathcal O(n^2)) 的,而 (mathtt{FFT}) 可以优化到 (mathcal O(nlog n))

1.2. 系数多项式 ( ightarrow) 点值多项式

将系数多项式按下标的奇偶分成两个函数(认为 (n)(2) 的正整数幂,如果不够可以补 (0) 系数):

[f_1(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{frac{n}{2}-1} ]

[f_2(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{frac{n}{2}-1} ]

那么有:

[f(x)=f_1(x^2)+x imes f_2(x^2) ]

傅里叶说:我们选择 (omega_n^i) 为点值多项式代入的 (x_i)

那么对于 (i<frac{n}{2}) 就有:

[f(omega_n^i)=f_1(omega_n^{2i})+omega_n^i imes f_2(omega_n^{2i}) ]

[=f_1(omega_{n/2}^{i})+omega_n^i imes f_2(omega_{n/2}^{i}) ]

[f(omega_n^{i+frac{n}{2}})=f_1(omega_n^{2i} imes omega_n^n)+omega_n^{i+frac{n}{2}} imes f_2(omega_n^{2i} imes omega_n^n) ]

[=f_1(omega_{n/2}^{i})-omega_n^{i} imes f_2(omega_{n/2}^{i}) ]

同时 (f_1,f_2) 又是新的 (f),这就是一个迭代的过程。

1.3. 点值多项式 ( ightarrow) 系数多项式

得到的 (y_i=sum_{j=0}^{n-1}a_j(omega_n^i)^j)

傅里叶说:我们选择 (omega_n^{-i}) 为点值多项式代入的 (x_i)(y_i) 为新的系数,再做一次。

[z_i=sum_{j=0}^{n-1}y_j(omega_n^{-i})^j ]

[=sum_{j=0}^{n-1}a_jsum_{k=0}^{n-1}(omega_n^{j-i})^k ]

后面那一坨是等比数列,如果公比不为 (1) 时易得为 (0),当公比为 (1)(i=j) 时为 (n)

所以有:

[z_i=a_i imes n ]

我们就转回去了。

1.4. 代码

说了这么多,终于来点实在的了。

由于我们递归地划分,原多项式和末多项式相同下标有奇怪的联系:二进制是相反的。

为了模拟递归,我们可以先求出末多项式再倒回来算。具体算法和解析戳这

具体看代码吧。

#include <cstdio>

#define rep(i,_l,_r) for(register signed i=(_l),_end=(_r);i<=_end;++i)
#define fep(i,_l,_r) for(register signed i=(_l),_end=(_r);i>=_end;--i)
#define erep(i,u) for(signed i=head[u],v=to[i];i;i=nxt[i],v=to[i])
#define efep(i,u) for(signed i=Head[u],v=to[i];i;i=nxt[i],v=to[i])
#define print(x,y) write(x),putchar(y)

template <class T> inline T read(const T sample) {
    T x=0; int f=1; char s;
    while((s=getchar())>'9'||s<'0') if(s=='-') f=-1;
    while(s>='0'&&s<='9') x=(x<<1)+(x<<3)+(s^48),s=getchar();
    return x*f;
}
template <class T> inline void write(const T x) {
    if(x<0) return (void) (putchar('-'),write(-x));
    if(x>9) write(x/10);
    putchar(x%10^48);
}
template <class T> inline T Max(const T x,const T y) {if(x>y) return x; return y;}
template <class T> inline T Min(const T x,const T y) {if(x<y) return x; return y;}
template <class T> inline T fab(const T x) {return x>0?x:-x;}
template <class T> inline T gcd(const T x,const T y) {return y?gcd(y,x%y):x;}
template <class T> inline T lcm(const T x,const T y) {return x/gcd(x,y)*y;}
template <class T> inline T Swap(T &x,T &y) {x^=y^=x^=y;}

#include <cmath>
#include <iostream>
using namespace std;

const int maxn=3e6+5;
const double Pi=acos(-1.0);

int n,m,lim=1,bit,rev[maxn];

struct cp {
	double x,y;
	cp operator + (const cp t) {
		return (cp) {x+t.x,y+t.y};
	}
	cp operator - (const cp t) {
		return (cp) {x-t.x,y-t.y};
	}
	cp operator * (const cp t) {
		return (cp) {x*t.x-y*t.y,y*t.x+x*t.y};
	}
	void Give(const double X,const double Y) {
		x=X,y=Y;
	}
} a[maxn],b[maxn],wn,w,tmp;

void init() {
	while(lim<=n+m) lim<<=1,++bit;
	// n+m 是最高项指数,由于 n 是 [0,n-1],所以这里是 <=
	rep(i,0,lim-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<bit-1);
}

void FFT(cp *t,int op) {
	rep(i,0,lim-1) if(i<rev[i]) swap(t[i],t[rev[i]]);
	// 求末多项式,要求 i<rev[i] 是为避免重复交换
	for(int mid=1;mid<lim;mid<<=1) {
		// mid 是枚举所在递归层 x 对应 f 的长度,其中递归层 x 已经被计算,现在要计算 x 上一层(即从 f_1,f_2 贡献到 f)
		wn.Give(cos(Pi/mid),sin(Pi/mid)*op);
		// wn 是 w_{mid*2},就是 x 上一层的单位复根
		// *op 是求 wn 的共轭复数,即 wn^{-1}(转回去需要),显然把这两个用欧拉公式可以得到两个互为相反数的辐角,所以横坐标不变,纵坐标相反
		for(int i=0;i<lim;i+=(mid<<1)) {
			// i 是找 x 的上一层的每一段的开头
			w.Give(1,0);
			// w 用来模拟单位复根的几次幂,初始化幅角为 0(cos(0)=1,sin(0)=0)
			for(int j=0;j<mid;++j,w=w*wn) {
				// x 这一层的每一段用 j 来遍历
				/*
                0 1 2 3 4 5 6 7
                0 2 4 6|1 3 5 7  y
                0 4|2 6|1 5|3 7  x
                0|4|2|6|1|5|3|7
                
                这是一个递归
                y 层的 0,4 根据公式由 x 层的 0,2 贡献
                */
				tmp=w*t[i+j+mid];
				t[i+j+mid]=t[i+j]-tmp,t[i+j]=t[i+j]+tmp;
			}
		}
	}
}

int main() {
	n=read(9),m=read(9);
	rep(i,0,n) scanf("%lf",&a[i].x);
	rep(i,0,m) scanf("%lf",&b[i].x);
	init(); 
	FFT(a,1),FFT(b,1);
	rep(i,0,lim) a[i]=a[i]*b[i];
	FFT(a,-1);
	rep(i,0,n+m) printf("%d ",(int)(a[i].x/lim+0.5));
	return 0;
}

2. 拓展

拆系数 (mathtt{FFT}),欢迎资瓷!

原文地址:https://www.cnblogs.com/AWhiteWall/p/14363886.html