Skip to content

Introduction to Modular Arithmetic

In algorithm competitions, an important part of number theory is modular arithmetic, which involves performing various integer operations under a given modulus. In addition to basic arithmetic operations and exponentiation, it can also be conveniently used for operations such as taking logarithms, extracting roots, and computing factorials and binomial coefficients.

Modular arithmetic appears in various types of problems, not limited to number theory. The actual answer to many problems may be very large, exceeding the storage range of common integer types. To avoid introducing big integer arithmetic and outputting long numbers, problems often require outputting the answer modulo a given number. This requires proficiency in various modular arithmetic techniques.

Integer Division and Modulo Operations in C/C++

In C/C++, integer division and modulo operations differ from the mathematical conventions for division and modulo.

For all standard versions of C/C++, integer division is defined as follows:

  1. When the divisor is 0, the behavior is undefined;
  2. Otherwise, (a / b) * b + a % b evaluates to a.

That is, the sign of the modulo operation depends on how division rounds; and how division rounds is implementation-defined (determined by the compiler).

Starting from the C99 and C++11 standards, it is specified that division truncates toward zero (discarding the fractional part); the sign of the modulo result is the same as the sign of the dividend. Since then, the following assertions are guaranteed to be true:

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

Modular Integer Class

Modular arithmetic can be viewed as performing various operations on congruence classes under a given modulus. If a structure is used to represent a congruence class, and operations such as addition, subtraction, and multiplication between congruence classes are encapsulated as methods or operator overloads of the structure, then modular arithmetic can be naturally implemented as a modular integer class. Below is a simple example that supports addition, subtraction, multiplication, and fast exponentiation for 32-bit signed integers under a modulus \(M < 2^{30}\):

A Simple Modular Integer Class 
 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;
  }
};

This implementation intentionally reduces the number of modulo operations because modulo operations are typically much more time-consuming than ordinary addition, subtraction, multiplication, or comparison operations. The code comments provide equivalent and more straightforward implementations. The main idea behind these simple optimizations is that when adding or subtracting two integers in \([0,M)\), the result must lie in the interval \((-M,2M)\), so it can be adjusted back to \([0,M)\) with one addition or subtraction. The exponentiation operation in this implementation uses the fast exponentiation technique.

In addition to these basic operations, the following operations can be performed under various moduli:

These operations are usually easier under prime moduli. For composite moduli, extended versions of the corresponding algorithms and the Chinese Remainder Theorem are often required. Most of these modular operations can be viewed as solving certain congruence equations. For general methods of solving congruence equations, refer to the Congruence Equation page.

This section introduces several methods for optimizing modulo, multiplication, and fast exponentiation operations under a modulus. For most problems, the simple implementation provided earlier is efficient enough. However, when a problem has strict requirements on algorithm constants, these optimization methods can be useful by reducing unnecessary computations and modulo operations, further lowering the algorithm's time overhead.

Fast Multiplication

In primality testing and integer factorization, multiplication modulo a modulus within the long long range is often encountered. To avoid integer overflow during computation, this section introduces a "fast multiplication" method that can handle moduli within the long long range without using __int128 and has complexity \(O(1)\). This algorithm requires that the evaluation system represents long double as at least an 80-bit extended precision floating-point number1.

Assume \(0 \le a, b < m\), and we want to compute \(ab\bmod m\). Note that:

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

Using natural overflow of 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}. \]

As long as we can compute the quotient \(\left\lfloor\dfrac{ab}m\right\rfloor\), the multiplication and subtraction in the rightmost expression can be computed directly using unsigned long long.

Next, we only need to consider how to compute \(\left\lfloor\dfrac {ab}m\right\rfloor\). The solution is to first compute \(\dfrac am\) using long double and then multiply by \(b\). Since long double is used, there will inevitably be precision errors. Assuming long double is represented as an 80-bit extended precision floating-point number (i.e., 1 sign bit, 15 exponent bits, 64 mantissa bits), then long double can represent at most 64 bits of precision2. So \(\dfrac am\) has an error starting at the 65th bit at worst, with an error range3 of \(\left(-2^{-64},2^{-64}\right)\). Multiplying by \(b\), a 64-bit signed integer, gives an error range of \((-0.5,0.5)\). To simplify subsequent discussion, we can add \(0.5\) and then floor, resulting in a final error range of \(\{0,1\}\).

