快速傅里叶变换（FFT）学习笔记

搬运自洛谷博客文章写于2024-01-31

看了一万篇博客终于学会了FFT，以后还要继续努力！

1. DFT与IDFT

假设有两个 \(n\) 次多项式 \(F(x)\) 和 \(G(x)\)，它们卷积的结果为 \(H=F*G\)

则有 \(H[x^k]=\sum_{i=0}^{k}F[x^i]G[x^{k-i}]\)

暴力算需要 \(O(n^2)\) 的时间复杂度

考虑插值的思想，如果知道 \(2n\) 次多项式 \(H(x)\) 经过的 \(2n+1\) 个点 \((x_1,y_1),(x_2,y_2),\dots,(x_{2n+1},y_{2n+1})\) ，那么可以唯一的确定 \(H(x)\) 的表达式

因此可以通过取 \(x_1,x_2,\dots,x_{2n+1}\) 的方式，求出 \(F(x_1)G(x1),F(x_2)G(x_2),\dots,F(x_{2n+1})G(x_{2n+1})\) ，即可得到 \(H(x)\)

其中，把系数表达转换为点值表达的算法叫做 DFT，把点值表达转换为系数表达的算法叫做 IDFT

2. DFT的原理

先把 \(n\) 扩大成 \(2^k\)

考虑 \(2^k-1=n-1\) 次多项式 \(F(x)\) ，按照每项次数分成奇偶两部分

\[F(x)=(F[x^0]+F[x^2]x^2+\dots+F[x^{n-2}]x^{n-2})+(F[x^1]x+F[x^3]x^3+\dots+F[x^{n-1}]x^{n-1}) \]

定义

\[FL(x)=F[x^0]+F[x^2]x+\dots+F[x^{n-2}]x^{n/2-1}\\ FR(x)=F[x^1]+F[x^3]x+\dots+F[x^{n-1}]x^{n/2-1} \]

那么 \(F(x)=FL(x^2)+xFR(x^2)\)

记 \(n\) 个 \(n\) 次单位根为 \(\omega_1,\dots,\omega_n\) ，根据单位根的性质，可以表示为 \(\omega_n^0,\omega_n^1,\dots,\omega_n^{n-1}\)

将这 \(n\) 个单位根分别代入 \(F(x)\)，使用如下方式：

\(k\) 从 \(0\) 循环至 \(n/2-1\)

先代入 \(\omega_n^k\) ，

\[\begin{aligned} F(\omega_n^k)&=FL((\omega_n^k)^2)+\omega_n^kFR((\omega_n^k)^2)\\ &=FL(\omega_{n/2}^k)+\omega_n^kFR(\omega_{n/2}^k) \end{aligned} \]

再代入 \(\omega_n^{k+n/2}\)，

\[\begin{aligned} F(\omega_n^{k+n/2})&=FL((\omega_n^{k+n/2})^2)+\omega_n^{k+n/2}FR((\omega_n^{k+n/2})^2)\\ &=FL(\omega_n^{2k+n})+\omega_n^{k+n/2}FR(\omega_n^{2k+n})\\ &=FL(\omega_n^{2k})+\omega_n^{k+n/2}FR(\omega_n^{2k})\\ &=FL(\omega_{n/2}^k)+\omega_n^{k+n/2}FR(\omega_{n/2}^k)\\ &=FL(\omega_{n/2}^k)-\omega_n^{k}FR(\omega_{n/2}^k) \end{aligned} \]

因此，我们可以根据

\[F(\omega_n^k)=FL(\omega_{n/2}^k)+\omega_n^kFR(\omega_{n/2}^k)\\ F(\omega_n^{k+n/2})=FL(\omega_{n/2}^k)-\omega_n^{k}FR(\omega_{n/2}^k) \]

进行分治，时间复杂度 \(O(nlog_2n)\)

3. IDFT的原理

假设 \(n\) 个单位根代入 \(F(x)\) 后得到的值为 \(y_0,\dots,y_{n-1}\)

