Number Theory and Bit Tricks
Modular arithmetic without overflow, inverses, sieves, binomial coefficients modulo a prime, the Chinese remainder theorem and the C++20 bit functions, each checked against brute force.
Answers that do not fit in 64 bits
Counting problems end with a familiar line: "print the answer modulo 10^9 + 7". The count itself can be astronomically large: the number of ways to order 30 items is 30!, about 2.65·10^32, far beyond the largest long long (about 9.22·10^18). The problem asks for the remainder so the answer fits in a machine integer, and so that the judge can compare one number.
The trap is that the remainder must be taken while computing, not at the end. A product that overflows has already lost information; taking % kMod afterwards produces a confident wrong number. This program computes n! modulo 10^9 + 7 both ways. The "reduce at the end" version uses std::uint64_t, whose overflow is defined to wrap around modulo 2^64, so the program runs and quietly computes the wrong thing:
#include <cstdint>
#include <print>
using ll = long long;
constexpr ll kMod = 1'000'000'007;
int main() {
// Wrong: multiply everything, reduce once at the end.
// (std::uint64_t wraps silently, so this runs; it just computes garbage.)
// Right: reduce after every multiplication, so values stay below kMod.
std::uint64_t full = 1;
ll reduced = 1;
for (int n = 1; n <= 25; ++n) {
full *= std::uint64_t(n);
reduced = reduced * n % kMod;
if (n >= 19 && n <= 22)
std::println("{:>2}! mod p: reduce at end {:>10}, reduce each step {:>10}", n,
full % kMod, reduced);
}
// % truncates toward zero, so a negative left operand gives a negative result.
ll a = 3, b = 5, m = 7;
std::println("(3 - 5) % 7 = {}", (a - b) % m);
std::println("((3 - 5) % 7 + 7) % 7 = {}", ((a - b) % m + m) % m);
std::println("-7 / 3 = {}, -7 % 3 = {}", -7 / 3, -7 % 3);
}Press Run to compile and run this example.
19! mod p: reduce at end 557316307, reduce each step 557316307
20! mod p: reduce at end 146326063, reduce each step 146326063
21! mod p: reduce at end 908159300, reduce each step 72847302
22! mod p: reduce at end 662000402, reduce each step 602640637
(3 - 5) % 7 = -2
((3 - 5) % 7 + 7) % 7 = 5
-7 / 3 = -2, -7 % 3 = -120! is about 2.43·10^18 and still fits in 64 bits, so both columns agree. 21! does not fit, and from then on the first column is wrong. With int or long long instead of std::uint64_t it would be worse: signed overflow is undefined behavior, and the sanitizers used by this course stop the program.
The last three lines show the second trap. Since C++11, integer division truncates toward zero, and % is defined so that (a / b) * b + a % b == a. So -7 % 3 is -1, not 2, and a difference of two residues can be negative.
This lesson collects the number theory that contest and interview problems actually use: modular arithmetic, gcd and extended Euclid, fast powers and inverses, sieves, binomial coefficients, the Chinese remainder theorem, and the bit functions of C++20. Every piece is checked against a brute force, as in Fast I/O, Contest Templates and Stress Testing. The last lesson, The Playbook, then puts these tools and those of the earlier lessons to work on problems you have not seen.
Modular arithmetic: the rules that keep numbers small
Write a mod m for the remainder in the range [0, m). Addition, subtraction and multiplication are compatible with taking remainders:
(a + b) mod m = ((a mod m) + (b mod m)) mod m(a - b) mod m = ((a mod m) - (b mod m) + m) mod m, with the+ mmaking the C++ result non-negative(a · b) mod m = ((a mod m) · (b mod m)) mod m
So you may reduce after every operation, and you must, because that keeps every intermediate value below m. The question then is how large an intermediate can get before it is reduced. For a product of two residues it is (m - 1)^2. That fits in a signed 64-bit integer as long as m is below about 3.03·10^9, which covers the usual moduli 10^9 + 7 and 998244353. For a modulus near 10^18, the product needs about 120 bits: multiply in __int128, which GCC and Clang provide on 64-bit targets as an extension (it is not standard C++), as the hashing code in the previous lesson did.
Division is the odd one out. (a / b) mod m is not (a mod m) / (b mod m); for example (12 / 4) mod 7 = 3, but 12 mod 7 = 5 is not divisible by 4 at all. Instead, you multiply by the modular inverse of b: a number b' with b · b' mod m = 1. It exists exactly when gcd(b, m) = 1, which is one reason the standard moduli are primes: modulo a prime, every nonzero value has an inverse.
gcd, lcm and extended Euclid
std::gcd and std::lcm live in <numeric> since C++17. std::gcd(m, n) works on the absolute values and returns 0 when both arguments are 0. The behavior is undefined if an argument's absolute value does not fit in the common type, and for std::lcm, also if the result does not fit. So std::lcm of two values near 10^12 in long long is not "a wrong number": it is undefined. When the lcm might overflow, compute a / std::gcd(a, b) * b yourself, dividing first, and check the bound before multiplying.
Euclid's algorithm rests on gcd(a, b) = gcd(b, a mod b): any common divisor of a and b also divides a - q·b. Each two steps at least halve the larger number, so it takes O(log min(a, b)) steps.
The extended algorithm also finds integers x and y with a·x + b·y = gcd(a, b) (Bezout's identity). The recursion computes them for (b, a mod b) first and then rewrites a mod b as a - (a / b)·b:
#include <algorithm>
#include <cstdlib>
#include <numeric>
#include <print>
using ll = long long;
struct Bezout { ll g, x, y; }; // a*x + b*y == g
// Extended Euclid: the same steps as gcd, carrying the coefficients along.
Bezout extGcd(ll a, ll b) {
if (b == 0) return {a, 1, 0};
auto [g, x, y] = extGcd(b, a % b); // b*x + (a % b)*y == g
return {g, y, x - (a / b) * y}; // rewrite a % b as a - (a / b)*b
}
ll bruteGcd(ll a, ll b) {
for (ll d = std::max(a, b); d >= 1; --d)
if (a % d == 0 && b % d == 0) return d;
return 0; // only when a == b == 0
}
int main() {
auto [g, x, y] = extGcd(240, 46);
std::println("gcd(240, 46) = {} = 240*({}) + 46*({})", g, x, y);
std::println("std::gcd(240, 46) = {}, std::lcm(240, 46) = {}",
std::gcd(240, 46), std::lcm(240, 46));
long long checked = 0, wrong = 0, large = 0;
for (ll a = 1; a <= 300; ++a)
for (ll b = 1; b <= 300; ++b) {
auto r = extGcd(a, b);
++checked;
wrong += r.g != bruteGcd(a, b) || r.g != std::gcd(a, b);
wrong += a * r.x + b * r.y != r.g;
// the coefficients never outgrow the inputs: |x| <= b/g, |y| <= a/g
large += std::abs(r.x) * r.g > b || std::abs(r.y) * r.g > a;
}
std::println("checked {} pairs: {} wrong, {} with large coefficients",
checked, wrong, large);
}Press Run to compile and run this example.
gcd(240, 46) = 2 = 240*(-9) + 46*(47)
std::gcd(240, 46) = 2, std::lcm(240, 46) = 5520
checked 90000 pairs: 0 wrong, 0 with large coefficientsCheck the first line by hand: 240·(-9) = -2160 and 46·47 = 2162, which sum to 2. The last check matters in practice: on every pair tested, the coefficients stay within b / g and a / g, so extended Euclid on 64-bit inputs does not overflow on its own. The products you form from its results later still can.
Fast exponentiation and modular inverses
Computing b^e mod m by multiplying e times is hopeless for e = 10^18. Repeated squaring writes e in binary: b^13 = b^8 · b^4 · b^1 because 13 = 1101 in binary. The loop squares the base once per bit of e and multiplies it into the result when that bit is 1, so it needs at most two multiplications per bit of e, about 2·log2(e).
That gives the first way to invert. Fermat's little theorem says that for a prime p and a not divisible by p, a^(p-1) mod p = 1. So a^(p-2) is the inverse of a. The second way works for any modulus: if gcd(a, m) = 1, extended Euclid gives a·x + m·y = 1, and reducing modulo m leaves a·x mod m = 1. The program counts the multiplications of the fast power and stress tests both inverses against a direct search over every a in [-300, 300] and every modulus up to 300:
#include <print>
using ll = long long;
int multiplications = 0;
ll mulMod(ll a, ll b, ll m) {
++multiplications;
return a * b % m; // safe while m <= about 3*10^9: a*b < 2^63
}
// base^e mod m by repeated squaring: O(log e) multiplications.
ll powMod(ll base, ll e, ll m) {
ll result = 1 % m;
base = ((base % m) + m) % m; // a negative base would give negative results
while (e > 0) {
if (e & 1) result = mulMod(result, base, m);
base = mulMod(base, base, m);
e >>= 1;
}
return result;
}
struct Bezout { ll g, x, y; };
Bezout extGcd(ll a, ll b) {
if (b == 0) return {a, 1, 0};
auto [g, x, y] = extGcd(b, a % b);
return {g, y, x - (a / b) * y};
}
// Inverse of a modulo any m >= 2, or -1 if gcd(a, m) != 1.
ll invGeneral(ll a, ll m) {
auto [g, x, y] = extGcd(((a % m) + m) % m, m);
return g == 1 ? ((x % m) + m) % m : -1;
}
// Inverse modulo a prime p, for a not divisible by p (Fermat).
ll invPrime(ll a, ll p) { return powMod(a, p - 2, p); }
int main() {
constexpr ll kMod = 1'000'000'007;
std::println("{:>20} {:>14}", "exponent", "multiplications");
for (ll e : {10LL, 1'000LL, 1'000'000LL, 1'000'000'000LL, kMod - 2,
1'000'000'000'000'000'000LL}) {
multiplications = 0;
powMod(3, e, kMod);
std::println("{:>20} {:>14}", e, multiplications);
}
ll a = 123'456'789;
ll inv = invPrime(a, kMod);
std::println("inverse of {} mod p = {}, check: {}", a, inv, a * inv % kMod);
// Stress: every a in [-300, 300], every modulus m in [2, 300], brute force.
long long checked = 0, wrong = 0;
for (ll m = 2; m <= 300; ++m) {
bool prime = true;
for (ll d = 2; d * d <= m; ++d) prime = prime && m % d != 0;
for (ll x = -300; x <= 300; ++x) {
ll brute = -1;
for (ll c = 0; c < m && brute < 0; ++c)
if ((((x % m) + m) % m) * c % m == 1) brute = c;
++checked;
wrong += invGeneral(x, m) != brute;
if (prime && x % m != 0) wrong += invPrime(x, m) != brute;
}
}
std::println("checked {} (a, m) pairs, {} wrong", checked, wrong);
}Press Run to compile and run this example.
exponent multiplications
10 6
1000 16
1000000 27
1000000000 43
1000000005 45
1000000000000000000 84
inverse of 123456789 mod p = 18633540, check: 1
checked 179699 (a, m) pairs, 0 wrongAn exponent of 10^18 costs 84 multiplications; each factor of 1000 in the exponent adds about 10 squarings and up to 10 more multiplications. Computing one Fermat inverse modulo 10^9 + 7 costs 45. That is cheap once, but not inside a loop that runs a million times, which is why the next sections compute inverses in bulk.
The line base = ((base % m) + m) % m is there because the stress test includes negative a. With the shorter base %= m, the same test reports 17,872 wrong answers out of 179,699: a negative base makes every product negative.
Sieves: all primes up to n
Testing each number for primality by trial division up to its square root costs O(n·sqrt(n)) for all numbers up to n. The sieve of Eratosthenes turns it around: for each prime p, cross out its multiples. Two details make it fast. Only primes need to cross anything out, and crossing out can start at p·p, because every smaller multiple k·p with k below p has a smaller prime factor and was crossed out already. The total work is the sum of n / p over primes p, which is O(n log log n).
The linear sieve (Gries and Misra, 1978) crosses out every composite exactly once: the number i·p is marked only by its smallest prime factor p. As a by-product it stores that smallest prime factor, spf[x], for every x. The program counts the cross-outs of both:
#include <cmath>
#include <print>
#include <vector>
// Sieve of Eratosthenes; returns how many times a cell was crossed out.
long long eratosthenes(int n, std::vector<char>& composite) {
composite.assign(n + 1, 0);
long long marks = 0;
for (long long p = 2; p * p <= n; ++p) {
if (composite[p]) continue;
for (long long q = p * p; q <= n; q += p) { // smaller multiples are done
composite[q] = 1;
++marks;
}
}
return marks;
}
// Linear sieve: spf[x] = smallest prime factor of x. Each composite is set once.
long long linearSieve(int n, std::vector<int>& spf, std::vector<int>& primes) {
spf.assign(n + 1, 0);
primes.clear();
long long marks = 0;
for (int i = 2; i <= n; ++i) {
if (spf[i] == 0) { spf[i] = i; primes.push_back(i); }
for (int p : primes) {
if (p > spf[i] || 1LL * i * p > n) break;
spf[i * p] = p; // p is the smallest prime factor of i*p
++marks;
}
}
return marks;
}
int main() {
std::println("{:>9} {:>8} {:>11} {:>10} {:>9} {:>11}", "n", "primes",
"erat marks", "per n", "ln ln n", "linear");
for (int n : {1'000, 10'000, 100'000, 1'000'000, 10'000'000}) {
std::vector<char> composite;
std::vector<int> spf, primes;
long long e = eratosthenes(n, composite);
long long l = linearSieve(n, spf, primes);
std::println("{:>9} {:>8} {:>11} {:>10.3f} {:>9.3f} {:>11}", n, primes.size(), e,
double(e) / n, std::log(std::log(double(n))), l);
}
}Press Run to compile and run this example.
n primes erat marks per n ln ln n linear
1000 168 1411 1.411 1.933 831
10000 1229 16981 1.698 2.220 8770
100000 9592 193078 1.931 2.443 90407
1000000 78498 2122048 2.122 2.626 921501
10000000 664579 22850051 2.285 2.780 9335420Read the "per n" column: multiplying n by 10 four times raises the marks per number only from 1.4 to 2.3, and the increase tracks ln ln n (the gap between the two columns stays near 0.5). That is what O(n log log n) looks like; for any n you will meet, it behaves like a small constant times n. The linear sieve's count is exactly the number of composites up to n (for n = 1000: 1000 - 1 - 168 = 831).
Fewer marks do not automatically mean faster. The linear sieve stores an int per number instead of a byte and has a more complex inner loop. Use Eratosthenes when you need only "is x prime"; use the linear sieve when you need spf, which the next section uses.
Factorizing many numbers with smallest prime factors
With spf in hand, factorizing x is a loop: divide by spf[x] until it no longer divides, move on to the new spf[x]. Each division at least halves x, so a factorization costs O(log x), independent of how large the primes are.
An original problem to try it on. A puzzle generator draws numbers up to 10^6 and needs, for each one, its prime factorization and its number of divisors. If x = p1^e1 · p2^e2 · ... · pk^ek, a divisor picks each exponent independently from 0 to ei, so there are (e1 + 1)(e2 + 1)...(ek + 1) divisors. Input: the number of queries, then the queries.
4
1
360
999983
720720#include <iostream>
#include <print>
#include <sstream>
#include <utility>
#include <vector>
constexpr int kMaxX = 1'000'000;
std::vector<int> smallestPrimeFactors(int n) {
std::vector<int> spf(n + 1, 0), primes;
for (int i = 2; i <= n; ++i) {
if (spf[i] == 0) { spf[i] = i; primes.push_back(i); }
for (int p : primes) {
if (p > spf[i] || 1LL * i * p > n) break;
spf[i * p] = p;
}
}
return spf;
}
// (prime, exponent) pairs in increasing order: O(log x) divisions.
std::vector<std::pair<int, int>> factorize(int x, const std::vector<int>& spf) {
std::vector<std::pair<int, int>> f;
while (x > 1) {
int p = spf[x], e = 0;
while (x % p == 0) { x /= p; ++e; }
f.push_back({p, e});
}
return f;
}
std::vector<std::pair<int, int>> trialDivision(int x) {
std::vector<std::pair<int, int>> f;
for (int p = 2; 1LL * p * p <= x; ++p) {
int e = 0;
while (x % p == 0) { x /= p; ++e; }
if (e > 0) f.push_back({p, e});
}
if (x > 1) f.push_back({x, 1});
return f;
}
void solve(std::istream& in, std::ostream& out, const std::vector<int>& spf) {
int q;
in >> q;
while (q--) {
int x;
in >> x;
auto factors = factorize(x, spf);
long long divisors = 1;
out << x << ':';
if (factors.empty()) out << " (no prime factors)";
for (auto [p, e] : factors) {
out << ' ' << p << '^' << e;
divisors *= e + 1; // choose each exponent from 0 to e
}
out << " -> " << divisors << " divisors\n";
}
}
int main() {
auto spf = smallestPrimeFactors(kMaxX);
std::istringstream in("4\n1\n360\n999983\n720720\n");
solve(in, std::cout, spf);
int wrong = 0;
for (int x = 1; x <= 200'000; ++x) wrong += factorize(x, spf) != trialDivision(x);
std::println("x = 1..200000 against trial division: {} wrong", wrong);
}Press Run to compile and run this example.
1: (no prime factors) -> 1 divisors
360: 2^3 3^2 5^1 -> 24 divisors
999983: 999983^1 -> 2 divisors
720720: 2^4 3^2 5^1 7^1 11^1 13^1 -> 240 divisors
x = 1..200000 against trial division: 0 wrongThe input includes the two edge cases worth testing every time: 1, which has no prime factors, and a large prime (999983), where spf[x] == x. Building spf up to 10^6 takes O(n) time and 4 MB of memory; after that each query is O(log x). If the numbers go up to 10^12 instead, a sieve up to 10^12 is impossible, but a sieve of the primes up to 10^6 still lets you factor by trial division over primes only.
Counting with binomial coefficients modulo p
C(n, k), the number of ways to choose k items from n, appears in almost every counting problem. The formula n! / (k! (n - k)!) contains a division, so modulo a prime it becomes n! · inv(k!) · inv((n - k)!). Precompute the factorials up to the largest n, then the inverse factorials. Computing each inverse by Fermat would cost a log p factor per entry; there is a better way. Invert only the largest factorial, then walk down with inv((i - 1)!) = inv(i!) · i, because (i - 1)! = i! / i. After O(n + log p) preprocessing, every C(n, k) is two modular multiplications.
This needs n below p: if n reaches p, then n! is divisible by p, has no inverse, and the formula breaks. With p = 10^9 + 7 and n up to a few million, that never happens. For small primes and huge n, look up Lucas's theorem.
An original problem. A delivery drone flies over a city grid from (0, 0) to (a, b), one block east or one block north at a time. How many different routes are there, modulo 10^9 + 7? Every route has a + b moves, and choosing which a of them go east fixes the route, so the answer is C(a + b, a). Input: the number of queries, then a b per query, with a and b up to 5·10^5.
3
2 2
10 7
500000 500000#include <iostream>
#include <print>
#include <sstream>
#include <vector>
using ll = long long;
constexpr ll kMod = 1'000'000'007; // prime, and larger than any n we use
ll powMod(ll base, ll e) {
ll result = 1;
base %= kMod;
for (; e > 0; e >>= 1, base = base * base % kMod)
if (e & 1) result = result * base % kMod;
return result;
}
struct Binomial {
std::vector<ll> fact, invFact;
explicit Binomial(int n) : fact(n + 1), invFact(n + 1) {
fact[0] = 1;
for (int i = 1; i <= n; ++i) fact[i] = fact[i - 1] * i % kMod;
invFact[n] = powMod(fact[n], kMod - 2); // one Fermat inverse...
for (int i = n; i > 0; --i) // ...then walk down: 1/(i-1)! = i/i!
invFact[i - 1] = invFact[i] * i % kMod;
}
ll choose(int n, int k) const {
if (k < 0 || k > n) return 0;
return fact[n] * invFact[k] % kMod * invFact[n - k] % kMod;
}
};
// Routes from (0, 0) to (a, b) moving one block east or north: choose which
// a of the a + b moves go east.
void solve(std::istream& in, std::ostream& out, const Binomial& c) {
int q;
in >> q;
while (q--) {
int a, b;
in >> a >> b;
out << c.choose(a + b, a) << '\n';
}
}
int main() {
Binomial c(1'000'000);
std::istringstream in("3\n2 2\n10 7\n500000 500000\n");
solve(in, std::cout, c);
// Stress: Pascal's triangle, computed with additions only.
const int n = 300;
std::vector<std::vector<ll>> pascal(n + 1, std::vector<ll>(n + 1, 0));
int wrong = 0, checked = 0;
for (int i = 0; i <= n; ++i) {
pascal[i][0] = 1;
for (int k = 1; k <= i; ++k)
pascal[i][k] = (pascal[i - 1][k - 1] + pascal[i - 1][k]) % kMod;
for (int k = -1; k <= i + 1; ++k) {
ll expected = (k < 0 || k > i) ? 0 : pascal[i][k];
++checked;
wrong += c.choose(i, k) != expected;
}
}
std::println("checked {} values against Pascal's triangle, {} wrong", checked, wrong);
}Press Run to compile and run this example.
6
19448
996692777
checked 46053 values against Pascal's triangle, 0 wrongThe first answer is small enough to check by drawing: from (0, 0) to (2, 2) there are 6 routes. The reference here is Pascal's rule C(n, k) = C(n - 1, k - 1) + C(n - 1, k), which uses only additions and so cannot share a bug with the inverse-based code. The test also asks for k = -1 and k = n + 1, where the answer is 0; the bounds check in choose is what makes those safe, and those are exactly the calls a dynamic program makes at its edges.