FFT 小记

写在前面

$Q:$ 为什么会~~心血来潮~~去学 FFT

$A:$ 当本蒟蒻还在努力消化凸包时：。所以本蒟蒻也来看一下

等等 ~~摸头警告~~ 。思维已经废了

About FFT

FFT( $Fast\ Fourier\ Transformation$ )

中文名：快速傅里叶变换

~~Fast Fast TLE~~

作用：在 $O(n\log n)$ 内求出多项式卷积

多项式

一个形如 $A(x)=a_0+a_1x+\cdot\cdot\cdot+a_{n-2}x^{n-2}+a_{n-1}x^{n-1}=\sum_{i=0}^{n-1}a_ix^i$ 的柿子，称其为多项式

系数表示法

将 $n-1$ 次多项式看成一个 $n$ 维向量 $\vec a=(a_0,a_1,\cdot\cdot\cdot,a_{n-1})$ 即为多项式的系数表示
点值表示法

$n-1$ 次多项式 $A(x)$ 将 $n$ 个不同的 x 代入可得到 n 个点，能唯一确定多项式 $A(x)$
多项式乘法

$A(x)=\sum_{i=0}^{n-1}a_ix^i,B(x)=\sum_{i=0}^{n-1}b_ix^i$ 则 $C(x)=A(x)*B(x)=\sum_{i=0}^{2n-2}\sum_{j+k=i}a_jb_kx_i$
卷积

两个向量 $\vec a=(a_0,a_1,\cdot\cdot\cdot,a_{n-1}),\vec b=(b_0,b_1,\cdot\cdot\cdot,b_{n-1})$

有卷积 $\vec a \otimes \vec b=c(c_0,c_1,\cdot\cdot\cdot,c_{2n-2})$ ，其中 $c_k=\sum_{i+j=k}a_ib_j$

系数表示法时计算是 $O(n^2)$ 的

但是对于两个点值表达式的多项式，可以 $O(n)$ 的计算出多项式乘积

便得出了 FFT 的三个步骤

系数表示法转为点值表示法， $DFT,O(n\log n)$
点值表示法相乘， $O(n)$
点值表示法转为系数表示法， $IDFT,O(n\log n)$

复数 Complex

复数由实部和虚部构成，可用二元组 $(a,b)$ 表示复数 $a+bi,i^2=-1$ 可将其理解为一个点或向量

加法

$(a_1,b_1)+(a_2,b_2)=(a_1+a_2,b_1+b_2)$
减法

$(a_1,b_1)-(a_2,b_2)=(a_1-a_2,b_1-b_2)$
乘法

$(a+bi)(c+di)=(ac-bd)+(bc+ad)i=(ac-bd,bc+ad)$
除法
- 共轭复数：实部相同，虚部互为相反数的两个复数， $z(a,b)=a+bi$ 的共轭复数为 $\bar z(a,-b)=(a-bi)$
  
  ~~有趣的~~性质：对于一个 $z=(a,b)$ ， $|z|=|\bar z|,z\cdot\bar z=a^2+b^2$
$\dfrac{a+bi}{c+di}=\dfrac{(a+bi)(c-di)}{(c+di)(c-di)}=\dfrac{ac+bd}{c^2+d^2}+\dfrac{bc-ad}{c^2+d^2}i$
单位根

将单位圆（以原点为圆心，半径为 1）n 等分，得到 n 个模长为 1 的复数

将点从 0 开始标号，设第 0 个点为 $\omega_n^0$ 以此类推

以 (1,0) 为起点，则 $\omega_n^i=(\cos(\dfrac{i}{n}2\pi),\sin(\dfrac{i}{n}2\pi))$ （用弧度制）

把这些复数称为 n 次单位根

性质：$\omega_n^n=1,\omega_n^k=\omega_n^{2k},\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^k$

DFT

将系数表达式转化为点值表达式，就 DFT 来说，它分治地来求当 $x=\omega_n^k$ 的时候 $F(x)$ 的值

使 $n=2^m$ ，不够位就补，不会对答案有影响

对于多项式 $A(x)=\sum_{i=0}^{n-1}a_ix^i$ ，将其按下标奇偶性分类

$A(x)=(a_0+a_2x^2+\cdot\cdot\cdot+a_{n-2}x^{n-2})+(a_1x+a_3x^3+\cdot\cdot\cdot+a_{n-1}x^{n-1})$

