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

zoukankan html css js c++ java

FFT(快速傅里叶变换)学习笔记
che dan环节：

概述

快速傅里叶变换(Fast Fourier Transformation)，简称(FFT)，它可以在(O(n log n))的时间内算两个多项式的乘积的所有系数。

是一个听起来逼格很高，实际逼格更高的算法。

一些定义

下面我们来交代一些奇怪的东西。

一个(n)项(n-1)次多项式，(F(x) = a_0x^0 + a_1x^1 + dots + a_{n-1}x^{n-1})

则多项式(F(x))可以用 __系数表示法__表示成(F(x) = { a_0,a_1,dots,a_{n-1} })

方便的，我们可以把(F)看成一个数组，用(F[k])表示多项式(F)的(k)次项系数，即(F[k] = a_k)，如果(F)不存在(k)次项的话系数就是零即(F[k] = 0)。

那么多项式(F)((n)次)，(G)((m)次)的乘积会得到一个(n + m)次多项式(H)就满足(H[k] = sumlimits_{i = 0}^{k} F[i]G[k-i])，按这样算多项式乘法就是(O(n^2))的。

如果您是个明白人的话会发现多项式乘法就是一个加法卷积（当然如果不知道也没关系）。

多项式还有一个比较有用的表示法叫 点值表示法，就是把多项式看做函数。

记一个结论：用(n)个点能确定一个(n-1)次多项式 ~~意会，意会～~~。

即代(n)个(x)到(F(x))里，就可以确定这个多项式$ F(x)(的)n$个系数(不要忘了常数项)。

假设我们代入的是({ x_0,x_1, dots ,x_{n-1} })，令(y_i = F(x_i))

那么(F(x))就可以表示为(F(x) = { (x_0,y_0),(x_1,y_1), dots,(x_{n-1},y_{n-1}) })

总结

总结一下，多项式(F(x) = a_0x^0 + a_1x^1 + dots + a_{n-1}x^{n-1})有两种表示方法
- 系数表示法：(F(x) = { a_0,a_1,dots,a_{n-1} })
- 点值表示法：即把(F(x))看做一个函数，用(n)个(x)就可以确定(F(x)),(F(x) = { (x_0,y_0),(x_1,y_1), dots,(x_{n-1},y_{n-1}) })
那这些有什么用呢？我们用系数表示法做乘法的复杂度是铁定(n)方的。

那用点值表示法呢？

如果要计算(F(x) cdot G(x))（(F)有(n)次，(G)有(m)次）的话，我们先把同一组({ x_i : 0 le i le n+m })（注意是(n + m + 1)个(x)，实际上多代了也不亏。）代入(F,G)，那么设(H(x) = F(x) cdot G(x))，多项式(H)的点值表示法就是

[H = { (x_0, F(x_0) cdot G(x_0)),dots,(x_{n+m},F(x_{n+m})cdot G(x_{n+m})) } ]
这样求出(F)和(G)的点至表示法后我们就可以通过(O(m +n))的算出(H)的点值表示。

然后呢？然后啥啊。。。

我们会得到$n + m + 1 $个点，但我们要算各项的系数啊。

那咋算呢？(FFT)就是做的这个。

提前交代一下(FFT)的流程
- (1.) 先分别算出两个多项式的点值表示法（我们称之为绝活(DFT)）
- (2.) 把这些对应点分别乘起来，得到乘积多项式(H)的点值表示法
- (3.) 再把(H)的点值表示法转回系数表示法（即(DFT)的逆变换，我们叫它(IDFT)）
前置姿势

复数

定义

形如(z = a + bi)的数我们管他叫复数，其中(i = sqrt{-1})，我们称为虚数单位，实数(a)称为(z)的实部，实数(b)称为(z)的虚部。

一个复数(z = a + bi)可以用复平面上的一个点((a,b))表示。

复平面：

可以这么理解，把他当成直角坐标系的话，上面的一点((x,y))就对应了一个复数(x + iy)。

原来直角坐标系下的(x)轴在复平面下对应实轴（上面的数表示实部(a)），(y)轴上的数对应虚部(b)，我们叫虚轴。

模长：

