【瞎口胡】快速傅里叶变换 / FFT

快速傅里叶变换（FFT）是一种在 $O(n \log n)$ 时间复杂度内求出两个 $n$ 次多项式乘积的算法。

系数表示法和点值表示法#

对于 $n$ 次多项式 $f(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n$，如果我们知道了每一个 $a_i$，那么这个多项式就唯一确定。于是我们用系数序列 $a=\{a_0,a_1,a_2,\cdots,a_n\}$ 来表示这个多项式，这被称作系数表示法。

而我们也可以取这个多项式在 $n+1$ 个不同的 $x$ 处的取值来表示这个多项式。根据高斯消元法，这 $n+1$ 个 $x$ 以及 $f(x)$ 的取值可以确定这个多项式。这被称作点值表示法。

一个多项式从系数表示法转换为点值表示法的过程，被称作离散傅里叶变换（DFT）。反之则是离散傅里叶逆变换（IDFT）。

FFT 取一些特殊的 $x$，来加速 DFT 和 IDFT 的过程。这些 $x$ 甚至不在实数域，而是一些复数。

特殊的复数#

在复平面上，一个以原点为圆心，半径为 $1$ 的圆被称作单位圆。从实轴（$x$ 轴）正方向开始，逆时针作 $n$ 个向量将单位圆 $n$ 等分，则这些向量与单位圆相交形成的 $n$ 个交点被称作 $n$ 次复根，第一个辐角为正的向量与单位圆的交点被称作 $n$ 次单位复根，记作 $\omega_n$。

根据复数的乘法运算法则「模长相乘，辐角相加」，可知 $n$ 次复根的 $n$ 次方都是 $1$。而所有 $n$ 次复根都可以用 $n$ 次单位复根的幂表示。

我们知道，$2 \pi \operatorname{rad} = 360^\circ$。根据三角函数的知识，$\omega_n$ 的实部即为 $\cos(\dfrac{2\pi}{n})$，虚部为 $\sin(\dfrac{2\pi}{n})$。于是，$\omega_n=\cos(\dfrac{2\pi}{n})+\sin(\dfrac{2\pi}{n})i$。

单位复根有一些性质：

• $\omega_n^n=\omega_0^0=1$
• 对于 $n=2m$，$\omega_n^{k}=-\omega_{n}^{k+m}$
• $\omega_n^k=\omega_{2n}^{2k}$

这些奇妙的性质将会在接下来的环节中充分发挥作用。

快速傅里叶变换 / FFT#

在 FFT 中，我们将 $x=\omega_{n}^{k}(0 \leq k \leq n-1)$ 依次 $f(x)$ 带入求值，便得到了一个 $n-1$ 次多项式的点值表示法。

但这样不够快，我们考虑分治。

对于 $n=2^k(k \in \mathbb Z_+)$，设 $n-1$ 次多项式

$$f(x)=a_0+a_1x+a_2x^2+\cdots + a_{n-1}x^{n-1}$$

$$=(a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+a_5x^5+...a_{n-1}x^{n-1})$$

$$=(a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+a_5x^4+...a_{n-1}x^{n-2})$$

设

$$A(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n-2}{2}}$$

$$B(x)=a_1+a_3x+a_5x^2+...a_{n-1}x^{\frac{n-2}{2}}$$

则有

$$f(x)=A(x^2)+x\times B(x^2)$$

带入 $\omega_n^k(0 \leq k < \dfrac n2)$：

$$A((\omega_n^k)^2)+\omega_n^k\times B((\omega_n^k)^2)$$

$$=A(\omega_n^{2k})+\omega_n^k\times B(\omega_n^{2k})$$

$$=A(\omega_{\frac n2}^{k})+\omega_n^k\times B(\omega_{\frac n2}^{k})$$

带入 $\omega_n^{k+\frac n2}(0 \leq k < \dfrac n2)$：

$$A((\omega_n^{k+\frac n2})^2)+\omega_n^{k+\frac n2}\times B((\omega_n^{k+\frac n2})^2)$$

$$=A(\omega_{n}^{2k+n})-\omega_n^{k}\times B(\omega_{n}^{2k+n})$$

$$=A(\omega_{n}^{2k})-\omega_n^{k}\times B(\omega_{n}^{2k})$$

