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

FFT ( 快速傅里叶变换 ) 学习笔记

本文有错误之处请诸位大佬多多指正！

$FFT$ ：快速傅里叶变换的英文缩写，快速傅里叶变换是对离散傅里叶变换 $DFT$ 的优化,用来解决多项式上的操作如 卷积 等问题。

参考文章：

https://www.luogu.com.cn/blog/command-block/fft-xue-xi-bi-ji

https://blog.csdn.net/weixin_43346722/article/details/116353343

https://www.cnblogs.com/zwfymqz/p/8244902.html

多项式

系数表示法

一般在 初中 数学上，表示一个多项式我们用的是系数表示法
设 $A(x)$ 表示一个关于$x$ 的 $n$ 次的多项式，那么：
$$A(x)=\sum _{i=0}^n{a}_i\ast x^i$$
展开后就是：
$$A(x)=a_0+a_1\ast x+a_2\ast x^2+\cdot \cdot \cdot \cdot \cdot \cdot +a_n\ast x^n$$
然后两个多项式相乘，每一项都要和另一个多项式的每一项相乘，时间复杂度$O(n^2)$

点值表示法

我们找 $n+1$ 组不同的 $x$ 代入多项式，就可以得出 $n+1$ 个 $y$ 值，
我们把它们对应：$(x_0,y_0),(x_1,y_1)\cdot \cdot \cdot \cdot \cdot \cdot (x_n,y_n)$

其中 $y_i=\sum _{j=0}^{n-1}{a}_j\ast x_i^j$
听说这玩意叫做 $DFT$，怎么和我的名字缩写一样
当然，你会发现时间复杂度还是没有变化，选点 $O(n)$，计算 $O(n)$，总共 $O(n^2)$
然后好像没有啥卵用
最后你可悲的发现，好像时间复杂度在$O(n^2)$处达到了下限，没有什么优化的空间。
但是数学家们的智慧不是你能揣测的。

前方第一波高能！！！

复数

前置芝士

向量

具有大小又有方向的量，在物理上又称矢量，比如说力，加速度等都是向量。
在图上一般用一根箭头来表示。

弧度制

$1{}^{\circ }=\frac{\pi }{180}rad$
$180{}^{\circ }=1\pi rad$

定义

我们设 $a,b\in R$，$i^2=-1$,所有形如 $a+bi$ 的数通称为复数，其中 $i$ 被称为虚数单位，$bi$ 被称为复数的虚部，$a$ 被称为实部，如果一个复数只有 $bi$ 部分，那么这个复数被称为纯虚数。
复数域 $C$ 是目前最大的域，实数集是复数集的真子集。

如何表示复数？

我们建立一个平面直角坐标系，称之为复平面， $x$ 轴代表实轴，$y$ 轴代表虚轴 (除原点外)，那么 $x$ 轴上的点表示的就是实数（类比数轴），$y$ 轴除原点外表示的都是纯虚数，原点表示实数 $0$ 。那么对于每一个复数 $a+bi$，我们都可以用复平面上的点 $(a,b)$ 来表示这个复数，或者我们可以看做是起点为 $(0,0)$ 终点为 $(a,b)$ 的向量。

复数的模长：就是向量的大小，具体来说就是$\sqrt{a^2+b^2}$
幅角：假以逆时针为正方向，从 $x$ 轴正半轴到已知向量的转角的有向角叫做幅角.

运算法则

加法遵循向量的平行四边形法则，这个玩意长这样：（盗的别人的图）

$\overrightarrow{AB}+\overrightarrow{BC}=\overrightarrow{AD}+\overrightarrow{DC}=\overrightarrow{AC}$

减法：加法的时候将向量变成反过来的（终点变起点，起点变终点）

乘法：
几何定义：模长相乘，幅角相加（虽然我也不懂）