复数所对应的点即到点((0,0))的距离，记为(|z|)。

辐角：

复数(a + bi)的辐角即实轴按逆时针方向旋转到点((a,b))的旋转角度，类似直角坐标系下的极角。

记为(arg (z))

运算

加减法：

这个应该很好理解吧，实部虚部分别加减就好了，即((a + bi) pm (c+di) = (apm c) + (bpm d)i)，

乘法

恩，这个比较重要。

可以暴力的乘起来：((a + bi)(c + di) = ac + adi + cbi + bdi^2 = (ac - bd) + (ad + cb)i)（因为(i^2 = -1)）。

即点((a,b))与点((c,d))通过乘法得到点((ac - bd,ad + cb))

再记个结论：

设(z = z_0 imes z_1)，则(|z| = |z_0| cdot |z_1|,arg(z) = arg(z_0) + arg(z_1))

即 复数乘法，模长相乘，辐角相加（这玩意儿非常有用）

Code

这是复数运算的一些(Code)，并没有涉及到复数的除法，因为~~我不会~~(FFT)并不会用除法。
```
struct cpl{
	db x,y;
	cpl operator + (cpl k1)const{return (cpl){x+k1.x,y+k1.y};}
	cpl operator - (cpl k1)const{return (cpl){x-k1.x,y-k1.y};}
	cpl operator * (cpl k1)const{return (cpl){x*k1.x-y*k1.y,x*k1.y+y*k1.x};}
}
```
单位根

这个有点东西啊。

定义

可以理解为(x^n = 1)这个方程的根，我们称之为 (n)次单位根。

在实数域中，这个方程当(n)为奇数时有一个解(x = 1)，(n)为偶数时有两个解(x_0 = 1, x_1 =-1)

那扩展到复数域呢？

在复平面下，复数的乘法有一个非常(nb)的性质，叫 模长相乘，辐角相加（不知道的话百度一下吧！）

那么能成为(x^n = 1)这个方程的解的复数(x)，有一个前提条件叫(|x | = 1)（(x)的模长为(1)）。

那么单位根就只会在复平面的单位圆上（蛤？啥是单位圆）。

上图（单位圆）：

还有一点叫 辐角相加，那么(n)次单位根就一定是单位圆周上的(n)个(n)等分点。

因为这样的点辐角一定可以表示为(frac{2 pi}{n} cdot k)（弧度制），自乘(n)次后辐角为(2 pi k)，恰好绕圆周(k)圈后到点((1,0))

那么就表明(n)次单位根有且仅有(n)个，即单位圆周上的(n)等分点，一般的，我们把(n)次单位根记为(omega_{n}^0,omega_{n}^{1},dots,omega_{n}^{n-1})

举个例子，(n = 4)时，有(4)个单位根，如图

性质
- ((1)) (omega_{n}^{0} = 1)，显然。
- ((2)) (omega_{n}^{k} = left( omega_{n}^{1} ight) ^ k)更一般的((omega_{n}^{k}) ^ i = omega_{n}^{ik})，结合复数乘法的性质，辐角相加，也非常显然。
- ((3)) (omega_{n}^{i} cdot omega_{n}^{j} = omega_{n}^{i + j})
  
  由((2))知(omega_{n}^{i} = left( omega_{n}^{1} ight)^i, omega_{n}^{j} = left( omega_{n}^{1} ight)^j)，则(omega_{n}^{i} cdot omega_{n}^{j} = left ( omega_{n}^{1} ight)^{i+j} = omega_{n}^{i + j})
- ((4)) (omega_{n}^{k} = omega_{n}^{k mod n})，即绕了一圈又回来了。那么(omega_{n}^{-k} = omega_{n}^{n-k})。
- ((5)) 如果(n)是偶数，(omega_{n}^{k+n/2} = -omega_{n}^{k})，跑了半个圆周，刚好到相反方向。
- ((6)) (omega_{mn}^{mk} = omega_{n}^{k})，这个也挺显然的。
还有一些等用到了再说吧（其实鉴于复数乘法的优秀特性，这些性质都比较显然）。

好啦，前置姿势胡完了！

~~坐稳了！~~

DFT & IDFT

