FFT 小记

发布时间 2023-08-18 19:14:40作者: LQ636721

由于懒,所以没图。

写得时候有点抽风,可能有 typo,望指出。

复数

复数表述为 \(a+b\times i\),其中 \(i\) 是复数单位 \(\sqrt{-1}\),同时由此可得 \(i^2=-1\)

\(a\) 是实部(下文简称 real),\(b\) 是虚部(简称 imag)。对于一个实数,其 real 等于其值,imag 为 0。

任意一个复数都可一一映射为二维平面(称作复平面,由实轴和虚轴组成)上一个点 \((real,imag)\)

幅角是这个点与原点连线后,与实轴正方向形成的夹角;模长是这个点与原点连线的长度,即到原点的距离。

记幅角为 \(\theta\),模长为 \(u\),则复数还可表为等于 \((u\cos(\theta),iu\sin(\theta))\),证明就是在虚平面上把这个点和原点连起来。

复数的加减法相当于 real 和 imag 的分别相加减,乘法则是看作多项式相乘,除法相当于多项式除法(确信)。

复数共轭是把复数的 imag 取反,并且复数的共轭乘上自身可以得到一个实数。

For example:

\[(a+bi)+(c+di)=(a+c)+(b+d)i \\ (a+bi)-(c+di)=(a-c)+(b-d)i \\ (a+bi)\times(c+di)=c(a+bi)+di(a+bi)=ac+bci+adi+bdi^2=(ac-bd)+(bc+ad)i \\ (a+bi)\times(a-bi)=a(a+bi)-bi(a+bi)=a^2+abi-abi+b^2i^2=a^2-b^2 \\ \frac{a+bi}{c+di}=\frac{(a+bi)\times(c-di)}{(c+di)\times(c-di)}=\frac{(ac+bd)+(bd-ad)i}{c^2-d^2}=\frac{ac+bd}{c^2-d^2}+\frac{bd-ad}{c^2-d^2}i \]

最后一个除法用了共轭进行化简,把分母有理化了。

然后这套运算还有几何意义:加减法是点的加减,乘法是模长相乘,幅角相加,共轭是关于实轴对称。

复数还有指数形式的定义,即 \(e^{iu}=\cos(u)+i\sin(u)\)(是不是很眼熟,这个似乎叫欧拉公式),其中 \(e\) 即为自然对数。

具体为什么……我也不清楚。

单位复数根

\(n\) 次单位根(\(\omega_n\) )是满足 \(\omega_n^n=1\) 的复数。

根据我们滴数学啊,\(n\) 次单位根一共有 \(n\) 个,并且为 \(\omega_n\)\(0 \sim n-1\) 次幂(或者说 \(1 \sim n\) 次幂),称 \(\omega_n\) 为主 \(n\) 次单位根。

一个性质是:在复平面上以原点为圆心,画一个半径为 \(1\) 的圆(此时最好工整地在草稿本上画一个圆),这个圆上的点模长均为 \(1\);然后从与实轴正方向的交点开始,把圆周等分为 \(n\) 份,恰好可以得到 \(n\) 个单位根,\(\omega_n^0\) 恰为这个圆和实轴正方向的交点。

证明则考虑当模长小于 \(1\) 的复数幂次越高模长就越小,故不可能为单位根,大于 \(1\) 同理。

啊你说为什么是等分的?你发现模长相等,所以考虑幅角,若把圆周等分为 \(n\) 分,任意一份重复 \(n\) 次恰好就是一整个圆,而整个圆的模长又为 \(1\),所以显然符合(注意,我们这里任意一份对应的复数模长本来就为 \(1\),所以做乘法的时候仍然不考虑)。

而且我们可以得到 \(\omega_n^k=e^{2\pi ik/n}\),这个就是找交点然后揉进式子里。

那么这个 \(\omega_n\),有什么性质呢?

  1. \(\omega_{dn}^{dk}=\omega_n^k\)

    这个比较显然,也符合直觉。

    从代数上看,就是 \(e^{2\pi idk/dn}=e^{2\pi ik/n}\)

    几何上看,你从实轴正方向开始把连续 \(d\) 个部分合成,毫无问题。

    这个会被称作消去引理。

  2. \(\omega_n^{n/2}=\omega_2=1\)

    同样从原来的式子上看,\(e^{2\pi i\frac n 2/n}=e^{2\pi i/2}=\omega_2=-1\)

    你从几何上看,就是 \(\omega_n^0\) 关于虚轴的对称,两者都在实轴上且互为相反数。

  3. \(\left\{(\omega_n^i)^2\vert i\in \left[0,n-1\right]\right\}=\left\{\omega_{n/2}^i\vert i\in \left[0,n/2-1\right]\right\}\),要求 \(2\parallel n\)

    这个东西看这有点高级,它其实是想表述:\(\omega_n^{k}=-\omega_n^{k+n/2}\),同时 \((\omega_n^k)^2=\omega_{n/2}^k\)

    后一个可以由消去引理得到,前一个又是什么理由呢?就是上一条 \(\omega_n^{n/2}=-1\),然后带了系数。

    注意这个家伙非常滴重要,我们的 FFT 要用。

    这个会被称作折半引理。

  4. \(\sum\limits_{k=0}^{n-1}(\omega_n^d)^k=\frac{(\omega_n^d)^n-1}{\omega_n^d-1}\),要求 \(d\nparallel n\)

    等比数列求和……

    但是同样重要,我们的 FFT 也要用。

    注意那个限制,不满足的话分母会挂(\(\omega_n^d=\omega_{n/d}=1\)),此时和恰为 \(n-1\)