证明：
对于每一个复数，都可以表示成$r(\cos x+i\ast \sin x)$ 的形式，其中 $r$ 为复数的模长，$x$ 是幅角 (不理解的建议学学三角函数 勾股定理 )
对于两个复数 $r_1(\cos x+i\ast \sin x),r_2(\cos y+i\ast \sin y)$
相乘，得到的答案： ( $i^2=-1$ 所以第二项变成了减号 )
$$r_1\ast r_2\ast [(\cos x\cos y-\sin x\sin y+i\ast (\cos x\sin y+\cos y\sin x)]$$
根据三角恒等变换式 (不会的翻翻高中数学必修第一册最后一章)
$\cos (x+y)=\cos x\cos y-\sin x\sin y,\sin (x+y)=\sin x\cos y+\sin y\cos x$
再看看上面的式子，我们发现可以化简：
$$r_1{r}_2\ast (\cos (x+y)+i\ast \sin (x+y))$$
所以就有了模长相乘 (显然)，幅角相加 (还是显然)

代数定义：
$$(a+bi)\ast (c+di)=ac+adi+bci+bd{i}^2$$
因为$i^2=-1$ ，就有$$(ac-bd)+(ad+bc)i$$
搞定。

然后你会一头雾水的发现，这和多项式乘法有什么关系啊？
别急，真正的好东西来了。

单位根

注：以下提到的复数都是在单位圆上的复数。
在复平面上，以原点为圆心，$1$ 为半径作圆，所得的圆叫单位圆。
又去盗了一张图

以圆点为起点，圆的 $n$ 等分点为终点，做$n$个向量，设幅角为正且最小的向量对应的复数为 ${\omega }_n^1$，称为 $n$ 次单位根。
就是 $x^n=1$ 在复数域上的 $n$ 个解
根据复数乘法的运算，其余 $n-1$ 个复数为${\omega }_n^2,{\omega }_n^3,\dots ,{\omega }_n^n$
其中${\omega }_n^0={\omega }_n^n=1$（对应在 $x$ 轴上长度为 $1$ 的向量）
然后我们就有$${\omega }_n^k=\cos k\frac{2\pi }{n}+i\ast \sin k\frac{2\pi }{n}$$
单位根的几条性质：

$(1):{\omega }_n^k=({\omega }_n^1{)}^k$
证明：两个向量相乘，模长相乘还是 1 ，幅角相加，总共加了 $k$ 次 (也可以看做向逆时针方向旋转了 $\frac{k}{n}$ 个圆周)，所以就是${\omega }_n^k$ (显而易见)

(2)：${\omega }_n^i\ast {\omega }_n^j={\omega }_n^{i+j}$
证明：${\omega }_n^i=({\omega }_n^1{)}^i,{\omega }_n^j=({\omega }_n^1{)}^j,{\omega }_n^i\ast {\omega }_n^j=({\omega }_n^1{)}^i\ast ({\omega }_n^1{)}^j=({\omega }_n^1{)}^{i+j}={\omega }_n^{i+j}$

(3)：${\omega }_n^k={\omega }_{pn}^{pk}$
证明：把 ${\omega }_{pn}^{pk}$ 这玩意代入到单位根的表示的式子里，发现分子分母同时乘了一个 $p$ ，复数的值不变

(4)：${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$
证明：${\omega }_n^{\frac{n}{2}}=\cos \frac{n}{2}\cdot \frac{2\pi }{n}+i\ast \sin \frac{n}{2}\cdot \frac{2\pi }{n}=\cos \pi +i\ast \sin \pi$

众所周知 根据三角函数表，其中 $\cos \pi =-1,\sin \pi =0$

所以${\omega }_n^{\frac{n}{2}}=-1,{\omega }_n^{k+\frac{n}{2}}={\omega }_n^k\ast {\omega }_n^{\frac{n}{2}}=-{\omega }_n^k$

当然你可以看做把这个复数向逆时针方向旋转了$\frac{n}{2}$ 个圆周 (就是一个半圆)，然后这个新的复数和之前的复数就关于原点对称了 (建议画图理解)

(5)：${\omega }_n^k(k\ge n)={\omega }_n^{k\mathrm{m}\mathrm{o}\mathrm{d}n}$
证明：很明显这个复数围绕单位圆转了 $\frac{n}{n}=1$ 圈后回到了之前的位置，对答案没有影响。
用三角函数的角度来理解就是终边相同的角三角函数值相等

正文真正的开始

快速傅里叶变换：

前文说到，对于一个 $n$ 次多项式，我们可以用 $n+1$ 个点来确定
对于一个 $n-1$ 次多项式 $A(x)$，我们的超级神仙傅里叶 降下了一道救世之光，拯救了那些被多项式梦魇诅咒的善良之人 反手将这个多项式拍成了渣渣，然后按照奇偶性分成了两半，就有了这个东西：
$A(x)=(a_0+a_2{x}^2+\cdot \cdot \cdot +a_{n-2}{x}^{n-2})+(a_{1}x+a_3{x}^3+\cdot \cdot \cdot +a_{n-1}{x}^{n-1})$
这时候我们设
$A_1(x)=a_0+a_{2}x+\cdot \cdot \cdot +a_{n-2}{x}^{\frac{n}{2}-1}$
$A_2(x)=a_1+a_{3}x+\cdot \cdot \cdot +a_{n-1}{x}^{\frac{n}{2}-1}$
对应上面把$A(x)$切成两半的操作，我们 惊悚 惊喜地发现：
$$A(x)=A_1(x^2)+x{A}_2(x^2)$$
这时候我们把${\omega }_n^k$ ( $k\le \frac{n}{2}$ ) 代入这个式子：
$$A({\omega }_n^k)=A_1({\omega }_n^{2k})+{\omega }_n^k{A}_2({\omega }_n^{2k})$$
然后再把${\omega }_n^{k+\frac{n}{2}}$ 代入，得到：
$$A({\omega }_n^{k+\frac{n}{2}})=A_1({\omega }_n^{2k+n})+{\omega }_n^{k+\frac{n}{2}}{A}_2({\omega }_n^{2k+n})=A_1({\omega }_n^{2k})+{\omega }_n^{k+\frac{n}{2}}{A}_2({\omega }_n^{2k})$$
再根据 ${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$，得到
$$A({\omega }_n^{k+\frac{n}{2}})=A_1({\omega }_n^{2k})-{\omega }_n^k{A}_2({\omega }_n^{2k})$$
我们惊喜的发现，两个式子只有一个常数项不同，那么你每一次枚举一个式子的时候，你可以在 $O(1)$ 时间复杂度之内求出第二个式子的值，那么你的问题就缩小了一半！这时候你还会发现，对于第一个式子，继续拆分，它的问题规模又可以减半，那就可以递归实现了 （是不是很玄学）

不要以为这样就结束了，逼事还有一大堆呢

快速傅里叶逆变换 (IFFT)：

当你以为你搞定了前面的东西的时候，你发现你得到的是一堆点值表达式，好像并没有卵用
你悲哀的发现你只会系数表示法， 点值表示法你不会······
所以你就必须将这个点值表示法的多项式转换为系数表示法来表示。
所以神仙数学家们又搞了一个叫做快速傅里叶逆变换 ( $IFFT$ ) 的东西。

内容：

我们设$(y_0,y_1,y_2\cdot \cdot \cdot y_{n-1})$ 为多项式$(a_0,a_1,a_2\cdot \cdot \cdot a_{n-1})$ 的傅里叶变换
这时候我们再设 $(c_0,c_1,c_2\cdot \cdot \cdot c_{n-1})$ 是 $(y_0,y_1,y_2\cdot \cdot \cdot y_{n-1})$ 在 ${\omega }_n^0,{\omega }_n^{-1}\cdot \cdot \cdot ,{\omega }_n^{-(n-1)}$ 处的傅里叶变换
那么就有 $c_k=\sum _{i=0}^{n-1}{y}_i\ast ({\omega }_n^{-k}{)}^i$
推式子高能警告！！！
展开式子：
$c_k=\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^i{)}^j)\ast ({\omega }_n^{-k}{)}^i)$

