Modular arithmetic and fast exponentiation

Why this matters when you build with AI

Contest and production code often asks for a result modulo a prime such as 1000000007. An AI tool will write (a * b) % m with long long for a and b near 10^18. The product does not fit in 64 bits, so the program prints a number that looks like a remainder but is not one. Sometimes it is negative. We ran a version that wraps the product in 64 bits, on 999999999999999989 * 999999999999999989 % 1000000007, and got -326471392. The correct remainder is 1444. The program compiles, runs fast and prints a value, so the bug is silent.

A second common mistake is a negative remainder. In C++, -7 % 3 is -1, because the sign of the remainder follows the left side. An AI tool often assumes the answer is always between 0 and m - 1, and the output then carries a negative value into the next step.

This lesson covers four things. Section by section, you will learn how to keep remainders in range, how to multiply two large numbers under a modulus without overflow, how to compute a power with about 60 multiplications even when the exponent is 10^18, and how to sum a geometric series without a division.

The examples are C++17, compiled with g++ -std=c++17 -O2 -Wall -Wextra using GCC 15.2. The __int128 type used here is a GCC and Clang extension, not part of the C++ standard. It is available on 64-bit targets. You can check for it with the macro __SIZEOF_INT128__, which is defined when the type exists. A Python version is given at the end, because Python integers never overflow.

The idea: remainders and congruence

Write a mod m for the remainder of a when divided by m. Two numbers are congruent modulo m when they leave the same remainder. Mathematicians write this as a ≡ b (mod m), and it means that a - b is a multiple of m.

The useful fact is that addition, subtraction and multiplication respect congruence. If a and c leave the same remainder, and b and d leave the same remainder, then a + b and c + d leave the same remainder, and so do a × b and c × d. So you may reduce at every step instead of at the end. The results stay small.

Division does not respect congruence in this way. We cover modular division in the modular inverse lesson.

Keeping remainders in range

The remainder that C++ returns has the sign of the dividend. For -7 % 3 the result is -1, not 2. Most contest problems want the value in the range 0 to m - 1, so we add m when the result is negative. This is the norm function below.

-7 % 3 = -1       (C++ result, negative)
norm(-7, 3) = 2   (the value we want, in 0 to m - 1)

The modular functions below call norm on their inputs, so no negative value can leak through.

Overflow in products

The product of two numbers below m needs about twice as many bits as m. If m is 10^9, the product is up to 10^18, which fits in a signed 64-bit integer, whose top value is about 9.2 × 10^18. If m is 10^18, the product is up to 10^36, which does not fit. The result is undefined behaviour for signed types and a wrong value for unsigned ones.

There are three safe ways to multiply.

128-bit intermediate. Cast the operands to __int128, multiply, and take the remainder. This is the shortest version and is correct for any m below 2^63.
Double and add. Split one operand into its binary digits, add shifted copies of the other, and reduce after each addition. It takes about 60 steps for 64-bit values, and it needs no special type.
long double estimate. Estimate the quotient with a floating point number, then correct the remainder with 64-bit wrapping arithmetic. It is fast, but it depends on the width of long double. On x86 GCC and Clang it is an 80-bit type with a 64-bit mantissa, which works for moduli up to 10^18 in our tests. On some compilers, such as MSVC, long double is the same as double, and this version gives wrong answers.

Binary exponentiation

To compute a^e mod m, you could multiply a by itself e - 1 times. That takes e steps, which is hopeless when e is 10^18. Binary exponentiation uses the binary digits of e. Each step squares the base, and when the current binary digit of e is 1, the result is multiplied by the current base.

A small example: 3^13. The exponent 13 is 1101 in binary, which is 8 + 4 + 1. So 3^13 = 3^1 × 3^4 × 3^8. Start with result 1 and base 3. Bit 1 is 1, so the result becomes 3. The base becomes 9. Bit 0 is 0, so nothing changes. The base becomes 81. Bit 1 is 1, so the result becomes 3 × 81 = 243. The base becomes 6561. Bit 1 is 1, so the result becomes 243 × 6561 = 1594323, which is 3^13. Each bit causes one squaring, so the number of multiplications is about twice the number of bits. For e up to 10^18 that is about 120 multiplications, and each one takes a fixed amount of work.

Here is the code for the library. The functions are in one block, so the same block can be reused in every program in this lesson.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using u64 = unsigned long long;
using i128 = __int128;  // GCC and Clang extension, not in the C++ standard
 
// Returns a value in [0, m) that is congruent to a. Works for negative a. Needs m >= 1.
ll norm(ll a, ll m) {
    ll r = a % m;  // C++ keeps the sign of a, so r can be negative
    return r < 0 ? r + m : r;
}
 