Finally, when substituting into the above formula, we need to multiply by \(-m\), so the final error range is \(\{0,-m\}\). Since \(m\) is within the long long range, if the result \(r\in[0,m)\), return \(r\) directly; otherwise, return \(r+m\).

The code implementation is as follows:

Reference Implementation 
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;
}

Nowadays, most evaluation systems' C/C++ compilers support the __int128 type4, so we can also directly promote the multiplier type to __int128 and then compute the modulo:

Reference Implementation 
1
2
3
long long mul(long long a, long long b, long long m) {
  return (__int128)a * b % m;
}

Of course, the modulo operation for __int128 is not cheap. If further constant optimization is needed, consider the methods introduced in the next two sections.

Barrett Reduction

As mentioned earlier, division and modulo operations are usually more time-consuming than other arithmetic operations. To reduce the overhead of modulo operations, some algorithms can achieve the same result without directly performing modulo. The Barrett reduction algorithm introduced in this section is one of them.

Let \(m\) be a fixed modulus, and suppose we need to compute \(a\bmod m\) multiple times for different \(a > 0\). By the division algorithm:

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

The key is the computation of the quotient \(\left\lfloor\dfrac{a}{m}\right\rfloor\). Let \(R\) be some constant, then5

\[ \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. \]

If we choose \(R = 2^k\), then \(\left\lfloor\dfrac{R}{m}\right\rfloor\) can be preprocessed, and dividing by \(R\) can be done via bit shifting. So computing the quotient using the right-hand side requires only one multiplication and one shift operation. Substituting into the expression for \(a\bmod m\) gives an estimate \(z'\) for the remainder.

Now analyze the error. The floor function has the property: for \(x > y > 0\), \(\lfloor x\rfloor - \lfloor y\rfloor \le \lceil x - y\rceil\). So the error

\[ \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} \]

As long as \(a \le R\), the error \(\Delta\) does not exceed \(m\). Since \(z' \ge z\), the estimate \(z'\) can only be \(z\) or \(z + m\). After obtaining the estimate, if \(z' \ge m\), subtract the extra \(m\) to ensure the correct answer.

In the Barrett reduction computation, only two multiplications, one shift operation, and at most two subtractions are needed to perform integer modulo. However, the efficiency gain is not without cost. In fact, the intermediate variables involved in Barrett reduction are often longer than the input variables. It is easy to see that the longest intermediate variable in Barrett reduction is \(a\left\lfloor\dfrac{R}{m}\right\rfloor\). Let \(\ell(x)\) be the binary length of integer \(x\). Then,

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

Since \(R\) must satisfy \(a < R\), this length is at least \(2\ell(a) - \ell(m)\). However, when modulo is needed, we generally have \(\ell(m)\le\ell(a)\), so this intermediate variable may be longer than the input length \(\ell(a)\). For example, if we need to compute a 64-bit integer modulo a 32-bit integer, the intermediate variable actually requires \(64 \times 2 - 32 = 96\) bits.

One application of Barrett reduction is computing the remainder of a product \(ab\bmod m\). If one of the multiplicands is fixed, say \(b\) is fixed, we can estimate similarly via

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

as long as \(\left\lfloor\dfrac{bR}{m}\right\rfloor\) is preprocessed. This case where \(b\) is fixed is sometimes called Shoup modular multiplication6.

A more common scenario is when neither \(a\) nor \(b\) is fixed. In this case, we first compute \(ab\) and then use Barrett reduction to get \(ab\bmod m\). For example, when implementing modular multiplication, we need to compute \(ab\bmod m\) for \(0 \le a,b < m\). Here, the chosen \(r\) must satisfy \(ab < R\). According to the previous analysis, the longest intermediate variable involved has length \(2\ell(ab)-\ell(m)\). When \(\ell(a)\approx\ell(b)\approx\ell(m)\), this length is \(3\ell(m)\). That is, if we use Barrett reduction to implement modular multiplication for 32-bit integers, the intermediate variable requires 96 bits. This is also a limitation of Barrett reduction in practical algorithm competition applications.

As an example, a reference implementation of 32-bit signed integer modular multiplication using Barrett reduction is as follows:

Reference Implementation 
 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); }
};

