Modular inverse

The modular inverse of a number $a$ modulo $m$ is the number that undoes multiplying by $a$. It is how you divide when a problem asks for the answer modulo a prime such as 10^9 + 7. If the modulus is not prime, find the inverse with the extended Euclidean algorithm instead.

Explanation

Modulo 13, the number 24 leaves 11 and the number 6 leaves 6. 24 / 6 = 4, but the remainders alone don't give it: 11 / 6 isn't a whole number, so modular arithmetic has no division of its own. Look at 6 times 11 instead. It is 66, which leaves 1 modulo 13. So multiplying by 6 and then by 11 brings any remainder back to where it was. Instead of dividing 11 by 6, multiply by 11: 11 times 11 is 121, which leaves 4, and 4 is 24 / 6.

The number that undoes multiplying by 6 is its inverse. Modulo 13, every number from 1 to 12 has exactly one. 6 and 11 are each other's inverses, and so are 2 and 7, since 2 times 7 is 14, which leaves 1. 1 and 12 are their own.

Fig. 1. Each number from 1 to 12 has exactly one inverse modulo 13.

Modulo 10^9 + 7 a table like this would have a billion cells. What is the inverse of 2 modulo 10^9 + 7?

The natural way is to try 1, 2, 3 and so on, until 2 times the number leaves 1. Modulo 13 that takes at most 12 tries. Here the first number that works is 500000004: 2 times 500000004 is 1000000008, one more than the modulus. The search does find it, after 500000004 tries. By the rule of thumb of about 10^8 simple steps per second, that is about 5 seconds for a single division, and a counting problem may divide a million times.

So the question is how to reach the inverse without walking past every number before it. Binary exponentiation reaches any single power of a number in a few dozen multiplications. If the inverse of 2 were always a power of 2 at a known exponent, we could jump straight to it. Check it modulo 13: multiply by 2 again and again, keeping only the remainder, until it comes back to 1 (Fig. 2).

Start with an empty row. Cell k will hold 2 multiplied by itself k times, modulo 13.

The first power is 2.

2 · 2 = 4.

4 · 2 = 8.

8 · 2 = 16, which leaves 3.

3 · 2 = 6.

6 · 2 = 12.

12 · 2 = 24, which leaves 11.

11 · 2 = 22, which leaves 9.

9 · 2 = 18, which leaves 5.

5 · 2 = 10.

10 · 2 = 20, which leaves 7.

7 · 2 = 14, which leaves 1. The 12th power is back at 1, so 7, the 11th power, is the inverse of 2.

1/13

Fig. 2. Step through the powers of 2 modulo 13 with the buttons or the arrow keys.

The remainder comes back to 1 at the 12th power, so the 11th power, 7, is the inverse of 2: 2 times 7 is 14, which leaves 1. And 12 is one less than 13. If every number gets back to 1 at the power one less than the prime, the inverse is always the power two less than the prime. The powers of 3 modulo 13 are 3, 9, 1, 3, 9, 1, … They come back to 1 already after 3 steps, so they are at 1 again after 12, and the 11th power, 9, is the inverse of 3 (27 leaves 1).

This is Fermat's little theorem: for a prime $p$ and any number $a$ not divisible by $p$, $a$ to the power p - 1 leaves 1 modulo $p$.

For the inverse of 2 modulo 10^9 + 7 the power is 10^9 + 5. That number has 30 binary digits, 15 of them ones, so binary exponentiation needs 30 squarings and 15 more multiplications, 45 in all, against 500000004 tries for the search. The result is 500000004.

Computing the modular inverse a^{p-2} \bmod p takes O(\log p) multiplications. Modulo a prime $p$, every $a$ that is not a multiple of $p$ has an inverse.

When the modulus is not prime

Fermat's little theorem needs a prime modulus. Modulo 25, 2 to the power 23 leaves 8, and 2 · 8 = 16, not 1; the inverse of 2 is 13. For any modulus $m$, the extended Euclidean algorithm finds the inverse of $a$ whenever their greatest common divisor \gcd(a, m) = 1, in O(\log m) steps: the inverse of 7 modulo 25 is 18 (7 · 18 = 126). When \gcd(a, m) \ne 1, there is no inverse at all. For example, 6 has no inverse modulo 9: 6 times any number is a multiple of 3, and so is every multiple of 9, so the product never leaves 1.

