Bỏ qua

Giới thiệu số học module

算法竞赛中，数论部分的一个重要组成部分就是 模算术（modular arithmetic），也就是在某一模数下进行各种整数运算．除了基础的四则运算和求幂外，还可以方便地进行取对数、开各次方、求阶乘和组合数等运算．

模算术常见于各类问题中，而不仅仅局限于数论部分．很多问题的实际答案可能非常大，超过了常见的整型变量的存储范围．此时，为了避免引入大整数运算和输出长数字，题目常常要求对答案取模后输出．这就要求熟练掌握各类模算术技巧．

C/C++ 的整数除法和取模运算

在 C/C++ 中，整数除法和取模运算，与数学上习惯的取模和除法不一致．

对于所有标准版本的 C/C++，规定在整数除法中：

  1. 当除数为 0 时，行为未定义；
  2. 否则 (a / b) * b + a % b 的运算结果与 a 相等．

也就是说，取模运算的符号取决于除法如何取整；而除法如何取整，这是实现定义的（由编译器决定）．

从 C99 和 C++11 标准版本起，规定 商向零取整（舍弃小数部分）；取模的符号就与被除数相同．从此，以下断言保证为真：

1
2
3
4
assert(5 % 3 == 2);
assert(5 % -3 == 2);
assert(-5 % 3 == -2);
assert(-5 % -3 == -2);

模整数类

模算术可以看做是对某模数下的 同余类 进行各种运算．如果用一个结构体来表示一个同余类，并且将同余类之间的加、减、乘等运算封装为结构体的方法或运算符重载，那么模算术就可以自然地实现为一个模整数类．下面给出一个简单的例子，它支持模数 \(M < 2^{30}\) 下 \(32\) 位带符号整数的加法、减法、乘法以及快速幂运算：

一个简单的模整数类 
 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
// A simple ModInt implementation.
template <int M>
struct ModInt {
  struct skip_mod {};

  ModInt(int v, skip_mod) : v(v) {}

  int v;

  ModInt() : v(0) {}

  // Initialization: find remainder.
  // Equivalent to: v = int((x % M + M) % M)
  ModInt(long long x) {
    x %= M;
    if (x < 0) x += M;
    v = int(x);
  }

  // Addition.
  // Equivalent to: ModInt((l.v + r.v) % M)
  friend ModInt operator+(ModInt l, ModInt r) {
    int res = l.v + r.v;
    if (res >= M) res -= M;
    return ModInt(res, skip_mod{});
  }

  // Subtraction.
  // Equivalent to: ModInt((l.v - r.v + M) % M)
  friend ModInt operator-(ModInt l, ModInt r) {
    int res = l.v - r.v;
    if (res < 0) res += M;
    return ModInt(res, skip_mod{});
  }

  // Multiplication.
  friend ModInt operator*(ModInt l, ModInt r) {
    return ModInt(1LL * l.v * r.v % M, skip_mod{});
  }

  // Exponentiation.
  ModInt pow(long long b) const {
    ModInt res{1}, po{*this};
    for (; b; b >>= 1) {
      if (b & 1) res = res * po;
      po = po * po;
    }
    return res;
  }
};

这个实现有意地减少了取模运算的次数，因为取模操作通常比普通的加、减、乘或比较运算要耗时得多．代码注释中提供了等价且更为直接的实现方式．这些简单优化的主要思路是，两个 \([0,M)\) 内的整数做加减运算时，结果一定落在区间 \((-M,2M)\) 内，因此可以通过一次加减法调整回区间 \([0,M)\)．该实现中的取幂运算用到了 快速幂 的技巧．

除了这些基础运算外，还可以在各种模数下做如下运算：

这些运算通常在素数模数下比较容易．对于合数模数，往往需要用到对应算法的扩展版本和 中国剩余定理．这些模意义下的运算大多可以看作求解某种同余方程．对于求解同余方程的一般方法，可以参考 同余方程 页面．

相关算法

