Skip to content

Modular Multiplicative Inverse

This article introduces the multiplicative inverse under modular arithmetic and discusses common methods for its computation.

Basic Concepts

The multiplicative inverse of a non-zero real number \(a\in\mathbf R\) is its reciprocal \(a^{-1}\). Similarly, in number theory, one can define the inverse of an integer \(a\) modulo \(m\), denoted as \(a^{-1}\bmod m\), or simply \(a^{-1}\). This is called the modular multiplicative inverse, also known as the number-theoretic reciprocal.

Inverse

For non-zero integers \(a,m\), if there exists \(b\) such that \(ab\equiv 1\pmod m\), then \(b\) is called the inverse of \(a\) modulo \(m\).

This is equivalent to saying that \(b\) is a solution to the linear congruence equation \(ax\equiv 1\pmod m\). According to the properties of linear congruence equations, the inverse \(a^{-1}\bmod m\) exists if and only if \(\gcd(a,m)=1\), i.e., when \(a\) and \(m\) are coprime, and it is unique modulo \(m\).

Computing a Single Inverse

Using the extended Euclidean algorithm or the fast exponentiation method, the inverse of a single integer can be computed in \(O(\log m)\) time.

Extended Euclidean Algorithm

Computing the inverse is equivalent to solving a linear congruence equation. Therefore, the extended Euclidean algorithm can be used to compute the inverse in \(O(\log\min\{a,m\})\) time. Additionally, since the linear equation corresponding to the inverse is special, the steps can be appropriately simplified.

Reference Implementation 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
// Extended Euclidean algorithm.
void ex_gcd(int a, int b, int& x, int& y) {
  if (!b) {
    x = 1;
    y = 0;
  } else {
    ex_gcd(b, a % b, y, x);
    y -= a / b * x;
  }
}