// Inputs are reduced first, so the sum is below 2m and fits when m <= 1e18.
ll add_mod(ll a, ll b, ll m) { return norm(norm(a, m) + norm(b, m), m); }
ll sub_mod(ll a, ll b, ll m) { return norm(norm(a, m) - norm(b, m), m); }
 
// (a * b) mod m with a 128-bit intermediate. Safe for any m below 2^63.
ll mul_mod(ll a, ll b, ll m) {
    return (ll)((i128)norm(a, m) * norm(b, m) % m);
}
 
// Portable version with no 128-bit type: double and add. Takes O(log b) steps.
ll mul_mod_slow(ll a, ll b, ll m) {
    a = norm(a, m);
    b = norm(b, m);
    ll r = 0;
    while (b > 0) {
        if (b & 1) r = (r + a) % m;  // r + a < 2m, fine for m <= 1e18
        a = (a + a) % m;
        b >>= 1;
    }
    return r;
}
 
// Portable version using long double (80-bit on x86 GCC and Clang). Valid for m <= 1e18.
ll mul_mod_ld(ll a, ll b, ll m) {
    a = norm(a, m);
    b = norm(b, m);
    ll q = (ll)((long double)a * b / m);             // estimate of floor(a*b/m)
    ll r = (ll)((u64)a * (u64)b - (u64)q * (u64)m);  // wraps mod 2^64; true value is small
    r %= m;
    if (r < 0) r += m;
    return r;
}
 
// a^e mod m by binary exponentiation. Needs e >= 0. Takes O(log e) multiplications.
ll pow_mod(ll a, ll e, ll m) {
    ll result = 1 % m;  // 1 % 1 is 0, so m = 1 gives 0
    a = norm(a, m);
    while (e > 0) {
        if (e & 1) result = mul_mod(result, a, m);
        a = mul_mod(a, a, m);
        e >>= 1;
    }
    return result;
}
 
// Returns {sum of r^i for 0 <= i < n, mod m}, {r^n mod m}. Needs n >= 0 and m >= 1.
// Doubling: a sum over 2k terms is the sum over k terms times (1 + r^k).
// No division is used, so this also works when m is not prime.
pair<ll, ll> geo_pair(ll r, ll n, ll m) {
    if (n == 0) return {0, 1 % m};
    if (n & 1) {
        auto [s, p] = geo_pair(r, n - 1, m);  // s covers n-1 terms, p = r^(n-1)
        return {add_mod(s, p, m), mul_mod(p, r, m)};
    }
    auto [s, p] = geo_pair(r, n / 2, m);  // s covers n/2 terms, p = r^(n/2)
    return {mul_mod(s, add_mod(1, p, m), m), mul_mod(p, p, m)};
}
 
ll geo_sum(ll r, ll n, ll m) { return geo_pair(r, n, m).first; }

The pow_mod function starts with 1 % m, not 1, so that the answer is 0 when m is 1. It also reduces the base once at the start. Inside the loop, every product goes through mul_mod, so no intermediate value grows beyond m.

Powers of powers and huge exponents

Exponents often come as expressions. The rules below are exact, so use them when they apply.

(a^b)^c = a^(b*c). If b * c fits in 64 bits, compute pow_mod(a, b * c, m) directly. Our exercise uses b and c up to 10^9, so the product is at most 10^18 and fits.
a^(b^c) needs the exact value of b^c as the exponent. If b^c is small enough to store, compute it exactly and use it. This is a checked power that stops once the value passes a limit.
If b^c is too large to store, the exponent cannot be computed directly. The fix is to reduce the exponent first, using a fact about remainders that lesson 8 proves. We do not use that fact here.

The capped power below returns the exact b^c if it is at most a limit, and -1 if not. It stops early once the value goes past the limit, so it never overflows. The loop runs at most about 60 times when b is at least 2, because b^60 is already far above 10^18.

// b^c if it is at most limit, otherwise -1. Needs b >= 0 and c >= 0.
ll capped_pow(ll b, ll c, ll limit) {
    if (c == 0) return 1;
    if (b == 0) return 0;
    if (b == 1) return 1;
    ll r = 1;
    for (ll i = 0; i < c; i++) {  // b >= 2, so this stops after about 60 steps
        if (r > limit / b) return -1;
        r *= b;
    }
    return r;
}

The check r > limit / b avoids the multiplication when it would pass the limit. Integer division rounds down, so r > limit / b holds exactly when r * b would be above limit.

Geometric series without division

