Modular inverse and division

Why this matters when you build with AI

Many contest answers are fractions. A problem asks for the expected value of a random process, or the number of ways to choose items, and the result is a fraction p/q. The judge wants p times q^-1 modulo a prime, which means a whole number. An AI tool often writes ans = total / count % MOD and the program prints a number that is wrong. Integer division in C++ rounds toward zero, and % after a division does not undo it. The program runs, prints a value, and passes the small sample. It fails on the large tests.

We ran two small programs under the modulus 1000000007. The expression (14 / 2) % MOD printed 7, which is right because 2 divides 14. The expression (6 / 7) % MOD printed 0. The correct value is 6 times the inverse of 7. The inverse of 7 is 142857144, because 7 times 142857144 is 1000000008, which leaves remainder 1. So the correct answer is 857142864. Nothing crashed. The output looked like a plausible remainder, and the error only showed on inputs where the division was not exact.

This lesson covers the fix. Division under a modulus means multiplying by an inverse. The inverse of a number a modulo m is a number x such that a times x leaves remainder 1 when divided by m. We cover three ways to find it, a way to find many inverses at once, and a worked exercise that answers many division queries quickly.

The examples are C++17, compiled with g++ -std=c++17 -O2 -Wall -Wextra using GCC 15.2. The Python version at the end uses pow with a negative exponent, which needs Python 3.8 or later. The change that added this was CPython issue 36027. Results printed below come from our runs on the machine described in each section. Your output should match, but if a line depends on the machine, we say so.

The idea: division is multiplication by an inverse

Write a^-1 for the inverse of a modulo m. Then a / b modulo m is a times b^-1 modulo m. The rule uses the same idea as ordinary division: dividing by 4 is the same as multiplying by one quarter. The inverse replaces the fraction.

A small example with m = 13. We want 10 / 5 modulo 13. The inverse of 5 is 8, because 5 times 8 is 40, and 40 leaves remainder 1 when divided by 13 (39 plus 1). So 10 / 5 modulo 13 is 10 times 8, which is 80. Dividing 80 by 13 leaves 2, because 13 times 6 is 78. The answer is 2. Check it: 2 times 5 is 10, so 2 is the right quotient.

Not every number has an inverse. The inverse of a modulo m exists exactly when gcd(a, m) = 1, which means a and m share no prime factor. If the modulus is a prime p, every number from 1 to p - 1 has an inverse. If the modulus is 10, the number 4 has no inverse, because 4 times any x leaves remainder 0, 2, 4, 6 or 8, never 1. An AI tool often assumes that every nonzero number has an inverse. That holds for primes only.

Method 1: extended Euclid

The extended Euclidean algorithm finds integers x and y with a times x plus b times y equal to gcd(a, b). Lesson 1 covers the plain Euclidean algorithm for gcd. The extended version keeps track of the coefficients as it goes.

Take a = a0 and m as the modulus, with gcd(a0, m) = 1. The extended algorithm gives x and y with a0 times x plus m times y equal to 1. Reduce both sides modulo m. The term m times y leaves remainder 0, so a0 times x leaves remainder 1. The number x is the inverse, after we move it into the range 0 to m - 1.

Here is a hand trace for the inverse of 3 modulo 7. The call ext_gcd(3, 7) asks for x and y with 3x + 7y = 1. It calls ext_gcd(7, 3), which calls ext_gcd(3, 1), which calls ext_gcd(1, 0). The last call returns x = 1, y = 0, for 1 times 1 plus 0 times 0 equal to 1. Going back up, each level sets the new x to the old y, and the new y to the old x minus (a divided by b) times the old y.

ext_gcd(1, 0)  returns x = 1, y = 0
ext_gcd(3, 1)  x = 0, y = 1 - 3*0 = 1          check: 3*0 + 1*1 = 1
ext_gcd(7, 3)  x = 1, y = 0 - 2*1 = -2         check: 7*1 + 3*(-2) = 1
ext_gcd(3, 7)  x = -2, y = 1 - 0*(-2) = 1      check: 3*(-2) + 7*1 = 1

