Convolution

概念

卷积 (Convolution), 是透过两个函数 \(f\) 和 \(g\) 生成第三个函数的一种数学算子. ----Wikipedia

上面是卷积的数学定义, 讨论的是连续函数的卷积, 在计算机科学中我们常用的一般的卷积就是对多项式做乘法, 属于离散卷积.

假设我们有两个 \(n\) 项的多项式, \(f(x) = \sum_{i = 0}^{n - 1}a_ix^i\), \(g(x) = \sum_{i = 0}^{n - 1}b_ix^i\). 求两个多项式的差卷积 \(\sum_{i = 0}^{2n - 2}\sum_{j = 0}^{i} a_jb_{i - j}x_i\).

常用的离散卷积的计算方法有三种: 直接计算, 分段卷积和快速傅里叶变换.

前两种方法在两个多项式项数相差非常悬殊的时候效率较高, 但是因为最后一种方法在各种情况下都不低, 所以我们主要讨论最后一种.

Fourier transform

傅里叶变换 (Fourier transform) 属于傅里叶分析的领域, 可以将一个函数 \(f\) 生成另一个函数 \(\hat f\), 记作 \(\hat f = \mathcal{F}[f]\). 这个函数的函数值是复数, \(\hat f(x)\) 实部虚部分别表示复合信号 \(f\) 中, 频率为 \(x\) 的信号的振幅和相位.

傅里叶变换是用积分定义的:

它的逆变换是这样的:

对于卷积运算来说, 傅里叶变换的意义在于:

这样就允许我们把多项式的系数乘法转化为点值乘法, 因为求一个 \(n\) 项的多项式的系数, 最少需要 \(n\) 个点值. 先求出 \(f\), \(g\) 的 \(O(n)\) 个点值, 然后用 \(O(n)\) 的时间算出 \(f * g\) 的 \(O(n)\) 个点值, 进行傅里叶逆变换即可求出 \(f * g\) 的系数表达.

接下来需要解决的就是对 \(f\), \(g\) 进行傅里叶变换得到 \(\mathcal{F}[f]\), \(\mathcal{F}[g]\) 的点值, 和对求出的 \(\mathcal{F}[f * g]\) 的点值进行傅里叶逆变换 \(f * g\) 的过程.

Discrete-time Fourier Transform

离散时间傅里叶变换 (Discrete-time Fourier Transform, DTFT), 指的是在原函数 \(f\) 上对时间离散取样 \(x_i\), 根据这些 \(f(x_i)\), 求出连续函数 \(\mathcal{F}[f]\) 的操作.

Discrete Fourier Transform

离散傅里叶变换 (Discrete Fourier Transform, DFT), 则是在离散时间傅里叶变换的基础上, 只求出 \(\mathcal{F}[f]\) 的一些点值而不是一个连续的函数.

相比于傅里叶变换中无限的时域, DFT 的时域是 \([0, n)\), 表示离散周期信号的一个周期. \(f(j)\) 便是在时间取 \([0, n)\) 中任意整数值的时候, 原信号 \(f\) 的点值.

DFT 在傅里叶变换的基础上, 把积分变成了求和. 其中频域也是离散的, 为 \([0, n)\) 内所有的整数, \(\mathcal{F}[f] (k)\) 中 \(k\) 表示一个周期为单位时间内, 频率为 \(k\) 的分量的振幅相位.

式子中 \(e\) 和 \(i\) 分别是自然对数的底数和复数单位, 它们在逆变换的定义式中的意义与之相同.

DFT 的逆变换 IDFT 则被用来根据 \(\mathcal{F}[f]\) 的点值求出 \(f\) 的点值.

但是这样求一个点值就需要 \(O(n)\) 的时间, 总复杂度就需要 \(O(n^2)\). 如何快速求多项式的 DFT 和其逆变换便是我们求多项式乘法的关键.

Euler's formula

这里的欧拉公式是最优美的那个: \(e^{ix} = cos(x) + isin(x)\), 可以通过泰勒级数验证.

这个公式可以理解为用 \(e\) 的幂和辐角为 \(x\) 的单位复数的关系.

单位根

定义 \(n\) 次单位根的 \(n\) 次方等于 \(1\), 记作 \(w_n^k\), 共 \(n\) 个. \(w_n^k\) 是一个复数, 它的模长为 \(1\), 辐角为 \(\frac{2\pi k}{n}\), 也就是说, 它是一个单位复数. 用欧拉公式可以表示为:

代入 DFT 的定义式:

DFT 和点值

我们前面讨论的 DFT 求的是 \(\mathcal{F}[f]\) 的点值, 如果想求 \(f\) 的点值, 就需要选择特定的取样点.

如果说我们认为一个多项式 \(F(x)\) 是这样的:

这时我们取 \(x = w_n^{-t}\), 带入这个式子, 发现:

\(F(w_n^{-t})\) 的点值不就是 \(f(k)\) 的 DFT 吗? 这便是 DFT 和求点值的联系所在. 又因为 \(n\) 次单位根的周期性, 我们知道 \(w_n^{-t} = w_n^{n - t}\).

这样就可以直接把 DFT 得到的序列作为 \(f\) 的点值序列了.

IDFT 和插值

可以用 DFT 的逆变换 (IDFT) 来解决知道 \(F\) 在 \(w_n^t (t \in Z \cup [0, n))\) 的 \(n\) 个点值, 求 \(F\) 的系数表示的插值问题.

重新审视 IDFT 的定义式:

因为前面推出 \(F(w_n^{n - t}) = \mathcal{F}[f] (t) (t \in Z \cup [0, n))\), 代入定义式:

可以用 \(F(w_n^k)\) 直接 IDFT 得到 \(nf(j)\), 也就是 \(F\) 的系数序列的 \(n\) 倍.

Fast Fourier Transform

快速傅里叶变换 (Fast Fourier Transform, FFT), 是用来快速计算多项式的离散傅里叶变换和其逆变换的方法.

这里主要研究的是库利-图基 (Cooley-Tukey) 算法. 我们假设 \(n\) 是 \(2\) 的整数幂, 如果问题中不是, 那么后面的项可以认为是 \(0\), 将式子补齐. 这样做不会影响算法复杂度和正确性.

如果这时我们把 \(f\) 的偶数项变成 \(g(j) = f(2j) (j \in Z \cup [0, \frac n2))\), 奇数项变成 \(h(j) = f(2j + 1) (j \in Z \cup [0, \frac n2))\). 那么 \(\mathcal{F}[f]\) 就可以用 \(\mathcal{F}[g]\) 和 \(\mathcal{F}[h]\) 来表示.

当然, \(g\) 和 \(h\) 的频域可能不包含 \(k\), 但是因为离散周期信号是按周期循环的, 所以我们这里的 \(k\) 也可以是 \(k - \frac n2\).

对于只有一个项的序列 \(f\), 它的 DFT 可以 \(O(1)\) 求出:

那么我们需要做的就是对于 \(n > 1\) 的情况, 递归求解奇数偶数项的 DFT, 然后将它们合并. 第 \(x\) 层递归, 有 \(2^{\log n - x}\) 个不同的 \(k\) 值, 有 \(2 ^ x\) 组不同的系数序列. 每次除递归之外的时间复杂度是 \(O(1)\), 因此每层的时间复杂度为 \(O(n)\). 从一开始的第 \(0\) 层开始, 一共是 \(\log n + 1\) 层. 因此总复杂度是 \(O(n \log n)\).

对于 IDFT, 其原理也是一样的:

因此其实现和 DFT 是相同的, 只是把 \(-w_n\) 变成了 \(w_n\), 可以写成一个函数. 复杂度仍然是 \(O(n \log n)\). 边界条件:

因此我们把 \(f\) 序列进行 DFT 可以得到 \(\mathcal{F}[f]\) 序列, 把 \(\mathcal{F}[f]\) 序列进行 DFT 可以得到 \(nf\) 序列.

Decimation-in-time

递归毕竟是大常数的写法, 所以实践中尝试用更加方便高效的非递归写法.

定义运算 \(Rev_{x}(a)\) 表示把小于 \(2^x\) 的数字 \(a\) 当成长度为 \(x\) 的 01 串, 把这个串镜像过来之后得到的数值大小.

我们知道在 DFT 过程中, 递归时第 \(x\) 层有 \(2^x\) 个子问题. 回溯时我们需要把第 \(x\) 层的第 \(j\) 个子问题和第 \(j \oplus 2^{x - 1}\) 个子问题合并成第 \(x - 1\) 层的第 \(j \And (2^{x - 1} - 1)\) 个子问题. 其中, 两个子问题的第 \(k\) 项进行操作可以生成新问题的第 \(k\) 项和第 \(k + 2^{\log n - x}\) 项.