由 DFT 可得

\[y_k=\sum_{i=0}^{n-1}(\omega_n^k)^iF[x^i] \]

先给出 IDFT 的结论：

\[nF[x^k]=\sum_{i=0}^{n-1}(\omega_n^{-k})^iy_i \]

证明：

\[\begin{aligned} \sum_{i=0}^{n-1}(\omega_n^{-k})^iy_i&=\sum_{i=0}^{n-1}\omega_n^{-ki}(\sum_{j=0}^{n-1}\omega_n^{ij}F[x^j])\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}\omega_n^{ij-ik}F[x^j]\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}\omega_n^{i(j-k)}F[x^j]\\ &=\sum_{i=0}^{n-1}\omega_n^{i(k-k)}F[x^k]+\sum_{i=0}^{n-1}\sum_{\substack{j=0\\j\neq k}}^{n-1}\omega_n^{i(j-k)}F[x^j]\\ &=nF[x^k]+\sum_{\substack{j=0\\j\neq k}}^{n-1}F[x^j]\sum_{i=0}^{n-1}\omega_n^{i(j-k)}\\ &=nF[x^k]+\sum_{\substack{j=0\\j\neq k}}^{n-1}F[x^j]\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\\ &=nF[x^k] \end{aligned} \]

于是

\[F[x^k]=\frac{1}{n}{\sum_{i=0}^{n-1}(\omega_n^{-k})^iy_i} \]

接下来与 IDFT 的实现原理类似，不多赘述

3. 优化

如果仅通过上述原理进行实现，由于常数过大还是会导致超时

位逆序置换

以 \(8\) 项多项式为例，模拟拆分的过程：

初始序列为 \(\{x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7\}\)

一次二分之后 \(\{x_0, x_2, x_4, x_6\},\{x_1, x_3, x_5, x_7 \}\)

两次二分之后 \(\{x_0,x_4\} \{x_2, x_6\},\{x_1, x_5\},\{x_3, x_7 \}\)

三次二分之后 \(\{x_0\}\{x_4\}\{x_2\}\{x_6\}\{x_1\}\{x_5\}\{x_3\}\{x_7 \}\)

规律：其实就是原来的那个序列，每个数用二进制表示，然后把二进制翻转对称一下，就是最终那个位置的下标。比如 \(x_1\) 是 001，翻转是 100，也就是 4，而且最后那个位置确实是 4。我们称这个变换为位逆序置换（bit-reversal permutation）。

可以通过数组线性预处理每个数二进制反转后得到的结果

for(ll i=0;i<n;i++){
	rev[i]=(rev[i>>1ll]>>1ll)|((i&1ll)?n>>1ll:0);
}

蝶形运算优化

已知 \(G(\omega_{n/2}^k)\) 和 \(H(\omega_{n/2}^k)\) 后，需要使用下面两个式子求出 \(f(\omega_n^k)\) 和 \(f(\omega_n^{k+n/2})\)：

\[\begin{aligned} f(\omega_n^k) & = G(\omega_{n/2}^k) + \omega_n^k \times H(\omega_{n/2}^k) \\ f(\omega_n^{k+n/2}) & = G(\omega_{n/2}^k) - \omega_n^k \times H(\omega_{n/2}^k) \end{aligned} \]
使用位逆序置换后，对于给定的 \(n, k\)：

\(G(\omega_{n/2}^k)\) 的值存储在数组下标为 \(k\) 的位置，\(H(\omega_{n/2}^k)\) 的值存储在数组下标为 \(k + \dfrac{n}{2}\) 的位置。

\(f(\omega_n^k)\) 的值将存储在数组下标为 \(k\) 的位置，\(f(\omega_n^{k+n/2})\) 的值将存储在数组下标为 \(k + \dfrac{n}{2}\) 的位置。

因此可以直接在数组下标为 \(k\) 和 \(k + \frac{n}{2}\) 的位置进行覆写，而不用开额外的数组保存值。此方法即称为 蝶形运算，或更准确的，基 - 2 蝶形运算。

