模算术简介

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

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

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

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

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

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

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

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

5 % 3 == 2;
5 % -3 == 2;
-5 % 3 == -2;
-5 % -3 == -2;

模整数类

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

一个简单的模整数类

// 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, 2 M)$ 内，因此可以通过一次加减法调整回区间 $[0, M)$ ．该实现中的取幂运算用到了快速幂的技巧．

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

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

相关算法

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

快速乘

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

假设 $0 \leq a, b < m$ ，要计算 $a b mod m$ ．注意到：

a b mod m = a b - ⌊ \frac{a b}{m} ⌋ m .

利用 unsigned long long 的自然溢出：

a b mod m = a b - ⌊ \frac{a b}{m} ⌋ m = (a b - ⌊ \frac{a b}{m} ⌋ m) mod 2^{64} .

只要能算出商 $⌊ \frac{a b}{m} ⌋$ ，最右侧表达式中的乘法和减法运算都可以使用 unsigned long long 直接计算．

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

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

代码实现如下：

参考实现

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 类型⁴，因此也可以直接将乘数类型提升至 __int128 后取模计算：

参考实现

long long mul(long long a, long long b, long long m) {
  return (__int128)a * b % m;
}

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

Barrett 约减

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

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

z = a mod m = a - ⌊ \frac{a}{m} ⌋ m .

关键在于商数 $⌊ \frac{a}{m} ⌋$ 的计算．设 $R$ 是某个常数，就有⁵

⌊ \frac{a}{m} ⌋ = ⌊ a \frac{R}{m} / R ⌋ \approx ⌊ a ⌊ \frac{R}{m} ⌋ / R ⌋ .

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

现在分析这样做的误差．因为对于 $x > y > 0$ 都有 $⌊ x ⌋ - ⌊ y ⌋ \leq ⌈ x - y ⌉$ ，所以误差

