其他分享
首页 > 其他分享> > FFT 学习笔记

FFT 学习笔记

作者：互联网

目录

FFT 学习笔记

前置知识

FFT 原理初探

以下涉及到的运算如无特殊说明均为复数域上的运算。

记 \(h(x)=f(x)\cdot g(x)\)。如果 \(f,g\) 用的是同一 \(x\) 序列得到的点值表示法，那么直接将对应点值相乘即可得到 \(h\) 的点值表示，这样做的复杂度是 \(O(n)\) 的。

这启发我们先把 \(f,g\) 由系数表达转化为点值表达，得到 \(h\) 的点值表达后再转为系数表达。系数表达转点值表达，随便选一组 \(x\) 代入是 \(n^2\) 的，不能接受。

还记得前置知识里的 \(n\) 次单位根吗？它拥有很良好的性质，我们不妨尝试把所有的 \(n\) 次单位根即 \(w_n^0,w_n^1,…,w_n^{n-1}\) 代入 \(f\) 中求出点值表达，这个过程叫作 DFT。在此我们假定 \(n\) 是 \(2\) 的整次幂，如果 \(n\) 不足就补到最近的幂，还有个小细节是如果 \(n\) 本身就是 2 的整次幂，仍需把 \(n\) 补到 \(2n\)，因为 \(n\) 次函数需要 \(n+1\) 个点才能唯一确定。

DFT 中，点值表达的第 \(k\) 项 \(y_k=\sum_{i=0}^{n-1}a_i(\omega_n^k)^i\)，按照 \(i\) 的奇偶性拆开，\(y_k=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}(\omega_n^k)^{2i}+\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}(\omega_n^k)^{2i+1}\) 。注意到 \((\omega_n^k)^{2i}=(\omega_n^{2k})^i=(\omega_\frac{n}{2}^k)^i,(\omega_n^k)^{2i+1} = \omega_n^k\cdot(\omega_n^{2k})^i = \omega_n^k(\omega_\frac{n}{2}^k)^i\)，则 \(y_k=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}(\omega_{\frac{n}{2}}^{k})^i+w_{n}^{k}\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}(\omega_{\frac{n}{2}}^{k})^{i}\)。如果把 \(a_{2i}\) 看做一个次数为 \(n\) 的多项式 \(F\) 的系数表达，\(a_{2i+1}\) 看做另一个次数为 \(n\) 的多项式 \(G\) 的系数表达，前一个和式其实就是 DFT(F) 的第 \(k\% \frac{n}{2}\) 项，后一个和式则为 DFT(G) 的第 \(k\%\frac{n}{2}\) 项。所以我们只需要递归求出 DFT(F) 及 DFT(G) 就可以 \(O(n)\) 求出 DFT(F)。这实际上是一个分治的结构，在递归到 \(n=1\) 时，DFT(f) 就等于 \(f\) 系数表达中唯一的一项 \(a_0\)（这是由于 \(w_n^0=1\)），直接返回，然后自底向上合并，总时间复杂度显然为 \(O(n\log n)\)。

现在我们能够在 \(O(n\log n)\) 的时间内求出 DFT(h)，而由 DFT(h) 反推出 \(h\) 的过程叫作 idft。这里有个高明的结论是 h 的系数表达等于把 dft(h) 的结果当做另一个多项式的系数表达再做一遍 dft 然后把得到的序列除以 \(n\)，只不过这次 dft 代入的是 \(w_n^{-i}\) 而非 \(\omega_n^i\)，即

\[n\times h_i=\sum_{j=0}^{n-1}(\omega_n^{-i})^j\cdot dft_j(h)=\sum_{j=0}^{n-1}(\omega_n^{-i})^j \sum_{k=0}^{n-1}h_k(\omega_n^j)^k=\sum_{0\le j,k<n}h_k\omega_n^{j(k-i)} \]

对于每个 \(k\)，如果 \(k\not=i\)，那么 \(\omega_n^{j(k-i)}\) 构成了一个公比为 \(w_n^{k-i}\) 的等比数列，其首项为 \(\omega_n^0=1\)，末项为 \((\omega_n^{k-i})^{n-1}\)，由等比数列求和公式知 \(\sum_{0\le j<n}(\omega_n^{k-i})^j=\frac{(\omega_n^{k-i})^n-1}{\omega_n^{k-i}-1}=0\)。而如果 \(k=i\)，其贡献为 \(nh_k=nh_i\)，恰等于左式，证毕。所以我们可以把 dft 和 idft 放在一个函数里写：