The implementation requires 128-bit integers4.

Montgomery Modular Multiplication

The Montgomery modular multiplication algorithm is very similar in function to the Barrett algorithm; it also reduces the overhead of modulo operations during modular integer arithmetic. Unlike the previous two algorithms that approximate the quotient, Montgomery modular multiplication maps all integers to Montgomery space, where operations are relatively easier, thereby reducing overall computational cost.

Let the modulus \(m\) be odd, and choose \(R = 2^k > m\). Then, the Montgomery form corresponding to the congruence class \(a \bmod m\) is

\[ aR\bmod m. \]

Since \(R\perp m\), there is a bijection between the congruence class \(a \bmod m\) and its Montgomery form \(aR\bmod m\). Therefore, we can convert integers to Montgomery form, perform several modulo \(m\) operations, and then convert the resulting Montgomery form back to an integer, and the result will always be correct.

Many modular integer operations are convenient in Montgomery form. As already mentioned, comparing whether two congruence classes are equal only requires comparing their Montgomery forms. Also, since

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

addition and subtraction of congruence classes correspond to addition and subtraction of their Montgomery forms. However, to compute the multiplication of congruence classes, we cannot directly multiply two Montgomery forms. Because

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

when computing the product of two Montgomery forms, we need to perform the following Montgomery reduction operation on their product \(x\):

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

Using this operation, the Montgomery form of the product \(ab\) is \(\operatorname{REDC}((aR\bmod m)(bR\bmod m))\). The Montgomery reduction operation is the core of Montgomery modular multiplication:

  • Converting \(a\) to its Montgomery form is \(\operatorname{REDC}((a\bmod m)(R^2\bmod m))\).
  • Converting \(a\)'s Montgomery form back to \(a\bmod m\) is \(\operatorname{REDC}(aR\bmod m)\).
  • The Montgomery form corresponding to the modular inverse \(a^{-1}\bmod m\) is \(\operatorname{REDC}((aR\bmod m)^{-1}(R^3\bmod m))\).

Now discuss the implementation of the Montgomery reduction operation \(\operatorname{REDC}\). When computing \(\operatorname{REDC}(x)\), we always assume \(0 \le x < m^2\), which holds for the above cases. Since \(R\perp m\), by Bézout's identity, there exist integers \(R^{-1},m'\) such that

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

So, let \(q=\lfloor xm' / R\rfloor\), then

\[ \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} \]

Since \(0 \le x < m^2 < mR\) and \(0 \le xm'\bmod R < R\), we have

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

That is, this quotient differs from \(xR^{-1}\bmod m\) by at most \(m\). As long as we add \(m\) when the quotient is negative, we obtain \(\operatorname{REDC}(x)\). Computing this quotient requires only two integer multiplications, one integer subtraction, and two bit operations (modulo \(R=2^k\) and division, respectively). Therefore, Montgomery reduction can be performed efficiently.

To perform Montgomery modular multiplication, a series of constants need to be preprocessed. First, \(m' = m^{-1}\bmod R\) used in Montgomery reduction can be computed via the Newton–Hensel method introduced below. Second, when reducing different operations to Montgomery reduction, constants such as \(R^2\bmod m\) are involved. To obtain it, we need to compute \(R\bmod m\) once, and adding it to itself gives \(2R\bmod m\). Then, treating it as the Montgomery form of \(2\), directly computing fast exponentiation yields \(2^kR\bmod m = R^2\bmod m\).

As an example, the implementation of Montgomery modular multiplication for 32-bit signed integers is as follows:

Reference Implementation 
 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); }
};

