Introduction to Modular Arithmetic

In algorithmic competitions, an important part of number theory is modular arithmetic (modular arithmetic), that is, performing various integer operations modulo a fixed modulus. Besides basic arithmetic and exponentiation, one can also conveniently compute discrete logarithms, roots of various degrees, factorials, binomial coefficients, and more.

Modular arithmetic appears in all kinds of problems, not only in 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 printing long numbers, problems often ask for the answer modulo some number. This makes it important to master various modular-arithmetic techniques.

Integer Division and Remainder in C/C++

In C/C++, integer division and remainder do not follow the usual mathematical conventions for modular arithmetic and division.

For every standard version of C/C++, integer division is specified as follows:

1. When the divisor is zero, the behavior is undefined.
2. Otherwise, the result of (a / b) * b + a % b is equal to a.

In other words, the sign of the remainder depends on how division rounds, and the way division rounds is implementation-defined (determined by the compiler).

Starting with the C99 and C++11 standards, the quotient is specified to round toward zero (discarding the fractional part), so the remainder has the same sign as the dividend. Therefore, the following assertions are guaranteed to be true:

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

ModInt Class

Modular arithmetic can be viewed as performing operations on congruence classes modulo a fixed number. If a structure is used to represent a congruence class, and addition, subtraction, multiplication, and other operations between congruence classes are encapsulated as structure methods or overloaded operators, modular arithmetic can naturally be implemented as a ModInt class. The following is a simple example supporting addition, subtraction, multiplication, and fast exponentiation of $M < 2^{30}$-bit signed integers for moduli $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;
  }
};

This implementation deliberately reduces the number of remainder operations, because remainder operations are usually much slower than ordinary addition, subtraction, multiplication, or comparison. The code comments provide equivalent, more direct implementations. The main idea of these simple optimizations is that when two integers in $[0,M)$ are added or subtracted, the result always lies in $(-M,2M)$, so one addition or subtraction can bring it back into $[0,M)$. The exponentiation implementation uses the binary exponentiation technique.

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

• Multiplicative inverse
• Division
• Factorial
• Binomial coefficients
• Square roots
• Discrete logarithms
• Roots

These operations are usually easier for prime moduli. For composite moduli, extended versions of the corresponding algorithms and the Chinese remainder theorem are often needed. Most modular operations can be viewed as solving some kind of congruence equation. For general methods of solving congruence equations, see the congruence equations page.

Related Algorithms

This section introduces several methods for optimizing remainder, multiplication, and fast-exponentiation operations modulo a number. For most problems, the simple implementation given above is already efficient enough. However, when a problem imposes strict requirements on constant factors, these optimizations can reduce unnecessary computation and remainder operations, further lowering the running time.

Fast Multiplication

In primality testing and integer factorization, one often encounters modular multiplication with moduli in the long long range. To avoid integer overflow during the computation, this section introduces a “fast multiplication” method that handles moduli in the long long range without using __int128, in $O(1)$ time. This algorithm requires the judging system to represent long double with at least $80$-bit extended precision floating point¹.

Assume $0 \le a, b < m$. We want to compute $ab\bmod m$. Notice that:

Using the natural overflow of unsigned long long:

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

It remains to consider how to compute $\left\lfloor\dfrac {ab}m\right\rfloor$. The solution is to first use long double to compute $\dfrac am$ and then multiply by $b$. Since long double is being used, there will inevitably be precision error. Suppose that long double is an $80$-bit extended-precision floating-point type (with a $1$-bit sign, a $15$-bit exponent, and a $64$-bit significand). Then long double can represent at most $64$ significant bits exactly². Therefore, $\dfrac am$ may first become wrong at bit $65$, with an error range³ of $\left(-2^{-64},2^{-64}\right)$. After multiplying by the $b$-bit signed integer $64$, the error range is $(-0.5,0.5)$. To simplify the discussion, we can first add $0.5$ and then round; the final error range is $\{0,1\}$.

Finally, substituting into the expression above requires multiplication by $-m$, so the final error range is $\{0,-m\}$. Since $m$ is in the long long range, when the result lies in $r\in[0,m)$, return $r$ directly; otherwise return $r+m$.

The implementation is as follows:

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

Most judging systems now provide C/C++ compilers supporting the __int128 type⁴, so the multiplicand type can also be promoted to __int128 before taking the remainder:

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

Of course, taking the remainder of __int128 values is not particularly fast. If further constant-factor optimization is needed, the methods in the next two sections can be considered.

Barrett Reduction

As mentioned above, division and remainder operations are usually slower than the other four arithmetic operations. To reduce the cost of remainder operations, some algorithms can obtain the same result without directly taking a remainder. 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 > 0$ for different positive values $a\bmod m$. By the division algorithm:

The key is computing the quotient $\left\lfloor\dfrac{a}{m}\right\rfloor$. Let $R$ be a constant. Then⁵

If we choose $R = 2^k$, then in the right-hand side, $\left\lfloor\dfrac{R}{m}\right\rfloor$ can be precomputed, and division by $R$ can be performed by a shift operation. Thus, computing the quotient with the right-hand side requires only one multiplication and one shift. Substituting it into the expression for $a\bmod m$ gives the desired estimate $z'$ of the remainder.