A geometric series is the sum of powers of a fixed number r. The sum of r^i for i from 0 to n - 1 is (r^n - 1) / (r - 1) when r is not 1. The formula has a division, and under a modulus a division needs a modular inverse. Then it fails for three cases: r = 1, where the denominator is zero, a composite modulus where the inverse may not exist, and a modulus that divides r - 1.

The doubling method avoids the division. Let S(n) be the sum of the first n terms. Two facts give the recursion. For an even count, S(2k) = S(k) + r^k × S(k), because the second half of the terms is the first half multiplied by r^k. For an odd count, S(n) = S(n - 1) + r^(n-1). The function geo_pair returns both S(n) and r^n, so each step knows the power it needs. The recursion depth is about twice the number of bits of n.

For r = 1 the sum is simply n, and the doubling recursion computes that correctly, which is why it is the safer choice.

Worked example: geometric series queries

The problem. You get Q queries, with Q at most 100000. Each query has three integers r, n and m. The values satisfy |r| at most 10^18, 0 at most n at most 10^18, and 1 at most m at most 10^18. For each query print the sum of r^i for i from 0 to n - 1, taken modulo m, in the range 0 to m - 1.

The approach. Use geo_sum. It makes about 2 log2(n) recursive calls, and each call does a constant number of mul_mod operations.

The program. It is the library above, plus this main.

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int q;
    cin >> q;
    while (q--) {
        ll r, n, m;
        cin >> r >> n >> m;
        cout << geo_sum(r, n, m) << '\n';
    }
}

Here is the sample input we ran.

4
2 10 1000
-3 5 7
1 1000000000000000000 1000000007
10 0 1

The real output of our run was:

The first line is the sum 1 + 2 + ... + 512 = 1023, which is 23 modulo 1000. The second is 1 - 3 + 9 - 27 + 81 = 61, which is 5 modulo 7. This case has a negative ratio, and the answer is still in range. The third has r = 1, so the sum is just the number of terms, 10^18, taken modulo 1000000007, which gives 49. The fourth has m = 1, so the answer is 0 whatever the other values are.

How we checked it. The test program in C++ checked the norm output against exact 128-bit arithmetic for 200000 random values. It checked mul_mod, mul_mod_slow and mul_mod_ld against each other and against exact 128-bit arithmetic, for 300000 random triples with moduli up to 10^18. It checked pow_mod against repeated multiplication for 20000 random cases, and geo_sum against adding the terms one at a time, for 20000 random cases with r equal to 1 or -1 in some of them. All checks passed. In a separate Python script, we checked the three sample queries with Python's exact integers, and we checked 300 random prime-modulus cases with the closed formula (r^n - 1) × inverse of (r - 1). All of them matched.

These checks cover the code we ran. They do not prove it correct for every input. The long double version is only tested up to moduli of 10^18.

Exercise: patterns a^(b^c) and (a^b)^c

The problem. The first line has Q, the number of queries. Each query has four integers a, b, c and m, with 0 at most a, b, c at most 10^9 and 1 at most m at most 10^9. For each query print two numbers separated by a space. The first is (a^b)^c mod m. The second is a^(b^c) mod m if b^c is at most 10^18, and -1 if it is larger.

Your plan. For the first number, use the identity (a^b)^c = a^(b*c) and pow_mod. For the second, use capped_pow to decide whether b^c is small enough to store, then use the exact value as the exponent.

The program. It uses the library above, plus capped_pow and this main.

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    const ll LIMIT = 1000000000000000000LL;  // 1e18
    int q;
    cin >> q;
    while (q--) {
        ll a, b, c, m;
        cin >> a >> b >> c >> m;
        // (a^b)^c = a^(b*c). b and c are at most 1e9 here, so b*c fits in 64 bits.
        ll x = pow_mod(a, b * c, m);
        // a^(b^c): compute the exponent exactly when it is small enough.
        ll e = capped_pow(b, c, LIMIT);
        ll y = (e < 0) ? -1 : pow_mod(a, e, m);
        cout << x << ' ' << y << '\n';
    }
}

The sample input and real output of our run:

3
2 3 2 1000000007
2 3 4 1000
7 2 3 1

64 512
96 352
0 0

For the first query, (2^3)^2 = 64 and 3^2 = 9, so 2^9 = 512. For the second, (2^3)^4 = 4096, which is 96 modulo 1000, and 3^4 = 81, so 2^81 modulo 1000 is 352. For the third, a = 7 and m = 1, so both answers are 0.

We checked this program against Python with exact integers on 200 random inputs with a up to 30, b and c up to 6, and m up to 10^6. All answers matched. For the second number, Python computed the exponent exactly and then reduced, so the check did not rely on the same code.

Complexity