下面终于要讲一些有用的东西了。。。

DFT

上文也说了，(DFT)可以快速求多项式(F(x) = a_0x_0 + a_1x_1 + dots + a_{n-1}x^{n-1})（下面我们特别的规定(n)能表示成(2^k)）的一组点值表示。

上古时期，有一位名叫傅里叶的大佬，他发现把(omega_{n}^{0},dots,omega_{n}^{n-1})代入(F(x))，鉴于单位根的一些性质，可以很快的完成(DFT)求出一组(F(x))的点值表达。

为了方便，我们令(DFT(F(x)) = { (x_0,y_0),(x_1,y_1), dots,(x_{n-1},y_{n-1}) })，其中$x_i = omega_{n}^{i},y_i = F(x_i) $

我们也把(F(x))的一个(DFT)看做一个数组，并设(G = DFT(F))，有(G[k] = y_k)

傅里叶先把(F(x))按次数奇偶分组，得到两个小一点大小均为(n/2)的多项式

[f(x) = F[0]x_0 + F[2]x^1 + cdots + F[n-2]x^{n/2-1} ]
[g(x) = F[1]x_0 + F[3]x^1 + dots + F[n-1]x^{n/2-1} ]
有(F(x) = f(x^2) + xcdot g(x^2))。

发现如果(x = omega_{n}^{k} (k < n/2))的话

[egin{aligned}F(x) &= f((omega_{n}^{k}) ^ 2) + omega_{n}^{k} g((omega_{n}^{k}) ^ 2) \&= f(omega_{n/2}^{k}) + omega_{n}^{k} g(omega_{n/2}^{k})end{aligned} ]
我们设(f' = DFT(f), g' = DFT(g))，那么有(G[k] = f'[k] + omega_{n}^{k}g'[k] (k<n/2))

那如果代入(x = omega_{n}^{k + n/2}(k < n/2))呢？

[egin{aligned}F(x) &= f((omega_{n}^{k+n/2}) ^ 2) + omega_{n}^{k+n/2} g((omega_{n}^{k+n/2}) ^ 2) \&= f(omega_{n/2}^{k+n/2}) + omega_{n}^{k+n/2} g(omega_{n/2}^{k+n/2}) \&= f(omega_{n/2}^{k}) - omega_{n}^{k}g(omega_{n/2}^{k})end{aligned} ]
换句话说(G[k+n/2] = f'[k] - omega_{n}^{k}g'[k] (k < n/2))。

那么我们就有了两个柿子(k < n/2)：

[G[k] = f'[k] + omega_{n}^{k}g'[k] ]
[G[k+n/2] = f'[k] - omega_{n}^{k}g'[k] ]
我们只需要求出(f'[0...n/2-1])和(g'[0...n/2-1])就可以算出(G[0...n])了。

这样我们就可以分治算(G)就可以了。

啊，忘了。可以用三角函数来算出(omega_{n}^{1})，然后自乘(k)次就可以得到(omega_{n}^{k})了。

Code

所以代码：
```
#include <bits/stdc++.h>
using namespace std;
typedef double db;
const db PI=acos(-1.0);
const int N=4e6+10;
struct cpl{
	db x,y;
	cpl operator + (cpl k1)const{return (cpl){x+k1.x,y+k1.y};}
	cpl operator - (cpl k1)const{return (cpl){x-k1.x,y-k1.y};}
	cpl operator * (cpl k1)const{return (cpl){x*k1.x-y*k1.y,x*k1.y+y*k1.x};}
}buf[N];
void dft(cpl *f,int n){
	if (n==1)return;
	for (int i=0;i<n;i++) buf[i]=f[i]; // 别忘了保存下来。
	cpl *L=f,*R=f+(n>>1); // 直接用f的空间。
	for (int i=0;i<(n>>1);i++) L[i]=buf[i<<1],R[i]=buf[i<<1|1];
	dft(L,n>>1),dft(R,n>>1);
	cpl w=(cpl){1,0},w1=(cpl){cos(2*PI/n),sin(2*PI/n)}; // w是不是很像$omega$?
	for (int i=0;i<(n>>1);i++){ // 上文中的k
		buf[i]=L[i]+w*R[i];
		buf[i+(n>>1)]=L[i]-w*R[i];
		w=w*w1;		
	}
	for (int i=0;i<n;i++) f[i]=buf[i]; // 复制回来。
}
cpl f[N];
int main(){
	int n; scanf("%d",&n); // n次n-1项柿。
	for (int i=0;i<n;i++) scanf("%lf",&f[i].x); // 各项系数。
	int limit=1; while (limit<n) limit<<=1; // 注意补成n=2^k的形式
	dft(f,limit);
	for (int i=0;i<n;i++) printf("(%.3f,%.3f)
",f[i].x,f[i].y);
	return 0;
}
```
IDFT