The x returned is -2. It is negative, so we move it into the range 0 to 6 by adding 7, which gives 5. Check: 3 times 5 is 15, and 15 leaves remainder 1 when divided by 7. So the inverse of 3 is 5.

The most common error is forgetting the normalising step. C++ keeps the sign of the dividend, so x % m can be negative, and the function must add m when it is. The code below does this in inv_ext.

The next block is the library for this lesson. Every function in it is used in the rest of the lesson. It is the code we compiled and ran.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
 
// Returns g = gcd(a, b) and sets x, y so that a*x + b*y = g. Needs a, b >= 0.
ll ext_gcd(ll a, ll b, ll& x, ll& y) {
    if (b == 0) { x = 1; y = 0; return a; }
    ll x1, y1;
    ll g = ext_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}
 
// Inverse of a modulo m, in [0, m). Returns -1 if gcd(a, m) != 1. Works for any m >= 1.
ll inv_ext(ll a, ll m) {
    ll x, y;
    ll a0 = ((a % m) + m) % m;
    ll g = ext_gcd(a0, m, x, y);
    if (g != 1) return -1;
    return ((x % m) + m) % m;
}
 
ll pow_mod(ll a, ll e, ll m) {
    ll r = 1 % m;
    a %= m;
    while (e > 0) {
        if (e & 1) r = r * a % m;   // m <= 1e9 here, so products fit in 64 bits
        a = a * a % m;
        e >>= 1;
    }
    return r;
}
 
// Inverse by Fermat's little theorem. Needs p prime and p does not divide a.
ll inv_fermat(ll a, ll p) { return pow_mod(a, p - 2, p); }
 
// inv[i] for 1 <= i <= n, modulo a prime p with n < p. O(n).
vector<ll> linear_inverses(int n, ll p) {
    vector<ll> inv(n + 1, 0);
    inv[1] = 1;
    for (int i = 2; i <= n; i++)
        inv[i] = (p - (p / i) * inv[p % i] % p) % p;
    return inv;
}
 
// Inverses of all values in v, modulo a prime p. Each v[i] must be in [1, p - 1].
// One pow_mod call plus about 3n multiplications.
vector<ll> batch_inverse(const vector<ll>& v, ll p) {
    int n = v.size();
    vector<ll> pre(n + 1, 1), out(n);
    for (int i = 0; i < n; i++) pre[i + 1] = pre[i] * v[i] % p;   // pre[i] = v[0]*...*v[i-1]
    ll cur = pow_mod(pre[n], p - 2, p);                           // inverse of the full product
    for (int i = n - 1; i >= 0; i--) {
        out[i] = cur * pre[i] % p;                               // inverse of v[i]
        cur = cur * v[i] % p;                                    // drop v[i] from the running inverse
    }
    return out;
}
 
int main() {
    ll p = 13;
    cout << "ext_gcd(240, 46): ";
    ll x, y;
    ll g = ext_gcd(240, 46, x, y);
    cout << "g=" << g << " x=" << x << " y=" << y << " check=" << 240 * x + 46 * y << "\n";
    cout << "inv_ext(3, 7) = " << inv_ext(3, 7) << "\n";
    cout << "inv_ext(2, 4) = " << inv_ext(2, 4) << "\n";
    cout << "inv_fermat(5, 13) = " << inv_fermat(5, p) << "\n";
    vector<ll> li = linear_inverses(12, p);
    cout << "linear 1..12 mod 13:";
    for (int i = 1; i <= 12; i++) cout << ' ' << li[i];
    cout << "\n";
    vector<ll> v = {2, 5, 7, 12};
    vector<ll> bi = batch_inverse(v, p);
    cout << "batch inverses mod 13:";
    for (ll t : bi) cout << ' ' << t;
    cout << "\n";
    return 0;
}

The real output of our run was:

ext_gcd(240, 46): g=2 x=-9 y=47 check=2
inv_ext(3, 7) = 5
inv_ext(2, 4) = -1
inv_fermat(5, 13) = 8
linear 1..12 mod 13: 1 7 9 10 8 11 2 5 3 4 6 12
batch inverses mod 13: 7 8 2 12

The line ext_gcd(240, 46) shows that the coefficients can be negative and large. The check column confirms that 240 times -9 plus 46 times 47 equals 2, the gcd. The line inv_ext(2, 4) returns -1, because 2 and 4 share the factor 2.

Method 2: Fermat's little theorem

Fermat's little theorem says that if p is a prime and a is not a multiple of p, then a^(p-1) leaves remainder 1 when divided by p. Dividing both sides by a gives a^(p-2) as the inverse. So for a prime p, the inverse is pow_mod(a, p - 2, p).

This method is short, and it costs one exponentiation, about 30 squarings and up to 30 more multiplications for p near 10^9. It is the method most contest code uses, and many contest problems use the prime 1000000007. The limit is clear. Fermat needs p to be prime. For a composite modulus the result can be wrong without any warning. We tested the modulus 10 and the number 3. The exponent is 10 - 2 = 8, and 3^8 leaves remainder 1 when divided by 10. But 3 times 1 is 3, not 1, so the value is not an inverse. Our run printed the following.

fermat 3^(10-2) mod 10 = 1, times 3 = 3

The same case appears in the common mistakes section.

Method 2 cannot handle a prime p that divides a. For example, 0 has no inverse, and 14 has none modulo 7. The function returns 0 for those values, and 0 times anything is 0, which hides the problem. Check the input before you call it.

Method 3: linear inverses for 1 to n

Sometimes you need the inverse of every number from 1 to n, for a prime p larger than n. One exponentiation per number costs n times about 30 multiplications, which is fast enough for n near 10^6. There is a faster method that takes one pass.

The idea is an identity. Write p as q times i plus r, where q is p divided by i, and r is p modulo i, so r is less than i. Reduce the equation q times i plus r = p modulo p. The right side leaves remainder 0, so q times i plus r is congruent to 0. Multiply both sides by the inverses of i and r:

q * i + r ≡ 0 (mod p)
i^-1 ≡ -q * r^-1 (mod p)

The number r is smaller than i, so its inverse was already computed when we reach i. One loop gives all the inverses. For example, with p = 13 and i = 9, q = 1 and r = 4. So the inverse of 9 is minus 1 times the inverse of 4. The inverse of 4 is 10, because 4 times 10 is 40, which leaves 1. So the inverse of 9 is minus 10, which is 3 modulo 13. Check: 9 times 3 is 27, which leaves 1. The table printed above lists the same value for 9.

This method needs p to be prime, and it needs n to be less than p. If n is at least p, one of the numbers is a multiple of p and has no inverse.

Method 4: batch inversion

Often you need the inverses of many values, and each one costs an exponentiation. Batch inversion gets all of them from a single exponentiation, plus about three multiplications for each value. The idea comes from a trick for working with fractions. Multiply the values together, invert the product once, then recover each inverse by multiplying with the other values.

A worked example with the modulus 7 and the values 2 and 3. The product is 6. The inverse of 6 is 6, because 36 leaves remainder 1. The inverse of 2 is the inverse of 6 times 3, which is 18, and 18 leaves remainder 4 when divided by 7. The inverse of 3 is the inverse of 6 times 2, which is 12, and 12 leaves remainder 5. Check: 2 times 4 is 8, which leaves 1, and 3 times 5 is 15, which leaves 1.

The code makes two passes. The first computes prefix products, where pre[i] is the product of the first i values. The second goes from the end. It keeps cur, the inverse of the product of the values that are still in front of it. At position i, cur times pre[i] is the inverse of v[i]. Then it removes v[i] from cur by multiplying by v[i], so cur is ready for the next position.

Batch inversion needs every value to be nonzero modulo p. One zero makes the product 0, and 0 has no inverse. The result is then 0 for every value, not only for the zero. We ran a check with the values 2, 0 and 3 modulo 13. The product was 0, and all three outputs came out as 0, while the correct inverses of 2 and 3 are 7 and 9. The symptom is that many answers become 0 at once. If your input can contain 0, handle it before the batch.