Now analyze the error. The floor function has the property that for all $x > y > 0$, we have $\lfloor x\rfloor - \lfloor y\rfloor \le \lceil x - y\rceil$. Therefore, the error is

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, it is sufficient to subtract an extra $z' \ge m$ when $m$ to ensure the answer is correct.

The Barrett reduction computation uses only two multiplications, one shift, and at most two subtractions to perform an integer remainder. This efficiency improvement is not free: 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 bit length of the binary representation of integer $x$. Then:

Because the choice of $R$ must satisfy $a < R$, this length is at least $2\ell(a) - \ell(m)$. However, when a remainder is needed, we usually have $\ell(m)\le\ell(a)$, so this intermediate variable may be longer than the input length $\ell(a)$. For example, if a $64$-bit integer is reduced modulo a $32$-bit integer, the intermediate variable actually needs $64 \times 2 - 32 = 96$ bits.

One application of Barrett reduction is computing the remainder of the product $ab\bmod m$. If one multiplicand is fixed, for example when $b$ is fixed, we can use:

to make a similar estimate; it is only necessary to precompute the value of $\left\lfloor\dfrac{bR}{m}\right\rfloor$. This fixed-$b$ case is sometimes called Shoup modular multiplication⁶.

The more common case is that $a, b$ are both variable. First compute $ab$, then use Barrett reduction to obtain $ab\bmod m$. For example, when implementing modular multiplication, compute $0 \le a,b < m$ for $ab\bmod m$. In this case, the chosen $r$ must satisfy $ab < R$. By the analysis above, the longest intermediate variable has length $2\ell(ab)-\ell(m)$. When $\ell(a)\approx\ell(b)\approx\ell(m)$, this length is $3\ell(m)$. In other words, to implement modular multiplication for $32$-bit integers using Barrett reduction, the intermediate variable needs $96$ bits. This is one limitation of Barrett reduction in practical algorithm-contest applications.

As an example, a reference implementation of modular multiplication for signed 32-bit integers using Barrett reduction is:

// 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 uses 128-bit integers⁴.

Montgomery Modular Multiplication

Montgomery modular multiplication serves a purpose very similar to Barrett reduction: it also reduces the cost of remainder operations during modular arithmetic. Unlike the previous two algorithms, which approximate the quotient, Montgomery multiplication maps all integers into Montgomery space, where operations are relatively easy, thereby reducing the 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

Because $R\perp m$, there is a bijection between the congruence class $a \bmod m$ and its Montgomery form $aR\bmod m$. Therefore, integers can be converted to Montgomery form, several operations modulo $m$ can be performed, and the resulting Montgomery form can then be converted back to an integer; the result is always correct.

Montgomery form makes many modular operations convenient. As noted above, whether two congruence classes are equal can be checked by comparing their Montgomery forms. Also,

so addition and subtraction of congruence classes correspond to addition and subtraction of their Montgomery forms. However, multiplication of congruence classes cannot be performed by directly multiplying the two Montgomery forms. Indeed,

Thus, when multiplying two Montgomery forms, we need to apply the following Montgomery reduction operation to their product $x$:

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

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

Now discuss how to implement the Montgomery reduction operation $\operatorname{REDC}$. When computing $\operatorname{REDC}(x)$, we always assume $0 \le x < m^2$, which holds in all the cases above. Since $R\perp m$, Bézout’s identity gives integers $R^{-1},m'$ such that

Therefore, let $q=\lfloor xm' / R\rfloor$. Then

Since $0 \le x < m^2 < mR$ and $0 \le xm'\bmod R < R$,

That is, this quotient differs from $xR^{-1}\bmod m$ by at most one $m$. Adding $m$ when the quotient is negative gives $\operatorname{REDC}(x)$. Computing this quotient requires only two integer multiplications, one integer subtraction, and two bit operations (taking the remainder modulo and dividing by $R=2^k$). Thus, Montgomery reduction can be performed efficiently.

To perform Montgomery modular multiplication, a series of constants must be precomputed. First, Montgomery reduction uses $m' = m^{-1}\bmod R$, which can be computed with the Newton–Hensel method introduced below. In addition, converting different operations into Montgomery reductions involves constants such as $R^2\bmod m$. To obtain it, compute $R\bmod m$ once, then add it to itself to obtain $2R\bmod m$. Next, regard it as the Montgomery form of $2$ and perform fast exponentiation directly to obtain $2^kR\bmod m = R^2\bmod m$.

As an example, Montgomery modular multiplication for signed $32$-bit integers is implemented as follows:

// 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 with implementing modular multiplication using Barrett reduction, Montgomery multiplication involves several steps, including conversion, multiplication in Montgomery form, and inverse conversion. Therefore, conversion and inverse-conversion costs can be amortized only when enough modular operations are performed between them, leading to high overall efficiency. However, because the Montgomery implementation uses only intermediate variables of length $2\ell(m)$, it is more flexible. For example, multiplication of $32$-bit integers requires only an intermediate variable of $64$ bits. Thus, if a modular-integer class is needed for various number-theoretic computations, Montgomery multiplication is more suitable.

