FFT相关

FFT相关

FFT 数组记得开两倍！

参考博客

快速傅里叶变换（一）FFT

快速傅里叶变换（二）NTT

快速傅里叶变换（三）

超详细傅里叶变换

FFT自己卷自己

FFT

简介

用于求卷积（\(a,b\) 已知）：

\[\sum_{i=0}^n a_ib_{n-i} \]

或者多项式乘法（\(A(x),B(x)\) 已知）：

\[C(x)=A(x)B(x) \]

\(A(x)=\sum_{i=0}^{n} a_i x^i\\ B(x)=\sum_{i=0}^{m} b_i x^i\)

可见 \(C(x)\) 是 \(n+m\) 次多项式。

如果我们把卷积的 \(a_i,b_i\) 看成多项式的系数，卷积就变成求：

\[[x^n] C(x)=A(x)B(x) \]

其中 \([x^n]\) 表示 \(x^n\) 的系数。

求卷积或者多项式乘法的时间复杂度是 \(O(n^2)\) 的。使用 FFT 可以做到 \(O(n\log n)\)。

大体思路

设 \(C(x)\) 的项数为 \(n\)（不是次数），若 \(A,B\) 不足 \(n\) 项就补系数 \(0\)。

显然这个过程很难优化，我们从另一个角度去想。

对于一个多项式，求其在 \(x\) 处的值的时间复杂度是 \(O(n)\) 的，我们把这个操作叫做点值（DFT）。

\(n\) 个点可以唯一确定一个 \(n\) 项多项式（即 \(n-1\) 次多项式），证明不显然但是略过我不会。而据说拉格朗日插值法给出了一种 \(O(n^2)\) 的求出这个多项式（即求出它的系数）的方法我不会，这个由点值求系数的过程叫做插值（IDFT）。

因此求多项式乘法，可以变成求任意 \(n\) 个点在两个多项式 \(A,B\) 的点值，然后由 \(C(x)=A(x)B(x) \ O(n)\) 求出多项式 \(C\) 在这 \(n\) 个点的点值，然后做一次插值求出 \(C\) 的系数。当然这个也是 \(O(n^2)\) 的，但是聪明的傅里叶给出了一种基于单位根的特殊性质的分治方法求 DFT 和 IDFT，成为快速傅里叶变换（FFT）。

单位根

写作 \(\omega_n^k\)，读作 \(n\) 次单位根的 \(k\) 次方。

\(\omega_n^k\) 是在复数域上的向量，形如在复平面上分成 \(n\) 等分，其中 \(\omega_n^0=\omega_n^n=1\)。按逆时针分别为 \(\omega_n^0,\omega_n^1\dots \omega_n^{n-1}\)。

后面为了方便，我们会假设 \(n=2^k\)。

有几个重要的性质。

• \(\omega_n^n=1\) 正确性显然
• \(\omega_{an}^{ak}=\omega_n^k\) 在平面上想象一下，显然正确
• \(\omega_n^{k-\frac{n}{2}}=-\omega_n^k\) 相当于把向量转 \(180^。\)。

点值

基于这些性质，我们求 \(A,B\) 在 \(\omega_n^0,\omega_n^1\dots \omega_n^{n-1}\) 处的点值。

以求 \(A(x)\) 为例，我们要求所有 \(A(\omega_n^k),0\le k<n\)。

我们把 \(A\) 按奇偶分为两部分：

设

\[A_1(x)=a_0+a_2 x+a_4 x^2+\dots a_{n-2} x^{\frac{n-2}{2}}\\ A_2(x)=a_1+a_3x^2+a_5x^4+\dots a_{n-1}x^{\frac{n-2}{2}}\]

因此有：

\[A(x)=A_1(x^2)+xA_2(x^2) \]

代入 \(x=\omega_n^k\)，有：

\[A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^k A_2(\omega_n^{2k})\\ A(\omega_n^k)=A_1(\omega_{\frac{n}{2}}^k)+\omega_n^k A_2(\omega_{\frac{n}{2}}^k)\]

