快速傅里叶变换(FFT)
FFT 是之前学的,现在过了比较久的时间,终于打算在回顾的时候系统地整理一篇笔记,有写错的部分请指出来啊 qwq。
更好的阅读体验:
https://www.cnblogs.com/Tenshi/p/15434004.html
ACWING 的编辑器的矩阵怎么渲染不出来啊(哭
upd:感谢 @垫底抽风 qwq,它可以了!acwing 编辑器写矩阵和 typora 不太一样啊awa。
卷积
卷积、旋积或褶积(英语:Convolution)是通过两个函数 $f$ 和 $g$ 生成第三个函数的一种数学算子。
定义
设 $f,g$ 在 $R1$ 上可积,那么 $h(x) = \int_{-∞}^∞f(\tau)g(x-\tau)d\tau$ 称为 $f$ 与 $g$ 的卷积。
对于整系数多项式域,$n-1$ 次多项式 $A,B$ 的相乘可以得到 $h(x) = A(x)B(x) = \sum_{i=0}^{2n-2} \sum_{j+k=i}a_j b_k = \sum_{i=0}^{2n-2} \sum_{j=0}^i a_j b_{i-j} x^i$,对应的卷积为 $C_i = \sum A_j B_{i-j}$。
系数表示法
即用多项式各项系数来刻画这个多项式,例如 $n-1$ 次多项式就可以写成这样:$A(x) = a_0 + a_1x+…+a_{n-1}x^{n-1}$
点值表示法
我们知道,$n$ 个不同的点可以确定一个 $n-1$ 次的多项式,所以我们可以使用 $n$ 个(不同)点来刻画一个 $n-1$ 次多项式。
这样做会有什么方便呢?
例如 $f(x) = (x_0, f(x_0)),…(x_n,f(x_n)),g(x) = (x_0, g(x_0)),…(x_n, g(x_n))$,那么它们的卷积 $h(x) = (x_0, f(x_0)g(x_0)),…(x_n, f(x_n)g(x_n))$。
这意味着在系数表示法中需要 $O(n^2)$ 次的乘法运算在点值表示法中只需要 $O(n)$ 次。
系数表示法转点值表示法(DFT)
下面考虑如何将 $n-1$ 次多项式从系数表示法转为点值表示法。
因为用普通的方法选取 $n$ 个点然后将系数表示法转为点值表示法的复杂度为 $O(n^2)$(因为需要选 $n$ 个点,然后对于每个点 $x$ 需要计算共 $n$ 项的结果),我们考虑如何优化这一步。
注意到满足 $w^n=1$ 的单位根 $w$ 有 $n$ 个,故从这里入手。
我们记方程 $w^n = 1$ 的第 $k$ 个单位根为 $w_n^k$。
方便起见,设 $n$ 为 $2$ 的幂(就算不是也可以看作是高次项的系数为 $0$)。
将 $A(x)$ 按照次数的奇偶性分别分成两组 $F(x),G(x)$,并表示为 $A(x) = F(x^2) + xG(x^2)$
例如 $A(x) = a_0+a_1x+a_2x^2+a_3x^3$,那么 $F(x) = a_0+a_2x,G(x)=a_1+a_3x$。
将 $x=w_n^k$ 代入 $A(x)$,由复数的性质,
$A(w_n^k) = F(w_{\frac{n}{2}}^k) + w_n^k G(w_{\frac{n}{2}}^k)$ ,类似地 $A(w_{n}^{k+\frac{n}{2}}) = F(w_{\frac{n}{2}}^k) - w_n^k G(w_{\frac{n}{2}}^k)$。
推导:
$A(w_{n}^{k+\frac{n}{2}}) = F(w_n^{2k+n}) + w_n^{k+\frac{n}{2}}G(w_n^{2k+n}) \\\= F(w_{\frac{n}{2}}^k) + w_n^{k+\frac{n}{2}} G(w_{\frac{n}{2}}^k) \\\= F(w_{\frac{n}{2}}^k) - w_n^k G(w_{\frac{n}{2}}^k)$
可以发现对于两个相应的单位根 $w_n^k,w_n^{\frac{n}{2}+k}$,可以用对应的 $F,G$ 算出(可以递归地实现这个过程),而且计算的范围折半,所以一共需要计算 $O(logN)$ 层,每一层执行 $O(n)$ 次运算,所以复杂度为 $O(NlogN)$。
点值表示法转系数表示法(IDFT)
下面考虑如何将 $n-1$ 次多项式从点值表示法转为系数表示法。
因为对于每个点值 $y_i = \sum_{k=0}^{n-1}w_n^{ki}$,其中 $i\in[0, n-1]$,我们可以写出等式:
$ \left\[ \begin{matrix} 1 & 1 & 1 & \cdots & 1\\\ 1 & w_n^1 & w_n^2 & \cdots & w_n^{n-1}\\\ 1 & w_n^2 & w_n^4 & \cdots & w_n^{2(n-1)}\\\ \vdots & \vdots & \vdots & \ddots & \vdots \\\ 1 & w_n^{n-1} & w_n^{2(n-1)} & \cdots & w_n^{(n-1)(n-1)} \end{matrix} \right\] \left\[ \begin{matrix} a_0\\\ a_1\\\ a_2\\\ \vdots\\\ a_{n-1} \end{matrix} \right\]=\left[ \begin{matrix} y_0\\\ y_1\\\ y_2\\\ \vdots\\\ y_{n-1} \end{matrix} \right] $
现在我们已经有向量 $y$ 了(就是右式),因此,如果要得到向量 $a$,只需要两边乘上 $w$ 矩阵的逆即可。
这里的 $w$ 矩阵正是著名的范德蒙矩阵,它的逆正好是每一项都取倒数,然后除以 $n$。
因此有 $a_i = \frac{1}{n}\sum_{k=0}^{n-1}w_n^{-ki}$,其中 $i\in[0, n-1]$。
有没有发现 $a_i,y_i$ 的形式非常接近?据此,我们可以在实现的时候在同一个函数中写出逆变换和正变换,然后在得到的结果 res
中除以 $n$ 就可以了。(参照下面的代码)
至此,FFT 的基本原理讲述完毕,下面是优化。
位逆序置换
按照上文的讲述,如果不看下面的代码,那么编写出来的是递归版本,但是这个版本的常数太大了,因此运行起来的效果不好,故使用位逆序置换来降低常数。
我们看看递归过程是什么样的,以 $n=8$ 为例:
这里就有一个非常神奇的规律:在最后一行中,原下标所对应的二进制数翻转正好是在最后一行的序数。例如 $x_6$ 的下标是 $6 = 110_{(2)}$,那么它的序数正好是 $011_{(2)} = 3$。
据此,可以处理出 rev
数组,它记录的正是最后一行所有元素对应的下标。
简单地说,递归形式是自上而下地做 FFT,而利用位逆序置换我们可以自下而上地做 FFT,它们在实际运行中有着常数上的区别。
模板题及代码
https://www.luogu.com.cn/problem/P3803
给定一个 $n$ 次多项式 $F(x)$,和一个 $m$ 次多项式 $G(x)$。
请求出 $F(x)$ 和 $G(x)$ 的卷积。
#include<bits/stdc++.h>
using namespace std;
const int N=3e5+5;
const double pi=acos(-1);
int n, m;
// 复数类
struct Complex{
double x, y;
Complex operator + (const Complex &o)const { return {x+o.x, y+o.y}; }
Complex operator - (const Complex &o)const { return {x-o.x, y-o.y}; }
Complex operator * (const Complex &o)const { return {x*o.x-y*o.y, x*o.y+y*o.x}; }
};
Complex a[N], b[N];
int res[N];
int rev[N], bit, tot;
void fft(Complex a[], int inv){ // inv 指示正变换、逆变换。
for(int i=0; i<tot; i++) if(i<rev[i]) swap(a[i], a[rev[i]]);
for(int mid=1; mid<tot; mid<<=1){
auto w1=Complex({cos(pi/mid), inv*sin(pi/mid)});
for(int i=0; i<tot; i+=mid*2){
auto wk=Complex({1, 0});
for(int j=0; j<mid; j++, wk=wk*w1){
auto x=a[i+j], y=wk*a[i+j+mid];
a[i+j]=x+y, a[i+j+mid]=x-y;
}
}
}
}
int main(){
cin>>n>>m;
for(int i=0; i<=n; i++) cin>>a[i].x;
for(int i=0; i<=m; i++) cin>>b[i].x;
while((1<<bit)<n+m+1) bit++; // 结果次数分布在 [0, n+m] 内,一共有 n+m+1 位。
tot=1<<bit; // 得到上文所说的 n
for(int i=0; i<tot; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
fft(a, 1), fft(b, 1); // 正变换 DFT
for(int i=0; i<tot; i++) a[i]=a[i]*b[i];
fft(a, -1); // 逆变换 IDFT
for(int i=0; i<=n+m; i++) res[i]=(int)(a[i].x/tot+0.5), printf("%d ", res[i]);
return 0;
}
位逆序最高位预处理可以用
__lg(n + m)
实现学到了!翎羽 Orz
天子姐姐!!
天子姐姐!!
qwq
acwing的矩阵挺奇怪的,具体可以看看这里{:target=”_blank”}。
谢谢抽风dalao,成功让它显示出来了