// Returns the modular inverse of a modulo m.
// Assumes that gcd(a, m) = 1, so the inverse exists.
int inverse(int a, int m) {
  int x, y;
  ex_gcd(a, m, x, y);
  return (x % m + m) % m;
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
# Extended Euclidean algorithm.
def ex_gcd(a, b):
    if b == 0:
        return 1, 0
    else:
        x1, y1 = ex_gcd(b, a % b)
        x = y1
        y = x1 - (a // b) * y1
        return x, y


# Returns the modular inverse of a modulo m.
# Assumes that gcd(a, m) = 1, so the inverse exists.
def inverse(a, m):
    x, y = ex_gcd(a, m)
    return (x % m + m) % m

This algorithm works for all cases where the inverse exists.

Fast Exponentiation Method

This method is primarily applicable when the modulus is a prime \(p\). In this case, by Fermat's little theorem, for any \(a\perp p\), we have

\[ a\cdot a^{p-2} = a^{p-1} \equiv 1 \pmod p. \]

Due to the uniqueness of the inverse, the inverse \(a^{-1}\bmod p\) equals \(a^{p-2}\bmod p\). Therefore, it can be directly computed using fast exponentiation in \(O(\log p)\) time:

Reference Implementation 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
// Binary exponentiation.
int pow(int a, int b, int m) {
  long long res = 1, po = a;
  for (; b; b >>= 1) {
    if (b & 1) res = res * po % m;
    po = po * po % m;
  }
  return res;
}

// Returns the modular inverse of a prime modulo p.
int inverse(int a, int p) { return pow(a, p - 2, p); }

1
2
3
4
# Returns the modular inverse of a prime modulo p.
# Use built-in pow function.
def inverse(a, p):
    return pow(a, p - 2, p)

In theory, this method can be generalized to arbitrary modulus \(m\) using Euler's theorem, i.e., computing the inverse as \(a^{\varphi(m)-1}\bmod m\). However, computing the Euler totient function \(\varphi(m)\) once is not trivial, so this algorithm is generally inefficient.

Computing Multiple Inverses

In some scenarios, it is necessary to quickly compute the inverses of multiple integers \(a_1,a_2,\cdots,a_n\) modulo \(m\). Computing them individually would require \(O(n\log m)\) time. However, if they are processed together, all inverses can be computed in \(O(n+\log m)\) time.

Consider the prefix products of the sequence \(\{a_i\}\):

\[ S_0 = 1,~ S_i = a_iS_{i-1},~ i=1,2,\cdots,n. \]

As long as each \(a_i\) is coprime with \(m\), their product \(S_n\) is also coprime with \(m\). Therefore, \(S_n^{-1}\bmod m\) can be computed using the methods described earlier. Since the inverse of a product is the product of inverses, starting from \(S_n^{-1}\) and traversing the sequence in reverse yields the inverses of each \(S_i\):

\[ S_{i-1}^{-1} = a_iS_i^{-1} \bmod m,~ i = n,n-1,\cdots,1. \]

Thus, the inverse of a single \(a_i\) can be computed as:

\[ a_i^{-1} = S_{i-1}S_i^{-1} \bmod m,~ i = 1,2,\cdots,n. \]

A reference implementation 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
// Returns the modular inverses for each x in a modulo m.
// Assume x mod m exists for each x in a.
std::vector<int> batch_inverse(const std::vector<int>& a, int m) {
  int n = a.size();
  std::vector<int> prod(n);
  long long s = 1;
  for (int i = 0; i < n; ++i) {
    // prod[i] = product of a[0...i-1]; prod[0] = 1.
    prod[i] = s;
    s = s * a[i] % m;
  }
  // s = product of all elements in a.
  s = inverse(s, m);
  std::vector<int> res(n);
  for (int i = n - 1; i >= 0; --i) {
    res[i] = s * prod[i] % m;
    s = s * a[i] % m;
  }
  return res;
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
# Returns the modular inverses for each x in a modulo m.
# Assume x mod m exists for each x in a.
def batch_inverse(a, m):
    n = len(a)
    prod = [0] * n
    s = 1
    for i in range(n):
        # prod[i] = product of a[0...i-1]; prod[0] = 1.
        prod[i] = s
        s = s * a[i] % m
    # s = product of all elements in a.
    s = inverse(s, m)
    res = [0] * n
    for i in reversed(range(n)):
        res[i] = s * prod[i] % m
        s = s * a[i] % m
    return res

In this algorithm, only one single inverse is computed, so the total time complexity is \(O(n+\log m)\).

Linear-Time Preprocessing of Inverses

If one needs to precompute the inverses of the first \(n\) positive integers modulo a prime \(p\), a recurrence relation discussed in this section can be used to compute them in \(O(n)\) time. This method is often used for preprocessing the reciprocals of factorials of the first \(n\) positive integers in combinatorial computations.

For a positive integer \(i\) with \(1< i < p\), consider the division with remainder:

\[ p = \left\lfloor \dfrac{p}{i} \right\rfloor i + (p\bmod i). \]

Taking this equation modulo the prime \(p\) yields

\[ 0 \equiv \left\lfloor \dfrac{p}{i} \right\rfloor i + (p\bmod i) \pmod p. \]

Multiplying both sides by \(i^{-1}(p\bmod i)^{-1}\) gives

\[ i^{-1} \equiv - \left\lfloor \dfrac{p}{i} \right\rfloor (p\bmod i)^{-1} \pmod p. \]

This is the formula for linear-time recurrence computation of inverses. Since \(p\bmod i < i\), this formula reduces the problem of computing \(i^{-1}\bmod p\) to the smaller problem of computing \((p\bmod i)^{-1}\bmod p\). Therefore, starting from \(1^{-1}\bmod p=1\) and applying this formula sequentially for each \(i\), the inverses of the first \(n\) integers can be obtained in \(O(n)\) time.

A reference implementation is as follows:

Reference Implementation 
1
2
3
4
5
6
7
8
9
// Precomputes modular inverses of all integers from 1 to n modulo prime p.
std::vector<int> precompute_inverses(int n, int p) {
  std::vector<int> res(n + 1);
  res[1] = 1;
  for (int i = 2; i <= n; ++i) {
    res[i] = (long long)(p - p / i) * res[p % i] % p;
  }
  return res;
}

1
2
3
4
5
6
7
# Precomputes modular inverses of all integers from 1 to n modulo prime p.
def precompute_inverses(n, p):
    res = [0] * (n + 1)
    res[1] = 1
    for i in range(2, n + 1):
        res[i] = (p - p // i) * res[p % i] % p
    return res

This algorithm only works when the modulus is prime. For a non-prime modulus \(m\), it cannot be guaranteed that \(m\bmod i\) in the recurrence formula remains coprime with \(m\), so the required \((m\bmod i)^{-1}\) for the recurrence might not exist. An example is \(m=8,i=3\). Here, \(m\bmod i = 2\), which has no inverse modulo \(m\).

Additionally, after obtaining this recurrence formula, a natural idea is to directly recursively compute the inverse of any number \(a\). Each recursion uses the recurrence formula to transform it into the inverse of the smaller remainder \(p\bmod a\), stopping when the remainder becomes \(1\). The complexity of this approach is currently unclear1, so it is recommended to use the conventional methods described earlier.

Exercises

References and Notes

  1. riteme's answer on Zhihu points out that the theoretical upper bound for the complexity of this approach is known to be \(O(p^{1/3+\varepsilon})\), while its performance on random data is close to \(O(\log p)\). ↩

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