如果想要避免递归, 一个很重要的目标是实现原址操作. 假设一个下标 \(j\), 前 \(x\) 位从右往左读是 \(j_{x, 1}\), 后 \(\log n - x\) 位从左往右读是 \(j_{x, 2}\) (先读的是高位, 后读的是低位). 如果在第 \(\log n\) 层, 使得第 \(j\) 位存储第 \(j_{x, 1}\) 个子问题的唯一的一项. 如果每一层都能保证第 \(j\) 位存储的是第 \(j_{x, 1}\) 个子问题的第 \(j_{x, 2}\) 项, 并且保证回溯对每两个数进行计算时不会影响其它位置, 就可以通过下标快速计算一些所需的变量, 从而方便地原址求 DFT 了.

利用归纳法, 假设我们在第 \(x\) 层满足第 \(j\) 位存储的是第 \(j_{x, 1}\) 个子问题的第 \(j_{x, 2}\) 项, 回溯到第 \(x - 1\) 层需要第 \(k_1\) 和第 \(k_1 \oplus 2^{x - 1}\) 个子问题各自的第 \(k_2\) 项相互计算出 \(k_1 \And (2^{x - 1} - 1)\) 的第 \(k_2\) 项和第 \(k_2 + 2^{\log n - x}\) 项.

根据一开始的规定, 第 \(x\) 层的第 \(k_1\) 个子问题的第 \(k_2\) 项的下标是 \(2^{\log n - x}Rev_{x}(k_1) + k_2\), 第 \(x\) 层的第 \(k_1 \oplus 2^{x - 1}\) 个子问题的第 \(k_2\) 项的下标是 \(2^{\log n - x}Rev_{x}(k_1 \oplus 2^{x - 1}) + k_2\). 变形整理第二个下标:

因为 \(k_2\) 是第 \(x\) 层的子问题中的项, 所以一定小于 \(2^{\log n - x}\). 因此:

继续变形下标:

因此这两个下标是计划中第 \(x - 1\) 层的第 \(k_1 \And (2^{x - 1} - 1)\) 个子问题的第 \(k_2\) 和 \(k_2 + 2^{\log n - x}\) 项, 因此可以保持原址操作.

对于长度为 \(8\) 的序列, 其过程如图所示, 图中 \(f(a, b, c)\) 表示第 \(a\) 层回溯时, 第 \(b\) 个子问题的第 \(c\) 项:

因为这个过程的输入是信号, 是在时域的每个点取样, 所以又叫时域抽取法 (Decimation-in-time, DIT). 使用 DIT 时需要先把变换序列的第 \(j\) 项赋值为原信号的第 \(Rev_{\log n}(j)\) 项, 最后变换得到的结果序列中 \(j\) 位置存储的则是频谱中的第 \(j\) 项.

下面是 DIT 实现的库利-图基算法的代码, 其中 Cplx(x) 表示 \(w_n\) 或 \(-w_n\) 的 \(x\) 次方, 如果是 DFT 则是 \(-w_n\), 否则是 \(w_n\).

inline void DIT(Cplx* f) {
  for (unsigned i(Lgn), j(1); i; --i, j <<= 1) {
    for (unsigned k(0); k < n; ++k) if (!(k & j)) {
      Cplx Tma(f[k]), Tmb(f[k + j]), W(Cplx((k & ((j << 1) - 1)) << (i - 1)) * Tmb);
      f[k] = Tma + W;
      f[k + j] = Tma - W;
    }
  }
}

Decimation-in-frequency

DIT 需要把信号以二进制镜像的下标存储, 我们如果仅使用 DIT, 那么一次多项式乘法就需要进行 \(3\) 次转置 (输入的两个序列的转置和点值乘法后的转置), 这无疑是没有必要的. 如果考虑逆过程, 直接把频谱扔进一个函数, 得到的是转置后的信号, 和 DIT 同时使用就可以完全避免转置.

与 DIT 相对的, 频域抽取法 (Decimation-in-frequency, DIF), 是前者的逆过程, 它模拟的是 DIT 的逆过程.

已知

如果已知 \(\mathcal{F}[f]\), 求 \(\mathcal{F}[g]\) 和 \(\mathcal{F}[h]\).

因此直接把频谱喂给 DIF 实现的库利-图基算法, 就可以得到 \(n\) 倍的原信号转置后的序列. 因为输入是频谱, 定义在频域上, 所以称为频域抽取法下面是代码.

