彼岸弗弗欢
彼岸弗弗欢
其实是傅里叶变换。
因为发现自己之前学得像阿罗拉臭臭泥,所以再学一遍,写写之前没注意到的。
会有很多毫无意义的瞎扯,以及大量的对引用文章的解释,配合我东拼西凑的理解。
这东西就当个乐子看,学的话还是看文末好文共赏的那几篇,卡住了可以过来找找灵感。
过两天可能还会写份详细版,这个东西就一个雏形。
一些不严谨的约定
- 序列 \(f\) 的长度为 \(n\),则意为序列 \(f\) 由 \(\left\{f_0,f_1,\cdot,f_{n-1}\right\}\) 组成。
- \(F(x)\) 表示一形式幂级数 \(\sum_{i\ge0}{f_ix^i}\),用 \(F_i\) 表示 \(F(x)\) 的第 \(i\) 次项系数。
- \(F(c)\) 表示 \(F(x)\) 在 \(x=c\) 时的值。
- \(\operatorname{DFT}(F)\) 表示对 \(F(x)\) 进行傅里叶变换之后的序列,相应的 \(\operatorname{IDFT}(f)\) 表示对 \(f\) 进行傅立叶逆变换后的形式,即 \(F(x)\) 本身。
- 一般不区分 FFT 和 NTT,即不考虑数域。
- 卷积指代加法卷积,用 \(f*g\) 表示 \(f\) 与 \(g\) 进行卷积。
- 单位根即本元单位根,既可指单位根,也可指原根,事实上这俩玩意儿是一个东西,只是数域不同。
- 多项式和系数序列会突然转化,标志(大概)是字母的大小写改变。
李三傅里叶变换
其实是离散傅里叶变换。
我们回到 FFT 的起点,我们的思路是利用高效的多项式点值乘法进行加速,于是需要对输入多项式先进行求值,再把结果插值得到我们想要的系数表示(因为系数表示更常用,比如可以看作形式幂级数)。
其中求值部分被称为 DFT,即离散傅里叶变换,插值则为 IDFT。若不考虑 FFT 的加速,二者的原始做法应该是这样:
- \(f_k=\operatorname{DFT}(F)_k=\sum\limits_{i=0}^{n-1}{F_i(\omega_n^k)^i}\);
- \(F_k=\operatorname{IDFT}(f)_k=\dfrac 1 n\sum\limits_{i=0}^{n-1}{f_i(\omega_n^k)^{-i}}\).
可以用单位根反演证明二者为互逆操作,而 DFT 后乘法的正确性可以简单地展开式子来说明。
不难发现 \(f_k=F(\omega_n^k)\),也就是说 DFT 可以看作是 \(F(x)\) 对序列 \(\left\{\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}\right\}\) 求值的结果。
IDFT 可以看作将 \(f\) 视作多项式后对序列 \(\left\{\omega_n^{-0},\omega_n^{-1},\cdots,\omega_n^{-(n-1)}\right\}\) 求值再除以 \(n\) 的结果,当然这个说法不严谨。
快速彼岸欢
其实是快速傅里叶变换。
基于求值这个想法,我们可以推出另一种 FFT 的实现\(^{\left[1\right]}\),当然依旧假定 \(n=2^k\)。
先提点前置知识。
也就是多项式取模的一点性质,首先模完了是常数这点显然,然后写开来是 \(F(c)=F(x)-P(x)\times(x-c)\) 考虑把 \(x\) 这个占位符换成 \(c\) 就得到。
同时还有:
这个证明就是单纯写开了。
其实 \(\left(1\right)\) 和 \(\left(2\right)\) 在多项式多点求值的传统方法里面也用上了。
考虑我们求 \(F(\omega_n^{0}),F(\omega_n^{1})\cdots,F(\omega_n^{n-1})\),也就可以变成求 \(F(x)\bmod(x-\omega_n^{0}),F(x)\bmod(x-\omega_n^{1}),\cdots,F(x)\bmod(x-\omega_n^{n})\)。
那我们直接用上多点求值那一套,做成分治树的形式,使用多项式取模那一套做到 \(\mathcal{O}(n\lg^2n)\)。
考虑单位根的好性质,我们可以尝试将单位根进行一定的分组凑出 \((x-\omega_n^k),(x-\omega_n^{k+\frac{n}{2}})=(x+\omega_n^k)\),然后它们的乘积就是 \((x^2-{\omega_n^k}^2)=(x^2-\omega_n^{2k})=(x^2-\omega_{n/2}^k)\),由于单位根的好性质一定都能凑出来,故不赘述。
然后我们这么一路合上去,每层的多项式数量都能折半,而顶层就是 \(x^n-1\),然后不知道有什么用。
不过我们考虑倒过来做,假定我们目前得到了 \(F(x)\bmod(x^{2n}-r^2)\),目的是求 \(F(x)\bmod(x^n-r)\) 与 \(F(x)\bmod(x^n+r)\),根据 \(\left(2\right)\) 中的 Lemma 这样做是合法的。
以 \(F(x)\bmod(x^n-r)\) 为例,考虑 \(F(x)=\sum_{i\ge0}{F_ix^i}\),对其前 \(n\) 项取模无影响,而 \(\forall k\in\left[0,n\right),F_{k+n}x^{k+n}\bmod(x^n-r)=F_{k+n}x^n\times x^k\bmod(x^n-r)=(F_{k+n}(x^n-r)+F_{k+n}r)\times x^k\bmod(x^n-r)=F_{k+n}(x^n-r)\times x^k+F_{k+n}r\times x^k\bmod(x^n-r)=F_{k+n}r\times x^k\),也就是第 \(k+n\) 项的系数贡献到了第 \(k\) 项。
同理对于 \(F(x)\bmod(x^n+r)\),\(\forall k\in\left[0,n\right),F_{k+n}x^{k+n}\bmod(x^n+r)=-F_{k+n}r\times x^k\)。
做法于是呼之欲出了。
考虑我们沿用分治的策略,当前层有 \(F(x)\bmod (x^{2n}-r^2)\) 的结果,而我们每次用 \(0,1,\cdots.n-1\) 项的位置存放 \(F(x)\bmod(x^n-r)\) 的结果,而 \(n,n+1,\cdots,2n-1\) 的位置存放 \(F(x)\bmod (x^n+r)\) 的结果,然后两部分显然互不影响,可以直接分别丢下去分治。
啥你疑惑为什么 \(F(x)\bmod(x^n+r)\) 跟我们 \(F(x)\bmod(x^{2n}-r^2)\) 可以是一个形式,emm……虽然在 \(\mathbb{Z}\) 上确实不一定找不到一个 \(t\) 使得 \(t^2=-r\)(当然 \(r<0\) 时可以,但此时找不到 \(x^n-r\) 对应的 \(t^2=r\)),但是我们可以扩域到 \(\mathbb{C}\) 或者 \(\mathbb{F_p}\)(前者可以开虚根,后者考虑解二次剩余)甚至一些奇奇怪怪的数域,总之能搞。
在最顶层模 \(x^n-1\),这样在底就可以得到一系列点值,并且可以发现分别对应了某一个 \(\omega_n^k\),而 \(k\) 次单位根对应到的位置恰好是按位逆序置换后的位置。
大概原理是这样,没看懂的话去原文看图,我懒得画了虽然链接那篇文章也是用的原论文的图。
实质上这个做法又称 DIF-FFT,而我们以前的写法又称 DIT-FFT,它们在信号分析中有不同的意义(时域和频域),但是我也不懂。
另外有讲义说 DIF 和 DIT 的做法互为转置,但我不懂转置原理,所以只提一嘴。
值得指出的是考虑 DIT-FFT 的迭代版是先位逆序置换,然后按 \(0,1,\cdots,n-1\) 顺序输出,而恰好 DIF-FFT 是按 \(0,1,\cdots,n-1\) 顺序输入,然后位逆序置换输出,结合二者本质相同的特点,我们可以用 DIF-FFT 作 DFT,然后直接用 DIT-FFT 作 IDFT,最后翻转一次序列(因为位逆序置换过的原因,但注意 \(0\) 位置是不动的),这样便能省去手动的位逆序置换,减小常数。
同时我们还能对单位根的预处理进行调整,使得缓存友好来减小常数,然而暂时不想写。
闲行卷积
其实是线性卷积。
即考虑两数列 \(f\) 和 \(g\),令 \(h=f*g\),则 \(h_k=\sum\limits_{i,j}{\left[i+j=k\right]f_i\times g_j}=\sum\limits_i{f_i\times g_{k-i}}\)。
容易发现这和多项式乘法一个形式,使用 FFT 加速,于是线性卷积是好做的。
卷积寻欢
其实是循环卷积。
仍即考虑两数列 \(f\) 和 \(g\),它们的长度均为 \(n\),令 \(h=f*g\),则 \(h_k=\sum\limits_{i,j}{\left[i+j\bmod n=k\right]f_i\times g_j}\)。
容易发现这和线性卷积差不多形式,做完合并贡献,于是循环卷积是好做的。
但是有更好的方法,考虑我们单位根的性质,会发现 \(n\) 个 \(n\) 次单位根是均匀分布且循环的,同时 \(\omega_n^{k}=\omega_n^{n+k}\)。
而我们之前做卷积时的方法是:对于长度分别为 \(n\) 和 \(m\) 的 \(f\) 和 \(g\),我们找到一个 \(l=2^t\) 使得 \(l\ge n+m-1\),然后点乘后再转回系数。
那如果我们的 \(\max\left\{n,m\right\}<l<n+m-1\) 会怎么样?考虑我们单位根的循环性质,不难猜到会出现重复贡献。
具体的,考虑 \(h=f*g\),假定 \(f\) 和 \(g\) 长度均不超过 \(n\) 且卷积使用 \(n=2^t\) 长度,然后把这个按照原始定义写开:
式子很长,但干的事情很简单,就是把 \(H\) 写成了 \(F*G\) 的形式,并且把 DFT 展开。
不难发现当 \(n\) 的长度不够时,\(k_1+k_2\) 就会贡献到 \((k_1+k_2)\bmod n\) 而非 \(k_1+k_2\),但由于我们此前一直保证了 \(n\) 的长度足够,因此我们并未发现自己做的是循环卷积。
当然我们一般的循环卷积仍能用线性卷积代替,但问题是一些乱七八糟的科技比如循环卷积快速幂里面需要直接卷积这个性质(主要是能求逆)来套回线性卷积的 \(\mathcal{O}(n\lg n)\) 做法,不然只能暴力快速幂俩 \(\lg\)。
现在我们知道直接用 FFT 可以做 \(n=2^t\) 的循环卷积,但对于任意长度的循环卷积,我们暂时无法解决。
于是请出我们的 Bluestein 或者说 Chirp Z 或者说 CZT,总之就是任意长度循环卷积。
这个东西的得到比较暴力,就直接单位根反演硬写式子,最后一步人类智慧给他换成线性卷积的样子。
啥你说最后还有一部卷积那它能干嘛?事实上尽管我们最后仍然用了 FFT 做线性卷积,但是我们因为构造了一摊系数于是最后线性卷积得到的值会是原序列循环卷积的值。
不过这个系数我还没推过,所以先鸽了,所以这部分配合好文共赏 \(\left[2\right]\sim\left[3\right]\)。
肘积好文共赏
其实肘鸡是一位硕人。
\(\left[1\right]\):https://loj.ac/d/3165.
\(\left[2\right]\):https://www.luogu.com.cn/blog/Troverld/fft-li-xing-xia-che.
\(\left[3\right]\):https://www.luogu.com.cn/blog/105254/qian-tan-fft-zong-ft-dao-fft.
\(\left[4\right]\): https://charleswu.site/archives/3065.
注意 \(\left[4\right]\) 中提到的几篇文章都值得看一看。