Compared to Barrett reduction for implementing modular multiplication, Montgomery modular multiplication involves multiple steps: conversion, multiplication in Montgomery form, inverse conversion, etc. Therefore, only when the number of modular operations between conversion and inverse conversion is sufficiently large can the costs of conversion and inverse conversion be amortized, leading to higher overall efficiency. However, since the implementation of Montgomery modular multiplication only involves intermediate variables of length \(2\ell(m)\), it is more flexible. For example, modular multiplication for 32-bit integers only requires 64-bit integer intermediate variables. Thus, if a modular integer class is needed for various number theory computations, Montgomery modular multiplication is more suitable.

Modular Integer Class for Powers of Two

This section discusses the implementation of a modular integer class when the modulus is a power of two. In this special case, division and modulo operations can be implemented via bit operations, resulting in high computational efficiency. Both Barrett reduction and Montgomery modular multiplication exploit this property to accelerate operations when \(2^e\) is used as the divisor or modulus. In particular, when the modulus is exactly \(2^{32}\), \(2^{64}\), etc., the modular integer class can be implemented using unsigned integers of the corresponding bit length combined with natural overflow, without any explicit modulo operations. Even if the modulus is not exactly such a number, it can be transformed into these special modulus cases. For example, when the modulus is \(2^{58}\), intermediate computations can be performed under modulus \(2^{64}\), and the final result can be taken modulo \(2^{58}\). Besides convenient modulo operations, other operations for modular \(2^e\) integer classes also have special implementations. This section focuses on the implementation of modular inverse and exponentiation operations.

First, the modular inverse operation: given an odd integer \(a\) and modulus \(m=2^e~(e > 2)\), we need to compute \(a^{-1}\bmod m\). Common methods for computing modular inverses include the extended Euclidean algorithm and fast exponentiation. The extended Euclidean algorithm involves modulo operations with a general modulus; ordinary fast exponentiation requires computing \(a^{\varphi(m)-1}\bmod{m}\), which requires \(\Theta(e)\) integer multiplications. A more efficient method is the Newton–Hensel method. Specifically, consider applying the following result:7

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

According to this expression, starting from \(x = 1\) and repeatedly applying \(x \gets x(2-mx)\), we can obtain \(m^{-1}\bmod R\) after \(\lceil\log_2 e\rceil\) iterations.

As an example, a reference implementation for modular inverse under \(2^{32}\) is as follows:

Reference Implementation