但只有(DFT)还是不够啊。。。我们还要把它再搞回系数表示。

于是就有了(DFT)的逆变换$IDFT $。

~~记个结论：~~

(DFT)那里我们已经知道了(G = DFT(F))，且(G[k] = sumlimits_{i = 0}^{n-1} (omega_{n}^{k})^i F[i])。

结论就是(n cdot F[k] = sum_{i=0}^{n-1} (omega_{n}^{-k})^iG[i])，证明还是要有滴~

等式右边

[egin{aligned}sum_{i=0}^{n-1} (omega_{n}^{-k})^iG[i]&= sumlimits_{i = 0} ^ {n - 1} omega_{n}^{-ki} sumlimits_{j = 0}^{n-1} (omega_{n}^{i}) ^ j F[j] \&= sumlimits_{i = 0} ^ {n - 1} omega_{n}^{-ki} sumlimits_{j = 0}^{n-1} omega_{n}^{ji} F[j] \&= sumlimits_{j = 0} ^ {n - 1} F[j] sumlimits_{i=0}^{n-1} (omega_{n}^{j-k})^iend{aligned} ]
分类讨论一下

(j = k)时，因为(omega_{n}^{0} = 1)，所以那一项的贡献就是

[F[j] sumlimits_{i=0}^{n-1} (omega_{n}^{j-k})^i = F[j] cdot n cdot 1 = n F[k] ]
(j ot = k)呢？观察到第二个(sum)是个公比为(omega_{n}^{j-k})的等比数列。

求和一下

[sumlimits_{i=0}^{n-1} (omega_{n}^{j-k})^i = frac{(omega_{n}^{j-k})^n - (omega_{n}^{j-k})^0}{omega_{n}^{j-k} - 1} ]
由于((omega_{n}^{j-k})^n = omega_{n}^{n(j-k)} = omega_{n}^{0} = 1,(omega_{n}^{j-k})^0 = 1)，所以分子是(0)，这一项的贡献也是(0)。

因此(n cdot F[k] = sum_{i=0}^{n-1} (omega_{n}^{-k})^iG[i])。

这东西的本质好像是什么 ~~复读机~~单位根反演？不会。。。

这玩意儿直接在$DFT $的代码中改一个符号就好了。

Code
```
void idft(cpl *f,int n){
	if (n==1)return;
	for (int i=0;i<n;i++) buf[i]=f[i]; // 别忘了保存下来。
	cpl *L=f,*R=f+(n>>1); // 直接用f的空间。
	for (int i=0;i<(n>>1);i++) L[i]=buf[i<<1],R[i]=buf[i<<1|1];
	idft(L,n>>1),idft(R,n>>1);
	cpl w=(cpl){1,0},w1=(cpl){cos(2*PI/n), - sin(2*PI/n)}; // w是不是很像$omega$?
									// 就是这个负号！
    for (int i=0;i<(n>>1);i++){
		buf[i]=L[i]+w*R[i];
		buf[i+(n>>1)]=L[i]-w*R[i];
		w=w*w1;		
	}
	for (int i=0;i<n;i++) f[i]=buf[i]; // 复制回来。
}
```
main里面别忘了(IDFT)后把每一项的值除以一个(n)。

Code

