[多项式学习]快速数论变换

Posted on | In

快速数论变换将 FFT 的作用域移到有限环 $\Z_m$ 上来，更加符合实际需求（如无精度误差、可实时取模）。

算法思路

如同在 $\C$ 内寻找单位根 $\omega_n$ 一样，我们需要在 $\Z_m$ 上寻找与单位根性质类似的元素，在这里则是原根。

阶

设 $a\perp m$，则由 Euler 数论定理知，$a^{\varphi(m)}\equiv 1\pmod m$。

因此使 $a^n\equiv 1\pmod m$ 成立的最小正整数 $n$ 存在，称其为 $a$ 在模 $m$ 下的阶，记作 $\delta_m(a)$，有些地方也记作 $\operatorname{ord}_m(a)$。

在抽象代数中，阶的定义可以推广到任意群，下面讨论模 $m$ 简化剩余系上阶的性质大多数也可以推广到任意群。

性质 1

$a^1,a^2,\cdots,a^{\delta_m(a)}$ 模 $m$ 两两不同余。

证明 考虑反证法，设 $a^i\equiv a^j\pmod p$ 且 $i<j\le \delta_m(a)$，则 $a^{j-i}\equiv 1\pmod m$，而 $j-i<\delta_m(a)$，与 $\delta_m(a)$ 的最小性矛盾。$\square$

性质 2

若 $a^n\equiv 1\pmod m$，则 $\delta_m(a)\mid n$。

证明 对 $n$ 除以 $\delta_m(a)$ 做带余除法，$n=q\delta_m(a)+r$，$0\le r<\delta_m(a)$。若 $r\ne 0$，则 $a^n=a^{q\delta_m(a)}\cdot a^r=\left(a^{\delta_m(a)}\right)^q\cdot a^r\equiv a^r\equiv 1\pmod m$，而 $r<\delta_m(a)$，与 $\delta_m(a)$ 的最小性矛盾。$\square$

推论 1

若 $a^p\equiv a^q\pmod m$，则 $p\equiv q\pmod {\delta_m(a)}$。

性质 3

设 $a,b\perp m$，则 $\delta_m(ab)=\delta_m(a)\delta_m(b)$ 的充要条件是 $\delta_m(a)\perp \delta_m(b)$。

证明

必要性：由 $a^{\delta_m(a)}\equiv 1\pmod m$ 和 $b^{\delta_m(b)}\equiv 1\pmod m$ 可知：

由推论 2 知：

所以，

于是 $\delta_m(a)\perp \delta_m(b)$。

充分性：由 $(ab)^{\delta_m(ab)}\equiv 1\pmod m$ 可知：

由推论 2 知，$\delta_m(a)\mid \delta_m(ab)\delta_m(b)$。结合 $\delta_m(a)\perp \delta_m(b)$ 得

同理有

再根据 $\delta_m(a)\perp \delta_m(b)$ 有

另一方面，

因此

综合两点可得

$\square$

性质 4

设 $k\in \N$，$a\perp m$，则

证明 注意到

运用性质 2 有

进而

又由

知

综合得到

$\square$

原根

若 $a\perp m$，且 $\delta_m(a)=\varphi(m)$，则称 $a$ 为模 $m$ 的一个原根。

在抽象代数中，原根就是循环群的生成元。若模 $m$ 的简化剩余系与一个循环群同构，则说明 $m$ 有原根。

原根判定定理

设 $m\ge 3$，$a\perp m$，则 $a$ 是模 $m$ 的原根的充要条件是，对于 $\varphi(m)$ 的每一个素因子 $p$，$a^{\frac{\varphi(m)}p}\not\equiv 1\pmod m$。

证明 必要性显然，接下来证明充分性。

若 $a$ 不是模 $m$ 的原根，则 $\delta_m(a)<\varphi(m)$。通过 $a^{\varphi(m)}\equiv 1\pmod m$ 我们有 $\delta_m(a)\mid\varphi(m)$，因此 $\dfrac{\varphi(m)}{\delta_m(a)}>1$。容易知道 $\exist p:p\mid \dfrac{\varphi(m)}{\delta_m(a)}$，进而 $\delta_m(a)\mid\dfrac{\varphi(m)}p$，于是 $a^{\frac{\varphi(m)}p}\equiv 1\pmod m$，与条件矛盾。$\square$

原根个数

若 $m$ 有原根，则它原根个数为 $\varphi(\varphi(m))$。

证明 设 $g$ 为 $m$ 的一个原根，则运用性质 4：

观察到 $g^k$ 为原根当且仅当 $\gcd(\delta_m(g),k)=1$，而 $k$ 在 $1\sim \varphi(m)$ 中取值，根据定义知满足条件的 $k$ 共有 $\varphi(\varphi(m))$ 个。

注意到 $g^i,i\le \varphi(m)$ 已经取遍模 $m$ 的简化剩余系，因此不存在剩余的原根。$\square$

原根存在定理

$m$ 有原根当且仅当 $m=2,4,p^\alpha,2p^{\alpha}$，其中 $p$ 为奇素数，$\alpha\in\N^\ast$。

