专栏文章

你能在不会虚数的情况下通关 FFT 吗？

rroger_yrj
算法·理论参与者 1已保存评论 0

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
0 条
当前快照
1 份
快照标识符
@miqn1kx8
此快照首次捕获于
2025/12/04 07:28
3 个月前
此快照最后确认于
2025/12/04 07:28
3 个月前
查看原文

你能在不会虚数的情况下通关 FFT 吗？

0.前言

本文主要是为了萌新 OIer 快速入门 FFT 而准备的，所以大部分证明将会跳过，想看证明请去别的文章。
众所周知 OI 不考证明。如果你觉得能在考场写出 FFT 并正确使用就够了，并不想了解详细证明，那么这篇文章可能适合你。
前置知识：多项式与生成函数，插值。
（如果你不会生成函数就没必要学 FFT 了）

1.应用

FFT 主要是用来解决多项式乘法

2.原理

对于一个 nn 次多项式，如果你知道 n+1n+1 个点值，就能确定这个多项式。
一个 nn 次多项式 F(x)F(x) 和一个 mm 次多项式 G(x)G(x) 相乘的结果是一个 n+mn+m 次多项式 H(x)H(x)，如果我们知道 n+m+1n+m+1 个不同的点值 (xi,F(xi))(x_i,F(x_i)) 和 (xi,G(xi))(x_i,G(x_i))，那么我们就能求出 (xi,H(xi))=(xi,F(xi)×G(xi))(x_i,H(x_i))=(x_i,F(x_i)\times G(x_i))，插值得到 H(x)H(x)。乘法这一步是 O(n)O(n) 的，问题转化成了如何快速点值和插值。

3.实现

为了方便后面的分治操作，我们可以在 F(x),G(x),H(x)F(x),G(x),H(x) 后面加若干系数为 00 的项，使得他们都变成 N=2kN=2^k 项的多项式。（N−1N-1 次多项式）
为了快速点值和插值，如何选择 xix_i 是关键。
定义一个虚数 ωn\omega_n，使得 ωnn=1\omega_n^n=1。你可以理解为 ωn=11n\omega_n=1^{\frac{1}{n}}（看起来很奇怪但它都是虚数了就顺从它吧）
性质：ω2n2k=ωnk\omega_{2n}^{2k}=\omega_n^k。你可以理解为 ω2n2k=12k2n=1kn=ωnk\omega_{2n}^{2k}=1^{\frac{2k}{2n}}=1^{\frac{k}{n}}=\omega_n^k
值：ωn=cos⁡(2πn)+sin⁡(2πn)i\omega_{n}=\cos(\frac{2\pi}{n})+\sin(\frac{2\pi}{n})i（直接记就行）
我们的目标是求出 F(x)F(x) 在 ωN0,ωN1,…,ωNN−1\omega_N^0,\omega_N^1,\dots,\omega_N^{N-1} 上的值。G(x)G(x) 同理。
设：
F(x)=a0+a1x1+a2x2+⋯+aN−1xN−1F(x)=a_0+a_1x_1+a_2x^2+\dots+a_{N-1}x^{N-1}
定义：
F[0](x)=a0+a2x1+a4x2+⋯+aN−2xN−22F^{[0]}(x)=a_0+a_2x^1+a_4x^2+\dots+a_{N-2}x^{\frac{N-2}{2}}
F[1](x)=a1+a3x1+a5x2+⋯+aN−1xN−22F^{[1]}(x)=a_1+a_3x^1+a_5x^2+\dots+a_{N-1}x^{\frac{N-2}{2}}
∴F(x)=F[0](x2)+xF[1](x2)\therefore F(x)=F^{[0]}(x^2)+xF^{[1]}(x^2)
问题转化为求 F[0](x)F^{[0]}(x) 在 ωN20,ωN21,…,ωN2N−1\omega_{\frac{N}{2}}^0,\omega_{\frac{N}{2}}^1,\dots,\omega_{\frac{N}{2}}^{N-1} 的值。F[1]F^{[1]} 同理。
根据定义，我们发现 ωN20=ωN2N2,ωN21=ωN2N2+1,…,ωN2N2−1=ωN2N−1\omega_{\frac{N}{2}}^0=\omega_{\frac{N}{2}}^{\frac{N}{2}},\omega_{\frac{N}{2}}^1=\omega_{\frac{N}{2}}^{\frac{N}{2}+1},\dots,\omega_{\frac{N}{2}}^{\frac{N}{2}-1}=\omega_{\frac{N}{2}}^{N-1}，所以只要求出 F[0](x)F^{[0]}(x) 在 ωN20,ωN21,…,ωN2N2−1\omega_{\frac{N}{2}}^0,\omega_{\frac{N}{2}}^1,\dots,\omega_{\frac{N}{2}}^{\frac{N}{2}-1} 的值，就能求出 F[0](x)F^{[0]}(x) 在 ωN2N2,ωN2N2+1,…,ωN2N−1\omega_{\frac{N}{2}}^{\frac{N}{2}},\omega_{\frac{N}{2}}^{\frac{N}{2}+1},\dots,\omega_{\frac{N}{2}}^{N-1} 的值。
所以我们可以通过不断将 aa 分组，然后就可以递归求解了。
你问我怎么插值？把 ωn\omega_n 换成 cos⁡(2πn)−sin⁡(2πn)i\cos(\frac{2\pi}{n})-\sin(\frac{2\pi}{n})i，然后把结果除以 NN 就行了。（证明跳过）