Worked example: queries of prefix sums divided by a value

The problem. You get an array A of N numbers, with 1 at most N at most 200000, and each A[i] in the range 0 to p - 1, where p = 1000000007. Then you get Q queries, with Q at most 200000. Each query has three integers l, r and v, with 1 at most l at most r at most N and 1 at most v at most p - 1. For each query print the sum of A[l] through A[r], modulo p, divided by v, modulo p, as a number in the range 0 to p - 1.

The approach. The sum of a range is a prefix sum difference, which you can get in O(1) after an O(N) pass. The division is multiplication by the inverse of v. The queries give Q values of v, which may all differ, so computing each inverse by Fermat costs Q exponentiations, roughly ten million modular multiplications at the largest sizes. Batch inversion computes all Q inverses with one exponentiation. The whole program is O(N + Q) plus one exponentiation.

The program. It uses the same two ideas as the library above. We wrote it with its own copies of pow_mod and batch_inverse, using the constant MOD, so that it can be read alone.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const ll MOD = 1000000007LL;
 
ll pow_mod(ll a, ll e, ll m) {
    ll r = 1 % m;
    a %= m;
    while (e > 0) {
        if (e & 1) r = r * a % m;
        a = a * a % m;
        e >>= 1;
    }
    return r;
}
 
vector<ll> batch_inverse(const vector<ll>& v, ll p) {
    int n = v.size();
    vector<ll> pre(n + 1, 1), out(n);
    for (int i = 0; i < n; i++) pre[i + 1] = pre[i] * v[i] % p;
    ll cur = pow_mod(pre[n], p - 2, p);
    for (int i = n - 1; i >= 0; i--) {
        out[i] = cur * pre[i] % p;
        cur = cur * v[i] % p;
    }
    return out;
}
 
int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int n, q;
    cin >> n >> q;
    vector<ll> pre(n + 1, 0);
    for (int i = 1; i <= n; i++) {
        ll a;
        cin >> a;
        pre[i] = (pre[i - 1] + a) % MOD;
    }
    vector<int> L(q), R(q);
    vector<ll> V(q), sums(q);
    for (int i = 0; i < q; i++) {
        cin >> L[i] >> R[i] >> V[i];
        sums[i] = ((pre[R[i]] - pre[L[i] - 1]) % MOD + MOD) % MOD;
    }
    vector<ll> invV = batch_inverse(V, MOD);   // one exponentiation for all q values of v
    for (int i = 0; i < q; i++) cout << sums[i] * invV[i] % MOD << '\n';
}

Here is the sample input and the real output of our run.

5 3
3 1 4 1 5
1 5 2
2 4 7
3 3 1000000006

7
857142864
1000000003

For the first query the sum is 3 + 1 + 4 + 1 + 5 = 14, and 14 divided by 2 is 7, which is exact. For the second the sum of positions 2 to 4 is 6, and 6 divided by 7 is the number 857142864, because 857142864 times 7 is 6000000048, which leaves remainder 6 when divided by 1000000007. For the third the value v is p - 1, which is -1 modulo p. The sum is 4, and 4 divided by -1 is -4, which is 1000000003 in the range 0 to p - 1.

How we checked it. We wrote a Python program with exact integers. It generated 300 random cases, with N up to 30, Q up to 30, values of v chosen from 1, p - 1, small numbers and random numbers, and it computed each answer as the sum times pow(v, -1, MOD). The C++ program matched on all 300 cases. We also checked the library functions. inv_ext was compared with a brute-force search for every modulus from 1 to 60 and every value from -130 to 130. ext_gcd was checked on 200000 random pairs, with the identity a times x plus b times y equal to the gcd. inv_fermat and linear_inverses were checked against products for several primes, and batch_inverse was checked on 2000 random arrays. All checks passed. The checks cover the inputs we tried. They do not prove the code correct for every input.

The same ideas in Python

