浅谈FFT(快速博立叶变换)&学习笔记

2023-08-30 11:10:28

0XFF---FFT是啥?

FFT是一种DFT的高效算法，称为快速傅立叶变换（fast Fourier transform），它根据离散傅氏变换的奇、偶、虚、实等 特性，对离散傅立叶变换的算法进行改进获得的。 ---百度百科

对于两个多项式 \(F(x)\) 和 \(G(x)\) ，要求你将他们乘起来。

那还不简单？直接暴力相乘啊：

设 \(F(x)\) 的系数数列为 \(C\)。

\(F(x) \times G(x) = C_nx^nG(x) + C_{n-1}x^{n-1}G(x) + C_{n-2}x^{n-2}G(x) \cdots C_2x^2G(x) + C_1x^1G(x) + C_0G(x)\)

这样下来需要做做 \(n\) 次单项式乘多项式，每次的时间复杂度 \(O(n)\) ，则总复杂度高达 \(O(n^2)\)

基本上 \(n\) 上了\(4000\) 就会被卡吧......那怎么提速呢？

这就需要我们伟大而又神奇的神器：\(FFT\) (快速博立叶变换)

复杂度就只有 \(O(nlogn)\) 了。

0X1F---FFT的前置知识.

1.复数是什么?

我们把形如 \(z=a+bi\)（ \(a,b\) 均为实数）的数称为复数，其中 \(a\) 称为实部， \(b\) 称为虚部， \(i\) 称为虚数单位。当虚部等于零时，这个复数可以视为实数；当z的虚部不等于零时，实部等于零时，常称z为纯虚数。复数域是实数域的代数闭包，即任何复系数多项式在复数域中总有根。 复数是由意大利米兰学者卡当在十六世纪首次引入，经过达朗贝尔、棣莫弗、欧拉、高斯等人 的工作，此概念逐渐为数学家所接受。 ---百度百科

想必大家都知道实数是啥(不知道重读幼儿园吧......)，实数位于数轴上，就像下图这样：

浅谈FFT(快速博立叶变换)&学习笔记

我们稍微观察一下，\(1\) 是怎么变到 \(-1\) 的呢？

浅谈FFT(快速博立叶变换)&学习笔记

在数轴上转了 \(180°\)。

如果，是 \(90°\) 的话，会发生什么呢？

这个时候，会转到 \(0\) 上面的位置，但是那里，好像没有数啊!

不对，其实是有的，只不过这个数不在实数轴上，而是在虚数轴上！

虚数轴的单位是 \(i\) ，我们可以这么表示：

浅谈FFT(快速博立叶变换)&学习笔记

嗯，对。这显然是一个平面坐标系。现在我们的数仅限于数轴上，如果是这个平面坐标系上的一个点怎么表达呢？

对于下面的红色点：

浅谈FFT(快速博立叶变换)&学习笔记

这个点的坐标很容易的可以得到：\((2,i)\) ，也可以表示成 \(2+i\) .

你没猜错！这个就叫复数！

一个很重要的结论：复数相乘时，模长相乘，幅角相加！

2.点值表示法是什么?

我们用一个二维平面坐标系，在上面画 \(N+1\) 个点，最终可以解出一个 \(n\) 元的函数。证明略。

同样，我们可以用 \(N-1\) 个点来表达一个多项式。

因为点值相乘的复杂度只有 \(O(n)\) 显然优秀许多。

3.单位根是什么?

*n次单位根(n为正整数)是n次幂为1的复数!

*n次单位根(n为正整数)是n次幂为1的复数!

*n次单位根(n为正整数)是n次幂为1的复数!

我们先在复平面上画个点，就像这样：

浅谈FFT(快速博立叶变换)&学习笔记

Ta叫做单位圆。

圆边上的任意一点的模长都是 \(1\).

只有单位圆上的点表示的复数才有可能成为\(n\)次单位根!

单位根的基本符号：\(ω\)

一个单位圆，我们将它切成 \(n\) 份，从 \((1,0)\) 开始旋转，每次旋转 \(\frac{1}{n} \times 360\) 度，每次旋转后的点都记为 \(ω_{n}^{k}\)，特别的，\(ω_{n}^{0}\) 和 \(ω_{n}^{n}\) 都是 \((1,0)\) 点。

还有，当 \(k>=n\) 或者 \(k<0\) 时，\(ω_{n}^{k}\) 也是合法的。

单位根的性质：

\(1.\) 对于任意的 \(n\) , \(ω_{n}^{0}\) 都为 \((1,0)\) 点。

\(2.\) $ω_{n}^{a} \times ω_{n}^{b} = ω_{n}^{a+b} $

\(3.\) $ω_{an}^{ak} = ω_{n}^{k} $

