zoukankan      html  css  js  c++  java
  • [整理]FFT(快速傅里叶变换)随记

    0.前言

    本文对信息学竞赛中的 FFT 算法及其推导作出一点粗浅的探讨,希望能够帮到大家一点点。
    另:作者默认大家已经有了一些高中数学基础,对于课本内基础知识部分可能会涉及较少。建议初中生先学完有关三角函数复数向量的部分之后再来学习 FFT 。

    1.前置知识

    多项式

    形如(A(x)=sum_{i=0}^n a_ix^i)的式子称作多项式(Polynomial),多项式可以进行加减乘(除法在此不涉及)等运算。
    这种表示方法称作多项式的系数表示(Coefficient representation),多项式还有一种点值表示(Point-value representation)法:
    众所周知,平面上的(n)个点((x_1,A(x_1))cdots(x_n,A(x_n)))可以唯一确定一个(n-1)次多项式,于是我们就用这些点来表示一个多项式。
    点值表示的优点:将多项式乘法的复杂度由系数表示的(O(n^2))减少到了(O(n))(只需要将点值相乘)。
    缺点:系数表示与点值表示的互化还是需要(O(n^2)),不过它是可以优化的(这也正是 FFT 做的事情)。

    复数及其运算

    首先我们有(sqrt{-1}=i),并将形如(z=a+bi (a,binmathbb{R}))的数称为复数,并定义其共轭(ar{z}=a-bi)
    其中,(Re z=a)称为实部(Im z=b)称为虚部
    以上是复数的代数表示,其运算同多项式运算,运算时要注意(i^2=-1)。(例如((a+bi)(c+di)=(ac-bd)+(ad+bc)i)
    复数还可以与平面上的向量一一对应,我们把这种表示叫做复数的几何表示,该平面叫做复平面
    对于一个复数,它的几何表示由模长(向量的长度)和幅角(从(x)轴正半轴逆时针转到该向量的角度)组成。
    复数几何表示的加减法运算同向量运算,而乘法的规则是模长相乘、幅角相加
    欧拉公式:(e^{i heta}=cos heta+isin heta)

    单位根及其性质

    接下来车速加快,请仔细阅读!
    定义方程(z^n=1)的复数根为(n)单位根,记作(omega_n^k=e^{frac{2pi ki}{n}} (k=0,1,cdots,n-1))。容易发现,这(n)个单位根恰好将单位圆(n)等分。
    如图,这是方程(z^3=1)的三个复数根(1,-dfrac 12+dfrac{sqrt{3}}2i,-dfrac 12-dfrac{sqrt{3}}2i)
    1.0
    如果一个单位根的(0,cdots,n-1)次方恰好构成互不相同的所有单位根,我们就将其称为本原单位根
    显然(omega_n=e^{frac{2pi i}{n}})是一个(n)次本原单位根(为方便,下文中“本原单位根”即特指(omega_n))。
    单位根有一些美妙的性质,我们来看一下。(假设以下的(n)都是正偶数)
    ( ext{Theorem 1: }omega_n^{2k}=omega_{frac{n}2}^k)
    ( ext{Proof: })利用复数几何表示的相乘法则证明。
    ( ext{Theorem 2: }omega_n^{frac{n}2+k}=-omega_n^k)
    ( ext{Proof: })可以理解为多转了半圈之后恰与原来相反。
    习题一:写出方程(z^6=1)的所有复数根。

    2.Cooley-Tukey算法

    该算法分为两个过程: DFT(Discrete Fourier Transform)IDFT(Inverse Discrete Fourier Transform) ,主要思想是分治
    铺垫了这么多,我们要想,如何利用单位根的一些性质来减少运算量呢?
    两位巨佬 Cooley 和 Tukey 想出了一种方法:将每一项按照指数奇偶分类。
    我们要计算一个多项式的 DFT ,也就是计算其点值表示,需要将每一个单位根代入计算。
    (A(omega_n^k)=sum_{j=0}^{n-1}a_iomega_n^{kj}=sum_{j=0}^{frac{n}2-1}a_{2j}omega_n^{2kj}+omega_n^ksum_{j=0}^{frac{n}2-1}a_{2j+1}omega_n^{2kj})
    根据刚刚提到的( ext{Theorem 1}),我们可以将其改写成(sum_{j=0}^{frac{n}2-1}a_{2j}omega_{frac{n}2}^{kj}+omega_n^ksum_{j=0}^{frac{n}2-1}a_{2j+1}omega_{frac{n}2}^{kj})
    此时再根据( ext{Theorem 2}),我们可以写出(A(omega_n^k))(A(omega_n^{k+frac{n}2}))的表示:

    [egin{aligned} A(omega_n^k)&=sum_{j=0}^{frac{n}2-1}a_{2j}omega_{frac{n}2}^{kj}+omega_n^ksum_{j=0}^{frac{n}2-1}a_{2j+1}omega_{frac{n}2}^{kj}\ A(omega_n^{k+frac{n}2})&=sum_{j=0}^{frac{n}2-1}a_{2j}omega_{frac{n}2}^{kj}+omega_n^{k+frac{n}2}sum_{j=0}^{frac{n}2-1}a_{2j+1}omega_{frac{n}2}^{kj}\ &=sum_{j=0}^{frac{n}2-1}a_{2j}omega_{frac{n}2}^{kj}-omega_n^ksum_{j=0}^{frac{n}2-1}a_{2j+1}omega_{frac{n}2}^{kj} end{aligned} ]

    我们成功了!现在只需要分别求奇数项和偶数项的 DFT 就能求出来两个值了!运算量减少了一半!(可以证明这个过程的确是(O(nlog n))的)
    但是先别高兴得太早,我们还有一步:将运算完的点值表示转化回系数表示(毕竟大多数时候我们要的是乘积多项式的系数)。这时候应该怎么办呢?
    Notice:由于证明较复杂且超出了高中的知识范围,建议从此处开始跳过直接看结论。
    注意到我们求值的过程实际上就是计算了一个矩阵乘法:

    [egin{bmatrix} (omega_n^0)^0&(omega_n^0)^1&cdots&(omega_n^0)^{n-1}\ (omega_n^1)^0&(omega_n^1)^1&cdots&(omega_n^1)^{n-1}\ vdots&vdots&ddots&vdots\ (omega_n^{n-1})^0&(omega_n^{n-1})^1&cdots&(omega_n^{n-1})^{n-1} end{bmatrix} egin{bmatrix} a_0\a_1\vdots\a_{n-1} end{bmatrix}= egin{bmatrix} A(omega_n^0)\A(omega_n^1)\vdots\A(omega_n^{n-1}) end{bmatrix} ]

    而我们的任务就是求系数矩阵(假设叫做(V))的逆矩阵。
    注意到这个矩阵是一个 Vandermonde 矩阵,对此,我们可以很容易地求出它的逆矩阵。
    两个比较繁复的证明:
    利用行列式及伴随矩阵
    利用拉格朗日插值
    然而最终的结果非常优美:设(V={v_{ij}}),则(V^{-1}=frac 1n{v_{ij}^{-1}})
    Notice:现在你可以在此停止了,结论就在下方。
    所以,我们只需要将 DFT 中(omega_n^k)换成(omega_n^{-k})(体现为虚部取相反数),再将答案乘以(frac 1n)就行了!
    由此,我们就可以用一个函数来实现 DFT 和 IDFT 的过程了。

    3.递归实现 FFT

    一层层向下递归,求出 FFT :

    void FFT(C *f,int len,int opt){//f是要进行DFT的数组,len是长度,opt是要进行的操作DFT/IDFT
    	if(len==1)return;
    	C *f0=f,*f1=f+(len>>1);
    	for(rg int i=0;i<len;i++)tmp[i]=f[i];
    	for(rg int i=0;i<(len>>1);i++){//按奇偶性分成两个多项式
    		f0[i]=tmp[i<<1],f1[i]=tmp[i<<1|1];
    	}
    	FFT(f0,len>>1,opt),FFT(f1,len>>1,opt);//分别进行FFT
    	C wn=C(cos(2*Pi/len),opt*sin(2*Pi/len)),w=C(1,0);//wn是本原单位根,w是wn的幂
    	for(rg int i=0;i<(len>>1);i++){
    		tmp[i]=f0[i]+w*f1[i];//按照上面推出来的式子进行计算
    		tmp[i+(len>>1)]=f0[i]-w*f1[i];
    		w=w*wn;
    	}
    	for(rg int i=0;i<len;i++)f[i]=tmp[i];
    }
    

    C++ STL 为我们提供了复数类,但是据说常数巨大,而手写复数类也不难写,所以还是建议大家手写:

    struct C {
    	double re,im;
    	C(double _re=0.0,double _im=0.0){
    		re=_re,im=_im;
    	}
    }f[N],g[N],tmp[N];
    il C operator + (C a,C b){
    	return C(a.re+b.re,a.im+b.im);
    }
    il C operator - (C a,C b){
    	return C(a.re-b.re,a.im-b.im);
    }
    il C operator * (C a,C b){
    	return C(a.re*b.re-a.im*b.im,a.re*b.im+a.im*b.re);
    }
    

    而最终实现可以这么写:

    FFT(f,len,1),FFT(g,len,1);
    ...//进行计算
    FFT(f,len,-1);
    

    4.迭代实现 FFT

    注意到上面的实现方式需要进行大量递归,常数巨大(据说在洛谷会只有77分,但我实测能很宽松地通过),我们考虑如何优化常数。
    这时候就需要拿出这张著名的图片了(注意下面最后两个框顺序反了):
    4.0
    这是我们进行递归的顺序,大家观察一下二进制有什么发现?
    可以发现,我们其实就是按照反转二进制排序之后依次分组的!如(0000,1000,0100,1100)倒过来就是(0000,0001,0010,0011)
    所以,我们可以将要进行 FFT 的数组按照反转二进制排序,然后直接迭代即可!

    for(rg int i=1;i<len;i++){
    	r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));//求出i的反转二进制(其中l是二进制位数)
    }
    

    而我们不需要递归的 FFT 函数变成了这样:

    void FFT(C *f,int len,int opt){
    	for(rg int i=0;i<len;i++){
    		if(i<r[i])swap(f[i],f[r[i]]);//排序
    	}
    	for(rg int i=2;i<=len;i<<=1){
    		int mid=i>>1;
    		C wn=C(cos(2*Pi/i),opt*sin(2*Pi/i));
    		for(rg int j=0;j<len;j+=i){
    			C w=C(1,0);
    			for(rg int k=0;k<mid;k++){
    				C z=f[j+k+mid]*w;
    				f[j+k+mid]=f[j+k]-z,f[j+k]=f[j+k]+z;
    				w=w*wn;
    			}
    		}
    	}
    }
    

    用时和内存都有了优化。
    另外,关于 FFT 还有一个优化:三次变两次。实现方法是将多项式(g)放到多项式(f)的虚部,这时(f^2)的虚部除以二就是我们要的答案。这个方法也可以省去不少常数。

    5.应用

    FFT 的最大应用是优化多项式卷积的计算。所谓卷积其实就是多项式乘法,它的定义是:对于两个函数(f(n),g(n)),它们的卷积((f*g)(n)=sum_{i+j=n}f(i)g(j))(就是个听起来挺nb的名字)
    所以,当我们碰到一个恶心的式子时,可以尝试将其化为卷积的形式。
    例题:洛谷P3338 [ZJOI2014]力
    题目让我们求(E_j=sum_{i=1}^{j-1}dfrac{q_i}{(i-j)^2}-sum_{i=j+1}^{n}dfrac{q_i}{(i-j)^2})的值,我们想办法将它化一化。
    为了好看,我们设(f_i=q_i,g_i=dfrac{1}{i^2})。这时候原式变成了(sum_{i=1}^{j-1}f_ig_{j-i}-sum_{i=j+1}^nf_ig_{i-j})
    我们优化一下上下界,令(f_0=g_0=0),有(E_j=sum_{i=0}^jf_ig_{j-i}-sum_{i=j}^nf_ig_{i-j})
    发现第一项已经是卷积形式,我们来着重推一推第二项(sum_{i=j}^nf_ig_{i-j})
    运用和式中变量代换的技巧,这个式子就等于(sum_{i=0}^{n-j}f_{i+j}g_i)
    这里还有一个常用的技巧,就是引入反转函数(f'_i=f_{n-i})。此时原式变成了(sum_{i=0}^{n-j}f'_{n-j-i}g_i),再代换一次就是(sum_{i=0}^mf'_{m-i}g_i)
    我们把它化成了卷积形式,现在可以用 FFT 做了!
    核心代码(代码里h[i]就是我们说的(f'_i)):

    int main(){
    	Read(n);
    	for(rg int i=1;i<=n;i++){
    		cin>>f[i].re;h[n-i].re=f[i].re;
    		g[i].re=1.0/(1.0*i*i);
    	}
    	while(len<n+n+1)len<<=1,l++;
    	for(rg int i=1;i<len;i++){
    		r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    	}
    	FFT(f,len,1),FFT(g,len,1),FFT(h,len,1);
    	for(rg int i=0;i<=len;i++)f[i]=f[i]*g[i],h[i]=h[i]*g[i];
    	FFT(f,len,-1),FFT(h,len,-1);
    	for(rg int i=1;i<=n;i++){
    		printf("%.3lf
    ",(f[i].re-h[n-i].re)/len);
    	}
    	return 0;
    }
    

    6.总结

    FFT 是多项式的基础算法,也是许多新手学多项式要过的第一关。理解了此处的知识,才能为之后的学习打下良好基础。
    习题二:洛谷P3803 【模板】多项式乘法(FFT)
    习题三:洛谷P1919 【模板】A*B Problem升级版(FFT快速傅里叶)

    7.完结撒花(凑数)

  • 相关阅读:
    Spring Boot2 系列教程(二十)Spring Boot 整合JdbcTemplate 多数据源
    Spring Boot 如何给微信公众号返回消息
    Spring Boot2 系列教程(十九)Spring Boot 整合 JdbcTemplate
    Spring Boot2 系列教程(十八)Spring Boot 中自定义 SpringMVC 配置
    Spring Boot 开发微信公众号后台
    Spring Boot2 系列教程(十七)SpringBoot 整合 Swagger2
    Spring Boot2 系列教程(十六)定时任务的两种实现方式
    Spring Boot2 系列教程(十五)定义系统启动任务的两种方式
    Spring Boot2 系列教程(十四)CORS 解决跨域问题
    JavaScript二维数组
  • 原文地址:https://www.cnblogs.com/juruoajh/p/14250217.html
Copyright © 2011-2022 走看看