证明 分成 $m=2,4$，$m=p^\alpha$，$m=2p^\alpha$，$m\ne 2,4,p^\alpha,2p^\alpha$ 四个部分证明。

  • $m=2,4$：显然分别有原根 $2,3$。

  • $m=p^\alpha$：

    • $\alpha=1$：

      引理 1 设 $a\perp p$，$b\perp p$，则 $\exist c\in \Z:\delta_p(c)=\operatorname{lcm}(\delta_p(a),\delta_p(b))$

      证明 将 $\delta_m(a)$，$\delta_m(b)$ 表示成素因数分解的形式：

      然后按指数大小关系分解：

      其中

      则由阶的性质 4 可得：

      同理

      显然有 $Y\perp W$，$YW=\operatorname{lcm}(\delta_p(a),\delta_p(b))$，因此根据阶的性质 3 可得：

      令 $c=a^Xb^Z$ 则引理 1 得证。$\square$

      对 $1\sim (p-1)$ 依次两两使用引理 1，可知存在 $g\in \Z$ 使得

      这表明 $\delta_p(j)\mid \delta_p(g)(j=1,\cdots,p-1)$，因此 $j=1,\cdots,p-1$ 都是方程

      的解。

      Lagrange 定理 设 $p$ 为素数，则模 $p$ 意义下的整系数多项式

      的同余方程 $f(x)\equiv 0\pmod p$ 至多有 $n$ 个不同解。

      证明 对 $n$ 使用归纳法。当 $n=0$ 时，方程显然无解。

      假设命题对 $\deg f<n$ 的 $f$ 都成立，则假设存在一个 $\deg f=n$ 且拥有至少 $n+1$ 个不同的解 $x_0,\cdots,x_n$ 的 $f$。

      因为 $f(x)-f(x_0)=(x-x_0)g(x)$，其中 $g(x)$ 是一个最多 $n-1$ 次的多项式。

      根据假设容易知道对于 $i=1,\cdots,n$，有

      因而 $g$ 有至少 $n$ 个解 $x_1,\cdots,x_n$，矛盾。$\square$

      因为方程有至少 $p-1$ 个解，则说明方程的次数 $\delta_p(g)\ge p-1$，由 $\delta_p(g)\mid p-1$ 容易知道 $\delta_p(g)=p-1=\varphi(p)$，那么 $g$ 是 $p$ 的原根。$\square$

    • $\alpha>1$：考虑平移 $p$ 的原根 $g$。

      引理 2 存在模 $p$ 的原根 $g$ 使得 $g^{p-1}\not\equiv 1\pmod {p^2}$。

      证明 任取 $p$ 的原根 $g$，若 $g$ 不满足条件，则认定 $g+p$ 满足条件。

      我们有

      其中最后一步是因为 $p\nmid g^{p-2}$，因而 $p^2\nmid pg^{p-2}$。$\square$

      考虑证明：若 $g$ 是一个满足引理 2 所述条件的原根，则对于任意 $\alpha\in\N^\ast$，$g$ 是模 $p^\alpha$ 的原根。

      引理 3 $\forall \beta\in\N^\ast:\exist k\beta:g^{\varphi(p^\beta)}=1+p^\beta\cdot k\beta\land p\nmid k_\beta$。

      证明 当 $\beta=1$ 时根据 $g$ 的选取显然成立。假设上式对 $\beta$ 成立，则

      结合 $p\nmid k_\beta$ 知命题对 $\beta+1$ 成立。

      所以命题对任意 $\beta\in\N^\ast$ 都成立。$\square$

      记 $\delta=\delta_{p^\alpha}(g)$，则易知 $\delta\mid \varphi(p^\alpha)=p^{\alpha-1}(p-1)$。

      由于 $g$ 为模 $p$ 的原根，且 $g^\delta\equiv 1\pmod{p^\alpha}$，于是 $g^\delta\equiv 1\pmod p$，然后有 $\delta_p(g)=\varphi(p)=p-1\mid \delta$，所以可设 $\delta=p^{\beta-1}(p-1)$，其中 $1\le\beta\le\alpha$。根据引理 3，

      观察可得

      由 $g^\delta\equiv 1\pmod{p^\alpha}$ 知 $\beta+1>\alpha$，即 $\beta\ge \alpha$，结合 $\alpha\ge\beta$ 知 $\alpha=\beta$，即：

      从而，$g$ 是模 $p^\alpha$ 的原根。$\square$

  • $m=2p^\alpha$：

    设 $g_0$ 是模 $p^\alpha$ 的原根，则 $g_0+p^\alpha$ 也是模 $p^\alpha$ 的原根。在 $g_0$ 与 $g_0+p^\alpha$ 中有且仅有一个是奇数，设其为 $g$，则因为 $g\perp p^\alpha$，$g\perp 2$，因此 $g\perp 2p^\alpha$。

    容易知道，$\delta{2p^\alpha}(g)\mid \varphi(2p^\alpha)$。又因为 $g^{\delta{2p^\alpha}(g)}\equiv 1\pmod {2p^\alpha}$，进而

    于是 $\varphi(p^\alpha)=\delta{p^\alpha}(g)\mid \varphi(2p^\alpha)$。综合两者与 $\varphi(p^\alpha)=\varphi(2p^\alpha)$ 得到 $\delta{2p^\alpha}(g)=\varphi(2p^\alpha)$，说明 $g$ 是模 $2p^\alpha$ 的原根。$\square$

  • $m\ne 2,4,p^\alpha,2p^\alpha$：

    • $m=2^\alpha$，其中 $\alpha\in\N^\ast$ 且 $\alpha\ge 3$，则对于任意的奇数 $a=2k+1$，都有：

      进而

      因此 $m$ 无原根。$\square$

    • $m\ne 2^\alpha$，则设 $m=rt$，这里 $2<r<t$ 且 $r\perp t$。

      若 $a\perp m$，则 $a\perp r$，$a\perp t$，进而

      所以有

      注意到当 $n>2$ 时 $2\mid\varphi(n)$，所以

      因而 $m$ 不存在原根。$\square$