$$=A(\omega_{\frac n2}^{k})-\omega_n^k\times B(\omega_{\frac n2}^{k})$$

因此当我们求出了 $A(\omega_{\frac n2}^{k}),B(\omega_{\frac n2}^{k})$ 之后，就可以求出 $A(\omega_{n}^{k}),B(\omega_{n}^{k})$。该算法的复杂度为 $T(n)=2T(\dfrac n2)+O(n)=O(n \log n)$。

但是这样要递归，不够快！于是我们考虑优化。我们如果能求出每个系数最后到了哪个位置，就可以不断地合并这些系数，然后求解。

对于 $n=8$，我们来看看每次递归之后，系数是怎么被分类的：

$$\{a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7\}$$

$$\{a_0,a_2,a_4,a_6\},\{a_1,a_3,a_5,a_7\}$$

$$\{a_0,a_4\},\{a_2,a_6\},\{a_1,a_5\},\{a_3,a_7\}$$

$$\{a_0\},\{a_2\},\{a_4\},\{a_6\},\{a_1\},\{a_5\},\{a_3\},\{a_7\}$$

把下标用二进制表示

$$000,010,100,110,001,101,110,111$$

我们观察到，原来在第 $i$（从 $0$ 开始）个位置的系数，最后变到了二进制翻转之后的那个位置。举个例子，$a_3$ 在现在在第 $6$ 个位置。$(3)_{10}=(011)_2$，翻转之后就是 $(110)_2=(6)_{10}$。

记 $r_i$ 为 $i$ 二进制翻转之后的值。我们递推求 $r_i$。已知 $r_0=0$。当 $i>0$ 时，$r_{\left \lfloor\frac i2\right \rfloor}$ 已知。我们考虑它在二进制下与 $r_i$ 的关系：

$$r_{\left \lfloor\frac i2\right \rfloor}=A+0$$

$$r_{i}=B+A$$

其中 $A$ 为 $n-1$ 位 $01$ 串，$B$ 为 $1$ 位 $01$ 串，$+$ 表示连接两个 $01$ 串。

为什么是 $A+0$ 呢？观察到 $\left \lfloor\frac i2\right \rfloor$ 的最高位一定是 $0$，所以翻转之后的最低位一定是 $0$。

同时，观察到 $B$ 是翻转后的最高位，即 $i$ 的最低位。

于是，我们得到了 $r_i(i>0)$ 的递推式：

$$r_i =\left \lfloor \dfrac{r_{\left \lfloor\frac i2\right \rfloor}}{2} \right \rfloor + 2^{k-1} \times (i \mod 2)$$

其中 $k$ 满足 FFT 的长度 $n=2^k$。

这个操作叫做蝴蝶变换，也称位逆序变换。

inline void change(Poly &f,int len){
	for(int i=0;i<=len;++i){
		rev[i]=(rev[i>>1]>>1);
		if(i&1){
			rev[i]|=(len>>1);
		}
	}
	for(int i=0;i<=len;++i){ // 保证每一个对只会被翻转一次
		if(i<rev[i]){
			std::swap(f[i],f[rev[i]]); // 直接将系数扔到最后的位置 然后 FFT
		}
	}
	return;
}

点值相乘#

在求出多项式 $f,g$ 的 $n$ 组点值后，我们要求 $h=f\times g$ 的 $n$ 组点值。显然，$f(i)g(i)=h(i)$。于是我们对求出点值进行对位相乘的操作，就得到 $h$ 的 $n$ 组点值。

下文中，我们会讲解如何通过带入的 $n$ 组特殊点值还原出多项式本身。容易发现，对 $h$ 进行这个过程，还原出的多项式就是 $f\times g$。

快速傅里叶逆变换 / IFFT#

我们通过 FFT 求出了 $n$ 组点值。记它们为 $y_0,y_1,\cdots,y_{n-1}$，其中 $y_i=f(\omega_n^i)$。

设多项式

$$A(x)=y_0 +y_1x+y_2x^2+\cdots+y_{n-1}x^{n-1}$$

则取 $x=\omega_n^{-i}(0 \leq i \leq n-1)$ 对 $A$ 进行 FFT，得到的点值序列就是原来 $a$ 序列的 $n$ 倍。

接下来来证明一下：