递归代码

CPP
void FFT(int N,complex *a){
	if(limit==1)return;
	complex a1[N>>1],a2[N>>1];
	for(int i=0;2*i<=N;i++)a1[i]=a[i<<1],a2[i]=a[i<<1|1];
	FFT(N>>1,a1);
	FFT(N>>1,a2);
    complex wn=complex(cos(2.0*Pi/N),sin(2.0*Pi/N)),w=complex(1,0);
	for(int i=0;i<(N>>1);i++,w*=wn)a[i]=a1[i]+w*a2[i],a[i+(N>>1)]=a1[i]-w*a2[i];
}
但是这种方法需要动态开区间，导致常数不可避免的大。

4.优化

我们观察 aa 经过分组后会变成什么样：
23926.png (800×136)
(图片来自 FlashHu 大佬的博客)
观察二进制：
CPP
000 001 010 011 100 101 110 111
变为：
CPP
000 100 010 110 001 101 011 111
我们发现就是上面的二进制翻转后的结果。
所以我们可以 O(n)O(n) 求出翻转后的结果 rr：
CPP
for(int i=1;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(k-1));
接下来只需从下到上不断合并即可。具体实现可以看看代码。

5.完整代码

CPP
#include<bits/stdc++.h>
#define C complex<double>
using namespace std;
const int N=4e6+10;
const double PI=acos(-1);
int n,m,r[N],l;
C f[N],g[N];
void FFT(C *a,int op){
	for(int i=0;i<n;i++)if(i<r[i])swap(a[i],a[r[i]]);
	for(int i=1;i<n;i<<=1){
		C wn{cos(PI/i),sin(PI/i)*op};
		for(int j=0;j<n;j+=(i<<1)){
			C w{1,0};
			for(int k=0;k<i;k++,w*=wn){
				C x=a[j+k],y=a[j+k+i];
				a[j+k]=x+w*y;
				a[j+k+i]=x-w*y;
			}
		}
	}
}
int main(){
	cin>>n>>m;
	for(int i=0,x;i<=n;i++)cin>>x,f[i]=x;
	for(int i=0,x;i<=m;i++)cin>>x,g[i]=x;
	m+=n;n=1;l=-1;
	while(n<=m)n<<=1,l++;
	for(int i=1;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<l);
	FFT(f,1);FFT(g,1);
	for(int i=0;i<n;i++)f[i]*=g[i];
	FFT(f,-1);
	for(int i=0;i<=m;i++)printf("%.0lf ",fabs(f[i].real()/n));
	return 0;
}

评论

共 0 条评论，欢迎与作者交流。

正在加载评论...