本节将介绍几种在模意义下优化取模、乘法和快速幂运算的方法．对于绝大多数题目来说，前文提供的简单实现已经足够高效．然而，当题目对算法常数有严格要求时，这些优化方法就可以发挥作用，通过减少不必要的计算和取模操作，进一步降低算法的时间开销．

快速乘

在素性测试与质因数分解中，经常会遇到模数在 long long 范围内的乘法取模运算．为了避免运算中的整型溢出问题，本节介绍一种可以处理模数在 long long 范围内，不需要使用 __int128 且复杂度为 \(O(1)\) 的「快速乘」．本算法要求测评系统中，long double 至少表示为 \(80\) 位扩展精度浮点数1．

假设 \(0 \le a, b < m\)，要计算 \(ab\bmod m\)．注意到：

\[ ab\bmod m=ab-\left\lfloor \dfrac{ab}m \right\rfloor m. \]

利用 unsigned long long 的自然溢出：

\[ ab\bmod m=ab-\left\lfloor \dfrac{ab}m \right\rfloor m=\left(ab-\left\lfloor \dfrac{ab}m \right\rfloor m\right)\bmod 2^{64}. \]

只要能算出商 \(\left\lfloor\dfrac{ab}m\right\rfloor\)，最右侧表达式中的乘法和减法运算都可以使用 unsigned long long 直接计算．

接下来，只需要考虑如何计算 \(\left\lfloor\dfrac {ab}m\right\rfloor\)．解决方案是先使用 long double 算出 \(\dfrac am\) 再乘上 \(b\)．既然使用了 long double，就无疑会有精度误差．假设 long double 表示为 \(80\) 位扩展精度浮点数（即符号为 \(1\) 位，指数为 \(15\) 位，尾数为 \(64\) 位），那么 long double 最多能精确表示的有效位数为 \(64\)2．所以 \(\dfrac am\) 最差从第 \(65\) 位开始出错，误差范围3为 \(\left(-2^{-64},2^{-64}\right)\)．乘上 \(b\) 这个 \(64\) 位带符号整数，误差范围为 \((-0.5,0.5)\)．为了简化后续讨论，可以先加一个 \(0.5\) 再取整，最后的误差范围是 \(\{0,1\}\)．

最后，代入上式计算时，需要乘以 \(-m\)，所以最后的误差范围是 \(\{0,-m\}\)．因为 \(m\) 在 long long 范围内，所以当结果 \(r\in[0,m)\) 时，直接返回 \(r\)，否则返回 \(r+m\)．

代码实现如下：

参考实现 
1
2
3
4
5
long long mul(long long a, long long b, long long m) {
  long long c = (unsigned long long)a * b -
                (unsigned long long)((long double)a / m * b + 0.5L) * m;
  return c < 0 ? c + m : c;
}

如今，绝大多数测评系统所配备的 C/C++ 编译器已支持 __int128 类型4，因此也可以直接将乘数类型提升至 __int128 后取模计算：

参考实现 
1
2
3
long long mul(long long a, long long b, long long m) {
  return (__int128)a * b % m;
}

当然，__int128 的取模运算耗时并不少．如果需要进一步卡常，可以考虑接下来两节介绍的方法．

Barrett 约减

前文提到，除法和取模运算通常比其他四则运算更为耗时．为了减少取模运算的开销，有一些算法可以在不直接做取模的情况下得到相同的结果．本节要介绍的 Barrett 约减算法就是其中之一．

设 \(m\) 为固定模数，假设要对不同的 \(a > 0\) 多次计算 \(a\bmod m\)．由带余除法可知

\[ z = a\bmod m = a - \left\lfloor\dfrac{a}{m}\right\rfloor m. \]

关键在于商数 \(\left\lfloor\dfrac{a}{m}\right\rfloor\) 的计算．设 \(R\) 是某个常数，就有5

\[ \left\lfloor\dfrac{a}{m}\right\rfloor = \left\lfloor a\dfrac{R}{m} / R\right\rfloor \approx \left\lfloor a\left\lfloor\dfrac{R}{m}\right\rfloor/R\right\rfloor. \]