其实我们完全可以把(DFT & IDFT)写一起，这样我们就得到了板子题的代码。
```
#include <bits/stdc++.h>
using namespace std;
typedef double db;
const db PI=acos(-1.0);
const int N=4e6+10;
struct cpl{
	db x,y;
	cpl operator + (cpl k1)const{return (cpl){x+k1.x,y+k1.y};}
	cpl operator - (cpl k1)const{return (cpl){x-k1.x,y-k1.y};}
	cpl operator * (cpl k1)const{return (cpl){x*k1.x-y*k1.y,x*k1.y+y*k1.x};}
}buf[N];
void fft(cpl *f,int n,int k1){ // k1==1 ? DFT : IDFT
	if (n==1)return;
	for (int i=0;i<n;i++) buf[i]=f[i];
	cpl *L=f,*R=f+(n>>1);
	for (int i=0;i<(n>>1);i++) L[i]=buf[i<<1],R[i]=buf[i<<1|1];
	fft(L,n>>1,k1),fft(R,n>>1,k1);
	cpl w=(cpl){1,0},w1=(cpl){cos(2*PI/n),k1*sin(2*PI/n)};
	for (int i=0;i<(n>>1);i++){
		buf[i]=L[i]+w*R[i];
		buf[i+(n>>1)]=L[i]-w*R[i];
		w=w*w1;		
	}
	for (int i=0;i<n;i++) f[i]=buf[i];
}
cpl f[N],g[N];
int main(){
	int n,m; scanf("%d%d",&n,&m);
	for (int i=0;i<=n;i++) scanf("%lf",&f[i].x);
	for (int i=0;i<=m;i++) scanf("%lf",&g[i].x);
	int limit=1; while (limit<=n+m)limit<<=1; // 补项
	fft(f,limit,1),fft(g,limit,1); // DFT
	for (int i=0;i<limit;i++) f[i]=f[i]*g[i]; // 点值相乘，求出答案的点值表示。
	fft(f,limit,-1); // IDFT回系数表示
	for (int i=0;i<=m+n;i++) printf("%d ",(int)(f[i].x/limit+0.5)); // 别忘了除以n
	return 0;
}
```
然后愉快的交一发：23333

淦。。。我写了一个假的(FFT)......

常数优化

我(O(n log n))的那么优秀的复杂度啊！咋(T)了呢？

冷静分析.jpg......

注意到板子题里(n le 10^6)，再加上我们要做两边(DFT)和一边(IDFT)，每次复杂度是(n log n)的，还有系统栈的巨大常数。。。

好吧，下面我们来卡卡常。

首先观察我们(DFT & IDFT)的过程。
```
第一层：0 1 2 3 4 5 6 7  (f数组的下标)
第二层：0 2 4 6|1 3 5 7 
第三层：0 4|2 6|1 5|3 7 
第四层：0|4|2|6|1|5|3|7 
```
发现原来(f)数组中的第(i)个数会出现在最后一层中第((i)的二进制翻转)个数中。

举个例子，(1:(001)_2)出现在了(4: (100)_2)上，(3: (011)_2)出现在了(6:(110)_2)上。

这样我们就可以把最后一层的情况求出来，然后再递归分治。

那怎么求最后一层呢？也就是求(i)的二进制翻转。（这个好像叫蝴蝶变换）