最小原根问题

素数 $p$ 的最小原根实际上非常的小，王元和 Burgess 分别在 1959 年和 1962 年证明了 $g_p=O(p^{1/4+\epsilon})$，而 Shoup 于 1990 年证明了（在假设广义 Riemann 猜想成立的前提下），$g_p=O(\log^6p)$。

当然我们知道，使得 $\log^6p<p^{1/4}$ 的数需要非常大，至少得要 $e^{-24W(1/24)}\approx2.125\times10^{49}$ 级别，这还没算上那个 $\epsilon$（两个人的原论文都找不到了，找不到具体的叙述），因而我们可以在平时采用 $O(p^{1/4+\epsilon})$ 这个界，要计算时暴力枚举寻找原根。

快速数论变换

考虑寻找一个拥有原根 $g$ 的素数 $p$，其中 $p-1=q\times 2^k$，$k$ 足够大使得 $g^q$ 能够起到类似 FFT 中 $\omega_n$ 的效果。

常用的三个模数，最小原根皆为 $3$：

以 $p$ 为模数，以 $g^q$ 为 $\omega_n$ 的变换称为 快速数论变换 (Number-Therotic Transform, NTT)。

代码实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
/**
 * @file #108-3. 多项式乘法.cpp
 * @author Kuriyama Mirai (hermione_granger@foxmail.com)
 * @brief NTT版本
 * 不要再错同一个地方啊！
 * @date 2022-01-31
 * @link https://loj.ac/s/1368985
 * @link https://www.luogu.com.cn/record/68403193
 * 
 * @copyright Copyright (c) 2022
 * 
 */
#include <cstdio>
#include <algorithm>

namespace mirai {

constexpr int MAXN = 2200005;
constexpr int MOD = 998244353;
typedef long long ll;
ll a[MAXN], b[MAXN];
ll rev[MAXN];

ll pow(ll a, ll b) {
  ll res = 1;
  while (b) {
    if (b & 1) {
      res = res * a % MOD;
    }
    a = a * a % MOD;
    b >>= 1;
  }
  return res;
}
inline ll inv(ll const& a) {
  return pow(a, MOD - 2);
}

void ntt(ll arr[], int len, int on) {
  for (int i = 0; i < len; ++i) {
    if (i < rev[i]) {
      std::swap(arr[i], arr[rev[i]]);
    }
  }
  for (int n = 2; n <= len; n *= 2) {
    ll gq_n = pow(on == 1 ? 3 : inv(3), (MOD - 1) / n);
    for (int i = 0; i < len; i += n) {
      ll gq = 1;
      for (int j = i; j < i + n / 2; ++j) { // 注意这里不是从 0 到 n/2-1!
        ll u = arr[j], v = arr[j + n / 2];
        arr[j] = (u + gq * v) % MOD;
        arr[j + n / 2] = (u - gq * v % MOD + MOD) % MOD;
        gq = gq * gq_n % MOD;
      }
    }
  }
  if (on == -1) {
    ll invlen = inv(len);
    for (int i = 0; i < len; ++i) {
      arr[i] = arr[i] * invlen % MOD;
    }
  }
}

int main(int argc, char** argv) {
  int n, m;
  std::scanf("%d%d", &n, &m);
  for (int i = 0; i <= n; ++i) {
    std::scanf("%lld", &a[i]);
  }
  for (int i = 0; i <= m; ++i) {
    std::scanf("%lld", &b[i]);
  }
  int len = 1;
  while (len < n + m + 2) {
    len <<= 1;
  }
  for (int i = 1; i < len; ++i) {
    rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (len >> 1));
  }
  ntt(a, len, 1);
  ntt(b, len, 1);
  for (int i = 0; i < len; ++i) {
    a[i] = a[i] * b[i] % MOD;
  }
  ntt(a, len, -1);
  for (int i = 0; i <= n + m; ++i) {
    std::printf("%lld ", a[i]);
  }
  std::printf("\n");
  return 0;
}

}

int main(int argc, char** argv) {
  return mirai::main(argc, argv);
}