#define db double
const int N=3e6+5;
const db Pi=acos(-1);
struct cp{
	db x,y;//重载复数类
	cp(){}
	cp(db a,db b){x=a,y=b;} 
}w[N],tmp[N];
cp operator+(const cp &x,const cp &y){return cp(x.x+y.x,x.y+y.y);}
cp operator-(const cp &x,const cp &y){return cp(x.x-y.x,x.y-y.y);}
cp operator*(const cp &a,const cp &b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
cp f[N],g[N],h[N],a[N];
int n,m;

void dft(cp *f,int n,int flg){
	if(n==1) return;
	fo(i,0,n-1) tmp[i]=f[i];
	cp *L=f,*R=f+n/2;
	fo(i,0,n/2-1) L[i]=tmp[i<<1];
	fo(i,0,n/2-1) R[i]=tmp[i<<1|1];
	dft(L,n/2,flg);
	dft(R,n/2,flg);
	cp sum=cp(1,0),I=cp(cos(2*Pi/n),flg*sin(2*Pi/n));
	fo(i,0,n-1){
		w[i]=sum;
		sum=sum*I;//这里是精度误差的重灾区
	}
	fo(i,0,n/2-1){
		cp x=L[i],y=R[i]*w[i];
		L[i]=x+y,R[i]=x-y;
	}
}

signed main(){
	cin>>n>>m;int lim=1;
	while(lim<=n+m) lim<<=1;
	fo(i,0,n) f[i]=cp(read(),0);
	fo(i,0,m) g[i]=cp(read(),0);
	dft(f,lim,1);
	dft(g,lim,1);
	fo(i,0,lim-1) h[i]=f[i]*g[i];
	dft(h,lim,-1);
	fo(i,0,n+m) cout<<(int)(h[i].x/lim+0.49)<<' ';
    //FFT 精度丢失很严重，你必须四舍五入才能保证正确性……
	return 0;
}

常数优化

蝴蝶变换

上述代码在实际题目中的表现极差，因为它采用递归实现，并且进行了大量的内存拷贝。细想一下，我们不断向下递归其实只是为了把系数序列按照奇偶性划分成两边。而在回溯时我们只是把这两边合并起来，并没有进行内存的拷贝等操作。如果能预先处理出递归到最底层时序列被划分成的模样，是不是就可以避免递归和内存拷贝了？

直接抛结论：最终得到的序列的第 \(i\) 项等于原序列的第 \(rev_i\) 项，其中 \(rev_i\) 表示把 \(i\) 的二进制翻转得到的数，如 \(rev_6=3(rev_{110}=011)\)。为什么？考虑最普通的序列分治（如归并排序）其实是按二进制下的最高位分成两边，最终它得到的还是原序列，因为原序列本来就是按照二进制最高位的顺序排列的。而在 dft 中，按奇偶性分类实则是按二进制下的最低位分成两边，等价于按照最低位的大小排序，最终得到的序列就是按 \(rev_i\) 从小到大排序的结果！因此我们可以提前预处理出 \(rev\) 数组，然后把序列按 \(rev\) 数组变换为它在原本递归到最底层时的样子，自底向上合并，这个过程被称为蝴蝶变换（蝶形优化）。你发现这样做连递归都省了，代码也更为短小精悍：

void fft(cp *f,int n,int flg){
	fo(i,0,n-1) if(i<rev[i]) swap(f[i],f[rev[i]]);
	w[0]=cp(1,0);
	for(int j=2;j<=n;j<<=1){
		cp I=cp(cos(2*Pi/j),flg*sin(2*Pi/j));
		go(i,j-1,0){
			if(i&1) w[i]=w[i>>1]*I;
			else w[i]=w[i>>1];
		}
		for(int i=0;i<n;i+=j){
			fo(k,i,i+j/2-1){
				cp qwq=w[k-i]*f[k+j/2];
				f[k+j/2]=f[k]-qwq;
				f[k]=f[k]+qwq;
			}
		}
	}
}

signed main(){
	cin>>n>>m;lim=1;
	fo(i,0,n) f[i]=cp(read(),0);
	fo(i,0,m) g[i]=cp(read(),0);
	while(lim<=n+m) lim<<=1;
	fo(i,0,lim-1) rev[i]=(rev[i>>1]>>1)|(i&1?(lim>>1):0);//你品，你细品
	fft(f,lim,1);
	fft(g,lim,1);
	fo(i,0,lim-1) h[i]=f[i]*g[i];
	fft(h,lim,-1);
	fo(i,0,n+m) cout<<(int)(h[i].x/lim+0.49)<<' '; 
	return 0;
}

还有个小 trick 是把单位根也一块预处理出来放在同一个数组里，不再赘述

三次变两次

一个不容忽视的事实是 \(f,g,h\) 这三个多项式的系数都只有实部没有虚部，在 fft 时却不得不当做完整的复数来运算，这难免会造成一些计算上的冗余。

考虑把 \(g\) 塞到 \(f\) 的虚部上去，令 \(P(x)=f(x)+g(x)i\)，则 \(P^2(x)=f^2(x)-g^2(x)+2f(x)g(x)i\)，因此 \(P^2(x)\) 的虚部除以 \(2\) 就是 \(f(x)g(x)\)。所以我们先算出 \(dft(P)\)，然后用它乘自己得到 \(dft(P^2)\)，再 idft 回去即可。这样只需要做两次 \(dft\)。

	cin>>n>>m;lim=1;
	fo(i,0,n) f[i].x=read();
	fo(i,0,m) f[i].y=read();
	while(lim<=n+m) lim<<=1;
	fo(i,0,lim-1) rev[i]=(rev[i>>1]>>1)|(i&1?(lim>>1):0);
	//out(rev,0,lim-1);
	fft(f,lim,1);
	fo(i,0,lim-1) h[i]=f[i]*f[i];
	fft(h,lim,-1);
	fo(i,0,n+m) cout<<(int)(h[i].y/2/lim+0.49)<<' ';

标签：lim,sum,FFT,笔记,学习,cp,omega,2i,fo
来源： https://www.cnblogs.com/vectorwyx/p/15725887.html