Python integers never overflow, so the Python version needs no 128-bit care. From Python 3.8, pow(a, -1, m) returns the inverse directly and raises ValueError when no inverse exists. The program below uses that function as a check against the Euclid version, and it runs the batch method on random arrays. It is input-free, so the site runs it as written.

Python · editable

MOD = 1000000007

def inv_ext(a, m):
    # Python 3.8+ also has pow(a, -1, m); this shows the extended Euclid version.
    def eg(a, b):
        if b == 0:
            return a, 1, 0
        g, x, y = eg(b, a % b)
        return g, y, x - (a // b) * y
    g, x, _ = eg(a % m, m)
    if g != 1:
        return None
    return x % m

def batch_inverse(v, p):
    pre = [1] * (len(v) + 1)
    for i, t in enumerate(v):
        pre[i + 1] = pre[i] * t % p
    cur = pow(pre[-1], p - 2, p)
    out = [0] * len(v)
    for i in range(len(v) - 1, -1, -1):
        out[i] = cur * pre[i] % p
        cur = cur * v[i] % p
    return out

print("inv_ext(3, 7) =", inv_ext(3, 7))
print("inv_ext(2, 4) =", inv_ext(2, 4))
print("pow(3, -1, 7) =", pow(3, -1, 7))
vals = [2, 5, 7, 12]
print("batch inverses mod 13:", batch_inverse(vals, 13))
print("check:", [v * w % 13 for v, w in zip(vals, batch_inverse(vals, 13))])
print("batch inverses mod 1e9+7 of 3, 10^9, 123456789:", batch_inverse([3, 10**9, 123456789], MOD))
import random
random.seed(5)
ok = True
for _ in range(2000):
    arr = [random.randint(1, MOD - 1) for _ in range(random.randint(1, 40))]
    inv = batch_inverse(arr, MOD)
    ok &= all(a * b % MOD == 1 for a, b in zip(arr, inv))
    ok &= all(inv_ext(a, MOD) == pow(a, MOD - 2, MOD) for a in arr[:3])
print("random checks passed:", ok)

The real output of our run was:

inv_ext(3, 7) = 5
inv_ext(2, 4) = None
pow(3, -1, 7) = 5
batch inverses mod 13: [7, 8, 2, 12]
check: [1, 1, 1, 1]
batch inverses mod 1e9+7 of 3, 10^9, 123456789: [333333336, 857142863, 18633540]
random checks passed: True

The value for 10^9 is minus 7, whose inverse is minus the inverse of 7. That is 857142863 in the range 0 to p - 1.

Complexity

ext_gcd and inv_ext take O(log m) steps, the same as the Euclidean algorithm.
inv_fermat takes O(log p) multiplications, about 30 for p near 10^9.
linear_inverses(n, p) takes O(n) time and O(n) memory.
batch_inverse on n values takes O(n) multiplications and one exponentiation.
The query program takes O(N + Q) time, plus one exponentiation.

Common mistakes

Dividing with the integer division operator. The symptom is an answer that is correct on the sample and wrong on larger inputs, for example 0 printed for 6 / 7. Multiply by the inverse instead.
Forgetting the sign after extended Euclid. The symptom is a negative number in the output, such as -9 where the answer should lie between 0 and p - 1. Add m and take the remainder again, as inv_ext does.
Using Fermat with a composite modulus. The symptom is a value that is not an inverse. For modulus 10 and the number 3, pow_mod(3, 8, 10) returns 1, and 3 times 1 is 3, not 1. Use ext_gcd for any modulus.
Calling Fermat on a multiple of p. The symptom is the answer 0 for a value that should be an error. Check that v is not 0 modulo p before you invert it.
Passing a zero to batch inversion. The symptom is that many answers become 0 at once. Remove the zeros, or handle them separately.
Starting the linear table at index 0. The symptom is an out-of-range access or a zero for every entry, because 0 has no inverse. The loop starts at i = 2 with inv[1] = 1.

A prompt to check AI output

Give the prompt below to an AI tool with your inverse code. A good answer names the modulus it assumes, says whether the modulus must be prime, and gives an input that breaks the code.

Here is my C++ function that divides a by b modulo m, using pow(b, m - 2, m) for the
inverse. Do not rewrite it yet. Tell me whether this is correct for every m, and if not,
give one value of m and one value of b where it returns a number that is not the quotient.
Then say what the function returns for b = 0, and whether it prints a warning.

A correct answer says the function is correct only for a prime m, gives a composite example such as m = 10 and b = 3, and says that b = 0 returns 0 with no warning. If the AI says the code is correct without a condition, run the modulus 10 case yourself before you trust it.

Checklist

Before a division, check that the divisor is not a multiple of the modulus.
Use ext_gcd when the modulus might not be prime, and check that the gcd is 1.
Use Fermat only with a prime modulus.
Normalise every inverse into the range 0 to m - 1.
Batch the inverses when there are many divisors, and remove zeros first.
Test with a modulus of 1 and a divisor equal to the modulus minus 1.
Compare one batched answer with one answer from a single inverse on the same input.

Practice

Each practice problem has a hint. Try each one before you read the answer under the line "Show the answer".

Find the inverse of 7 modulo 26 by extended Euclid, by hand. Hint: run the Euclidean steps 26 = 3 times 7 + 5, then 7 = 1 times 5 + 2, then 5 = 2 times 2 + 1, and write 1 as a combination of 26 and 7.

Show the answer The inverse is 15, because 7 times 15 is 105, which is 4 times 26 plus 1. Working back: 1 = 5 - 2 times 2 = 5 - 2 times (7 - 5) = 3 times 5 - 2 times 7 = 3 times (26 - 3 times 7) - 2 times 7 = 3 times 26 - 11 times 7. So x = -11, and -11 leaves remainder 15 modulo 26. Python's pow(7, -1, 26) also returns 15.

Without a program, find the inverse of 2 modulo 1000000007. Hint: the answer is (p + 1) divided by 2, because 2 times that value is p + 1.

Show the answer The inverse is 500000004. Check: 2 times 500000004 is 1000000008, which is 1000000007 plus 1, so the remainder is 1. Python's pow(2, -1, 1000000007) gives the same value.

Use the linear table formula to find the inverse of 9 modulo 11. Hint: write 11 = 1 times 9 + 2, so q = 1 and r = 2, and look up the inverse of 2.

Show the answer The inverse of 2 modulo 11 is 6, because 12 leaves 1. By the formula, the inverse of 9 is minus 1 times 6, which is -6, and -6 leaves remainder 5 modulo 11. Check: 9 times 5 is 45, which is 4 times 11 plus 1.

You have the values 2 and 3 modulo 7. Compute both inverses with one exponentiation and the batch method. Hint: the product is 6, and 6 is its own inverse modulo 7.

Show the answer The inverse of the product 6 is 6. The inverse of 2 is 6 times 3, which is 18, and 18 leaves 4. The inverse of 3 is 6 times 2, which is 12, and 12 leaves 5. Check: 2 times 4 is 8, which leaves 1, and 3 times 5 is 15, which leaves 1.

Explain why inv_ext(4, 8) returns -1, and say what a program should do with that value. Hint: the gcd of 4 and 8 is 4, not 1.

Show the answer The number 4 shares the factor 4 with the modulus 8, so no x makes 4 times x leave remainder 1. The function returns -1 to signal that. The caller should treat it as an error, or check gcd(a, m) first. It must not use -1 as an inverse, because -1 modulo 8 is 7, and 4 times 7 is 28, which leaves 4.

Summary

Division under a modulus means multiplying by an inverse. Use extended Euclid for any modulus, and check that the gcd is 1. Use Fermat only for a prime modulus, and never for a divisor that is a multiple of the prime. For many inverses of one prime, use the linear table for 1 to n, or batch inversion for an arbitrary list, which needs one exponentiation in total. Test with a composite modulus, a zero divisor and a negative intermediate value, because those are the cases where the sample passes and the judge fails.