FFT 小记

目录

• 复数
• 单位复数根
• Poly
• FFT
• Quicklier?
• End
• 参考资料

由于懒，所以没图。

写得时候有点抽风，可能有 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：

$$\begin{gathered} (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+bd{i}^2=(ac-bd)+(bc+ad)i \\ (a+bi)\times (a-bi)=a(a+bi)-bi(a+bi)=a^2+abi-abi+b^2{i}^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 \end{gathered}$$

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

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

复数还有指数形式的定义，即 $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. $\{({\omega }_n^i{)}^2|i\in [0,n-1]\}=\{{\omega }_{n/2}^i|i\in [0,n/2-1]\}$，要求 $2\mid 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 _{k=0}^{n-1}({\omega }_n^d{)}^k=\frac{({\omega }_n^d{)}^n-1}{{\omega }_n^d-1}$，要求 $d\nmid n$

   等比数列求和……

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

   注意那个限制，不满足的话分母会挂（${\omega }_n^d={\omega }_{n/d}=1$），此时和恰为 $n$。

   这玩意儿好像还能导出单位根反演啥的。

上面的目前就够了。

Poly

对于一个多项式，一般有两种表述：

1. $F(x)=\sum _{i=0}^{n-1}{f}_i{x}^i$，这称作系数形式，同时这个多项式是 $n$ 次的。
2. $\{(x_i,y_i)|i\in [0,n-1],F(x_i)=y_i\}$，这称作点值形式。

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

For example：

$$\begin{gathered} F(x)+G(x)=\sum _{i=0}^{n-1}(f_i+g_i)x^i \\ F(x)-G(x)=\sum _{i=0}^{n-1}(f_i-g_i)x^i \\ F(x)\ast G(x)=\sum _{i=0}^{n-1}\sum _{j=0}^{m-1}{f}_i{g}_j{x}^{i+j} \end{gathered}$$

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

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

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

$$\begin{gathered} \{(x_0,y_0)\cdots (x_{n-1},y_{n-1})\}\text{plusmn}\{(x_0,y_0^{\prime })\cdots (x_{n-1},y_{n-1}‘)\}=\{(x_0,y_0\text{plusmn}{y}_0^{\prime })\cdots (x_{n-1},y_{n-1}\text{plusmn}{y}_{n-1}^{\prime })\} \\ \{(x_0,y_0)\cdots (x_{n-1},y_{n-1})\}\times \{(x_0,y_0^{\prime })\cdots (x_{n-1},y_{n-1}‘)\}=\{(x_0,y_0\times y_0^{\prime })\cdots (x_{n-1},y_{n-1}\times y_{n-1}^{\prime })\} \end{gathered}$$

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

但是点值形式一般不常用，而系数形式和点值形式的转化（求值和插值）直接做又同样是 $\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_{1}x+f_2{x}^2+\cdots +f_{n-1}{x}^{n-1}$ 拆做如下两部分：

$$\begin{gathered} F^{\left[0\right]}(x)=f_0+f_{2}x+f_4{x}^2+\cdots +f_{n-2}{x}^{n/2-1} \\ {F}^{\left[1\right]}(x)=f_1+f_{3}x+f_5{x}^2+\cdots +f_{n-1}{x}^{n/2-1} \end{gathered}$$

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

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

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

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

$$\begin{gathered} F^{\left[0\right]}(x):\{({\omega }_{n/2}^0,y_0)\cdots ({\omega }_{n/2}^{n/2-1},y_{n/2-1})\}=\{({\omega }_n^0,y_0)\cdots ({\omega }_n^{n-2},y_{n/2-1})\} \\ {F}^{\left[1\right]}(x):\{({\omega }_{n/2}^0,y_0^{\prime })\cdots ({\omega }_{n/2}^{n/2-1},y_{n/2-1^{\prime }})\}=\{({\omega }_n^0,y_0^{\prime })\cdots ({\omega }_n^{n-2},y_{n/2-1^{\prime }})\} \end{gathered}$$

直接对应位置做乘法，就可以得到 $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 系数矩阵（叫什么范德蒙德矩阵），然后说明它俩互逆，中途要用到求和引理。

简单口胡一下要点：就是你发现最后乘出来的东西就是求和引理，然后 $i\ne j$ 的时候这个位置就为 $0$，而 $i=j$ 的时候这个位置为 $n$（求和引理那里提到过，就是个等差数列），然后你除掉个 $n$ 才能把它构造成单位矩阵。

懒了，矩阵可以去别的地方看（以后也许会敲上来）。

所以结论是：

• DFT 的表述为 $y_i=\sum _{k=0}^{n-1}{a}_k{\omega }^{ki}$
• IDFT 的表述为 $f_i=\frac{1}{n}\sum _{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_();
}

其实这个做法就跟原来的 插值/求值 基本不沾边了，更贴近于 拆分+合并 的形式。

Quicklier?

考虑进行了 $\mathcal{O}(\lg n)$ 次拷贝，非常慢。

第二次拷贝回去很好处理，如下：

		Complex *g=f,*h=f+(n>>=1);
		///...
		for(int i=0;i<n;++i){
			//交换顺序即可
			//这里把 w*h[i] 存下来，因为复数乘法很慢
			Complex t=w*h[i];
			f[i+n]=g[i]-t;
			f[i]=g[i]+t;
			w*=step;
		}
		// n<<=1;
		// for(int i=0;i<n;++i) f[i]=tmp[i];

上面那个要用到叫做蝴蝶变换的技巧，就是划分到最底层的每个数恰好在它二进制翻转的位上。

如：$6(0110)\to 6(0110),7(0111)\to 14(1110)$。

这个规律手玩一下即可，打表也显然。

然后 $\mathcal{O}(n)$ 递推预处理一下即可：

	for(int i=0;i<x;++i) r[i]=(r[i>>1]>>1)|((i&1)?(x>>1):0);
	//即取出最低位放在最高位，其它位翻转然后平移
	for(int i=0;i<x;++i) if(i<r[i]) swap(f[i],f[r[i]]);
	//另外的同理

这时我的代码就从 3s 降到 1.7s 了。

还不够，我们可以把这 FFT 的部分由递归转为非递归，又降了 0.2s 左右。

    void FFT(Complex *f,int n,int type){
		for(int i=0;i<n;++i) if(i<r[i]) swap(f[i],f[r[i]]);
		for(int p=1,len=2;len<=n;p=len,len<<=1){
			Complex step(cos(PI/p),type*sin(PI/p));
			for(int i=0;i<n;i+=len){
				Complex *g=f+i,*h=f+i+p;
				Complex w(1,0);
				for(int k=0;k<p;++k){
					Complex t=w*h[k];
					g[k+p]=g[k]-t;
					g[k]=g[k]+t;
					//注意一下这里下标等等有微调
					w*=step;
				}
			}
		}
    }

完了吗？没有。

zyxawa 提到了一个叫做“三次变两次”的优化，这个优化的意思是只做 DFT 和 IDFT 各一次。

这个恰好也是算导的习题。

怎么实现呢？考虑这个式子：

$$(a+bi)(c+di)=ac-bd+(bc+ad)i$$

容易发现我们算了四次乘法，但如果只关心 imag，我们只需要做两次。

然后我们整出一个类似的复数域上的多项式：

$$\begin{array}{rl}F(x)& =G(x)+H(x)i\\ F(x{)}^2& =(G(x)+H(x)i{)}^2\\ F(x{)}^2& =G(x{)}^2-H(x{)}^2+2G(x)H(x)i\end{array}$$

发现这时 imag 正好对应了 $2G(x)H(x)$。

那么我们构造这么一个 $F(x)$（用 $G(x)$ 和 $H(x)$ 分别充作 real 和 imag），然后对其平方，就可以得到 $G(X)H(X)$ 了。

这个过程只做了两次 FFT，非常高效，让我跑进了 1s。

现在的 FFT 代码如下：

namespace MAIN{
    using QL::IO::lq;
    using QL::Base_Tools::max;
	using QL::Base_Tools::swap;
    constexpr int N=2<<20;
    int n,m;
    using Complex=std::complex<double>;
    const double PI=acos(-1);
    Complex f[N];
	int r[N];
    void FFT(Complex *f,int n,int type){
		for(int i=0;i<n;++i) if(i<r[i]) swap(f[i],f[r[i]]);
		for(int p=1,len=2;len<=n;p=len,len<<=1){
			Complex step(cos(PI/p),type*sin(PI/p));
			for(int i=0;i<n;i+=len){
				Complex *g=f+i,*h=f+i+p;
				Complex w(1,0);
				for(int k=0;k<p;++k){
					Complex t=w*h[k];
					g[k+p]=g[k]-t;
					g[k]=g[k]+t;
					w*=step;
				}
			}
		}
    }
    signed _main_(){
        lq.read(n,m);
        for(int i=0;i<=n;++i) f[i]=Complex(lq.read<int>(),0);
        for(int i=0;i<=m;++i) f[i]=Complex(f[i].real(),lq.read<int>());
        int x=1<<(32-__builtin_clz(max(n,m))+1);
		for(int i=0;i<x;++i) r[i]=(r[i>>1]>>1)|((i&1)?(x>>1):0);
        FFT(f,x,1);
        for(int i=0;i<x;++i) f[i]*=f[i];
        FFT(f,x,-1);
        for(int i=0;i<=(n+m);++i) lq.write(int(f[i].imag()/x/2+0.5),' ');
        return 0;
    }
}
signed main(){
    return MAIN::_main_();
}

End

上面的代码中提到了，复数的计算涉及到实数，存在误差（cmd-block 的博客对减少误差的方法有涉及，可自取）。

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

加之计数题多取模，所以不够实用。

所以在这里挖个坑：NTT。

To be continue...

参考资料

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

OI wiki，这个偏一般。

算法导论，非常给力。