我们先求出 \(k< \frac{n}{2}\) 的 \(A_1,A_2\) 的点值。这个范围是缩小了一半的。然后把它们相加就得到了 \(A\) 在 \(k< \frac{n}{2}\) 的点值。

然后我们求剩下一半的 \(A_1,A_2\) 的点值。仍然设 \(k< \frac{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(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k})-\omega_n^{k} A_2(\omega_n^{2k})\\ A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_{\frac{n}{2}}^k)-\omega_n^{k} A_2(\omega_{\frac{n}{2}}^k)\]

因此你发现，这俩十分地相似，因此你求出 \(k< \frac{n}{2}\) 的 \(A_1,A_2\) 的点值之后，可以直接 \(O(n)\) 求出 \(A\) 的 \(n\) 个点值了。然后就这样分治下去，点值时间复杂度为 \(O(n\log n)\)。

因为要一直对 \(n\) 除以 \(2\) 所以令 \(n=2^k\) 意义就在此。（位数不足高位补系数 \(0\)）

插值

求插值的过程是类似的，改几个参数就行了。我还没搞懂是怎么推出来的

然后我们就这么求两次点值，再求一次插值就做完了。

改进

code

然后你会发现过不了板子……

因为这个递归的过程常数很大，假设 FFT 的常数本身就大，然后就超时了。

考虑从下往上递推分治的过程。

\[\left(a_{0}, a_{1}, a_{2}, a_{3}, a_{4}, a_{5}, a_{6}, a_{7}\right)\\ \left(a_{0}, a_{2}, a_{4}, a_{6}\right)\left(a_{1}, a_{3}, a_{5}, a_{7}\right)\\ \left(a_{0}, a_{4}\right)\left(a_{2}, a_{6}\right)\left(a_{1}, a_{5}\right)\left(a_{3}, a_{7}\right)\\ \left(a_{0}\right)\left(a_{4}\right)\left(a_{2}\right)\left(a_{6}\right)\left(a_{1}\right)\left(a_{5}\right)\left(a_{3}\right)\left(a_{7}\right) \]

然后你惊奇地发现最下层的顺序是 \(000,100,010,110,001,101,011,111\)，刚好是 \(0\sim 7\) 的二进制的 reverse。

然后就有了如下板子：

Code

struct fushu {
	double x,y;
	fushu (double _x=0,double _y=0):x(_x),y(_y){}
}a[N],b[N];
fushu operator + (fushu a,fushu b) { return {a.x+b.x,a.y+b.y}; }
fushu operator - (fushu a,fushu b) { return {a.x-b.x,a.y-b.y}; }
fushu operator * (fushu a,fushu b) { return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x}; }
int len;
const double pi=acos(-1.0);
int re[N];
void FFT(fushu *c,int type) {
	rep(i,0,(1<<len)-1) {
		if(i<re[i]) swap(c[i],c[re[i]]);
	}
	for(int k=1;k<(1<<len);k<<=1) {
		fushu wn(cos(pi/k),type*sin(pi/k));
		for(int r=k<<1,j=0;j<(1<<len);j+=r) {
			fushu w(1,0);
			for(int i=0;i<k;i++,w=w*wn) {
				fushu x=c[j+i],y=w*c[j+k+i];
				c[j+i]=x+y;
				c[j+k+i]=x-y;
			}
		}
	}
}
int n,m;
int main(){
	sf("%d%d",&n,&m);
	rep(i,0,n) sf("%lf",&a[i].x);
	rep(i,0,m) sf("%lf",&b[i].x);
	while((1<<len)<=n+m) len++;
	rep(i,0,(1<<len)-1) {
		re[i]=(re[i>>1]>>1)|((i&1)<<(len-1));
	}
	FFT(a,1),FFT(b,1);
	rep(i,0,(1<<len)-1) a[i]=a[i]*b[i];
	FFT(a,-1);
	rep(i,0,n+m) pf("%d ",(int)(a[i].x/(1<<len)+0.5));
}

算法缺陷