```cpp // 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; }

$$

Next, discuss the exponentiation operation: given \(x,a,b\) and modulus \(m=2^e~(e > 2)\), we need to compute \(xa^b\bmod m\), where \(a\) is odd. According to the analysis of the multiplicative structure of integers modulo \(2^e\), \(a\) can always be written in the form \(\pm g^{\ell}\)8, and the negative sign appears if and only if \(a\equiv 3\pmod 4\). For this case, we can replace \(a\) with \(-a\) and multiply the final result by \((-1)^b\). Therefore, we can assume \(a\equiv 1\pmod 4\) holds. The core idea of the algorithm is to write \(a\) as \(g^{L(a)}\bmod m\) and then compute the desired power as \(xg^{bL(a)}\bmod m\).

Computing \(L(a)\) is computing the discrete logarithm \(\operatorname{ind}_ga\). Note that as long as \(a\equiv 1\pmod 4\), \(a\) can always be written in the form:

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

where \(1 < e_1 < e_2 < \cdots < e_s < e\). This is because expanding this product reveals that the second least significant bit equal to \(1\) in the binary representation of \(a\) is exactly the \(e_1\)-th bit (index starting from \(0\)), allowing recursive discovery of this representation. According to the properties of discrete logarithms, we have

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

Since the modulus of the discrete logarithm equals the order \(\delta_m(g)=2^{e-2}=m/4\), we multiply the entire congruence by \(4\) to ensure computation can be performed in the residue classes modulo \(m\). Thus, by preprocessing all \(4L(2^d+1)\) for \(1 < d < e\), we can quickly compute the value of \(4L(a)\).

Conversely, it is easy to obtain the value of \(g^a\bmod{m}\) from \(L(a)\). According to the binomial theorem, for \(1 < d < e\), we have

\[ \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} \]

so \(\delta_m(2^d+1) = 2^{e-d}\). By properties of order,

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

Thus, \(\gcd(\delta_m(g), \operatorname{ind}_g(2^d+1)) = 2^{d-2}\). This implies \(L(2^d+1) = \operatorname{ind}_g(2^d+1) = 2^{d-2}r\), where \(2\nmid r\). So the least significant bit equal to \(1\) in the binary representation of \(4L(2^d+1)\) is exactly the \(d\)-th bit (index starting from \(0\)). Therefore, we can also recursively decompose \(4L(a)\) into sums of the form \(4L(2^d+1)\) via binary representation. From this, we can obtain the value of \(a\).

In concrete implementation, there are further optimizations. First, when decomposing \(a\) into product form, division is still needed. It is more convenient to compute the decomposition of \(a^{-1}\), i.e., find \(1 < e_1 < e_2 < \cdots < e_s < e\) such that

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

holds. Again, we determine \(e_1\) by finding the second least significant bit equal to \(1\), but to eliminate the factor \(2^{e_1}+1\) from \(a^{-1}\), we only need to multiply \(a\) by \(2^{e_1}+1\), which can be done via bit operations. Since \(4L(a^{-1})=-4L(a)\), subtraction is used instead of addition when accumulating \(4L(a)\). Second, for specially chosen base \(g\), iteration need not proceed to \(d = e-1\), but only to \(d = \lceil e/2\rceil - 1\). For this, we need to choose \(g\) such that

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

For \(d \ge e / 2\), we have

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

So, by induction starting from \(d = \lceil e/2\rceil\), \(L(2^d+1)=2^d\) holds for all \(d \ge e/2\). Furthermore, as long as \(e/2 \le e_1 < e_2 < \cdots < e_s < e\), we have

\[ (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} \]

and

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

Therefore, after processing all binary bits with \(d < e/2\), the discrete logarithm for the remaining part can be obtained directly without bit-by-bit computation. Applying the first optimization, the entire exponentiation operation requires only \(O(e)\) additions/subtractions and bit operations and \(1\) multiplication operation; applying the second optimization saves about half the additions/subtractions and bit operations but requires an additional \(1\) multiplication operation.

As an example, a reference implementation for exponentiation under \(2^{32}\) is as follows:

Reference Implementation

```cpp // 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 xa 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); }

$$

Preprocessing of discrete logarithms can be done via the Pohlig–Hellman algorithm, and the base \(g\) can be chosen as

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

References and Notes

  1. This applies to GCC or Clang compilers on most 64-bit systems. ↩

  2. See Double-precision floating-point format - Wikipedia. ↩

  3. Here we use the condition \(a < m\), i.e., \(a / m \in [0,1)\). ↩

  4. In current mainstream compilation environments, only MSVC on Windows does not support the __int128 type. If cross-platform compatible code is needed, the MSVC compilation environment can be detected via the macro _MSC_VER, and under that condition, include the <intrin.h> header file and use its built-in functions (such as _umul128, etc.) to indirectly implement 128-bit integer operations (only available on 64-bit platforms). ↩↩

  5. Here \(\left\lfloor\dfrac{r}{m}\right\rfloor\) can also be replaced by other integer estimates of \(\dfrac{r}{m}\), such as the ceiling function \(\left\lceil\dfrac{r}{m}\right\rceil\) or rounding function \(\left\lfloor\dfrac{r}{m}\right\rceil\), etc., as long as the error correction step for the estimate is adjusted accordingly. ↩

  6. Shoup implemented this extension of Barrett reduction in his number theory library NTL, hence the name. ↩

  7. Direct verification: from \(mx \equiv 1 \pmod{2^e}\), let \(mx = 1 + \lambda 2^e\), then \(mx(2-mx) = (1+\lambda 2^e)(1-\lambda 2^e) = 1 - \lambda^2 2^{2e} \equiv 1\pmod{2^{2e}}\). ↩

  8. The referenced page only proves that \(g\) can be taken as \(5\). Actually, repeating the proof entirely shows that \(g\) can be taken as any integer congruent to \(5\) modulo \(8\). The selection of \(g\) will be discussed later. ↩

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用