Implementation

The program reads a, b and a prime MOD (the $p$ of Fermat's little theorem) and prints a / b modulo MOD. It uses the function power(x, n, m) from binary exponentiation, which returns x to the power n modulo m. There are three steps:

1. Reduce: inverse(b) takes b % MOD, moved into 0..MOD-1, because % in C++ keeps the sign of b. If it is 0, there is no inverse, and the function returns −1.

2. Raise: otherwise inverse(b) returns power(b, MOD - 2, MOD).

3. Multiply: main multiplies a by inverse(b) instead of dividing by b.

For 24 / 6 modulo 13 the input is:

24 6 13

#include <iostream>
using namespace std;

long long MOD;   // the prime modulus

// x to the power n modulo m, binary exponentiation
long long power(long long x, long long n, long long m) {
    long long result = 1 % m;
    x %= m;
    while (n > 0) {
        if (n % 2 == 1) result = result * x % m;
        x = x * x % m;
        n /= 2;
    }
    return result;
}

// the inverse of b modulo the prime MOD, or -1 if b is a multiple of MOD
long long inverse(long long b) {
    // reduce
    b = (b % MOD + MOD) % MOD;
    if (b == 0) return -1;
    // raise
    return power(b, MOD - 2, MOD);
}

int main() {
    long long a, b;
    cin >> a >> b >> MOD;
    long long ib = inverse(b);
    if (ib == -1) {
        cout << "no inverse\n";
        return 0;
    }
    // multiply: a / b modulo MOD is a times the inverse of b
    a = (a % MOD + MOD) % MOD;
    cout << a * ib % MOD << '\n';
}

The program makes O(\log \text{MOD}) multiplications and uses $O(1)$ memory.

Inverses of all numbers from 1 to n

When a problem divides by many different numbers up to $n$, calling inverse for each takes O(n \log \text{MOD}). A short loop builds all of them in $O(n)$, each from the inverse of a smaller number.

Take MOD = 13 and i = 5. Since 13 = 2 · 5 + 3 and 13 leaves 0, the numbers 3 and −2 · 5 leave the same remainder. Multiply both by inv[3] · inv[5]. On the left, 3 · inv[3] leaves 1, so what remains is inv[5]. On the right, 5 · inv[5] leaves 1, so what remains is −2 · inv[3] = −2 · 9. Adding 13 to −2 keeps the remainder, so inv[5] = 11 * 9 % 13, which is 99 % 13 = 8. Indeed 5 · 8 = 40 = 3 · 13 + 1.

In general, for each i from 2 to n:

1. Split: MOD = (MOD / i) * i + MOD % i.

2. Look up: inv[MOD % i] is already known, because MOD % i is less than i.

3. Combine: inv[i] = (MOD - MOD / i) * inv[MOD % i] % MOD.

Step 3 follows from step 1 as in the case i = 5. Write q = MOD / i and r = MOD % i. Then q * i + r leaves 0, so r and -q * i leave the same remainder. Multiplying both by inv[r] * inv[i] leaves inv[i] on the left and -q * inv[r] on the right, and adding MOD to -q gives (MOD - q) * inv[r]. Here r is never 0, because a prime MOD has no divisor between 2 and MOD - 1. Modulo 12, i = 2 gives r = 0, and inv[0] is not the inverse of anything.

inv[8] is then built from inv[5], since 13 = 1 · 8 + 5 (Fig. 3).

Fig. 3. The inverses of 1 to 12 modulo 13. Each arc points from i to 13 % i, the entry it is built from.

vector<long long> inv(n + 1);   // needs a prime MOD and n < MOD
inv[1] = 1;
for (int i = 2; i <= n; i++)
    inv[i] = (MOD - MOD / i) * inv[MOD % i] % MOD;   // split, look up, combine

Inverse factorials

Most counting problems divide by factorials, as in the binomial coefficient C(n, k) = \frac{n!}{k! \, (n - k)!}. Compute fact[i] for i up to n first. Then one call to power gives inv_fact[n], and the others follow downwards: dividing by (i - 1)! is the same as dividing by i! and multiplying by i. For MOD = 13 and n = 6, fact[6] leaves 5, so inv_fact[6] is 8, the inverse of 5. Then inv_fact[5] = 8 · 6 % 13 = 9, inv_fact[4] = 9 · 5 % 13 = 6, and so on down to inv_fact[0] = 1.

vector<long long> fact(n + 1), inv_fact(n + 1);   // needs a prime MOD and n < MOD
fact[0] = 1;
for (int i = 1; i <= n; i++)
    fact[i] = fact[i - 1] * i % MOD;
inv_fact[n] = power(fact[n], MOD - 2, MOD);
for (int i = n; i >= 1; i--)
    inv_fact[i - 1] = inv_fact[i] * i % MOD;

// C(n, k) for 0 <= k <= n
long long C(int n, int k) {
    return fact[n] * inv_fact[k] % MOD * inv_fact[n - k] % MOD;
}

The tables take O(n + \log \text{MOD}) time, and then each binomial coefficient takes $O(1)$.

Correctness proof

Take a prime $p$ and a number $a$ not divisible by $p$. Multiply each of 1, 2, \dots, p - 1 by $a$ and take the remainders modulo $p$. None of them is 0, because a prime that divides a product divides one of its factors, and $p$ divides neither $a$ nor a number below $p$. No two of them are equal: if a \cdot i and a \cdot j left the same remainder, $p$ would divide a \cdot (i - j), so it would divide i - j, which lies strictly between $-p$ and $p$ and is not 0. So the p - 1 remainders are the numbers 1, \dots, p - 1 in another order. For a = 2 and p = 13, Fig. 4 shows the twelve products.

Fig. 4. Multiplying 1 to 12 by 2 modulo 13 gives the same twelve numbers in another order.

Since the two lists hold the same numbers, their products leave the same remainder. The product of the second list is a^{p-1} \cdot (p - 1)!, the product of the first is (p - 1)!, so $p$ divides (p - 1)! \cdot (a^{p-1} - 1). It doesn't divide (p - 1)!, so it divides a^{p-1} - 1, and a^{p-1} leaves 1. For a = 2, p = 13: $12!$ leaves 12, and 2^{12} = 4096 leaves 1. Finally, a \cdot a^{p-2} = a^{p-1} leaves 1, so a^{p-2} is the inverse of $a$.

The same argument shows the inverse is unique: the remainder 1 appears exactly once in the second list.

Common mistakes

• Using power(b, MOD - 2, MOD) when MOD isn't prime. It returns a wrong number without any error, as with 2 modulo 25 above. Use the extended Euclidean algorithm instead.

• Inverting a number divisible by MOD. power(0, MOD - 2, MOD) returns 0 without complaint. With MOD = 7, C(8, 1) = 8 leaves 1, but fact[8] is 0, so fact[8] * inv_fact[1] % MOD * inv_fact[7] % MOD gives 0. The factorial tables work only for n < MOD, which always holds for 10^9 + 7 and n up to 10^6.

• Multiplying in int. Modulo 10^9 + 7, 500000004 · 500000004 = 250000004000000016 overflows int. Keep power, the tables and every product in long long. That is safe while the modulus is below about 3 \cdot 10^9, since (3 \cdot 10^9)^2 = 9 \cdot 10^{18} is still below the long long limit of about 9.2 \cdot 10^{18}.

Similar algorithms

The extended Euclidean algorithm is the fallback for a modulus that is not prime; it takes a few more lines of code.

When $n$ is at most a few thousand, or the modulus isn't prime, Pascal's triangle builds binomial coefficients with additions only and needs no division.

Euler's theorem extends Fermat's little theorem to any modulus: a^{\varphi(m) - 1} is the inverse of $a$ when \gcd(a, m) = 1. Computing $\phi (m)$ needs the prime factors of $m$, so the extended Euclidean algorithm is usually simpler.

Practice