由于复数域是用浮点数计算的，所以会存在掉精度问题。如果答案是对一个特别的指数取模，如著名的 \(998244353\)，可以使用原根代替单位根计算，在剩余系里计算而不是在复数域计算。详见 NTT。

NTT

简介

NTT 的原理是用原根代替单位根。质数 \(p\) 的原根在模 \(p\) 剩余系意义下具有我们利用的单位根的性质的相同性质。因此如果要求的多项式的系数只需要模 \(p\) 意义下的，可以用原根代替单位根。因为原根是整型，所以可以避免精度问题。常见质数 \(998244353\) 的原根是 \(3\)。

Code

const int N=4e6+7,mod=998244353,G=3,invG=332748118;
ll a[N],b[N];
int len;
int re[N];
ll ksm(ll a,ll b=mod-2) {
	ll s=1;
	while(b) {
		if(b&1) s=s*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return s;
}
void NTT(ll *c,int type) {
	rep(i,0,(1<<len)-1) {
		if(i<re[i]) swap(c[i],c[re[i]]);
	}
	for(int k=1;k<(1<<len);k<<=1) {
		ll wn=ksm(type==1?G:invG,(mod-1)/(k<<1));
		for(int r=k<<1,j=0;j<(1<<len);j+=r) {
			ll w=1;
			for(int i=0;i<k;i++,w=w*wn%mod) {
				ll x=c[j+i],y=w*c[j+k+i]%mod;
				c[j+i]=(x+y)%mod;
				c[j+k+i]=(x-y+mod)%mod;
			}
		}
	}
}
int n,m;
int main(){
	#ifdef LOCAL
	freopen("in.txt","r",stdin);
	freopen("my.out","w",stdout);
	#endif
	sf("%d%d",&n,&m);
	rep(i,0,n) sf("%lld",&a[i]),a[i]%=mod;
	rep(i,0,m) sf("%lld",&b[i]),b[i]%=mod;
	while((1<<len)<=n+m) len++;
	ll inv=ksm(1<<len);
	rep(i,0,(1<<len)-1) {
		re[i]=(re[i>>1]>>1)|((i&1)<<(len-1));
	}
	NTT(a,1),NTT(b,1);
	rep(i,0,(1<<len)-1) a[i]=a[i]*b[i]%mod;
	NTT(a,-1);
	rep(i,0,n+m) pf("%lld ",a[i]*inv%mod);
}

分治 FFT

Luogu 模板

形如 \(f_i=\sum_{j=1}^i f_{i-j}g_j\) 的卷积（\(f_0\) 需要初值）。

一般的卷积我们是知道卷起来的两个多项式的系数的，但是这个卷积显然我们只知道 \(g\) 的系数，却不知道 \(f\) 的系数，那么如何卷呢？可以分治 FFT 解决。

类似于 CDQ 分治，要求 $f_{1\sim n} $，先求 \(f_{1 \sim mid}\)，然后计算左边对右边的贡献，即算 \(f_{1\sim n}=\sum_{j=1}^{mid} f_{j}g_{i-j}\)，做一次 \(O(n\log n)\) 的 FTT，因为这个卷积不是标准卷积（指 \(j\) 不是枚举到 \(i\)），因此我们可以给 \(j>mid\) 的 \(f_j\) 当做 \(0\) 来做卷积，相当于求 \(\sum_{j=1}^i f_j g_{i-j}\)。然后我们再求 \(f_{mid+1\sim r}\)，因为左边已经对右边贡献过了，所以我们递归右边的时候，计算贡献不需要再带上左边。

也就是说，求 \(f_{l \sim r}\)，先算左边，然后算左边对右边的贡献，即 \(\sum_{j=1}^{r-l+1} f_{j+l-1}g_{(r-l+1)-j}\)，这是一个标准卷积， FFT 时间复杂度是 \(O(len \log len)\)。这样分治下去，直到长度为 \(1\) 时，不用卷积了，直接返回值。分治一共 \(\log\) 层，每一层没个区间都要做 \(O(len\log len)\) 的 FTT，总时间复杂度为 \(O(n\log^2 n)\)。