Re0 从零开始的多项式之旅

1|0Re0 从零开始的多项式之旅

1|1前の言

1|0温馨提示

1. 本文可能不适合 0 基础者观看，因为文章涉及一定的高中数学内容并未特殊说明（主要是懒）；
2. 本文可能存在节奏过快的问题，但未提及部分均为读者自推；
3. 本文还在更新中；
4. 本文可能存在 $LATEX$ 格式书写问题，希望能在评论区中指出批评。

1|0闲话

已经很久没有碰多项式了，一道突然的省选模拟题让死去的回忆复现，这使我充满了决心。当然，可能未曾在此文涉及题目，需要咕咕咕。

1|2快速傅里叶变换

1|0DFT 与 IDFT 的孽缘

多项式算法的最初研究目的就是解决多项式卷积。

但是，直接进行卷积操作不易优化（但可以通过分治——即 Karatsuba 算法，时间复杂度 $\mathcal{O}(n^{{\log }_{2}3})$），我们想到可以通过一个媒介转化。很显然的，点值表示法是不错的选择，它符合大众的认知（毕竟在初中时期我们就通过待定系数法用三点确定一条二次函数）。

至此，我们就正式引入了 $\mathtt{\text{DFT}}$（离散傅里叶变换，Discrete Fourier Transform）和 $\mathtt{\text{IDFT}}$（离散反傅里叶变换，Inverse Discrete Fourier Transform）。

1|0葬送的单位根

为快速转化，我们将选取一些特殊的点。

warning: 后文将涉及复数和原根内容，如有需要请移步 OI-Wiki

对于 $\mathtt{\text{FFT}}$ 而言，它所选取的是 $w_n^k=e^{\frac{2k\pi }{n}i}$，其中 $k$ 从 $0$ 取到 $n-1$，可表示最高次为 $n-1$ 的多项式。

对于 $\mathtt{\text{NTT}}$ 而言，它所选取的是 $w_n^k=g^{\frac{mod-1}{n}k}$，其中 $g$ 为原根， $k$ 从 $0$ 取到 $n-1$，可表示最高次为 $n-1$ 的多项式。

1|0单位根的性质

1. $w_n^{x+y}=w_n^x{w}_n^y$
2. $w_{2n}^{2k}=w_n^k$

有了如上基础，就可以来推式子了。

1|0DIT 和DIF 的超距作用

我们令 $F(k)=\sum _{i=0}^{n-1}{a}_i{w}_n^{ik}$。容易发现，$F(k)$ 所表示的就是第 $k+1$ 个我们所选取的点的函数值。

为方便分析，后文的 $n$ 为 $2$ 的幂次（实践中也需要先转成 $2$ 的幂次）。

1|0DIT： 按时域抽取（Decimation-In-Time）

$$\begin{array}{}\text{(1)}& F(k)& =\sum _{i=0}^{n-1}{a}_i{w}_n^{ik}\text{(2)}& & =\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i}{w}_n^{2ik}+\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{w}_n^{2ik+k}\text{(3)}& & =\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i}{w}_{\frac{n}{2}}^{ik}+w_n^k\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{w}_{\frac{n}{2}}^{ik}\end{array}$$

发现，当 $k<\frac{n}{2}$ 时，我们令 $k^{\prime }=k+\frac{n}{2}$，则有

$$\begin{array}{}\text{(4)}& F(k^{\prime })& =\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i}{w}_{\frac{n}{2}}^{i{k}^{\prime }}+w_n^{k^{\prime }}\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{w}_{\frac{n}{2}}^{i{k}^{\prime }}\text{(5)}& & =\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i}{w}_{\frac{n}{2}}^{ik}-w_n^k\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{w}_{\frac{n}{2}}^{ik}\end{array}$$

显然，由 $F(k)$ 和 $F(k^{\prime })$ 的相似性，我们可以直接递归求解，时间复杂度 $\mathcal{T}(n)=2\mathcal{T}(\frac{n}{2})+\mathcal{O}(n)$，即为 $\mathcal{O}(n\log n)$。

当然，除了递归，我们也可以使用迭代的方式求解。我们需要将 $a$ 序列变换位置，这被称为 蝴蝶变换（butterfly transform）。

观察发现，在上述式子中 $a_{2i}$ 和 $a_{2i+1}$ 分别被移到 $i$ 和 $i+\frac{n}{2}$ 位置递归。在二进制上，令 $2^l=n$，则 $a_i$ 最终变换到的位置是 $i$ 在 $l$ 位二进制上对称过来的二进制对应的值。例子如下：

13:01101 -> 10110:22

1|0DIF：按频域抽取（Decimation in Frequency）

$$\begin{array}{}\text{(6)}& F(k)& =\sum _{i=0}^{n-1}{a}_i{w}_n^{ik}\text{(7)}& & =\sum _{i=0}^{\frac{n}{2}-1}{a}_i{w}_n^{ik}+\sum _{i=0}^{\frac{n}{2}-1}{a}_{i+\frac{n}{2}}{w}_n^{(i+\frac{n}{2})k}\text{(8)}& & =\sum _{i=0}^{\frac{n}{2}-1}(a_i+(-1{)}^k{a}_{i+\frac{n}{2}})w_n^{ik}\end{array}$$