\(4.\) $(ω_{n}{x})y = (ω_{n}{y})x $

\(5.\) $ω_{n}^{k+n/2} = -ω_{n}^{k} $ if(n%2==0)

0X2F---FFT的求解过程.

  • 分治思想很重要！

我们将多项式 \(F(x)\) 按位置分成两块。

那么变成了(保证n是2的正整数次幂)：

\(F(x) = (C_0+C_2x^2+C_4x^4+ \cdots +C_{n-2}x^{n-2}) + (C_1x+C_3x^3+C_5x^5+ \cdots +C_{n-1}x^{n-1})\)

设两个多项式 \(F1(x),F2(x)\)。

\(F1(x) = C_0+C_2x+C_4x^2+ \cdots +C_{n-2}x^{n/2-1}\)

\(F2(x) = C_1x+C_3x+C_5x^2+ \cdots +C_{n-1}x^{n/2-1}\)

则我们可以得出：

\(F(x) = F1(x^2) + F2(x^2) \times x\)

设 \(k<n/2\) ， 将 \(ω_{n}^{k}\) 带入多项式 \(F(x)\).

\(F(ω_{n}^{k}) = F1((ω_{n}^{k})2) + F2((ω_{n}^{k})^2) \times ω_{n}^{k}\)

简化得： \(F(ω_{n}^{k}) = F1(ω_{n/2}^{k}) + F2(ω_{n/2}^{k}) \times ω_{n}^{k}\)

再假设 \(k<n/2\) ，将 \(ω_{n}^{k+n/2}\) 带入多项式 \(F(x)\).

\(F(ω_{n}^{k+n/2}) = F1((ω_{n}^{k+n/2})2) + F2((ω_{n}^{k+n/2})^2) \times ω_{n}^{k}\)

\(F(ω_{n}^{k+n/2}) = F1(ω_{n}^{2k+n}) + F2(ω_{n}^{2k+n}) \times ω_{n}^{k+n/2}\)

\(F(ω_{n}^{k+n/2}) = F1(ω_{n}^{2k}) + F2(ω_{n}^{2k}) \times ω_{n}^{k+n/2}\)

\(F(ω_{n}^{k+n/2}) = F1(ω_{n/2}^{k}) + F2(ω_{n/2}^{k}) \times ω_{n}^{k+n/2}\)

\(F(ω_{n}^{k+n/2}) = F1(ω_{n/2}^{k}) - F2(ω_{n/2}^{k}) \times ω_{n}^{k}\)

比较一下两个式子：

  • \(F(ω_{n}^{k}) = F1(ω_{n/2}^{k}) + F2(ω_{n/2}^{k}) \times ω_{n}^{k}\)
  • \(F(ω_{n}^{k+n/2}) = F1(ω_{n/2}^{k}) - F2(ω_{n/2}^{k}) \times ω_{n}^{k}\)

等式右边只有一个负号的差别!

这两个式子很关键！

0X3F---FFT的代码实现.

对于复数的使用

虽然 \(C++ STL\) 里面有复数 \((complex)\) 但是太慢不建议大家使用。

你可以自己手打 \(complex\)

  • 手打的 \(complex\) :
struct complex{complex(double a=0,double b=0){x=a,y=b;}double x,y;};

complex operator + (complex a,complex b){return complex(a.x+b.x,a.y+b.y);}

complex operator - (complex a,complex b){return complex(a.x-b.x,a.y-b.y);}