现在设 $A_1(x)=(a_0+a_2x+\cdot\cdot\cdot+a_{n-2}x^{\frac{n-2}{2}}),A_2(x)=(a_1+a_3x+\cdot\cdot\cdot+a_{n-1}x{\frac{n-2}{2}})$

则 $A(x)=A_1(x^2)+xA_2(x^2)$

对于 $k<\dfrac{n}{2}$ 有 $A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})\\ =A_1(\omega_{\frac{n}{2}}^k)+\omega_n^kA_2(\omega_{\frac{n}{2}}^k)$

对于 $k+\dfrac{n}{2}$ 有 $A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n})\\\ =A_1(\omega_{\frac{n}{2}}^k*\omega_n^n)+\omega_n^k*\omega_n^{\frac{n}{2}}A_2(\omega_{\frac{n}{2}}^k*\omega_n^n)$

因为 $\omega_n^n=(1,0),\omega_n^{\frac{n}{2}}=(-1,0)$ 所以 $A(\omega_n^{k+\frac{n}{2}})=A_1(w_{\frac{n}{2}}^k)-w_n^kA_2(w_{\frac{n}{2}}^k)$

于是问题被分解成了更小的子问题，递归求解即可

时间复杂度 $O(n\log n)$

IDFT

将点值表达式转回为系数表达式

以下内容摘自 OI-Wiki

已知 $y_i=F(\omega_n^i),i\in[0,n)$ 求 $\{a_0,a_1,\cdot\cdot\cdot,a_{n-1}\}$ 。构造式是：$A(x)=\sum_{i=0}^{n-1}y_ix^i$

相当于把 $\{y_0,y_1,\cdot\cdot\cdot,y_{n-1}\}$ 当作多项式 $A$ 的系数表示法，设 $b_i=\omega_n^{-i}$

则多项式 $A$ 在 $x=b_0,b_1,\cdots,b_{n-1}$ 处的点值表示法为 $\{A(b_0),A(b_1),\cdots\,A(b_{n-1})\}$

对 $A(x)$ 的定义式做一下变换，可以将 $A(b_k)$ 表示为

\[\begin{aligned} A(b_k) &= \sum_{i=0}^{n-1}F(\omega_n^i)\omega_n^{-ik} =\sum_{i=0}^{n-1}\omega_n^{-ik}\sum_{j=0}^{n-1}a_j(\omega_n^i)^j \\ &= \sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j\omega_n^{i(j-k)} =\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i \end{aligned} \]

记 $S(\omega_n^a)=\sum_{i=0}^{n-1}(\omega_n^a)^i$。当 $a=0$ 时，$S(\omega_n^a)=n$；当 $a\ne0$ 时，错位相减

\[\begin{aligned} S(\omega_n^a) &= \sum_{i=0}^{n-1}(\omega_n^a)^i \\ \omega_n^a S(\omega_n^a) &= \sum_{i=1}^n (\omega_n^a)^i\\ S(\omega_n^a) &= \dfrac{(\omega_n^a)^n-(\omega_n^a)^0}{\omega_n^a-1}=0 \end{aligned} \]