上面的目前就够了。

Poly

对于一个多项式,一般有两种表述:

  1. \(F(x)=\sum\limits_{i=0}^{n-1} f_i x^i\),这称作系数形式,同时这个多项式是 \(n\) 次的。
  2. \(\left\{\left(x_i,y_i\right)\vert i \in \left[0,n-1\right],F(x_i)=y_i \right\}\),这称作点值形式。

同时我们可以定义多项式的运算,就和竖式的运算类似(为了方便,后文用不同的大小写来区分多项式和多项式的各个系数)。

For example:

\[F(x)+G(x)=\sum\limits_{i=0}^{n-1} (f_i+g_i) x^i \\ F(x)-G(x)=\sum\limits_{i=0}^{n-1} (f_i-g_i) x^i \\ F(x)*G(x)=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{m-1} f_i g_j x^{i+j} \]

除法暂时不管,涉及到多项式求逆。

然后发现系数形式的多项式做加减法是 \(\mathcal{O}(n)\) 的,但乘法是 \(\mathcal{O}(n^2)\)

不过点值形式就会好做多:

\[\left\{(x_0,y_0)\cdots(x_{n-1},y_{n-1}) \right\}\plusmn\left\{(x_0,y_0’)\cdots(x_{n-1},y_{n-1}‘) \right\}=\left\{(x_0,y_0\plusmn y_0')\cdots(x_{n-1},y_{n-1}\plusmn y_{n-1}') \right\} \\ \left\{(x_0,y_0)\cdots(x_{n-1},y_{n-1}) \right\}\times \left\{(x_0,y_0’)\cdots(x_{n-1},y_{n-1}‘) \right\}=\left\{(x_0,y_0\times y_0')\cdots(x_{n-1},y_{n-1}\times y_{n-1}') \right\} \]

不过对于除法不成立(仍需用多项式求逆),可见《算法导论》习题。

但是点值形式一般不常用,而系数形式和点值形式的转化(求值和插值)直接做又同样是 \(\mathcal{O}(n^2)\)(插值要做到 \(\mathcal{O}(n^2)\) 还需套个拉格朗日插值)。

要快速进行多项式乘法,需要快速求值与插值,但求值和插值应当怎么优化呢?

FFT

正片开始。

傅里叶的做法是:带入 单位根 进行求值/插值。

有什么优势?

我们可以对于任意的 \(n\) 次多项式恰好带入 \(n\)\(n\) 次单位根,同时由折半引理, \(n\)\(n\) 次单位根的平方构成的集合恰好等价于 \(n/2\)\(n/2\) 次单位根勾成的集合,且同一对里面两个数互为相反数。

只要我们构造出两个部分,恰好取到这 \(n/2\) 个平方,就可以保证问题规模搞好缩小到一半。

为了方便,我们假定 \(n\) 都是 \(2^k,k\in \mathbf{N}\),不足的部分可以补 \(f_i=0\) 替代。

构造是容易的,我们可以把多项式 \(F(x)=f_0+f_1x+f_2x^2+\cdots+f_{n-1}x^{n-1}\) 拆做如下两部分:

\[F^{\left[0\right]}(x)=f_0+f_2x+f_4x^2+\cdots+f_{n-2}x^{n/2-1} \\ F^{\left[1\right]}(x)=f_1+f_3x+f_5x^2+\cdots+f_{n-1}x^{n/2-1} \\ \]

即按照下标的奇偶性划分成两个 \(n/2\) 次的多项式。

然后简单观察可以发现如下性质:

\[F(x)=F^{\left[0\right]}(x^2)+xF^{\left[1\right]}(x^2) \]

现在我们将右边的两个多项式计算出来,点值形式大概长成这个样子。

\[F^{\left[0\right]}(x):\left\{(\omega_{n/2}^0,y_0)\cdots(\omega_{n/2}^{n/2-1},y_{n/2-1}) \right\}=\left\{(\omega_{n}^0,y_0)\cdots(\omega_{n}^{n-2},y_{n/2-1}) \right\} \\ F^{\left[1\right]}(x):\left\{(\omega_{n/2}^0,y_0')\cdots(\omega_{n/2}^{n/2-1},y_{n/2-1'}) \right\}=\left\{(\omega_{n}^0,y_0')\cdots(\omega_{n}^{n-2},y_{n/2-1'}) \right\} \\ \]

直接对应位置做乘法,就可以得到 \(F(x)\)