如果选取 \(R = 2^k\)，那么，右式中 \(\left\lfloor\dfrac{R}{m}\right\rfloor\) 可以预处理，除以 \(R\) 的操作可以通过移位运算进行．所以，用右式计算商数，仅需要一次乘法和一次移位操作．再代入 \(a\bmod m\) 的表达式，就得到所求余数的估计 \(z'\)．

现在分析这样做的误差．取整函数 具有性质：对于 \(x > y > 0\) 都有 \(\lfloor x\rfloor - \lfloor y\rfloor \le \lceil x - y\rceil\)．所以，误差

\[ \begin{aligned} \Delta &= |z' - z| = m\left|\left\lfloor\dfrac{a}{m}\right\rfloor - \left\lfloor a\left\lfloor\dfrac{R}{m}\right\rfloor/R\right\rfloor\right|\\ &\le m\left\lceil a\left(\dfrac{R}{m} - \left\lfloor\dfrac{R}{m}\right\rfloor\right) /R\right\rceil \le m\left\lceil\dfrac{a}{R}\right\rceil. \end{aligned} \]

只要 \(a \le R\)，误差 \(\Delta\) 就不超过 \(m\)．由于 \(z' \ge z\)，所以估计值 \(z'\) 只能是 \(z\) 或 \(z + m\)．只要在得到估计值后，在 \(z' \ge m\) 时再减去多余的 \(m\)，就可以保证答案正确．

在 Barrett 约减的计算过程中，仅使用了两次乘法、一次移位操作和至多两次减法，就完成了整数取模．但效率的提升并非毫无成本，实际上，Barrett 约减涉及的中间变量长度往往长于输入变量长度．容易发现，Barrett 约减中涉及的最长中间变量为 \(a\left\lfloor\dfrac{R}{m}\right\rfloor\)．设 \(\ell(x)\) 为整数 \(x\) 的二进制表示长度．那么，有

\[ \ell\left(a\left\lfloor\dfrac{R}{m}\right\rfloor\right) \approx \ell(a) + \ell(R) - \ell(m). \]

由于 \(R\) 的选取需要满足条件 \(a < R\)，这一长度至少为 \(2\ell(a) - \ell(m)\)．但是，需要取模时，一般都有 \(\ell(m)\le\ell(a)\)，因此，这一中间变量的长度可能大于输入长度 \(\ell(a)\)．例如，如果需要将 \(64\) 位整数对 \(32\) 位整数取模，实际上中间变量需要 \(64 \times 2 - 32 = 96\) 位整数．

Barrett 约减的一个应用场景就是计算乘积的余数 \(ab\bmod m\)．如果其中一个乘数固定，比如 \(b\) 固定时，可以通过

\[ ab\bmod m = ab - \left\lfloor a\left\lfloor\dfrac{bR}{m}\right\rfloor/R\right\rfloor m \]

进行与上文类似的估计，只要预处理出 \(\left\lfloor\dfrac{bR}{m}\right\rfloor\) 的值即可．这种 \(b\) 固定的情形有时也称为 Shoup 模乘6．

更为常见的情形是 \(a, b\) 都不固定．此时，需要首先计算 \(ab\) 的值，再利用 Barrett 约减得到 \(ab\bmod m\)．例如，实现模意义下乘法时，需要对 \(0 \le a,b < m\) 计算 \(ab\bmod m\)．此时，选取的 \(r\) 需要满足 \(ab < R\)．根据前文分析，计算过程涉及的最长中间变量长度为 \(2\ell(ab)-\ell(m)\)．当 \(\ell(a)\approx\ell(b)\approx\ell(m)\) 时，该长度为 \(3\ell(m)\)．也就是说，如果要用 Barrett 约减实现 \(32\) 位整数的模乘，中间变量需要 \(96\) 位整数．这也是 Barrett 约减在算法竞赛中实际应用时的一个限制．

作为示例，利用 Barrett 约减实现 32 位有符号整数模乘的参考实现如下：

参考实现 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
// Modular multiplication of int32_t using Barrett reduction.
class Barrett {
  int32_t m;
  uint64_t r;

 public:
  Barrett(int32_t m) : m(m), r((uint64_t)(-m) / m + 1) {}

  // Barrett reduction: a % m.
  int32_t reduce(int64_t a) const {
    int64_t q = (__int128)a * r >> 64;
    a -= q * m;
    return a >= m ? a - m : a;
  }

  // Modular multiplication: (a * b) % m;
  // Assume that 0 <= a, b < m.
  int32_t mul(int32_t a, int32_t b) const { return reduce((int64_t)a * b); }
};

实现中需要用到 128 位整数4．

Montgomery 模乘

Montgomery 模乘算法的功能和 Barrett 算法十分类似，它同样可以减少模整数运算过程中取模运算的开销．与前两个算法都是在近似计算商数不同，Montgomery 模乘将所有整数都映射到 Montgomery 空间上，而 Montgomery 空间中的运算相对容易，进而降低了整体计算成本．

设模数 \(m\) 为奇数，并选取 \(R = 2^k > m\)．那么，同余类 \(a \bmod m\) 对应的 Montgomery 形式就是

\[ aR\bmod m. \]

因为 \(R\perp m\)，所以同余类 \(a \bmod m\) 与它的 Montgomery 形式 \(aR\bmod m\) 之间存在双射．因此，可以将整数转换为 Montgomery 形式后，进行若干模 \(m\) 的运算，再将得到的 Montgomery 形式转换回整数，结果总是正确的．

利用 Montgomery 形式可以方便地进行很多模整数的运算．刚刚已经说明，比较两个同余类是否相同，只要比较它们的 Montgomery 形式．又因为

\[ (a+b)R\bmod m = ((aR\bmod m)\pm(bR\bmod m)) \bmod{m}, \]

所以同余类的加法、减法就对应它们的 Montgomery 形式的加法、减法．但是，要计算同余类的乘法，并不能直接将两个 Montgomery 形式相乘．因为

\[ (ab)R\bmod m = ((aR\bmod m)(bR\bmod m)R^{-1}) \bmod{m}, \]

所以，计算两个 Montgomery 形式的乘法时，需要对它们的乘积 \(x\) 作如下 Montgomery 约减（Montgomery reduction）操作：

\[ \operatorname{REDC}: x \mapsto xR^{-1}\bmod m. \]

利用这一操作，乘积 \(ab\) 的 Montgomery 形式就是 \(\operatorname{REDC}((aR\bmod m)(bR\bmod m))\)．Montgomery 约减操作是 Montgomery 模乘的核心操作：

  • 将 \(a\) 转换为它的 Montgomery 形式就是 \(\operatorname{REDC}((a\bmod m)(R^2\bmod m))\)．
  • 将 \(a\) 的 Montgomery 形式转换回 \(a\bmod m\) 就是 \(\operatorname{REDC}(aR\bmod m)\)．
  • 模逆元 \(a^{-1}\bmod m\) 对应的 Montgomery 形式就是 \(\operatorname{REDC}((aR\bmod m)^{-1}(R^3\bmod m))\)．

现在讨论 Montgomery 约减操作 \(\operatorname{REDC}\) 的实现方法．在计算 \(\operatorname{REDC}(x)\) 时，总是假定 \(0 \le x < m^2\)，这对于以上情形都是成立的．因为 \(R\perp m\)，所以由 裴蜀定理，存在整数 \(R^{-1},m'\) 使得

\[ RR^{-1} + mm' = 1. \]

所以，设 \(q=\lfloor xm' / R\rfloor\)，就有

\[ \begin{aligned} xR^{-1} &= x\dfrac{1 - mm'}{R} \equiv \dfrac{x-xmm' + qmR}{R} = \dfrac{x - m(xm'\bmod R)}{R} \pmod{m}. \end{aligned} \]

因为 \(0 \le x < m^2 < mR\) 且 \(0 \le xm'\bmod R < R\)，所以

\[ -m < \dfrac{x - m(xm'\bmod R)}{R} < m. \]

也就是说，这个商和 \(xR^{-1}\bmod m\) 之间至多差一个 \(m\)．只要在商小于零时，再加上 \(m\) 就可以得到 \(\operatorname{REDC}(x)\)．计算这个商，只需要两次整数乘法、一次整数减法和两次位操作（分别是对 \(R=2^k\) 取模和做除法）．因此，Montgomery 约减操作可以高效进行．

为了进行 Montgomery 模乘操作，需要预处理出一系列常数．首先，Montgomery 约减中会用到 \(m' = m^{-1}\bmod R\)，可以通过 下文 介绍的 Newton–Hensel 方法计算．其次，将不同操作归约为 Montgomery 约减操作时，还涉及诸如 \(R^2\bmod m\) 这样的常数．为了得到它，需要计算一次 \(R\bmod m\)，将它与自身相加就得到 \(2R\bmod m\)．随后，将它看作 \(2\) 的 Montgomery 形式，直接计算快速幂，就可以得到 \(2^kR\bmod m = R^2\bmod m\)．

作为示例，\(32\) 位有符号整数的 Montgomery 模乘实现如下：

参考实现 
 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
// Montgomery modular multiplication.
// The modulus m must be odd. The constant r is 2^32.
class Montgomery {
  int32_t m;
  uint32_t mm, r2;

 public:
  Montgomery(int32_t m) : m(m), mm(1), r2(-m) {
    // Compute mm as inv(m) mod r.
    for (int i = 0; i < 5; ++i) {
      mm *= 2 - mm * m;
    }
    // Compute r2 as r * r mod m.
    // If allowed to use modular operation for uint64_t, simply use:
    //   r2 = (uint64_t)(-m) % m;
    r2 %= m;
    r2 <<= 1;
    if (r2 >= (uint32_t)m) r2 -= m;
    for (int i = 0; i < 5; ++i) {
      r2 = mul(r2, r2);
    }
  }

  // Montgomery reduction: x * inv(r) % m.
  // Also used to transform x from Montgomery space to the normal space.
  int32_t reduce(int64_t x) {
    uint32_t u = (uint32_t)x * mm;
    int32_t ans = (x - (int64_t)m * u) >> 32;
    return ans < 0 ? ans + m : ans;
  }

  // Multiplication in Montgomery space: x * y * inv(r) % m.
  int32_t mul(int32_t x, int32_t y) { return reduce((int64_t)x * y); }

  // Transform x from the normal space to Montgomery space.
  int32_t init(int32_t x) { return mul(x, r2); }
};

相对于 Barrett 约减实现模意义下乘法，Montgomery 模乘的计算涉及转换、Montgomery 形式的乘法、逆转换等多个步骤．因此，只有在转换和逆转换之间的模运算次数足够多时，转换和逆转换的成本才可以摊平，进而获得较高的整体效率．但是，由于 Montgomery 模乘的实现过程中只涉及长度为 \(2\ell(m)\) 的中间变量，所以实现起来更为灵活．例如，\(32\) 位整数的模乘仅需要 \(64\) 位整数的中间变量．所以，如果需要实现一个模整数类用于各种数论计算，Montgomery 模乘更为合适．

模 2 的幂次的整数类

本节讨论模数是 \(2\) 的幂次时，模整数类的实现．在这一特殊情形中，除法和取模运算可以通过位操作实现，计算效率很高．Barrett 约减和 Montgomery 模乘都是利用了 \(2^e\) 作为除数和模数时的这一特性来加速运算．特别地，当模数恰为 \(2^{32}\) 和 \(2^{64}\) 等特殊数字时，可以利用相应位长的无符号整数结合自然溢出实现模整数类，无需任何显式的取模运算．即使模数并非恰好如此，也可以转化为这些特殊模数的情形．例如模数为 \(2^{58}\) 时，可以在模数 \(2^{64}\) 下完成中间计算，最后再将结果对 \(2^{58}\) 取模．除了取模方便外，模 \(2^e\) 整数类的其他操作也有很多特殊实现．本节重点介绍逆元和取幂操作的实现方式．

首先是取逆操作：给定奇数 \(a\) 和模数 \(m=2^e~(e > 2)\)，需要求出 \(a^{-1}\bmod m\)．求逆元的常见方法包括扩展欧几里得算法和快速幂法．扩展欧几里得算法的过程涉及对一般模数取模；普通的快速幂法需要计算 \(a^{\varphi(m)-1}\bmod{m}\)，这需要 \(\Theta(e)\) 次整数乘法．更为高效的方法是 Newton–Hensel 方法．具体地，考虑应用如下结论：7

\[ mx \equiv 1 \pmod{2^e} \implies mx(2 - mx) \equiv 1\pmod{2^{2e}}. \]

根据这一表达式，只要从 \(x = 1\) 开始，反复应用 \(x \gets x(2-mx)\)，就可以在 \(\lceil\log_2 e\rceil\) 次迭代后得到 \(m^{-1}\bmod R\)．

作为示例，模 \(2^{32}\) 整数取逆操作参考实现如下：

参考实现 
1
2
3
4
5
6
7
8
// Compute 1/v mod 2^32 for odd v.
uint32_t inv_mod_2_32(uint32_t v) {
  uint32_t x = 1;
  for (int i = 0; i != 5; ++i) {
    x *= 2 - v * x;
  }
  return x;
}

接下来，讨论取幂操作：给定 \(x,a,b\) 和模数 \(m=2^e~(e > 2)\)，需要求出 \(xa^b\bmod m\)，其中，\(a\) 是奇数．根据对模 \(2^e\) 整数乘法结构的 分析 可知，\(a\) 总是可以写成 \(\pm g^{\ell}\) 的形式8，且负号出现且仅出现在 \(a\equiv 3\pmod 4\) 的情形．对于这种情况，可以将 \(a\) 替换成 \(-a\)，并将最终结果再乘上 \((-1)^b\)．因此，接下来不妨假设 \(a\equiv 1\pmod 4\) 成立．此时，算法的核心想法是，将 \(a\) 写成 \(g^{L(a)}\bmod m\) 的形式，然后用 \(xg^{bL(a)}\bmod m\) 计算所求的幂．

计算 \(L(a)\) 就是计算离散对数 \(\operatorname{ind}_ga\)．注意到，只要 \(a\equiv 1\pmod 4\)，那么 \(a\) 总能写成如下形式：

\[ a \equiv (2^{e_1}+1)(2^{e_2}+1)\cdots(2^{e_s}+1) \pmod{m}, \]

其中，\(1 < e_1 < e_2 < \cdots < e_s < e\)．这是因为直接将这一乘积展开可以发现，\(a\) 的二进制表示中等于 \(1\) 的次低位就是第 \(e_1\) 位（下标从 \(0\) 开始），由此就可以递归地找到这一表示．根据离散对数的 性质 可知，有

\[ 4L(a) \equiv 4L(2^{e_1}+1) + 4L(2^{e_2}+1) + \cdots + 4L(2^{e_s}+1) \pmod{m}. \]

由于离散对数的模数等于阶 \(\delta_m(g)=2^{e-2}=m/4\)，所以此处直接将整个同余式都乘以 \(4\)，以保证计算可以在模 \(m\) 剩余类中进行．由此，只需要对 \(1 < d < e\) 预处理出所有的 \(4L(2^d+1)\)，就可以快速计算 \(4L(a)\) 的值．

反过来，从 \(L(a)\) 也很容易得到 \(g^a\bmod{m}\) 的值．根据 二项式定理 可知，对于 \(1 < d < e\)，都有

\[ \begin{aligned} (2^d+1)^{2^{e-d}} \equiv 1 \pmod{m},\quad (2^d+1)^{2^{e-d-1}} \equiv 1 + 2^{e-1} \pmod{m}, \end{aligned} \]

所以，\(\delta_m(2^d+1) = 2^{e-d}\)．根据阶的性质可知

\[ \delta_m(2^d+1) = \dfrac{\delta_m(g)}{\gcd(\delta_m(g), \operatorname{ind}_g(2^d+1))}. \]

所以，\(\gcd(\delta_m(g), \operatorname{ind}_g(2^d+1)) = 2^{d-2}\)．这说明 \(L(2^d+1) = \operatorname{ind}_g(2^d+1) = 2^{d-2}r\)，其中，\(2\nmid r\)．所以，\(4L(2^d+1)\) 的二进制表示中等于 \(1\) 的最低位恰为第 \(d\) 位（下标从 \(0\) 开始）．因此，同样可以通过二进制表示递归地将 \(4L(a)\) 分解为形如 \(4L(2^d+1)\) 的和．由此，就可以得到 \(a\) 的值．

具体实现时，有一些可以进一步优化的点．首先，将 \(a\) 分解为乘积形式时，还是需要用到除法．更方便的是计算 \(a^{-1}\) 的分解，即寻找 \(1 < e_1 < e_2 < \cdots < e_s < e\) 使得

\[ a(2^{e_1}+1)(2^{e_2}+1)\cdots(2^{e_s}+1) \equiv 1 \pmod{m} \]

成立．同样是通过寻找等于 \(1\) 的次低位来确定 \(e_1\)，但是要在 \(a^{-1}\) 中消去 \(2^{e_1}+1\) 因子，只需要在 \(a\) 上乘以 \(2^{e_1}+1\) 即可，这可以通过位操作进行．又因为 \(4L(a^{-1})=-4L(a)\)，所以统计 \(4L(a)\) 时，需要用减法代替加法．其次，对于特殊选择的基底 \(g\)，迭代无需进行到 \(d = e-1\)，而只要进行到 \(d = \lceil e/2\rceil - 1\) 即可．为此，需要选择 \(g\) 使得

\[ 4L(2^{\lceil e/2\rceil} + 1) = 2^{\lceil e/2\rceil}. \]

对于 \(d \ge e / 2\)，都有

\[ (2^d+1)^2 = 2^{2d} + 2^{d+1} + 1 \equiv 2^{d+1} + 1 \pmod{m}. \]

所以，从 \(d = \lceil e/2\rceil\) 开始归纳可知，\(L(2^d+1)=2^d\) 对于所有 \(d \ge e/2\) 都成立．进而，只要 \(e/2 \le e_1 < e_2 < \cdots < e_s < e\)，就有

\[ (2^{e_1}+1)(2^{e_2}+1)\cdots(2^{e_s}+1) \equiv 1 + 2^{e_1} + 2^{e_2} + \cdots + 2^{e_s} \pmod{m} \]

以及

\[ 4L((2^{e_1}+1)(2^{e_2}+1)\cdots(2^{e_s}+1)) = 2^{e_1} + 2^{e_2} + \cdots + 2^{e_s}. \]

因此，处理完所有 \(d < e/2\) 的二进制位后，可以直接得到剩余部分的离散对数，而无需逐位计算．应用第一个优化后，整个取幂操作只需要 \(O(e)\) 次加减法和位操作和 \(1\) 次乘法操作；应用第二个优化后，可以省去约一半的加减法和位操作，但需要额外 \(1\) 次乘法操作．

作为示例，模 \(2^{32}\) 整数取幂操作参考实现如下：

参考实现 
 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
// Store 4L(a) for a = 2^d + 1, where L(a) is disc. log. base 388251981.
// The first two values are never used and thus set to zero.
// The base is chosen such that 4L(2^16+1) = 2^16.
constexpr uint32_t log_table[16] = {
    0x00000000, 0x00000000, 0xbba0267c, 0x49b9d1e8, 0xf0026f90, 0xd6e17e20,
    0xe78bf840, 0x039fe080, 0xaf7f8100, 0x60fe0200, 0xd1f80400, 0x23e00800,
    0x47801000, 0x8e002000, 0x18004000, 0x20008000,
};

// Compute 4L(v).
uint32_t log_mod_2_32(uint32_t x, uint32_t v) {
  for (int i = 2; i != 16; ++i) {
    if ((v >> i) & 1) {
      v += v << i;
      x -= log_table[i];
    }
  }
  x += v ^ 1;
  return x;
}

// Compute x*a for 4L(a) = v.
uint32_t exp_mod_2_32(uint32_t x, uint32_t v) {
  for (int i = 2; i != 16; ++i) {
    if ((v >> i) & 1) {
      x += x << i;
      v -= log_table[i];
    }
  }
  x *= v ^ 1;
  return x;
}

// Compute x*a^b for odd a.
uint32_t pow_odd_mod_2_32(uint32_t a, uint32_t b, uint32_t x) {
  if (a & 2) {
    a = -a;
    if (b & 1) {
      x = -x;
    }
  }
  return exp_mod_2_32(x, log_mod_2_32(0, a) * b);
}

// Compute x*a^b mod 2^32.
uint32_t pow_mod_2_32(uint32_t a, uint32_t b, uint32_t x = 1) {
  if (!a) return b == 0 ? x : 0;
  auto d = __builtin_ctz(a);
  if ((uint64_t)d * b >= 32) return 0;
  return pow_odd_mod_2_32(a >> d, b, x) << (d * b);
}

离散对数的预处理可以通过 Pohlig–Hellman 算法进行，基底 \(g\) 可以选择为

\[ 5^{\operatorname{ind}_5(2^{\lceil e/2\rceil})/2^{\lceil e/2\rceil - 2}}\bmod{2^e}. \]

参考资料与注释

  1. 这适用于大多数 64 位系统上的 GCC 或 Clang 编译器． ↩

  2. 参见 Double-precision floating-point format - Wikipedia． ↩

  3. 此处用到了条件 \(a < m\)，即 \(a / m \in [0,1)\)． ↩

  4. 在目前的主流编译环境中，只有 Windows 平台上的 MSVC 不支持 __int128 类型．若需要编写可在多平台上兼容的代码，可以通过宏 _MSC_VER 检测 MSVC 编译环境，并在该条件下包含 <intrin.h> 头文件，利用其提供的内建函数（如 _umul128 等）来间接实现 128 位整数运算（仅在 64 位平台上可用）． ↩↩

  5. 此处 \(\left\lfloor\dfrac{r}{m}\right\rfloor\) 也可以替换成 \(\dfrac{r}{m}\) 的其他整数估计，例如上取整函数 \(\left\lceil\dfrac{r}{m}\right\rceil\) 和四舍五入取整函数 \(\left\lfloor\dfrac{r}{m}\right\rceil\) 等，只要相应地调整对估计值的误差修正步骤． ↩

  6. Shoup 在他的数论计算库 NTL 中实现了 Barrett 约减的这一扩展，因此得名． ↩

  7. 直接验证：由 \(mx \equiv 1 \pmod{2^e}\)，可以设 \(mx = 1 + \lambda 2^e\)，那么就有 \(mx(2-mx) = (1+\lambda 2^e)(1-\lambda 2^e) = 1 - \lambda^2 2^{2e} \equiv 1\pmod{2^{2e}}\)． ↩

  8. 文中所引页面仅证明了 \(g\) 可以取 \(5\)．实际上，完全重复该证明，可以说明 \(g\) 可以取任何模 \(8\) 余 \(5\) 的整数．后文会讨论 \(g\) 的选取方法． ↩

Last updated on this page:, Update history
Found an error? Want to help improve? Edit this page on GitHub!
Contributors to this page:383494, buuzzing, c-forrest, cr4c1an, Emp7iness, Enter-tainer, Great-designer, HeRaNO, jifbt, Kaiser-Yang, Koishilll, ksyx, Marcythm, Qiu-Quanzhi, Saisyc, sshwy, StarryReverie, StudyingFather, Tiphereth-A, Xeonacid, xyf007
All content on this page is provided under the terms of the CC BY-SA 4.0 and SATA license, additional terms may apply