norm, add_mod and sub_mod take constant time.
mul_mod takes constant time with the 128-bit type. mul_mod_slow takes about log2(b) steps, which is about 60 for 64-bit values.
pow_mod(a, e, m) takes O(log e) calls to mul_mod, about 120 for e near 10^18.
geo_sum(r, n, m) takes O(log n) recursive calls, each with a constant number of mul_mod calls.
capped_pow(b, c, limit) takes at most about 60 steps when b is at least 2.

Common mistakes

A negative remainder in the output. The symptom is an answer such as -1 for a value that should be 2. Use norm on every result that might be negative. For -7 % 3, the remainder is -1.
Overflow in a * b % m. The symptom is a value that looks random, such as -326471392 where 1444 is correct. Use mul_mod. The naive product of two numbers near 10^18 does not fit in 64 bits.
Starting the power at 1 when m is 1. The symptom is an answer of 1 where 0 is correct. Start at 1 % m.
Using the geometric formula with a division. The symptom is a division by zero when r is 1, or a wrong answer for a composite modulus. Use the doubling method in geo_pair.
Passing a negative exponent to pow_mod. The loop does not run, and the function returns 1 without any error. Check the sign of the exponent in the caller.
Treating 0^0 as undefined. The convention used here, and in most libraries, is that 0^0 = 1. pow_mod(0, 0, m) returns 1 % m. Decide what the problem means and handle it explicitly if the problem defines it.
Computing b * c without checking it. For (a^b)^c = a^(b*c), the product must fit in 64 bits. Check the limits of the problem before you rely on it.

A prompt to check AI output

Give the prompt below to an AI tool with your modular code. A good answer shows the widest value that enters each multiplication and states whether it fits.

Here is my C++ function that computes a * b mod m for long long values, where m can be
up to 1e18. Do not rewrite it yet. Tell me the largest value the product can reach,
say whether it fits in a signed 64-bit integer, and give one input pair that breaks it.
Also say what the function returns for a = -7, b = 3, m = 3.

A correct answer names 10^36 as the largest product, says it does not fit, and gives an input such as two values near 10^18. It also says the answer for -7 and 3 modulo 3 is 0, because -7 times 3 is -21, which is a multiple of 3. If the AI skips the negative case, test it yourself with norm before you trust the function.

Checklist

Every value is passed through norm before it is used in a comparison or printed.
Every product of two values below m goes through mul_mod, not through * and %.
Powers start at 1 % m, so that m = 1 gives 0.
Exponents are checked for a negative sign before pow_mod is called.
Geometric sums use the doubling method, not a division by r - 1.
Products such as b × c in an exponent are checked against the 64-bit limit.
Tests include m = 1, a negative base, a zero exponent and a value near 10^18.

Practice

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

Compute 3^200 modulo 1000 by hand, using binary exponentiation with the bits of 200. Hint: 200 is 11001000 in binary.

Show the answer 3^200 modulo 1000 is 1. There are 7 squarings, one for each bit after the first, and the results are reduced modulo 1000 at each step. We checked it with Python's three-argument pow, which gave the same value.

Compute the sum of 2^i for i from 0 to 9, modulo 7, without adding all ten terms. Hint: use the formula (r^n - 1) / (r - 1) with an inverse, or add the terms in a loop and compare.

Show the answer The exact sum is 1023. The remainder of 1023 divided by 7 is 1, because 7 × 146 = 1022. A geometric formula gives the same result.

Find 10^18 modulo 7. Then find the sum of 1 for n terms modulo 7, with n = 10^18. Hint: a geometric series with r = 1 is a count.

Show the answer 10^18 modulo 7 is 1, and the sum of 10^18 ones modulo 7 is also 1. We checked both with Python's exact integers. The sum equals n, so the second value is n mod 7.

Compute (2^5)^3 modulo 1000 in two ways: as nested powers, and as a single power of 2 with the product of exponents. Hint: the two ways must agree.

Show the answer The nested value is (2^5)^3 = 32^3 = 32768, which is 768 modulo 1000. The single power is 2^(5 × 3) = 2^15 = 32768, which is also 768 modulo 1000.

Write norm(-5, 3) and give the value C++ returns for -5 % 3. Hint: the sign of a remainder in C++ follows the dividend.

Show the answer In C++, -5 % 3 is -2. The function norm(-5, 3) returns 1, because -5 modulo 3 is 1 in the range 0 to 2, and 1 is -2 + 3.

Summary

Keep every remainder in the range 0 to m - 1 with a single normalising function. Multiply under a modulus with a 128-bit intermediate, or with a portable method whose limits you have tested. Binary exponentiation computes huge powers in about 120 multiplications, and the doubling method sums a geometric series without a division. Check every product against the 64-bit limit before you trust it.