注意到 \(\omega_n^{k+n/2}=\omega_n^k\),所以实现时循环到 \(n/2\),另一部分只需把 \(F^{\left[1\right]}(x)\) 的系数取反即可。

上面的过程对多项式进行了求值,被称为 DFT(插值对应的叫做 IDFT),贴一个递归的代码实现:

	//Complex 可以手写,也可以借用 STL 中的 std::complex<double>
	void DFT(Complex *f,int n){
		if(n==1) return;//边界,此时本身就是点值形式
		for(int i=0;i<n;++i) tmp[i]=f[i];//临时数组记录f
		for(int i=0;i<n;++i)
			f[(i>>1)+((i&1)?(n>>1):0)]=tmp[i];//划分
		Complex *g=f,*h=f+(n>>1);
		DFT(g,n>>1),DFT(h,n>>1);//递归计算
		Complex w(1,0),step(cos(2*PI/n),sin(2*PI/n));//单位根公式
		for(int i=0;i<(n>>1);++i){
			tmp[i]=g[i]+w*h[i];
			tmp[i+(n>>1)]=g[i]-w*h[i];//上面说的
			w*=step;//维护单位根
		}
		for(int i=0;i<n;++i) f[i]=tmp[i];
	}

完事之后用点值形式对应相乘即可(返回的这个复数数组就是点值形式的 \(y\))。

然后还需要把点值形式变回系数形式,这时需要运用 IDFT。

你需要写个 too huge 的矩阵出来,然后说明这玩意儿是 DFT,然后又证明 IDFT 的系数是什么什么,中途要用到求和引理。

所以懒了,结论是:

  • DFT 的表述为 \(y_i=\sum\limits_{k=0}^{n-1}a_k\omega^{ki}\)
  • IDFT 的表述为 \(f_i=\frac{1}{n}\sum\limits_{k=0}^{n-1}y_k\omega_n^{-ki}\)

你发现这俩出奇地相似,所以可以套上面 DFT 的做法,写出这么一个函数:

	void IDFT(Complex *f,int n){
		if(n==1) return;
		for(int i=0;i<n;++i) tmp[i]=f[i];
		for(int i=0;i<n;++i)
			f[(i>>1)+(i&1?(n>>1):0)]=tmp[i];
		Complex *g=f,*h=f+(n>>1);
		IDFT(f,n>>1),IDFT(h,n>>1);
		Complex w(1,0),step(cos(2*PI/n),-sin(2*PI/n));//只有这一处不一样……
		for(int i=0;i<(n>>1);++i){
			tmp[i]=g[i]+w*h[i];
			tmp[i+(n>>1)]=g[i]-w*h[i];
			w*=step;
		}
		for(int i=0;i<n;++i) f[i]=tmp[i];
	}
	//主函数内:
	for(int i=0;i<n;++i) f[i]/=Complex(n,0);

然后这俩实在是太像了,为了降低代码长度,可以把他们合并了。

于是板子题的 AC 代码如下:

//主要部分
namespace MAIN{
	using QL::IO::lq;
	using QL::Base_Tools::max;
	constexpr int N=2<<20;
	int n,m;
	using Complex=std::complex<double>;
	const double PI=acos(-1);
	Complex f[N],g[N],h[N],tmp[N];
	void FFT(Complex *f,int n,int type){
		if(n==1) return;
		for(int i=0;i<n;++i) tmp[i]=f[i];
		for(int i=0;i<n;++i)
			f[(i>>1)+((i&1)?(n>>1):0)]=tmp[i];
		Complex *g=f,*h=f+(n>>=1);
		FFT(g,n,type),FFT(h,n,type);
		Complex w(1,0),step(cos(PI/n),type*sin(PI/n));
		for(int i=0;i<n;++i){
			tmp[i]=g[i]+w*h[i];
			tmp[i+n]=g[i]-w*h[i];
			w*=step;
		}
		n<<=1;
		for(int i=0;i<n;++i) f[i]=tmp[i];
	}
	signed _main_(){
		lq.read(n,m);
		for(int i=0,x;i<=n;++i) lq.read(x),f[i]=Complex(x,0);
		for(int i=0,x;i<=m;++i) lq.read(x),g[i]=Complex(x,0);
		int x=1<<(32-__builtin_clz(max(n,m))+1);
		FFT(f,x,1),FFT(g,x,1);
		for(int i=0;i<x;++i) h[i]=f[i]*g[i];
		FFT(h,x,-1);
		//注意浮点数有向下舍去的误差,这里要+0.5
		for(int i=0;i<=(n+m);++i) lq.write(int((h[i]/Complex(x,0)).real()+0.5),' ');
		return 0;
	}
}
signed main(){
	return MAIN::_main_();
}

End

上面的代码中提到了,复数的计算有误差。

而实际上,不仅有误差,它还跑得相当之慢……

所以在这里挖个坑:NTT。

To be continue...

参考资料

command-block 的博客,大佬写得很详细。

OI wiki,这个偏一般。

算法导论,非常给力。