complex operator * (complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
  • FFT:
complex a[N],b[N];

inline void FFT(complex *f,int len,short inv){

    if(!len)return;complex f1[len+1],f2[len+2];

    for(int k=0;k<len;++k)f1[k]=f[k<<1],f2[k]=f[k<<1|1];//按位置分

    FFT(f1,len>>1,inv);FFT(f2,len>>1,inv);//递归处理子问题

	complex tmp=complex(cos(PI/len),inv*sin(PI/len)),buf=complex(1,0);

	/*tmp:做一次平方后坐标的变换*/ /*buf:初始位置*/

	for(RI k=0;k<len;++k){

	    complex t=buf*f2[k];

	    f[k]=f1[k]+t,f[k+len]=f1[k]-t;buf=buf*tmp;//按照公式还原

	}return;

}

//注意，inv的作用是判断是 "系数转点值" 还是 "点值转系数"

\(Code\) 中提到的公式是这两项：

  • \(F(ω_{n}^{k}) = F1(ω_{n/2}^{k}) + F2(ω_{n/2}^{k}) \times ω_{n}^{k}\)
  • \(F(ω_{n}^{k+n/2}) = F1(ω_{n/2}^{k}) - F2(ω_{n/2}^{k}) \times ω_{n}^{k}\)

对于文中的"坐标的变换"：

我们依旧来看单位圆：

浅谈FFT(快速博立叶变换)&学习笔记

实际上，这个坐标的变换，直接用园中的三角形，运用三角函数就可以得出解了。

过程略.

最后我们得到的结果是：\(ω_{n}^{1} = (cos(\frac{2π}{n}),sin(\frac{2π}{n}))\)

求出 \(ω_{n}^{1}\) 后将它乘 \(n\) 次，可以得到：\({ω_{n}^{0},ω_{n}^{1},ω_{n}^{2},ω_{n}^{3}，ω_{n}^{4}，ω_{n}^{5} \cdots ω_{n}^{n-1}}\)

贴出最终的代码：

#include<cstdio>

#include<cmath>

#include<string>

#define ll long long

#define RI register int

#define inf 0x3f3f3f3f

#define PI 3.1415926535898

using namespace std;

const int N=6e4+2;

template <typename _Tp> inline _Tp max(const _Tp&x,const _Tp&y){return x>y?x:y;}

template <typename _Tp> inline _Tp min(const _Tp&x,const _Tp&y){return x<y?x:y;}

template <typename _Tp> inline void IN(_Tp&x){

    char ch;bool flag=0;x=0;

	while(ch=getchar(),!isdigit(ch))if(ch=='-')flag=1;

	while(isdigit(ch))x=x*10+ch-'0',ch=getchar();

	if(flag)x=-x;

}

struct complex{complex(double a=0,double b=0){x=a,y=b;}double x,y;};

complex operator + (complex a,complex b){return complex(a.x+b.x,a.y+b.y);}

complex operator - (complex a,complex b){return complex(a.x-b.x,a.y-b.y);}

complex operator * (complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

complex a[N],b[N];

inline void FFT(complex *f,int len,short inv){

    if(!len)return;complex f1[len+1],f2[len+2];

    for(int k=0;k<len;++k)f1[k]=f[k<<1],f2[k]=f[k<<1|1];

    FFT(f1,len>>1,inv);FFT(f2,len>>1,inv);

	complex tmp=complex(cos(PI/len),inv*sin(PI/len)),buf=complex(1,0);

	for(RI k=0;k<len;++k){

	    complex t=buf*f2[k];

	    f[k]=f1[k]+t,f[k+len]=f1[k]-t;buf=buf*tmp;

	}return;

}

int n,m;

int main(){

    scanf("%d%d",&n,&m);

    for(RI i=0;i<=n;++i)scanf("%lf",&a[i].x);

    for(RI i=0;i<=m;++i)scanf("%lf",&b[i].x);

    for(m+=n,n=1;n<=m;n<<=1);

    FFT(a,n>>1,1);FFT(b,n>>1,1);

    for(int i=0;i<n;++i)a[i]=a[i]*b[i];

    FFT(a,n>>1,-1);

    for(int i=0;i<=m;++i)printf("%.0f ",fabs(a[i].x)/n);

    putchar('\n');

    return 0;

}

听说可以优化，那啥的我还不会，就到这吧.

过了一会儿......

"原来FFT小优化这么简单啊!"

0X4F---FFT的一些小优化.

不用递归：

递归版(数组下标，先偶后奇，从0开始)：

0  1  2  3  4  5  6  7  --第1层

0  2  4  6 |1  3  5  7  --第2层

0  4 |2  6 |1  5 |3  7  --第3层

0 |4 |2 |6 |1 |5 |3| 7  --第4层

发现了什么吗？

最后的序列是原序列的二进制反转!

比如： \(6 = (110)_2\) 反过来变成了 \((011)_2 = 3\) ！

如何得到二进制翻转后的数列？递推即可!

for(RI i=0;i<n;++i)filp[i]=(filp[i>>1]>>1)|((i&1)?n>>1:0);

//filp[i] 即为 i 的二进制位翻转

Code:

#include<cstdio>

#include<cmath>

#include<string>

#define ll long long

#define RI register int

#define inf 0x3f3f3f3f

#define PI 3.1415926535898

using namespace std;

const int N=3e6+2;

int n,m,filp[N];

template <typename _Tp> inline _Tp max(const _Tp&x,const _Tp&y){return x>y?x:y;}

template <typename _Tp> inline _Tp min(const _Tp&x,const _Tp&y){return x<y?x:y;}

template <typename _Tp> inline void IN(_Tp&x){

    char ch;bool flag=0;x=0;

	while(ch=getchar(),!isdigit(ch))if(ch=='-')flag=1;

	while(isdigit(ch))x=x*10+ch-'0',ch=getchar();

	if(flag)x=-x;

}

struct complex{complex(double a=0,double b=0){x=a,y=b;}double x,y;};

complex operator + (complex a,complex b){return complex(a.x+b.x,a.y+b.y);}

complex operator - (complex a,complex b){return complex(a.x-b.x,a.y-b.y);}

complex operator * (complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

complex a[N],b[N];

inline void FFT(complex *f,short inv){

    for(RI i=0;i<n;++i)if(i<filp[i]){complex tmp=f[i];f[i]=f[filp[i]];f[filp[i]]=tmp;}

    /*换位置*/

    for(RI p=2;p<=n;p<<=1){//每局区间长度

        RI len=p/2;//合并子区间的长度(所以是p/2)

        complex tmp=complex(cos(PI/len),inv*sin(PI/len));

        for(RI k=0;k<n;k+=p){//每局左端点

            complex buf=complex(1,0);

            for(RI l=k;l<k+len;++l){//遍历区间

                complex t=buf*f[len+l];

				f[len+l]=f[l]-t,f[l]=f[l]+t,buf=buf*tmp;//赋值有微小的变化，注意顺序！

            }

        }

    }return;

}

/*主程序不变*/

int main(){

    scanf("%d%d",&n,&m);

    for(RI i=0;i<=n;++i)scanf("%lf",&a[i].x);

    for(RI i=0;i<=m;++i)scanf("%lf",&b[i].x);

    for(m+=n,n=1;n<=m;n<<=1);

    for(RI i=0;i<n;++i)filp[i]=(filp[i>>1]>>1)|((i&1)?n>>1:0);

    FFT(a,1);FFT(b,1);

    for(RI i=0;i<n;++i)a[i]=a[i]*b[i];

    FFT(a,-1);

    for(RI i=0;i<=m;++i)printf("%.0f ",fabs(a[i].x)/n);

    putchar('\n');

    return 0;

}

luogu上的题，递归的总是T最后一个点，改成非递归版的就A了？

emmmmmmmmmmmmmm

所有优化全开：

很作死，建议不要轻易尝试[滑稽]

#pragma GCC optimize(2)

#pragma GCC optimize(3)

#pragma GCC optimize("Ofast")

#pragma GCC optimize("inline")

#pragma GCC optimize("-fgcse")

#pragma GCC optimize("-fgcse-lm")

#pragma GCC optimize("-fipa-sra")

#pragma GCC optimize("-ftree-pre")

#pragma GCC optimize("-ftree-vrp")

#pragma GCC optimize("-fpeephole2")

#pragma GCC optimize("-ffast-math")

#pragma GCC optimize("-fsched-spec")

#pragma GCC optimize("unroll-loops")

#pragma GCC optimize("-falign-jumps")

#pragma GCC optimize("-falign-loops")

#pragma GCC optimize("-falign-labels")

#pragma GCC optimize("-fdevirtualize")

#pragma GCC optimize("-fcaller-saves")

#pragma GCC optimize("-fcrossjumping")

#pragma GCC optimize("-fthread-jumps")

#pragma GCC optimize("-funroll-loops")

#pragma GCC optimize("-fwhole-program")

#pragma GCC optimize("-freorder-blocks")

#pragma GCC optimize("-fschedule-insns")

#pragma GCC optimize("inline-functions")

#pragma GCC optimize("-ftree-tail-merge")

#pragma GCC optimize("-fschedule-insns2")

#pragma GCC optimize("-fstrict-aliasing")

#pragma GCC optimize("-fstrict-overflow")

#pragma GCC optimize("-falign-functions")

#pragma GCC optimize("-fcse-skip-blocks")

#pragma GCC optimize("-fcse-follow-jumps")

#pragma GCC optimize("-fsched-interblock")

#pragma GCC optimize("-fpartial-inlining")

#pragma GCC optimize("no-stack-protector")

#pragma GCC optimize("-freorder-functions")

#pragma GCC optimize("-findirect-inlining")

#pragma GCC optimize("-fhoist-adjacent-loads")

#pragma GCC optimize("-frerun-cse-after-loop")

#pragma GCC optimize("inline-small-functions")

#pragma GCC optimize("-finline-small-functions")

#pragma GCC optimize("-ftree-switch-conversion")

#pragma GCC optimize("-foptimize-sibling-calls")

#pragma GCC optimize("-fexpensive-optimizations")

#pragma GCC optimize("-funsafe-loop-optimizations")

#pragma GCC optimize("inline-functions-called-once")

#pragma GCC optimize("-fdelete-null-pointer-checks")

[滑稽][滑稽][滑稽][滑稽][滑稽][滑稽][滑稽][滑稽][滑稽][滑稽][滑稽]

最后，因为本人实在太弱了，太蒟了，所以实在写不出啥了。

\(by Qiuly\)

上一篇：Python办公自动化之Excel做表自动化：全网最全，看这一篇就够了！


下一篇：可视化设计，类Excel的快速开发平台