设补完项后的多项式长度为limit，令rev[i]表示(i)的二进制翻转，则rev[i]=rev[i>>1]>>1|((i&1)?limit>>1:0)，即不看(i)的最低位，把i>>1翻转过来后再把(i)的最低位拼上去。另外我们还能省几个数组。
```
#include <bits/stdc++.h>
using namespace std;
const int N=4e6+10;
typedef double db;
const db PI=acos(-1.0);
struct cpl{
	db x,y;
	cpl(db xx=0,db yy=0):x(xx),y(yy){}
	cpl operator + (cpl k1)const{return cpl(x+k1.x,y+k1.y);}
	cpl operator - (cpl k1)const{return cpl(x-k1.x,y-k1.y);}
	cpl operator * (cpl k1)const{return cpl(x*k1.x-y*k1.y,x*k1.y+y*k1.x);}
};
void fft(cpl *f,int n,int k1){
	if (n==1) return;
	fft(f,n>>1,k1),fft(f+(n>>1),n>>1,k1);
	cpl w(1,0),w1(cos(2*PI/n),sin(2*PI/n)*k1);
	for (int i=0;i<(n>>1);i++){
		cpl tmp=w*f[i+(n>>1)];
		f[i+(n>>1)]=f[i]-tmp; // 注意这里的赋值顺序！
		f[i]=f[i]+tmp;
		w=w*w1;
	}
}
int rev[N];
void change(cpl *f,int n){
	for (int i=0;i<n;i++) if (i<rev[i]) swap(f[i],f[rev[i]]);
}
cpl f[N],g[N];
int main(){
	int n,m; scanf("%d%d",&n,&m);
	for (int i=0;i<=n;i++) scanf("%lf",&f[i].x);
	for (int i=0;i<=m;i++) scanf("%lf",&g[i].x);
	int limit=1,len=0; while (limit<=n+m) limit<<=1,len++;
	for (int i=0;i<limit;i++) rev[i]=rev[i>>1]>>1|((i&1)?limit>>1:0);
	change(f,len),change(g,len);
	fft(f,limit,1),fft(g,limit,1);
	for (int i=0;i<limit;i++) f[i]=f[i]*g[i];
	change(f,len),fft(f,limit,-1);
	for (int i=0;i<=n+m;i++) printf("%d ",(int)(f[i].x/limit+0.5));
	return 0;
}
```
欸！那我们为什么不能直接自底向上合并呢？这样还省去了系统栈的常数。

最终Code
```
#include <bits/stdc++.h>
using namespace std;
const int N=4e6+10;
typedef double db;
const db PI=acos(-1.0);
struct cpl{
	db x,y;
	cpl operator + (cpl k1)const{return (cpl){x+k1.x,y+k1.y};}
	cpl operator - (cpl k1)const{return (cpl){x-k1.x,y-k1.y};}
	cpl operator * (cpl k1)const{return (cpl){x*k1.x-y*k1.y,x*k1.y+y*k1.x};}
};
int rev[N];
void fft(cpl *f,int n,int k1){
	for (int i=0;i<n;i++) if (rev[i]>i)swap(f[rev[i]],f[i]);
	for (int len=2;len<=n;len<<=1){
		cpl w1=(cpl){cos(2*PI/len),sin(2*PI/len)*k1};
		for (int i=0;i<n;i+=len){
			cpl w=(cpl){1,0};
			for (int j=i;j<i+(len>>1);j++){
				cpl tmp=w*f[j+(len>>1)];
				f[j+(len>>1)]=f[j]-tmp;
				f[j]=f[j]+tmp;
				w=w*w1;
			}
		}		
	}
}
cpl f[N],g[N];
int main(){
	int n,m; scanf("%d%d",&n,&m);
	for (int i=0;i<=n;i++) scanf("%lf",&f[i].x);
	for (int i=0;i<=m;i++) scanf("%lf",&g[i].x);
	int limit=1; for(;limit<=n+m;limit<<=1);
	for (int i=0;i<limit;i++) rev[i]=rev[i>>1]>>1|((i&1)?limit>>1:0);
	fft(f,limit,1),fft(g,limit,1);
	for (int i=0;i<limit;i++) f[i]=f[i]*g[i];
	fft(f,limit,-1);
	for (int i=0;i<=n+m;i++) printf("%d ",(int)(f[i].x/limit+0.5));
	return 0;
}
```
这样我们就能过模板题了。
查看全文

相关阅读:
js canvas登陆验证
 媒体查询
 js读取excel中日期格式转换问题
 jquery获取元素对应高度
 js引用类型的赋值
 asp.net core mvc视频A：笔记2-4.ActionResult（动作结果，即返回值）
asp.net core mvc视频A：笔记2-3.高级数据绑定
 asp.net core mvc视频A：笔记2-2.接收数据
 asp.net core mvc视频A：笔记2-1.控制器定义
 asp.net core mvc视频A：笔记1.基本概念介绍

原文地址：https://www.cnblogs.com/wxq1229/p/12232247.html

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

che dan环节：

概述

一些定义

总结

前置姿势

复数

定义

复平面：

模长：

辐角：

运算

加减法：

乘法

Code

单位根

定义

性质

DFT & IDFT

DFT

Code

IDFT

Code

Code

常数优化

最终Code