$$A(\omega_n^{-k})=\sum \limits_{i=0}^{n-1} y_i(\omega_n^{-k})^i$$

$$=\sum \limits_{i=0}^{n-1} \sum \limits_{j=0}^{n-1}a_j(\omega_n^i)^j(\omega_n^{-k})^i$$

$$=\sum \limits_{j=0}^{n-1}a_j \sum \limits_{i=0}^{n-1} \omega_n^{i(j-k)}$$

$$=\sum \limits_{j=0}^{n-1}a_j \sum \limits_{i=0}^{n-1} {(\omega_n^{j-k})}^{i}$$

后面的和式在 $\omega_n^{j-k}=1$ 时值为 $n$，此时 $n \mid (j-k)$，由这两个和式的范围可得 $j=k$。

在 $\omega_n^{j-k} \neq 1$（此时 $j \neq k$）时，$\sum \limits_{i=0}^{n-1} {(\omega_n^{j-k})}^{i}$ 是一个等比数列求和

$$\sum \limits_{i=0}^{n-1} {(\omega_n^{j-k})}^{i} = \dfrac{(\omega_n^{j-k})^{n} - (\omega_n^{j-k})^0}{\omega_n^{j-k}-1}$$

由单位复根的性质得该式值为 $0$。

则我们继续推导，

$$A(\omega_n^{-k})=\sum \limits_{j=0}^{n-1}a_j \sum \limits_{i=0}^{n-1} {(\omega_n^{j-k})}^{i}$$

$$=a_k \times n$$

于是我们只需要在 FFT 时将单位根变为 $\omega_n^{-1}$，再进行 FFT，就完成了 IFFT。

# include <bits/stdc++.h>

const int N=4000010;

struct Complex{
	double x,y;
	Complex(double _x=0.0,double _y=0.0){
		x=_x,y=_y;
		return;
	}
	Complex operator + (const Complex &b) const{
		return Complex(x+b.x,y+b.y);
	}
	Complex operator - (const Complex &b) const{
		return Complex(x-b.x,y-b.y);
	}
	Complex operator * (const Complex &b) const{
		return Complex(x*b.x-y*b.y,x*b.y+y*b.x);
	}
};

typedef std::vector <Complex> Poly;

const double PI=acos(-1.0);

int rev[N];

Poly F,G;

int n,m;

inline int read(void){
	int res,f=1;
	char c;
	while((c=getchar())<'0'||c>'9')
		if(c=='-')f=-1;
	res=c-48;
	while((c=getchar())>='0'&&c<='9')
		res=res*10+c-48;
	return res*f;
}
inline void change(Poly &f,int len){ // 蝴蝶变换
	for(int i=0;i<=len;++i){
		rev[i]=(rev[i>>1]>>1);
		if(i&1){
			rev[i]|=(len>>1);
		}
	}
	for(int i=0;i<=len;++i){
		if(i<rev[i]){
			std::swap(f[i],f[rev[i]]);
		}
	}
	return;
}
inline void fft(Poly &f,int len,double op){
	change(f,len);
	for(int h=2;h<=len;h<<=1){
		Complex wn(cos(2*PI/h),sin(op*2*PI/h));
		for(int j=0;j<len;j+=h){
			Complex w(1,0);
			for(int k=j;k<j+h/2;++k){
				Complex u=f[k],t=w*f[k+h/2];
				f[k]=u+t,f[k+h/2]=u-t;
				w=w*wn;
			}
		}
	}
	return; 
}
int main(void){
	n=read(),m=read();
	
	int maxlen=1;
	while(maxlen<=n+m){ // 注意是 <= 而非 <
		maxlen<<=1;
	}
	F.resize(maxlen+5),G.resize(maxlen+5);
	for(int i=0;i<=n;++i){
		F[i]=Complex(read(),0);
	}
	for(int i=0;i<=m;++i){
		G[i]=Complex(read(),0);
	}
	fft(F,maxlen,1),fft(G,maxlen,1);
	for(int i=0;i<=maxlen;++i){
		F[i]=F[i]*G[i];
	}
	fft(F,maxlen,-1);
	for(int i=0;i<=n+m;++i){
		printf("%d ",(int)(F[i].x/maxlen+0.5));
	}
	return 0;
}