也就是说 $S(\omega_n^a)=\left\{\begin{aligned}n,a=0\\0,a\ne0\end{aligned}\right.$，那么代回原式$A(b_k)=\sum_{j=0}^{n-1}a_jS(\omega_n^{j-k})=a_k\cdot n$

也就是说给定 $b_i=\omega_n^{-i}，则 $ $A$ 的点值表示法为

$\{(b_0,A(b_0)),(b_1,A(b_1)),\cdots,(b_{n-1},A(b_{n-1}))\}$

$=\{(b_0,a_0\cdot n),(b_1,a_1\cdot n,\cdots,(b_{n-1},a_{n-1}\cdot n\}$

综上所述，我们取单位根为其倒数，对 $\{y_0,y_1,\cdots,y_{n-1}\}$ 跑一遍 FFT ，然后除以 n 即可得到 $F(x)$ 的系数表示法。

迭代实现

然而，如果直接用递归来打，本蒟蒻会因为自带的 $O(\infin)$ 大常数被卡死。

可以研究以下序列变换过程

0 1 2 3 4 5 6 7

0 2 4 6 | 1 3 5 7

0 4 | 2 6 | 1 5 | 3 7

这时候发现：最终序列为原顺序二进制下的反转，直接模拟从下往上合并

可以用递推来预处理

$\lfloor\dfrac{x}{2}\rfloor$ 的翻转值是已知的，而这个值右移一位就是 x 除了二进制个位的翻转值

若个位是 0 ，反转后最高位就是 0，否则最高位就是 1

Code

Luogu P3803

#include<bits/stdc++.h>
using namespace std;
const int N=2100000;
const double PI=acos(-1.0);
struct Comp {
	double x,y;
	inline Comp(double p=0.0,double q=0.0):x(p),y(q) { }
	inline Comp operator +(Comp o) { return Comp(x+o.x,y+o.y); }
	inline Comp operator -(Comp o) { return Comp(x-o.x,y-o.y); }
	inline Comp operator *(Comp o) { return Comp(x*o.x-y*o.y,x*o.y+y*o.x); }
}x[N],y[N];
int r[N];
inline void FFT(Comp *a,int n,int on) {
	register int i,j,k;
	register Comp wn,w,x,y;
	for(i=1;i<=n;i++)
		if(i<r[i])swap(a[i],a[r[i]]);
	for(i=1;i<n;i<<=1) {
		wn=Comp(cos(PI/i),on*sin(PI/i));
		for(j=0;j<n;j+=(i<<1)) {
			w=Comp(1,0);
			for(k=0;k<i;k++,w=w*wn) {
				x=a[j+k],y=w*a[i+j+k];
				a[j+k]=x+y,a[i+j+k]=x-y;
			}
		}
	}
	if(on==-1)
	for(i=0;i<=n;i++)
	a[i].x=floor(a[i].x/n+0.5);
}
int n,m,lim,L;
int main() {
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++)scanf("%lf",&x[i].x);
	for(int i=0;i<=m;i++)scanf("%lf",&y[i].x);
	for(lim=1,L=-1;lim<=n+m;lim<<=1,++L);
	for(int i=0;i<lim;i++)r[i]=(r[i>>1]>>1)|((i&1)<<L);
	FFT(x,lim,1),FFT(y,lim,1);
	for(int i=0;i<lim;i++)x[i]=x[i]*y[i];
	FFT(x,lim,-1); 
	for(int i=0;i<=n+m;i++)printf("%d ",(int)x[i].x);
}

三步变两步优化

原方法中求了三次 FFT ，由于本蒟蒻常数较大，三次不是特别快

我们可以把 $B(x)$ 放到 $A(x)$ 的虚部上去，求出 $A(x)^2$，然后把 $A(x)$ 的虚部取出来除 2 就是答案了
正确性的证明：$(a+bi)^2=(a^2-b^2)+(2abi)$

这样的话效率是原来的 $\dfrac{2}{3}$

#include<bits/stdc++.h>
using namespace std;
const int N=2100000;
const double PI=acos(-1.0);
struct Comp {
	double x,y;
	inline Comp(double p=0.0,double q=0.0):x(p),y(q) { }
	inline Comp operator +(Comp o) { return Comp(x+o.x,y+o.y); }
	inline Comp operator -(Comp o) { return Comp(x-o.x,y-o.y); }
	inline Comp operator *(Comp o) { return Comp(x*o.x-y*o.y,x*o.y+y*o.x); }
}x[N];
int r[N];
inline void FFT(Comp *a,int n,int on) {
	register int i,j,k;
	register Comp wn,w,x,y;
	for(i=1;i<=n;i++)
		if(i<r[i])swap(a[i],a[r[i]]);
	for(i=1;i<n;i<<=1) {
		wn=Comp(cos(PI/i),on*sin(PI/i));
		for(j=0;j<n;j+=(i<<1)) {
			w=Comp(1,0);
			for(k=0;k<i;k++,w=w*wn) {
				x=a[j+k],y=w*a[i+j+k];
				a[j+k]=x+y,a[i+j+k]=x-y;
			}
		}
	}
	if(on==-1)
	for(i=0;i<=n;i++)
	a[i].x=floor(a[i].x/n+0.5),
	a[i].y=floor(a[i].y/n/2+0.5);
}
int n,m,lim,L;
int main() {
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++)scanf("%lf",&x[i].x);
	for(int i=0;i<=m;i++)scanf("%lf",&x[i].y);
	for(lim=1,L=-1;lim<=n+m;lim<<=1,++L);
	for(int i=0;i<lim;i++)r[i]=(r[i>>1]>>1)|((i&1)<<L);
	FFT(x,lim,1);
	for(int i=0;i<lim;i++)x[i]=x[i]*x[i];
	FFT(x,lim,-1); 
	for(int i=0;i<=n+m;i++)printf("%d ",(int)x[i].y);
}