inline void DIF(Cplx* f) {
  for (unsigned i(1), j(n >> 1); i <= Lgn; ++i, j >>= 1) {
    for (unsigned k(0); k < n; ++k) if (!(k & j)) {
      Cplx Tma(f[k]), Tmb(f[k + j]);
      f[k] = Tma + Tmb;
      f[k + j] = (Tma - Tmb) * Cplx((k & (j - 1)) << (i - 1));
    }
  }
}

代码的剩余部分

由于是复数操作, 所以我们首先定义一个复数类.

double Cp[2100000][2];
unsigned m, n, nn, Lgn(0);
unsigned A, B, C, D, t;
unsigned Cnt(0), Ans(0), Tmp(1);
char Inv(0);
struct Cplx {
  double Rel, Img;
  inline Cplx() {}
  inline Cplx(unsigned x) {
    Rel = Cp[x][0]; Img = Inv ? Cp[x][1] : -Cp[x][1];
  }
  inline Cplx operator + (const Cplx& x) {
    Cplx Rt(x);
    Rt.Rel += Rel, Rt.Img += Img;
    return Rt;
  }
  inline Cplx operator + (const double& x) {
    Cplx Rt(*this);
    Rt.Rel += x, Rt.Img;
    return Rt;
  }
  inline Cplx operator - (const Cplx& x) {
    Cplx Rt(x);
    Rt.Rel = Rel - Rt.Rel, Rt.Img = Img - Rt.Img;
    return Rt;
  }
  inline Cplx operator - (const double& x) {
    Cplx Rt(*this);
    Rt.Rel -= x, Rt.Img;
    return Rt;
  }
  inline Cplx operator * (const Cplx& x) {
    Cplx Rt;
    Rt.Rel = Rel * x.Rel - Img * x.Img, Rt.Img = Img * x.Rel + Rel * x.Img;
    return Rt;
  }
  inline Cplx operator * (const double& x) {
    Cplx Rt(*this);
    Rt.Rel *= x, Rt.Img *= x;
    return Rt;
  }
  inline Cplx operator / (const Cplx& x) {
    Cplx Rt;
    double TmpDe(x.Rel * x.Rel + x.Img * x.Img);
    Rt.Rel = (Rel * x.Rel + Img * x.Img) / TmpDe;
    Rt.Img = (Img * x.Rel - Rel * x.Img) / TmpDe;
    return Rt;
  }
  inline Cplx operator / (const double& x) {
    Cplx Rt(*this);
    Rt.Rel /= x, Rt.Img /= x;
    return Rt;
  }
}a[2100000], b[2100000];

然后是主函数. 我们用 DIF 求出两个输入信号的频谱的转置, 然后用 DIT 求出转置的频谱相乘得到的结果的原信号增强 \(n\) 倍的结果.

signed main() {
  nn = RD() + 1, m = RD() + 1, n = nn + m - 1;
  while (Tmp < n) Tmp <<= 1, ++Lgn; n = Tmp;
  for (unsigned i(0); i <= n; ++i) {//预处理 n 次单位根
    double Arc(Pi * 2 * i / n);
    Cp[i][0] = cos(Arc), Cp[i][1] = sin(Arc);
  }
  for (unsigned i(0); i < nn; ++i) a[i].Rel = RD();
  for (unsigned i(0); i < m; ++i) b[i].Rel = RD();
  DIF(a), DIF(b);//正变换
  for (unsigned i(0); i < n; ++i) a[i] = a[i] * b[i];
  Inv = 1, DIT(a);//逆变换
  for (unsigned i(0); i < n; ++i) a[i] = a[i] / n;//逆变换让值增加到原来的 n 倍
  for (unsigned i(0); i < m + nn - 1; ++i) printf("%u ", (unsigned)(a[i].Rel + 0.5));putchar(0x0A);
  return Wild_Donkey;
}

Prime Root

假设 \(a\), \(m\) 互质, 我们知道 \(a^d \equiv a^{d \% \phi(m)} \pmod m\), 因此任何整数 \(a\) 在模 \(m\) 意义下的整数幂的结果, 只有 \(\phi(m)\) 种.

如果对于正整数 \(a\), 使得 \(a^d \equiv 1 \pmod m\) 成立的最小的正整数 \(d\) 等于 \(\phi(m)\). 则称 \(a\) 是模 \(m\) 的一个原根 (Prime Root).

如果 \(m\) 可以用正整数 \(n\) 和奇质数 \(p\) 表示为 \(p^n\) 或 \(2p^n\), 又或是 \(m = 2\) 或 \(m = 4\), 则存在模 \(m\) 的原根.