$=\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^j{)}^i)\ast ({\omega }_n^{-k}{)}^i)$

$=\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^j{)}^i\ast ({\omega }_n^{-k}{)}^i)$

$=\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^{j-k}{)}^i)$

$=\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^{j-k}{)}^i)$

$=\sum _{j=0}^{n-1}{a}_j(\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i)$ 注意 $j,i$ 位置的调换

这时候我们设 $S(x)=\sum _{i=0}^{n-1}{x}_i$
代入 ${\omega }_n^k$ ，有
$$S({\omega }_n^k)=1+{\omega }_n^k+({\omega }_n^k{)}^2+\cdot \cdot \cdot +({\omega }_n^k{)}^{n-1}$$
当 $k!=0$ 时，这时候 ${\omega }_n^k!=1$
将这个式子乘上 ${\omega }_n^k$ ,就有了：
$${\omega }_n^{k}S({\omega }_n^k)={\omega }_n^k+({\omega }_n^k{)}^2+({\omega }_n^k{)}^3+\cdot \cdot \cdot +({\omega }_n^k{)}^n$$
后式减去前式，得到：

${\omega }_n^{k}S({\omega }_n^k)-S({\omega }_n^k)=({\omega }_n^k{)}^n-1$

$S({\omega }_n^k)=\frac{({\omega }_n^k{)}^n-1}{{\omega }_n^k-1}$