Integer Classes Modulo Powers of 2

This section discusses implementing a modular-integer class when the modulus is a power of $2$. In this special case, division and remainder operations can be implemented with bit operations, making computation very efficient. Barrett reduction and Montgomery multiplication both exploit this property when $2^e$ is used as divisor and modulus. In particular, when the modulus is exactly $2^{32}$ or $2^{64}$, an unsigned integer of the corresponding bit width combined with natural overflow can implement a modular-integer class without any explicit remainder operations. Even when the modulus is not exactly one of these values, it can be converted to one of these cases. For example, when the modulus is $2^{58}$, intermediate computations can be performed modulo $2^{64}$, and the result can finally be reduced modulo $2^{58}$. Besides the convenience of taking remainders, modular-integer classes modulo $2^e$ have many other special implementations. This section focuses on implementing inverses and exponentiation.

First, consider inversion: given odd $a$ and modulus $m=2^e~(e > 2)$, compute $a^{-1}\bmod m$. Common methods include the extended Euclidean algorithm and fast exponentiation. The extended Euclidean algorithm involves taking remainders modulo a general modulus; ordinary fast exponentiation requires computing $a^{\varphi(m)-1}\bmod{m}$, which takes $\Theta(e)$ integer multiplications. A more efficient method is the Newton–Hensel method. Specifically, apply the following fact:⁷

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

As an example, a reference implementation for inversion modulo $2^{32}$ is:

// 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, consider exponentiation: given $x,a,b$ and modulus $m=2^e~(e > 2)$, compute $xa^b\bmod m$, where $a$ is odd. From the analysis of the multiplicative structure modulo $2^e$, we know that $a$ can always be written in the form $\pm g^{\ell}$⁸, and the minus sign appears exactly when $a\equiv 3\pmod 4$. In this case, replace $a$ with $-a$ and multiply the final result by $(-1)^b$. Thus, assume from now on that $a\equiv 1\pmod 4$. The core idea of the algorithm is to write $a$ as $g^{L(a)}\bmod m$, then compute the desired power as $xg^{bL(a)}\bmod m$.

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

where $1 < e_1 < e_2 < \cdots < e_s < e$. This is because expanding the product directly shows that the binary representation of $a$ has $1$ as its second-lowest set bit, at position $e_1$ (with bit indices starting from $0$), so this representation can be found recursively. By the properties of discrete logarithms,

Because the modulus of the discrete logarithm is the order $\delta_m(g)=2^{e-2}=m/4$, multiply the entire congruence by $4$ here so that the computation can be carried out in the residue classes modulo $m$. Therefore, precomputing, for every $1 < d < e$, the quantity $4L(2^d+1)$ makes it possible to compute $4L(a)$ quickly.

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

Therefore, $\delta_m(2^d+1) = 2^{e-d}$. By the properties of orders,

Thus, $\gcd(\delta_m(g), \operatorname{ind}_g(2^d+1)) = 2^{d-2}$. This shows that $L(2^d+1) = \operatorname{ind}_g(2^d+1) = 2^{d-2}r$, where $2\nmid r$. Therefore, the least significant bit equal to $4L(2^d+1)$ in the binary representation of $1$ is exactly bit $d$ (with bit indices starting from $0$). Hence, the same binary-representation recursion can decompose $4L(a)$ into a sum of terms of the form $4L(2^d+1)$. This gives the value of $a$.

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

holds. Again, determine $1$ by finding the second-lowest bit equal to $e_1$, but to remove the factor $a^{-1}$ from $2^{e_1}+1$, it is enough to multiply $a$ by $2^{e_1}+1$, which can be done with bit operations. Since $4L(a^{-1})=-4L(a)$, subtraction must replace addition when accumulating $4L(a)$. Second, for a specially chosen base $g$, iteration need not continue to $d = e-1$; it only needs to reach $d = \lceil e/2\rceil - 1$. For this, choose $g$ such that

For $d \ge e / 2$, we have

Therefore, induction starting from $d = \lceil e/2\rceil$ shows that $L(2^d+1)=2^d$ for every $d \ge e/2$. Consequently, as long as $e/2 \le e_1 < e_2 < \cdots < e_s < e$, we have

and

Therefore, after processing all binary bits with $d < e/2$, the remaining discrete logarithm can be obtained directly without computing bit by bit. With the first optimization, the entire exponentiation requires $O(e)$ additions/subtractions and bit operations and $1$ multiplication; with the second optimization, about half of the additions/subtractions and bit operations can be omitted, at the cost of an additional $1$ multiplication.

As an example, exponentiation modulo $2^{32}$ is implemented as follows:

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

The discrete-logarithm preprocessing can be done with the Pohlig–Hellman algorithm. The base $g$ can be chosen as

References and Notes

• Fast Modular Multiplication by orz - Codeforces
• Barrett Reduction - Wikipedia
• Barrett Reduction - A41
• Barrett Reduction Principle and Correctness Proof by Chen - Zhihu
• Montgomery Multiplication - CP Algorithms
• Montgomery Modular Multiplication by Chen - Zhihu
• 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.