FFT | 1 | 基础傅立叶变换知识
引入问题
什么是多项式卷积
假设我们现在有多项式 $f(x)$ 和 $g(x)$ ，它可以被表示为
$$
f(x)=\sum_{i=0}^{n-1} a_i\cdot x^i\\
g(x)=\sum_{i=0}^{m-1} b_i\cdot x^i
$$
其中 $a$ 和 $b$ 为系数数组， $n$ 和 $m$ 分别为两个多项式的长度
那么它们的卷积为
$$
f(x)\bigotimes g(x)=\sum_{i=0}^{n-1} \sum_{j=0}^{m-1} a_i\cdot b_j\cdot x^{i+j}
$$
也可以表示成
$$
c_k=\sum_{i=0}^ka_i\cdot b_{k-i}
$$
其实就是简单的两个多项式相乘
如何快速求卷积？
为了方便表述，下文我们将假定两个多项式次数相等，即 $n=m$
显然，如果我们暴力枚举
1 | for (int i=0;i<n;i++) |
这样的复杂度是 $O(n^2)$ 的。
有没有什么突破口呢？
点值表达式
我们发现，如果我们仅仅是使用系数表达式是无法在时间上取得突破的，因此我们考虑使用点值表达式来表达一个多项式
多项式 $f(x)$ 和 $g(x)$ 的点值表达式可以用下面的式子来描述
$$
A={[x_0,f(x_0)],[x_1,f(x_1)],\cdots ,[x_{n-1},f(x_{n-1})]}\
B={[x_0,g(x_0)],[x_1,g(x_1)],\cdots ,[x_{n-1},g(x_{n-1})]}
$$
其中 $x$ 是任意的元素不重复的数组
众所周知，一个长度为 $n$ 的多项式可以由 $n$ 个多项式上取得的点来确定；那么有了这个 $A$ ，我们就可以唯一地确定多项式 $f(x)$ 了
我们发现，如果我们将两个多项式相乘，实际上就是将他们点值表达式中的 $y$ 值相乘
也就是说，如果
$$
h(x)=f(x)\bigotimes g(x)
$$
那么我们能得到 $h(x)$ 的点值表达式 $C$ 是
$$
C={[x_0,f(x_0)\cdot g(x_0)],[x_1,f(x_1)\cdot g(x_1)],\cdots ,[x_{n-1},f(x_{n-1})\cdot g(x_{n-1})]}
$$
这时我们就可以 $O(n)$ 实现点值多项式的卷积了！
然而问题并不简单
我们发现，尽管点值表示法的多项式卷积时间复杂度很低，但是我们无法在同样的时间复杂度内将一个系数表示的多项式转换为点值表示、或者是将点值表示的多项式转换为系数表示——实际上，朴素的系数转点值和点值转系数算法都是 $O(n^2)$ 的
那么我们现在的问题就转换为了
如何快速系数转点值/点值转系数？
$FFT$ 快速傅立叶变换
终于轮到今天的主角—— $FFT$ 出场了！
首先我们的傅立叶大神提出了单位根的概念：
$$
\omega_n^k=cos\frac{2\pi k}{n}+sin\frac{2\pi k}{n}\cdot i
$$
形象一点理解， $\omega_n^k$ 就是将复数平面上的单位圆平均分为 $n$ 份，然后从 $1+0i$ 开始逆时针旋转 $k$ 个部分得到的点
下图展示了 $n=8$ 时单位根 $\omega_8^k$ 的分布情况
有人问
这样不还是有 $n$ 个点需要分别代入计算吗？时间复杂度有减少吗？？？
傅立叶曰
等我秀一波操作
首先傅立叶展示了一堆易证的关于单位根的性质：
$$
(\omega_n^k)^t=w_n^{tk}\tag{1}
$$
$$
\begin{align}Proof:(\omega_n^k)^t&=(1,\frac{2\pi k}{n})^t\\
&=(1^t,t\cdot \frac{2\pi k}{n})\\
&=(1,\frac{2\pi (tk)}{n})\\
&=\omega_n^{tk}
\end{align}
$$
$$
\omega_n^k=\omega_{n/2}^{k/2}\tag{2}
$$
$$
\begin{align}
Proof:\omega_n^k&=(1,\frac{2\pi k}{n})\\
&=(1,\frac{2\pi \cdot(k/2)}{(n/2)})\\
&=\omega_{n/2}^{k/2}
\end{align}
$$
$$
\omega_n^{j+k}=\omega_n^j\cdot\omega_n^k\tag{3}
$$
$$
\begin{align}
Proof:\omega_n^{j+k}&=(1,\frac{2\pi (j+k)}{n})\\
&=(1\cdot 1,\frac{2\pi j}{n}+\frac{2\pi k}{n})\\
&=(1,\frac{2\pi j}{n})\cdot(1,\frac{2\pi k}{n})\\
&=\omega_n^j\cdot \omega_n^k
\end{align}
$$
$$
\omega_n^0=\omega_n^n=1+0i\\
\omega_n^{n/2}=-1+0i\tag{4}
$$
$$
w_n^{k+n/2}=-\omega_n^k\tag{5}
$$
$$
\begin{align}
Proof:\omega_n^{k+n/2}&=\omega_n^k\cdot \omega_n^{n/2}\\
&=-\omega_n^k
\end{align}
$$
然后傅立叶展现了他的操作（下文均要求 $n=2^k,k\in \mathbb{Z}$）
$$
\begin{align}
f(x)&=a_0+a_1x+a_2x^2+\cdots+a_nx^n\\
&=(a_0+a_2x^2+a_4x^4+\cdots+a_{n-1}x^{n-1})+(a_1x+a_3x^3+a_5x^5+\cdots+a_nx^n)
\end{align}
$$
我们令
$$
f1(x)=a_0+a_2x+a_4x+\cdots+a_nx^{(n+1)/2}\\
f2(x)=a_1+a_3x+a_5x+\cdots+a_{n-1}x^{(n+1)/2}
$$
那么易得
$$
f(x)=f1(x^2)+x\cdot f2(x^2)
$$
我们把之前的单位根 $\omega_n^k$ 代入这个式子中
$$
\begin{align}
f(\omega_n^k)&=f1[(\omega_n^k)^2]+\omega_n^k\cdot f2[(\omega_n^k)^2]\\
&=f1(\omega_n^{2k})+\omega_n^k\cdot f2(\omega_n^{2k})\\
&=f1(\omega_{n/2}^k)+\omega_n^k\cdot f2(\omega_{n/2}^k)
\end{align}
$$
傅立叶端详这个式子，随后又提笔
$$
\begin{align}
f(\omega_n^{k+n/2})&=f1[(\omega_n^{k+n/2})^2]+\omega_n^{k+n/2}\cdot f2[(\omega_n^{k+n/2})^2]\\
&=f1(\omega_n^{2k+n})+\omega_n^{k+n/2}\cdot f2(\omega_n^{2k+n})\\
&=f1(\omega_n^{2k})-\omega_n^k\cdot f2(\omega_n^{2k})\\
&=f1(\omega_{n/2}^k)-\omega_n^k\cdot f2(\omega_{n/2}^k)
\end{align}
$$
啊哈！出现了——分治！
两个式子都被展开成了 $f1(\omega_{n/2}^k)\pm \omega_n^k\cdot f2(\omega_{n/2}^k)$ 的形式
由于求 $f1(x)$ 、求 $f2(x)$ 和求 $f(x)$ 这三个问题的性质是相同的，并且在分治得到 $f1(\omega_{n/2}^k)$ 和 $f2(\omega_{n/2}^k)$ 后我们可以同时得出 $f(\omega_n^k)$ 和 $f(\omega_n^{k+n/2})$ ，那么系数转点值的问题就可以用 $O(n\cdot \log(n))$ 的复杂度来解决了！
等等：那么点值转系数呢？
$IFFT$ 快速傅立叶逆变换
与 $FFT$ 相反， $IFFT$ 是用来求系数转点值的算法
我们令 $A$ 为 $f(x)$ 的系数数组
我们假设我们现在已经得到了 $f(x)$ 的点值表达式 $D$ ：
$$
D={[\omega_n^0,f(\omega_n^0)],[\omega_n^1,f(\omega_n^1)],\cdots,[\omega_n^{n-1},f(\omega_n^{n-1})]}
$$
我们令 $Q$ 为 $D$ 的简单形式：
$$
Q={f(\omega_n^0),f(\omega_n^1),\cdots,f(\omega_n^{n-1})}
$$
我们可以得出
$$
\left[\begin{matrix}
(\omega_n^0)^0&(\omega_n^0)^1&(\omega_n^0)^2&\cdots&(\omega_n^0)^{n-1}\
(\omega_n^1)^0&(\omega_n^1)^1&(\omega_n^1)^2&\cdots&(\omega_n^1)^{n-1}\
\vdots&\vdots&\vdots&\ddots&\vdots\
(\omega_n^{n-1})^0&(\omega_n^{n-1})^1&(\omega_n^{n-1})^2&\cdots&(\omega_n^{n-1})^{n-1}
\end{matrix}\right]
\cdot
\left[\begin{matrix}
A_0\A_1\ \vdots\A_{n-1}
\end{matrix}\right]
\left[\begin{matrix}
Q_0\Q_1\ \vdots\Q_{n-1}
\end{matrix}\right]
$$
也就是说
$$
\left[\begin{matrix}
A_0\A_1\ \vdots\A_{n-1}
\end{matrix}\right]
\left[\begin{matrix}
Q_0\Q_1\ \vdots\Q_{n-1}
\end{matrix}\right]
\div
\left[\begin{matrix}
(\omega_n^0)^0&(\omega_n^0)^1&(\omega_n^0)^2&\cdots&(\omega_n^0)^{n-1}\
(\omega_n^1)^0&(\omega_n^1)^1&(\omega_n^1)^2&\cdots&(\omega_n^1)^{n-1}\
\vdots&\vdots&\vdots&\ddots&\vdots\
(\omega_n^{n-1})^0&(\omega_n^{n-1})^1&(\omega_n^{n-1})^2&\cdots&(\omega_n^{n-1})^{n-1}
\end{matrix}\right]\tag{6}
$$
众所周知，除以一个矩阵等于乘上这个矩阵的逆矩阵
那么怎么求出右边那个矩阵（姑且称其为 $ W$ ）的逆矩阵呢？我也不知道
最后得出来的结论就是
$$
W^{-1}=
\left[\begin{matrix}
1/(n\cdot W_{0,0})&1/(n\cdot W_{0,1})&\cdots&1/(n\cdot W_{0,n-1})\
1/(n\cdot W_{1,0})&1/(n\cdot W_{1,1})&\cdots&1/(n\cdot W_{1,n-1})\
\vdots&\vdots&\ddots&\vdots\
1/(n\cdot W_{n-1,0})&1/(n\cdot W_{n-1,1})&\cdots&1/(n\cdot W_{n-1,n-1})
\end{matrix}\right]
$$
就是把 $W$ 的每一个元素取倒数然后乘上 $1/n$
又因为
$$
\begin{align}
W_{i,j}&=(\omega_n^i)^j\
&=\omega_n^{ij}
\end{align}
$$
所以
$$
\begin{align}
W_{i,j}^{-1}&=(\omega_n^{ij})^{-1}/n\
&=\omega_n^{-ij}/n
\end{align}
$$
最后我们代回 $(6)$ 式，再进行一些变换，得到
$$
\left[\begin{matrix}
A_0\A_1\ \vdots\A_{n-1}
\end{matrix}\right]
\left[\begin{matrix}
Q_0\Q_1\ \vdots\Q_{n-1}
\end{matrix}\right]
\cdot
\left[\begin{matrix}
(\omega_n^0)^0&(\omega_n^0)^1&(\omega_n^0)^2&\cdots&(\omega_n^0)^{n-1}\
(\omega_n^{-1})^0&(\omega_n^{-1})^1&(\omega_n^{-1})^2&\cdots&(\omega_n^{-1})^{n-1}\
\vdots&\vdots&\vdots&\ddots&\vdots\
(\omega_n^{-(n-1)})^0&(\omega_n^{-(n-1)})^1&(\omega_n^{-(n-1)})^2&\cdots&(\omega_n^{-(n-1)})^{n-1}
\end{matrix}\right]\cdot \frac{1}{n}
$$
发现什么了吗？
我们只需要对 $Q$ 再做一次伪·$FFT$：也就是将 $FFT$ 过程中的所有的 $\omega_n^k$ 都替换为 $\omega_n^{-k}$ ，然后把结果数组所有元素都除以 $n$ 就好了！
说了这么多，让我们来看看代码吧？
代码实现
给定 $n$ 次多项式 $F(n)$ 和 $m$ 次项多项式 $G(n)$，求两个多项式的卷积
$(n, m\le10^6)$
首先是一个简单易懂的虚数结构
1 | struct Complex { |
然后就是我们期待已久的 $FFT$ 啦！
1 | Complex tmp[nmax]; |
然后是主程序部分
1 | cin>>n>>m; |
好的！让我们把这道模板题A掉吧！
…等等？
？？$O(n\cdot \log(n))$时间复杂度为何过不了？
因为递归常数太大了
蝴蝶变换优化
啊哈！玄学优化
请看下图
这张图展示了对一个长度为8的多项式进行 $FFT/IFFT$ 的过程
我们只需要能预先求得下面的后序列就可以用迭代的方式
由此我们发现新的序列其实是把原来的序列按位反转了，因此我们只需要处理出这个 $rev$ 数组就好了
通过观察，我们不难得出转移公式
$$
rev[i]=\left{\begin{align}0&&i=0\rev[i>>1]>>1&&i%2=0\(rev[i>>1]>>1)\text|(len>>1)&&otherwise.\end{align}\right.
$$
然后我们只需要枚举当前的 $mid$ 就好了
1 |
|
例题
P1919 - A*B Problem升级版
给定两个位数为 $n$ 的数 $A$ 和 $B$ ，求出这两个数的乘积
$n\le 60000$
这道题其实就是裸的 $FFT$
可以把 $A$ 和 $B$ 理解为两个多项式 $f(x)$ 和 $g(x)$ 在 $x=10$ 时的值，那么 $f(x)$ 和 $g(x)$ 从低到高的系数就分别对应 $A$ 和 $B$ 从低位到高位的数字
举个例子，假定我们有 $A=19260817$ ，那么
$$
f(x)=1+9x+2x^2+6x^3+0x^4+8x^5+1x^6+7x^7
$$
我们只需要对 $f(x)$ 和 $g(x)$ 求个卷积后输出就好了
提示
- 由于读入时是先读入高位的，所以程序处理中也应该从高位到序赋值
- 最后求出来的结果系数可能会出现大于10的情况，要记得取余
代码
实际上我们只需要将上一题的主程序部分改一下
1 | int main() { |
伪·A+B Problem
啊哈！这道题并不是一般的 A+B Problem
给定一个长度为 $n$ 的数组 $a$ ，要求出有多少个有序三元数对 $(i,j,k)$ 使得 $a_i+a_j=a_k$
$n\le 200000, -50000\le a_i\le 50000$
由于可能有重复的数字出现，且数据范围有限，我们可以开一个数组 $q$ ， $q_i$ 代表数字 $i$ 在数组 $a$ 中出现的次数
再开一个数组 $t$ ，$t_k$ 表示满足 $a_i+a_j=a_k$ 的方案数
因此我们可以明显得到
$$
t_k=\sum_{i=-50000}^{50000}q_i\cdot q_{k-i}
$$
再来看卷积的定义
$$
c_k=\sum_{i=0}^ka_i\cdot b_{k-i}
$$
这不就是一个卷积的形式吗？$FFT$ 切了就好
提示
- 考虑到数组下标会出现负数，所以我们读入时统一加上50000就好了
- 卷积中不会排除 $i=j$ 的情况，所以我们需要从得出的答案中手动减出去
P3338 - [ZJOI2014]力
给定长度为 $n$ 的数组 $q$ ，令
$$
E_i=\sum_{j<i}\frac{q_j}{(j-i)^2}-\sum_{j>i}\frac{q_j}{(j-i)^2}
$$求出数组 $E$
$n\le 100000$
我们来推一下这个式子（数组下标从0开始）
$$
\begin{align}
E_i&=\sum_{j<i}\frac{q_j}{(j-i)^2}-\sum_{j>i}\frac{q_j}{(j-i)^2}\\
&=\sum_{j<i}\left(q_j\cdot\frac{1}{(j-i)^2}\right)-\sum_{j>i}\left(q_j\cdot\frac{1}{(j-i)^2}\right)
\end{align}
$$
令
$$
w_i=i^{-2}\
q’i=q{n-i-1}
$$
那么原式变为
$$
\begin{align}
E_i&=\sum_{j<i}q_j\cdot w_{i-j}-\sum_{j>i}q_j\cdot w_{j-i}\\
&=\sum_{j=0}^{i-1}q_j\cdot w_{i-j}-\sum_{j=1}^{n-j-1}q_{i+j}\cdot w_j\\
&=\sum_{j=0}^{i-1}q_j\cdot w_{i-j}-\sum_{j=1}^{n-j-1}q’_{(n-i-1)-j}\cdot w_j
\end{align}
$$
好的！原式被完全转换为了卷积的形式，可以用 $FFT$ 来求了QwQ
提示
- 尽管我们最后并不会用到 $E$ 的 $n$ 或者以后的下标，但我们在决定 $len$ 的时候还是应该将 $len$ 定位 $2n$ ，因为 $FFT$ 在计算时会用到那后面的东西
（还是写一下代码吧
代码实现
1 |
|
更多例题
暑假培训Day31 C 逃
FFT | 1 | 基础傅立叶变换知识