$S({\omega }_n^k)=\frac{({\omega }_n^n{)}^k-1}{{\omega }_n^k-1}$

因为 ${\omega }_n^n=1$，所以这个分数的分子为 $0$，分数值为 $0$
当 $k=0$ 的时候，$S({\omega }_n^k)=1+1+1+\cdot \cdot \cdot \cdot \cdot =n$
这时候我们看回刚刚这个式子：
$c_k=\sum _{j=0}^{n-1}{a}_j(\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i)$
$=\sum _{j=0}^{n-1}{a}_{j}S({\omega }_n^{j-k})$
当 $j-k!=0$ 时，贡献为 $0$
当$j=k，j-k=0$ 时，贡献为 $n$
所以就有 $c_k=n\ast a_k$
最后$$a_k=\frac{c_k}{n}$$

这样我们就得到了点值表示法和系数表示法之间的关系啦。
这个过程和上一个过程一样，都可以分治递归来实现

迭代优化

你慢悠悠的打出了一个递归，交上洛谷，然后只有 $77$ 分
没错，递归常数太大，在 $10^6$ 的常数面前，它没了······ (但这并不表示这个算法是不可行的)
所以我们把递归扔掉，使用迭代
去盗了一个大佬的图

有没有看出什么显而易见的性质？
没错，我们需要求的序列实际是原序列下标的二进制反转 (真的玄学)
因此没有必要对数列的奇偶性进行分类
然后我就不会证了
来看看这位 command_block 大佬的博客
原文地址：https://www.luogu.com.cn/blog/command-block/fft-xue-xi-bi-ji

最后迭代就是从下往上合并

悄悄话

累死我了，终于打完了······
(热泪盈眶)

代码：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<queue>
#include<cstring>
#include<cmath>
#define r register
#define rep(i,x,y) for(r ll i=x;i<=y;++i)
#define per(i,x,y) for(r ll i=x;i>=y;--i)
using namespace std;
typedef long long ll;
const ll V=(1e7+10)/2;
const double pi=acos(-1.0);
ll in()
{
	ll res=0,f=1;
	char ch;
	while((ch=getchar())<'0'||ch>'9')
	 if(ch=='-') f=-1;
	res=res*10+ch-48;
	while((ch=getchar())>='0'&&ch<='9')
	 res=res*10+ch-48;
	return res*f;
}
ll n,m,f[V];
ll now,k;
struct complex 
{
	double x,y;
	complex (double xx=0,double yy=0) { x=xx,y=yy; }
}a[V],b[V];
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); }
void FFT(complex *a,ll v)
{
	rep(i,0,now-1)
	 if(i<f[i]) swap(a[i],a[f[i]]);
	for(r ll mid=1;mid<now;mid<<=1)
	{
		complex Wn(cos(pi/mid),v*sin(pi/mid));
		ll R=mid<<1;
		for(r ll j=0;j<now;j+=R)
		{
			complex val(1,0);
			for(r ll k=0;k<mid;++k,val=val*Wn)
			{
				complex x=a[j+k],y=val*a[j+mid+k];
				a[j+k]=x+y;
				a[j+mid+k]=x-y;
			}
		}
	}
}
int main()
{
	scanf("%lld%lld",&n,&m);
	rep(i,0,n) a[i].x=in();
	rep(i,0,m) b[i].x=in();
	now=1;//单位根中的k值，默认为2的自然数次幂 
	while(now<=n+m) now<<=1,++k;
	rep(i,0,now-1)
	 f[i]=(f[i>>1]>>1)|((i&1)<<(k-1));
	FFT(a,1);
	FFT(b,1);
	rep(i,0,now) a[i]=a[i]*b[i];
	FFT(a,-1);
	rep(i,0,m+n)
	 printf("%lld ",(ll)(a[i].x/now+0.49));
	return 0;
}