再详细说明一下如何借助蝶形运算完成所有段长度为 \(\frac{n}{2}\) 的合并操作：

令段长度为 \(s = \frac{n}{2}\)；

同时枚举序列 \(\{G(\omega_{n/2}^k)\}\) 的左端点 \(l_g = 0, 2s, 4s, \cdots, N-2s\) 和序列 \(\{H(\omega_{n/2}^k)\}\) 的左端点 \(l_h = s, 3s, 5s, \cdots, N-s\)；

合并两个段时，枚举 \(k = 0, 1, 2, \cdots, s-1\)，此时 \(G(\omega_{n/2}^k)\) 存储在数组下标为 \(l_g + k\) 的位置，\(H(\omega_{n/2}^k)\) 存储在数组下标为 \(l_h + k\) 的位置；

使用蝶形运算求出 \(f(\omega_n^k)\) 和 \(f(\omega_n^{k+n/2})\)，然后直接在原位置覆写。

除此之外，可以将 fft 函数迭代实现，避免函数递归浪费时间。

4. 代码实现

STL中提供了复数模板，但由于常数较大，还是使用了手写复数类。

计算 \(\pi\) 可以通过 \(arccos(-1)\) 得到。

加入位逆序置换和蝶形优化后，代码如下：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef long double ld; 
const ll N=3000009; 
const ld PI=acos(-1);
struct Complex{
	ld a,b;
	Complex(ld x=0,ld y=0):a(x),b(y){}
	Complex operator+(Complex x){
		return Complex(a+x.a,b+x.b);
	}
	Complex operator-(Complex x){
		return Complex(a-x.a,b-x.b);
	}
	Complex operator*(Complex x){
		return Complex(a*x.a-b*x.b,a*x.b+b*x.a);
	}
} f[N],g[N];
ll n,m;
inline ll read(){
	ll c=0,f=1; char ch=getchar();
	while(ch<'0'||ch>'9'){
		if(ch=='-') f=-1; ch=getchar();
	}
	while(ch>='0'&&ch<='9'){
		c=(c<<1ll)+(c<<3ll)+(ch^48ll); ch=getchar();
	}
	return c*f;
}
inline void write(ll x){
	if(x<0) x=-x,putchar('-');
	if(x>9) write(x/10);
	putchar('0'+x%10);
}
ll rev[N];
inline void fft(Complex *f,bool flag){
	for(ll i=0;i<n;i++){
		if(i<rev[i]){
			swap(f[i],f[rev[i]]);
		}
	}
	for(ll p=2;p<=n;p<<=1ll){
		ll len=p>>1ll;
		Complex wn(cos(2*PI/p),sin(2*PI/p));
		if(!flag) wn.b*=-1;
		for(ll k=0;k<n;k+=p){
			Complex w(1,0);
			for(ll i=k;i<k+len;i++){
				Complex now=w*f[len+i];
				f[len+i]=f[i]-now;
				f[i]=f[i]+now;
				w=w*wn;
			}
		}
	}
}
int main(){
	n=read(),m=read();
	for(ll i=0;i<=n;i++){
		f[i].a=read();
	}
	for(ll i=0;i<=m;i++){
		g[i].a=read();
	}
	for(m+=n,n=1;n<=m;n<<=1ll);
	for(ll i=0;i<n;i++){
		rev[i]=(rev[i>>1ll]>>1ll)|((i&1ll)?n>>1ll:0);
	}
	fft(f,1);
	fft(g,1);
	for(ll i=0;i<n;i++){
		f[i]=f[i]*g[i];
	}
	fft(f,0);
	for(ll i=0;i<=m;i++){
		write((ll)(f[i].a/n+0.5));
		putchar(' ');
	}
	puts("");
	return 0;
}

请全文朗读背诵并默写。

posted @ 2025-04-20 21:25 QWQcoding 阅读(260) 评论(0) 收藏举报

刷新页面返回顶部

QWQcoding