\begin{aligned} Δ & = | z^{'} - z | = m | ⌊ \frac{a}{m} ⌋ - ⌊ a ⌊ \frac{R}{m} ⌋ / R ⌋ | \\ \leq m ⌈ a (\frac{R}{m} - ⌊ \frac{R}{m} ⌋) / R ⌉ \leq m ⌈ \frac{a}{R} ⌉ . \end{aligned}

只要 $a \leq R$ ，误差 $Δ$ 就不超过 $m$ ．由于 $z^{'} \geq z$ ，所以估计值 $z^{'}$ 只能是 $z$ 或 $z + m$ ．只要在得到估计值后，在 $z^{'} \geq m$ 时再减去多余的 $m$ ，就可以保证答案正确．

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

ℓ (a ⌊ \frac{R}{m} ⌋) \approx ℓ (a) + ℓ (R) - ℓ (m) .

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

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

a b mod m = a b - ⌊ a ⌊ \frac{b R}{m} ⌋ / R ⌋ m

进行与上文类似的估计，只要预处理出 $⌊ \frac{b R}{m} ⌋$ 的值即可．这种 $b$ 固定的情形有时也称为 Shoup 模乘⁶．

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

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

参考实现

// 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 位整数⁴．

Montgomery 模乘

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

设模数 $m$ 为奇数，并选取 $R = 2^{k} > m$ ．那么，同余类 $a mod m$ 对应的 Montgomery 形式就是

a R mod m .

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

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

(a + b) R mod m = ((a R mod m) \pm (b R mod m)) mod m,

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

(a b) R mod m = ((a R mod m) (b R mod m) R^{- 1}) mod m,

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

REDC : x \mapsto x R^{- 1} mod m .

利用这一操作，乘积 $a b$ 的 Montgomery 形式就是 $REDC ((a R mod m) (b R mod m))$ ．Montgomery 约减操作是 Montgomery 模乘的核心操作：

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

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

R R^{- 1} + m m^{'} = 1.

所以，设 $q = ⌊ x m^{'} / R ⌋$ ，就有

\begin{aligned} x R^{- 1} & = x \frac{1 - m m^{'}}{R} \equiv \frac{x - x m m^{'} + q m R}{R} = \frac{x - m (x m^{'} mod R)}{R} (\mod m) . \end{aligned}

因为 $0 \leq x < m^{2} < m R$ 且 $0 \leq x m^{'} mod R < R$ ，所以

- m < \frac{x - m (x m^{'} mod R)}{R} < m .

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

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

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

参考实现

// 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 ℓ (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} mod m$ ．求逆元的常见方法包括扩展欧几里得算法和快速幂法．扩展欧几里得算法的过程涉及对一般模数取模；普通的快速幂法需要计算 $a^{φ (m) - 1} mod m$ ，这需要 $Θ (e)$ 次整数乘法．更为高效的方法是 Newton–Hensel 方法．具体地，考虑应用如下结论：⁷

m x \equiv 1 (\mod 2^{e}) ⟹ m x (2 - m x) \equiv 1 (\mod 2^{2 e}) .

根据这一表达式，只要从 $x = 1$ 开始，反复应用 $x \leftarrow x (2 - m x)$ ，就可以在 $⌈ \log_{2} e ⌉$ 次迭代后得到 $m^{- 1} mod R$ ．

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

参考实现

// 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)$ ，需要求出 $x a^{b} mod m$ ，其中， $a$ 是奇数．根据对模 $2^{e}$ 整数乘法结构的分析可知， $a$ 总是可以写成 $\pm g^{ℓ}$ 的形式⁸，且负号出现且仅出现在 $a \equiv 3 (\mod 4)$ 的情形．对于这种情况，可以将 $a$ 替换成 $- a$ ，并将最终结果再乘上 $(- 1)^{b}$ ．因此，接下来不妨假设 $a \equiv 1 (\mod 4)$ 成立．此时，算法的核心想法是，将 $a$ 写成 $g^{L (a)} mod m$ 的形式，然后用 $x g^{b L (a)} mod m$ 计算所求的幂．

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

a \equiv (2^{e_{1}} + 1) (2^{e_{2}} + 1) \dots (2^{e_{s}} + 1) (\mod m),

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

4 L (a) \equiv 4 L (2^{e_{1}} + 1) + 4 L (2^{e_{2}} + 1) + \dots + 4 L (2^{e_{s}} + 1) (\mod m) .

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

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

\begin{array}{r} (2^{d} + 1)^{2^{e - d}} \equiv 1 (\mod m), (2^{d} + 1)^{2^{e - d - 1}} \equiv 1 + 2^{e - 1} (\mod m), \end{array}

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

δ_{m} (2^{d} + 1) = \frac{δ_{m} (g)}{gcd (δ_{m} (g), {ind}_{g} (2^{d} + 1))} .

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

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

a (2^{e_{1}} + 1) (2^{e_{2}} + 1) \dots (2^{e_{s}} + 1) \equiv 1 (\mod m)

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

4 L (2^{⌈ e / 2 ⌉} + 1) = 2^{⌈ e / 2 ⌉} .

对于 $d \geq e / 2$ ，都有

(2^{d} + 1)^{2} = 2^{2 d} + 2^{d + 1} + 1 \equiv 2^{d + 1} + 1 (\mod m) .

所以，从 $d = ⌈ e / 2 ⌉$ 开始归纳可知， $L (2^{d} + 1) = 2^{d}$ 对于所有 $d \geq e / 2$ 都成立．进而，只要 $e / 2 \leq e_{1} < e_{2} < \dots < e_{s} < e$ ，就有

(2^{e_{1}} + 1) (2^{e_{2}} + 1) \dots (2^{e_{s}} + 1) \equiv 1 + 2^{e_{1}} + 2^{e_{2}} + \dots + 2^{e_{s}} (\mod m)

以及

4 L ((2^{e_{1}} + 1) (2^{e_{2}} + 1) \dots (2^{e_{s}} + 1)) = 2^{e_{1}} + 2^{e_{2}} + \dots + 2^{e_{s}} .

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

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

参考实现

// 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^{{ind}_{5} (2^{⌈ e / 2 ⌉}) / 2^{⌈ e / 2 ⌉ - 2}} mod 2^{e} .

参考资料与注释

Fast modular multiplication by orz - Codeforces
Barrett Reduction - Wikipedia
Barrett Reduction - A41
Barrett 约减原理及正确性证明 by Chen - 知乎专栏
Montgomery Multiplication - CP Algorithms
Montgomery 模乘 by Chen - 知乎专栏
Binary Exponentiation by Factoring - CP Algorithms
Barrett, Paul. "Implementing the Rivest Shamir and Adleman public key encryption algorithm on a standard digital signal processor." In Conference on the Theory and Application of Cryptographic Techniques, pp. 311-323. Berlin, Heidelberg: Springer Berlin Heidelberg, 1986.
Becker, Hanno, Vincent Hwang, Matthias J. Kannwischer, Bo-Yin Yang, and Shang-Yi Yang. "Neon NTT: Faster Dilithium, Kyber, and Saber on Cortex-A72 and Apple M1." IACR Transactions on Cryptographic Hardware and Embedded Systems (2022): 221-244.
Montgomery, Peter L. "Modular multiplication without trial division." Mathematics of computation 44, no. 170 (1985): 519-521.

这适用于大多数 64 位系统上的 GCC 或 Clang 编译器． ↩
参见 Double-precision floating-point format - Wikipedia． ↩
此处用到了条件 $a < m$ ，即 $a / m \in [0, 1)$ ． ↩
在目前的主流编译环境中，只有 Windows 平台上的 MSVC 不支持 __int128 类型．若需要编写可在多平台上兼容的代码，可以通过宏 _MSC_VER 检测 MSVC 编译环境，并在该条件下包含 <intrin.h> 头文件，利用其提供的内建函数（如 _umul128 等）来间接实现 128 位整数运算（仅在 64 位平台上可用）． ↩↩
此处 $⌊ \frac{r}{m} ⌋$ 也可以替换成 $\frac{r}{m}$ 的其他整数估计，例如上取整函数 $⌈ \frac{r}{m} ⌉$ 和四舍五入取整函数 $⌊ \frac{r}{m} ⌉$ 等，只要相应地调整对估计值的误差修正步骤． ↩
Shoup 在他的数论计算库 NTL 中实现了 Barrett 约减的这一扩展，因此得名． ↩
直接验证：由 $m x \equiv 1 (\mod 2^{e})$ ，可以设 $m x = 1 + λ 2^{e}$ ，那么就有 $m x (2 - m x) = (1 + λ 2^{e}) (1 - λ 2^{e}) = 1 - λ^{2} 2^{2 e} \equiv 1 (\mod 2^{2 e})$ ． ↩
文中所引页面仅证明了 $g$ 可以取 $5$ 。实际上，完全重复该证明，可以说明 $g$ 可以取任何模 $8$ 余 $5$ 的整数。后文会讨论 $g$ 的选取方法。 ↩

本页面最近更新：2025/10/30 13:44:01，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：c-forrest, 383494, buuzzing, cr4c1an, Emp7iness, Enter-tainer, Great-designer, HeRaNO, jifbt, Kaiser-Yang, Koishilll, ksyx, Marcythm, Qiu-Quanzhi, Saisyc, sshwy, StarryReverie, StudyingFather, Tiphereth-A, Xeonacid, xyf007
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用