同理，当 $k$ 为偶数时，我们令 $k^{\prime }=k+1$，则有

$$\begin{array}{}\text{(9)}& F(k)& =\sum _{i=0}^{\frac{n}{2}-1}(a_i+a_{i+\frac{n}{2}})w_{\frac{n}{2}}^{\frac{k}{2}i}\text{(10)}& F(k^{\prime })& =\sum _{i=0}^{\frac{n}{2}-1}(a_i-a_{i+\frac{n}{2}})w_n^i{w}_{\frac{n}{2}}^{\frac{k}{2}i}\end{array}$$

一样的，可以递归，也可以迭代。直接对原序列迭代，最后得到的 $F$ 是蝴蝶变换后的 $F$，再变换一次即可。

$\mathtt{\text{DIT}}$ 和 $\mathtt{\text{DIF}}$ 的代码之后给出（只有迭代式）。

1|0反方向的钟 —— IDFT 的奇妙冒险

虽然存在一种用矩阵理解的方式，但对于写博客的人来说太折磨了，所以我们用一种简单粗暴的方式来理解并证明它。

事实上，如果我们令 ${\mathrm{DFT}}_k(a,n)$ 表示当序列为 $a$，长度为 $n$ 是 $F(k)$ 的值。同理定义 ${\mathrm{IDFT}}_k(a,n)$，则有

$${\mathrm{IDFT}}_k(a,n)={\mathrm{DFT}}_k^{-1}(a,n)={\displaystyle \frac{1}{n}}\sum _{i=0}^{n-1}{a}_i{w}_n^{-ik}$$

事实上这不好理解，也不易证明，但我还是竭尽码力来书写一番。

1|0证明

即证，$IDF{T}_k(\{\sum _{i=0}^{n-1}{a}_i{w}_n^i,\sum _{i=0}^{n-1}{a}_i{w}_n^{2i},\dots ,\sum _{i=0}^{n-1}{a}_i{w}_n^{(n-1)i}\},n)=a_k$。

直接展开，得

$${\displaystyle \frac{1}{n}}\sum _{i=0}^{n-1}\left(\sum _{j=0}^{n-1}{a}_j{w}_n^{ij}\right)w_n^{-ki}={\displaystyle \frac{1}{n}}\sum _{j=0}^{n-1}{a}_j\sum _{i=0}^{n-1}{w}_n^{i(j-k)}$$

显然，当 $j\ne k$ 时，$\sum _{i=0}^{n-1}{w}_n^{i(j-k)}=0$；而当 $j=k$ 时，$\sum _{i=0}^{n-1}{w}_n^{i(j-k)}=\sum _{i=0}^{n-1}1=n$。

综上，原式等于 $a_k$，即证。

1|0我的代码实现大抵没问题 —— Coding

事实上，当我们从真正意义上将 $\mathtt{\text{DFT}}$ 与 $\mathtt{\text{IDFT}}$ 联系起来，就可以打出一份完整的代码了：（下为用 $\mathtt{\text{DIT}}$ 书写的板子）

// 数组实现版本
using comp = complex<double>;
const double PI = acos(-1);

void FFT(comp *const s, int op) {
    for (int i = 0; i < l; ++i) if (i < r[i]) swap(s[i], s[r[i]]); // 做蝴蝶变换
    for (int k = 1; k < l; k <<= 1) {
        const comp W(cos(PI / k), op * sin(PI / k)); // 单位根
        for (int i = 0, K = k << 1; i < l; i += K) {
            comp w(1);
            for (int j = 0; j < k; ++j, w *= W) {
                comp x = s[i + j], y = w * s[i + j + k];
                s[i + j] = x + y;
                s[i + j + k] = x - y;
            }
        }
    }

    if (op == -1) for (int i = 0; i < l; ++i) s[i] /= l;
}

// vector 实现版本（个人推荐）
using comp = complex<double>;
const double PI = acos(-1);

void FFT(vector<comp> &s, int op) {
    for (int i = 0; i < l; ++i) if (i < r[i]) swap(s[i], s[r[i]]); // 做蝴蝶变换
    for (int k = 1; k < l; k <<= 1) {
        const comp W(cos(PI / k), op * sin(PI / k)); // 单位根
        for (int i = 0, K = k << 1; i < l; i += K) {
            comp w(1);
            for (int j = 0; j < k; ++j, w *= W) {
                comp x = s[i + j], y = w * s[i + j + k];
                s[i + j] = x + y;
                s[i + j + k] = x - y;
            }
        }
    }

    if (op == -1) for (int i = 0; i < l; ++i) s[i] /= l;
}

然后是蝴蝶变换的递推式（以方便 $\mathcal{O}(n)$ 求出变换值）：

for (int i = 0; i < l; ++i)
    r[i] = (r[i >> 1] >> 1) | (i & 1 ? l >> 1 : 0);
// 读者易推

1|3参考博文

再探 FFT – DIT 与 DIF，另种推导和优化

__EOF__

本文作者：TrueFalse
本文链接：https://www.cnblogs.com/cqbzljh/p/18002240.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！