通过有关群论的知识我们知道, 如果一个数 \(m\) 有原根, 那么它的不同的原根数量为 \(\phi(\phi(m))\).

Number-Theoretic Transform

数论变换 (Number-Theoretic Transform, NTT), 是用原根代替 \(n\) 次单位根以避免复数运算的处理整数离散周期信号的算法.

仍然是把原来的序列加 \(0\) 补齐为长度为 \(2\) 的整数幂的序列, 长度为 \(n\). 我们选择一个质数 \(m\) 作为模数, 需要满足 \(n | m - 1\), 找出模 \(m\) 的一个原根 \(\alpha\), 把 \(\alpha^{\frac {m - 1}n} \pmod m\) 记作 \(w_n\).

NTT 的定义式和 DFT 的定义式形式上十分相似, 如果没有说明, 我们下面提到的所有运算都是模 \(m\) 意义下的, 用 \(Inv(x)\) 表示 \(x\) 模 \(m\) 意义下的乘法逆元.

先来看正变换:

然后是逆变换:

发现式子可以直接使用库利-图基算法求, 其流程和复数实现的 DFT 是一样的.

const unsigned long long Mod(998244353);
unsigned long long W, a[2100000], b[2100000];
unsigned m, n, nn, Lgn(0);
unsigned A, B, C, D, t;
unsigned Cnt(0), Ans(0), Tmp(1);
char Inv(0);
inline unsigned long long Pow(unsigned long long x, unsigned y) {
  unsigned long long PR(1);
  while (y) {
    if (y & 1) PR = PR * x % Mod;
    x = x * x % Mod, y >>= 1;
  }
  return PR;
}
inline void DIT(unsigned long long* f) {
  for (unsigned i(1), j(Lgn - 1); ~j; --j, i <<= 1) {
    unsigned long long Pri(Pow(Inv ? W : Pow(W, Mod - 2), 1 << j)), Po(1);
    for (unsigned k(0); k < n; ++k, Po = Po * Pri % Mod) if (!(k & i)) {
      unsigned long long Tma(f[k]), Tmb(f[k + i] * Po % Mod);
      f[k] = Tma + Tmb;
      f[k + i] = Mod + Tma - Tmb;
      if (f[k] >= Mod) f[k] -= Mod;
      if (f[k + i] >= Mod) f[k + i] -= Mod;
    }
  }
}
inline void DIF(unsigned long long* f) {
  for (unsigned i(n >> 1), j(0); i; ++j, i >>= 1) {
    unsigned long long Pri(Pow(Inv ? W : Pow(W, Mod - 2), 1 << j)), Po(1);
    for (unsigned k(0); k < n; ++k, Po = Po * Pri % Mod) if (!(k & i)) {
      unsigned long long Tma(f[k]), Tmb(f[k + i]);
      f[k] = Tma + Tmb;
      if (f[k] >= Mod) f[k] -= Mod;
      f[k + i] = (Mod + Tma - Tmb) * Po % Mod;
    }
  }
}
signed main() {
  nn = RD() + 1, m = RD() + 1, n = nn + m - 1;
  while (Tmp < n) Tmp <<= 1, ++Lgn; n = Tmp;
  W = Pow(3, (Mod - 1) / n);
  for (unsigned i(0); i < nn; ++i) a[i] = RD();
  for (unsigned i(0); i < m; ++i) b[i] = RD();
  DIF(a), DIF(b);
  for (unsigned i(0); i < n; ++i) a[i] = a[i] * b[i] % Mod;
  Inv = 1, DIT(a), W = Pow(n, Mod - 2);
  for (unsigned i(0); i < n; ++i) a[i] = a[i] * W % Mod;
  for (unsigned i(0); i < m + nn - 1; ++i) printf("%llu ", a[i]); putchar(0x0A);
  system("pause");
  return Wild_Donkey;
}

总结

综合前面两种计算卷积的方法的共同点, 只要一个数 \(\alpha\) 符合下面的条件, 都可以用来作为 DFT 计算卷积的旋转因子, 而 \(-\alpha\) 则被作为是 IDFT 过程的旋转因子:

• \(\alpha^{2^k} = 1\)

• \(\alpha^{j + 2^{k - 1}} = -\alpha^{j}\)

值得注意的是, 只有 \(-w_n\) 可以用来将信号转化为频谱, 其它的 \(\alpha\) 都只能用于计算卷积.