Why this matters when you build with AI
Ask an AI tool to compute 2^(10^18) mod 1000000007 and it may write a loop that multiplies 10^18 times, which never finishes. Ask for a shortcut and it may write pow(a, e % (m - 1), m). That line is right when m is a prime and wrong for many other moduli, and nothing in the code tells you which case you are in. The wrong answer looks like a normal number.
The second common problem is a silent edge case. A base that is a multiple of the modulus, such as 7 modulo 7, breaks the usual shortcut, and an AI tool will often not mention it.
Two terms first. Two numbers are coprime when their greatest common divisor is 1. The multiplicative order of a modulo m is the smallest k at least 1 such that a^k leaves remainder 1 when divided by m. A primitive root modulo a prime p is a number g whose order is p - 1, which means its powers hit every non-zero remainder before they repeat.
This lesson uses the totient function from the previous lesson, which counts the numbers from 1 to m that are coprime with m. Write it as phi(m). The binary exponentiation and mulmod functions come from the modular arithmetic lesson, and the code below uses a short copy of them so that it compiles on its own.
The idea: Fermat and Euler
Fermat's little theorem
It says: if p is a prime and a is not a multiple of p, then a^(p - 1) leaves remainder 1 when divided by p.
For a small example, take p = 7 and a = 2. The powers of 2 modulo 7 are 2, 4, 1, 2, 4, 1, and so on. The sixth power is 64 = 9 * 7 + 1, so the remainder is 1, as the theorem says. The powers cycle with period 3, which is a divisor of 6.
Euler's theorem
Euler's theorem is the same result for any modulus m. If a and m are coprime, then a^phi(m) leaves remainder 1 when divided by m. For a prime p, phi(p) = p - 1, so Euler's theorem gives Fermat's theorem as a special case.
Here is a short argument for why it holds. List the residues from 1 to m that are coprime with m. There are phi(m) of them. Multiplying each one by a, where a is coprime with m, produces the same set of residues in a different order, because no two of them can become equal. So the product of the list, call it P, is the same before and after multiplying by a. That gives a^phi(m) * P = P modulo m. P is coprime with m, so it can be cancelled, and the result is a^phi(m) = 1 modulo m.
Worked example: reducing an exponent
Take a = 3, m = 10 and e = 1000. Euler's theorem gives phi(10) = 4, so 3^4 = 81 = 1 modulo 10. The exponent can be reduced: 1000 = 250 * 4 + 0, so 3^1000 = (3^4)^250 = 1 modulo 10. Now take e = 1002. Then 1002 mod 4 = 2, so 3^1002 = 3^2 = 9 modulo 10. The rule is
a^e = a^(e mod phi(m)) modulo m, when a and m are coprime and e is at least 0.
This works because e = q * phi(m) + r with 0 <= r < phi(m), and a^e = (a^phi(m))^q * a^r = a^r modulo m. Reducing the exponent turns a huge number into a small one.
When the shortcut breaks
The rule above needs a and m to be coprime. If they share a factor, it can give the wrong answer.
Take m = 4 and a = 2, which share the factor 2. Here phi(4) = 2. Take e = 2. The direct answer is 2^2 = 4, which is 0 modulo 4. The reduced exponent is 2 mod 2 = 0, which gives 2^0 = 1 modulo 4. The two answers differ.
Take a multiple of the modulus as the base. For 7^6 modulo 7, the answer is 0, since 7 divides 7^6. The Fermat reduction gives 6 mod 6 = 0, which then gives 7^0 = 1. Again the answers differ.
A correct version for any base is this. Let K be the largest exponent of a prime in m, which is at most 62 when m fits in 63 bits. For every a and every e at least K,
a^e = a^(e mod phi(m) + phi(m)) modulo m.
Here is a proof sketch you can follow with pencil and paper. Look at one prime power p^k that divides m, with k at most K. If p does not divide a, the first formula works for p^k, because a^phi(p^k) = 1 modulo p^k and phi(p^k) divides phi(m). Both exponents are congruent modulo phi(p^k), so they give the same remainder modulo p^k. If p divides a, then a^e is a multiple of p^e, and both exponents are at least K, so both sides are 0 modulo p^k. Each prime power agrees, and the Chinese remainder theorem combines them into agreement modulo m.
The code below uses the threshold 64 for any modulus up to 2^63, because 64 is at least K for every such modulus. The extra phi(m) keeps the exponent large, so the shortcut is valid even when a shares a factor with m.
Huge exponents: a^(b^c) modulo m
A common contest pattern is a power whose exponent is itself a power. The exact exponent b^c can have billions of digits, so it cannot be stored. The reduction still works, because only b^c modulo phi(m) is needed. The one extra fact is whether b^c is at least 64, so the reduced exponent with phi(m) added can be used. Compute that test by multiplying b repeatedly and stopping once the value passes 64.
Binary exponentiation and mulmod
Fast power multiplies the base by itself and halves the exponent each step, so a^e costs about log2(e) multiplications. Products of two numbers below m can exceed 64 bits when m is large. The code uses __int128 for the product, which holds up to about 1.7 * 10^38. The modular arithmetic lesson explains this in more detail.
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
// (a * b) mod m for m up to about 1e18, without overflow.
ll mulmod(ll a, ll b, ll m) {
return (ll)((__int128)a * b % m);
}
// a^e mod m for e >= 0, by binary exponentiation.
ll powmod(ll a, ll e, ll m) {
ll result = 1 % m;
a %= m;
if (a < 0) a += m;
while (e > 0) {
if (e & 1) result = mulmod(result, a, m);
a = mulmod(a, a, m);
e >>= 1;
}
return result;
}Multiplicative order
The order of a modulo m, for a coprime with m, is the smallest k at least 1 with a^k = 1 modulo m. Euler's theorem says the order divides phi(m), because a^phi(m) = 1. So start with k = phi(m) and try to shrink it. For each distinct prime q of k, keep dividing k by q as long as a^(k / q) is still 1 modulo m. When no prime q can be removed, k is the order.
Take a = 10 and m = 13. Here phi(13) = 12. The prime factors of 12 are 2 and 3. For q = 2: 10^6 = 1 modulo 13, so k becomes 6. Then 10^3 = 12 modulo 13, which is not 1, so the loop stops for q = 2. For q = 3: 10^2 = 9 modulo 13, which is not 1, so k stays 6. The order is 6.
Each prime factor of phi(m) costs a few powers, and the factoring here uses trial division.
// The distinct prime factors of n, by trial division.
vector<ll> distinct_primes(ll n) {
vector<ll> primes;
for (ll p = 2; p * p <= n; ++p) {
if (n % p == 0) {
primes.push_back(p);
while (n % p == 0) n /= p;
}
}
if (n > 1) primes.push_back(n);
return primes;
}
// phi(n) by trial division, O(sqrt n).
ll phi_one(ll n) {
ll result = n;
for (ll p = 2; p * p <= n; ++p) {
if (n % p == 0) {
while (n % p == 0) n /= p;
result = result / p * (p - 1);
}
}
if (n > 1) result = result / n * (n - 1);
return result;
}
// The multiplicative order of a modulo m: the smallest k >= 1 with a^k = 1 (mod m).
// Needs gcd(a, m) = 1. The order divides phi(m), so start there and remove factors.
ll order_mod(ll a, ll m) {
ll k = phi_one(m);
for (ll q : distinct_primes(k)) {
while (k % q == 0 && powmod(a, k / q, m) == 1) k /= q;
}
return k;
}Primitive roots of a prime
A primitive root modulo a prime p is a number g with order p - 1. Its powers g^1, g^2, ..., g^(p - 1) are all the non-zero residues in some order. Such a g exists for every prime p, and there are phi(p - 1) of them. For p = 13 there are phi(12) = 4 of them, namely 2, 6, 7 and 11.
To test whether g is a primitive root, check that g^((p - 1) / q) is not 1 for every distinct prime q of p - 1. If some g^((p - 1) / q) were 1, the order would divide (p - 1) / q, which is smaller than p - 1. The code searches g = 2, 3, and so on, and returns the first one that passes.
// The smallest primitive root of a prime p >= 3.
ll primitive_root(ll p) {
vector<ll> qs = distinct_primes(p - 1);
for (ll g = 2; g < p; ++g) {
bool ok = true;
for (ll q : qs)
if (powmod(g, (p - 1) / q, p) == 1) { ok = false; break; }
if (ok) return g;
}
return -1;
}Take p = 7. Then p - 1 = 6 has the prime factors 2 and 3. For g = 2: 2^3 = 1 modulo 7, so 2 fails. For g = 3: 3^3 = 27 = 6 modulo 7, which is not 1, and 3^2 = 2 modulo 7, which is not 1. So 3 is a primitive root, and the search stops at 3.
Python version
Python has pow(a, e, m) built in, which does binary exponentiation for any size of integer, so the power functions need no code. The reduced-exponent rule is short enough to write in Python too. It uses phi_one from the lesson on totients, which is short and is shown there.
def pow_exact_exponent(a, e, m):
# a^e mod m for an exact integer e
if e < 64:
return pow(a, e, m)
phi = phi_one(m) # from the totient lesson
return pow(a, e % phi + phi, m)The Python version is fine for a few queries. For a batch of a million orders, use C++.
Worked example: order and primitive root
The problem below was written for this lesson.
Input: one integer t, then t lines. Each line holds an odd prime p with 3 <= p <= 10^9 and an integer a with 1 <= a < p. For each line print two numbers: the multiplicative order of a modulo p, then the smallest primitive root of p.
The solution uses order_mod and primitive_root, which were shown above.
// Input: t, then t lines "p a". p is an odd prime up to 1e9, 1 <= a < p.
// Output per line: the order of a modulo p, then the smallest primitive root of p.
int main() {
int t;
cin >> t;
while (t--) {
ll p, a;
cin >> p >> a;
cout << order_mod(a, p) << " " << primitive_root(p) << "\n";
}
}The program was compiled with g++ -O2 -std=c++17 -Wall -Wextra using g++ 15.2.0 on Linux. Here is its output on a three-line input.
3
7 2
13 10
998244353 33 3
6 2
998244352 3The order of 2 modulo 7 is 3, and 3 is the smallest primitive root of 7. The order of 10 modulo 13 is 6, and 2 is the smallest primitive root of 13. The order of 3 modulo 998244353 is 998244352, which is p - 1, so 3 is a primitive root of that prime. This matches the value of p - 1 and the earlier computation of the smallest primitive root of 998244353, which is also 3.
The program was tested on 2141 queries with primes up to 3000, and every order and primitive root matched a plain loop that tries every power. It was also tested on 400 random primes between 10^8 and 10^9, where the order and root were checked with Python's pow instead of a loop. All checks passed. The same tests also covered the power functions: 3000 cases of a^(b^c) modulo m, with b up to 10^9 and c up to 5, were compared with Python's pow(a, b ** c, m), which takes the exact huge exponent, and all matched. The reduction was also compared with direct multiplication for every base up to 210, every modulus up to 200 and every exponent up to 400. All matched.
Hands-on exercise
Work through these steps in a scratch folder on your own machine.
First, write order_mod by the plain loop: start at a and multiply by a until you reach 1, and count the steps. This is slow but obviously correct. Use it for every prime below 3000.
Second, write the fast version from this lesson. Compare the two on every a from 1 to p - 1 for each prime p up to 300. Then compare them on 200 random pairs with p up to 10^6. Keep the random seed fixed.
Third, test the reduction. Pick a = 2 and m = 15. Compute 2^20 modulo 15 directly, and then with the reduced exponent 20 mod phi(15) = 20 mod 8 = 4, and then with 20 mod 14, which uses the wrong modulus. The first two agree and the third does not.
Fourth, time 1000 random primes near 10^9 with primitive_root and order_mod. Record the number, which will depend on your machine.
A good answer gives the same order and root as the plain loop on every test, and matches 3 and 2 for the examples in this lesson. Check the order by two facts: the power a^order must be 1 modulo p, and for each prime q of the order, the power a^(order / q) must not be 1.
Prompt to try with an AI tool
Use this prompt, then check the answer with the steps above:
Write a C++17 function that computes a^(b^c) mod m for a, b, c, m up to 10^9, with
c up to 10^18. Use Euler's theorem to reduce the exponent and handle a that shares a
factor with m. Explain in two sentences why the reduction needs the phi(m) added back.
Then give the output for a = 2, b = 3, c = 2, m = 15.A good answer computes phi(m), checks whether b^c is at least 64, and uses the reduced exponent with phi(m) added. For a = 2, b = 3, c = 2, m = 15 the exponent is 3^2 = 9, so the answer is 2^9 modulo 15, which is 2. Check the code by comparing with Python's pow(2, 3 ** 2, 15), which prints 2, and check that it does not call pow with m - 1.
Complexity
Fast power takes O(log e) multiplications. The order takes O(log m) work per prime factor of phi(m) plus the cost of factoring phi(m) by trial division, which is O(sqrt m). The primitive root search takes O(sqrt p) for the factoring and a few powers per candidate g, and a small candidate g is found quickly in practice. The power a^(b^c) takes O(sqrt m + log phi(m)) work, mostly for the totient.
Common mistakes
These mistakes are listed with the output each one gave in the tests for this lesson.
The first mistake is reducing the exponent modulo m - 1 for a composite m. For m = 15 the value phi(m) is 8, not 14. The line powmod(2, 20 % 14, 15) printed 4, but the correct value of 2^20 modulo 15 is 1, which powmod(2, 20, 15) printed. The symptom is a wrong answer for some exponents and right answers for small ones, so it passes a quick test that uses small exponents only.
The second mistake is reducing when the base is a multiple of the modulus. For 7^6 modulo 7, the reduced line powmod(7, 6 % 6, 7) printed 1, but the direct value is 0. The same happens for 2^2 modulo 4: the reduced exponent printed 1, the direct value is 0. The fix is the version with phi(m) added back, which printed 0 in both cases.
The third mistake is computing the order of a number that is not coprime with m. For a = 2 and m = 4, order_mod(2, 4) printed 2, and the check powmod(2, 2, 4) printed 0, so 2 is not an order. The function does not check coprimality. Add a call to __gcd(a, m) == 1 before you use it.
The fourth mistake is overflow in a product of two residues. Signed 64-bit integers hold values up to 9223372036854775807, and the product of two residues below m fits only when m is at most 3037000500. Use __int128 in mulmod, as the code above does, for any larger modulus.
The fifth mistake is applying Fermat's theorem to a composite modulus. Fermat needs a prime. For m = 15 and a = 2, the reduction with m - 1 = 14 gives the wrong answer shown above, while the reduction with phi(15) = 8 is correct.
A prompt to check AI output
When an AI tool gives you a power function, ask it to show its own checks:
For the power function you wrote, give a^e mod m for these inputs and say which rule you
used for each: (2, 20, 15), (7, 6, 7), (2, 2, 4), (3, 1000, 10). Then list which of these
inputs the rule with m - 1 would get wrong, and why.A good answer gives 1, 0, 0 and 1 for the four inputs, names the phi-based rule for each, and says that the m - 1 rule is wrong for (2, 20, 15) and (7, 6, 7). Check each value with Python's pow, which is exact for any integers.
Checklist
- Use Fermat's theorem only for a prime modulus and a base that is not a multiple of it.
- Use Euler's theorem with phi(m) for any modulus, and only when the base is coprime with m.
- For a base that may share a factor with m, reduce with phi(m) added back, once the exponent is at least 64.
- Compute the order from phi(m), and check that the base is coprime with m first.
- A primitive root of a prime p has order p - 1. Test g against each prime factor of p - 1.
- Use
__int128for products when the modulus can be above about 3 * 10^9. - Compare every power and order function with Python's
powon random inputs before you trust it.
Practice
These take a few minutes each. Try them on paper or in your scratch folder before opening the answers.
Practice 1: the order of 10 modulo 13
Find the order of 10 modulo 13 by hand, without a computer.
Hint: Compute the powers 10, 10^2 and 10^3 modulo 13, and look for the first that gives 1 or -1.
Show the answer
10^1 = 10, 10^2 = 100 = 9, and 10^3 = 90 = 12 = -1 modulo 13. Then 10^6 = (-1)^2 = 1, so the order is 6, which is a divisor of phi(13) = 12 as expected.
Practice 2: all primitive roots of 13
List every primitive root of 13.
Hint: There are phi(12) = 4 of them. Test the numbers 2 through 12 with the factors 2 and 3 of 12.
Show the answer
The primitive roots of 13 are 2, 6, 7 and 11. The count 4 matches phi(12). The program printed the same four numbers in the test for this lesson.
Practice 3: the value of 7^(3^100) modulo 1000
Compute 7^(3^100) modulo 1000 without a loop that multiplies 3^100 times.
Hint: phi(1000) = 400, and 3^100 is at least 64. Find 3^100 modulo 400 first.
Show the answer
Since 3^100 is far larger than 64, the rule gives 7^(3^100 mod 400 + 400) modulo 1000. The program prints 7 for this value, and Python's pow(7, 3 ** 100, 1000) also gives 7.
Practice 4: why m - 1 is wrong for m = 15
Explain in plain words why the rule a^e = a^(e mod (m - 1)) gives the wrong answer for a = 2, m = 15 and e = 20.
Hint: For a composite m, the correct period of powers comes from phi(m), not from m - 1.
Show the answer
The value of phi(15) is 8, and 14 is not the period of powers of 2 modulo 15. The correct value 2^20 modulo 15 is 1, because 2^4 = 16 = 1 modulo 15. The rule with 14 reduces 20 to 6 and gives 2^6 = 64 = 4 modulo 15, which is wrong. The rule with phi(15) = 8 reduces 20 to 4, and 2^4 = 1, which is right.
Practice 5: the order of 2 modulo 15
Find the order of 2 modulo 15, and say why it divides phi(15).
Hint: Compute 2, 4, 8 and 16 modulo 15.
Show the answer
The powers of 2 modulo 15 are 2, 4, 8 and 1, so the order is 4. Euler's theorem gives 2^8 = 1 modulo 15, and the order must divide 8. Since 4 divides 8 and is the first power that gives 1, the order is 4. The code order_mod(2, 15) prints 4 in the test.
Summary
Fermat's theorem and Euler's theorem let you replace a huge exponent with a small one, but only when the base shares no factor with the modulus, or when phi(m) is added back for a large exponent. The order of a number divides phi(m) and can be found by removing prime factors while the power stays 1. A primitive root of a prime is a number whose order is p - 1, and it can be found by testing each candidate against the prime factors of p - 1. Compare every function with Python's pow on random inputs before you trust it.
Sources and scope
Fermat's little theorem and Euler's theorem are standard results of elementary number theory. Euler's theorem extends Fermat's to any modulus. The proof sketch in this lesson is a short argument written for this text, and the proof of the reduction for a base that shares a factor with the modulus is also written here. The claims about the programs come from runs with g++ 15.2.0 on Linux, compared with a plain loop and with Python's pow. The timing numbers depend on the machine and will differ on yours. The primitive root of 998244353 is shown as 3 by the program, and the same value is also checked by the order test above.