Number-Theory Toolkit

Intuition, Vocabulary and Foundations

A number-theory toolkit is a box of small machines with explicit input and output contracts. The pieces become dependable when they agree on conventions: gcd(a, b) is nonnegative, residues are canonical, failed inverses and incompatible congruences are represented explicitly, and factorization output is sorted. For example, factoring $(84)$ as $(2^2\cdot3\cdot7)$ should feed directly into divisor generation, while the compatible system $(x\equiv2\pmod6)$, $(x\equiv5\pmod9)$ should merge to $(x\equiv14\pmod{18})$.

Think of the API as a set of electrical adapters. Each algorithm may be correct alone, yet composition fails if one returns negative residues, another assumes a prime modulus, or a third uses zero as an error sentinel. The proof obligation therefore has two layers: prove each mathematical transformation and prove that every public contract preserves the next routine's preconditions.

The useful mental model is “normalize at the boundary, preserve an invariant inside, and return a typed success or failure.” Mathematical equivalence is not enough for a library: the representation, domain restrictions, and failure behavior are part of correctness.

Operation	Valid example and result	Invalid or misleading case	Contract and question answered
GCD	`gcd(-30, 18) = 6`	Returning `-6` makes downstream normalization unstable	Return the unique nonnegative greatest common divisor; `gcd(0, 0) = 0` by API convention
Modular inverse	$(4^{-1}\equiv4\pmod5)$	Asking for $(6^{-1}\pmod9)$, where the gcd is $(3)$	Require modulus (m>1) and gcd equal to one; report failure rather than fabricate a value
Modular power	`pow_mod(3, 5, 7) = 5`	Reducing a negative exponent without first finding an inverse	Maintain accumulator-times-base-to-exponent congruent to the original power
SPF sieve	`spf[49] = 7`	Calling `spf[1]` a prime factor	Store the smallest prime factor for every integer from $(2)$ through the bound
Factorization	$(84\mapsto[(2,2),(3,1),(7,1)])$	Unsorted factors or zero exponents	Return increasing prime bases with positive exponents whose product is the input
Divisors	$(12\mapsto[1,2,3,4,6,12])$	Treating zero as having a finite divisor list	Accept positive integers and return each positive divisor exactly once in sorted order
Generalized CRT	$((2\bmod6),(5\bmod9)\mapsto14\bmod18)$	$(x\equiv0\pmod2)$ and $(x\equiv1\pmod2)$	Merge iff residue differences are divisible by the moduli's gcd; otherwise return no solution

Step-by-Step Derivations and Algorithms

Euclid starts from the invariant $(gcd(a,b)=gcd(b,a\bmod b))$; extended Euclid carries coefficients $(s,t)$ so that $(sa+tb=gcd(a,b))$. An inverse of $(a)$ modulo $(m)$ exists exactly when that gcd is one, and the coefficient of $(a)$, reduced modulo $(m)$, is the inverse. Binary exponentiation writes the exponent in base two and repeatedly squares, using only $(O(\log e))$ modular multiplications.

For bulk work, a smallest-prime-factor sieve pays $(O(N\log\log N))$ preprocessing once. Repeated division by spf[n] then factors an integer in $(O(\log n))$ steps, and a Cartesian expansion of prime powers emits its divisors in time proportional to the output size. Generalized CRT composes congruences incrementally. Given $(x\equiv a\pmod m)$ and $(x\equiv b\pmod n)$, write $(x=a+mt)$. The reduced equation $((m/g)t\equiv(b-a)/g\pmod{n/g})$, where $(g=gcd(m,n))$, is solvable precisely when $(g)$ divides (b-a); its solution produces one residue modulo $(operatorname{lcm}(m,n))$.

Number-Theory Toolkit algorithm map

flowchart TD
    A[Validate and normalize public inputs] --> B{Requested operation}
    B -->|GCD or inverse| C[Run Euclid or extended Euclid]
    B -->|Power| D[Run binary modular exponentiation]
    B -->|Bulk factors| E[Build smallest-prime-factor table]
    E --> F[Extract prime powers]
    F --> G[Generate divisors or statistics]
    B -->|Congruences| H[Merge one CRT pair]
    H --> I{Compatible?}
    I -->|No| J[Return explicit failure]
    I -->|Yes| K[Normalize residue and continue]

Every routine maintains a checkable equality or congruence, and every composition validates the next routine's domain first. With an SPF table through $(N)$, preprocessing costs $(O(N\log\log N))$ time and $(O(N))$ space; gcd, inverse, power, and each CRT merge use logarithmically many arithmetic steps, while divisor generation necessarily costs $(\Theta(\tau(n)))$ output time.

Implementation, Complexity and Engineering

The reference implementation below uses Python integers so the mathematical contracts are visible without machine-overflow noise. A fixed-width port must widen products or use overflow-safe modular multiplication, especially in CRT reconstruction and modular exponentiation. Public functions reject invalid domains with ValueError; ordinary mathematical nonexistence uses None, allowing callers to distinguish malformed input from a valid but unsatisfiable request.

Test composition as well as examples. Useful oracles include Python's math.gcd, direct multiplication of returned prime powers, brute-force divisor enumeration for small inputs, and substitution of every CRT answer into all original congruences. Property tests should include signs, zero, unit inputs, incompatible systems, and values at the precomputation boundary.

from dataclasses import dataclass
from math import gcd as oracle_gcd


@dataclass(frozen=True)
class CRTSolution:
    residue: int
    modulus: int


def gcd(a: int, b: int) -> int:
    a, b = abs(a), abs(b)
    while b:
        a, b = b, a % b
    return a


def extended_gcd(a: int, b: int) -> tuple[int, int, int]:
    old_r, r = abs(a), abs(b)
    old_s, s = 1, 0
    old_t, t = 0, 1
    while r:
        q = old_r // r
        old_r, r = r, old_r - q * r
        old_s, s = s, old_s - q * s
        old_t, t = t, old_t - q * t
    return old_r, old_s * (-1 if a < 0 else 1), old_t * (-1 if b < 0 else 1)


def inverse_mod(a: int, modulus: int) -> int | None:
    if modulus <= 1:
        raise ValueError("modulus must exceed one")
    g, x, _ = extended_gcd(a, modulus)
    return x % modulus if g == 1 else None


def pow_mod(base: int, exponent: int, modulus: int) -> int | None:
    if modulus <= 1:
        raise ValueError("modulus must exceed one")
    if exponent < 0:
        base = inverse_mod(base, modulus)
        if base is None:
            return None
        exponent = -exponent
    result = 1 % modulus
    base %= modulus
    while exponent:
        if exponent & 1:
            result = result * base % modulus
        base = base * base % modulus
        exponent >>= 1
    return result


def smallest_prime_factors(limit: int) -> list[int]:
    if limit < 1:
        raise ValueError("limit must be positive")
    spf = list(range(limit + 1))
    spf[0] = spf[1] = 0
    p = 2
    while p * p <= limit:
        if spf[p] == p:
            for multiple in range(p * p, limit + 1, p):
                if spf[multiple] == multiple:
                    spf[multiple] = p
        p += 1
    return spf


def factorize(value: int, spf: list[int]) -> list[tuple[int, int]]:
    if value < 1 or value >= len(spf):
        raise ValueError("value must be in the SPF table's positive range")
    factors = []
    while value > 1:
        prime = spf[value]
        exponent = 0
        while value % prime == 0:
            value //= prime
            exponent += 1
        factors.append((prime, exponent))
    return factors


def divisors_from_factors(factors: list[tuple[int, int]]) -> list[int]:
    divisors = [1]
    previous = 1
    for prime, exponent in factors:
        if prime <= previous or exponent <= 0:
            raise ValueError("factors must use increasing bases and positive exponents")
        powers = [1]
        for _ in range(exponent):
            powers.append(powers[-1] * prime)
        divisors = [d * power for d in divisors for power in powers]
        previous = prime
    return sorted(divisors)


def divisors(value: int, spf: list[int]) -> list[int]:
    return divisors_from_factors(factorize(value, spf))


def crt_merge(left: CRTSolution, right: CRTSolution) -> CRTSolution | None:
    m, n = left.modulus, right.modulus
    if m <= 0 or n <= 0:
        raise ValueError("CRT moduli must be positive")
    a, b = left.residue % m, right.residue % n
    g, coefficient, _ = extended_gcd(m, n)
    difference = b - a
    if difference % g:
        return None
    reduced_modulus = n // g
    step = (difference // g * coefficient) % reduced_modulus
    combined_modulus = m * reduced_modulus
    return CRTSolution((a + m * step) % combined_modulus, combined_modulus)


def crt(congruences: list[tuple[int, int]]) -> CRTSolution | None:
    answer = CRTSolution(0, 1)
    for residue, modulus in congruences:
        answer = crt_merge(answer, CRTSolution(residue, modulus))
        if answer is None:
            return None
    return answer


for a in range(-30, 31):
    for b in range(-30, 31):
        assert gcd(a, b) == oracle_gcd(a, b)
        g, x, y = extended_gcd(a, b)
        assert g == gcd(a, b) and a * x + b * y == g

assert inverse_mod(4, 5) == 4
assert inverse_mod(6, 9) is None
assert pow_mod(3, 5, 7) == 5
assert pow_mod(3, -1, 7) == 5

spf = smallest_prime_factors(200)
assert factorize(84, spf) == [(2, 2), (3, 1), (7, 1)]
for value in range(1, 201):
    expected = [d for d in range(1, value + 1) if value % d == 0]
    assert divisors(value, spf) == expected

assert crt([(2, 6), (5, 9)]) == CRTSolution(14, 18)
assert crt([(0, 2), (1, 2)]) is None
assert crt([]) == CRTSolution(0, 1)
try:
    crt([(0, 0)])
    raise AssertionError("zero modulus was accepted")
except ValueError:
    pass

print("trace:", factorize(84, spf), divisors(12, spf), crt([(2, 6), (5, 9)]))
print("all tests passed")

Choice	Cost and suitable workload	Numeric or semantic risk	Engineering mitigation
Euclid per call	$(O(\log\min(\lvert a\rvert,\lvert b\rvert)))$; ideal for isolated gcds	Sign inconsistency at API boundaries	Normalize inputs and guarantee a nonnegative result
Binary modular power	$(O(\log e))$ multiplications, $(O(1))$ auxiliary space	Product overflow in fixed-width languages	Use a widened type, Montgomery reduction, or checked multiplication
SPF table	$(O(N\log\log N))$ build and $(O(N))$ memory; strong for many bounded queries	Out-of-range indexing and excessive memory	Validate bounds; use trial division when query volume is small
Trial division	$(O(\sqrt n))$ worst case and constant auxiliary space	Too slow for many large inputs	Switch to a sieve or stronger factorization backend at a measured threshold
Divisor materialization	$(\Theta(\tau(n)))$ output work	Output can be large; duplicates arise from malformed factors	Validate factor tuples and sort once at the boundary
Incremental generalized CRT	One extended-gcd merge per congruence	Growing modulus may overflow; systems may contradict	Widen intermediates, reduce after each merge, and return `None` on incompatibility
Exceptions plus optional results	Invalid domains raise; mathematical nonexistence returns `None`	Conflating malformed input with “no solution”	Document the distinction and test both paths

Applications, Case Studies and Exercises

An approachable use case is a contest CLI that reads a bound, builds one SPF table, and answers factor, divisors, inverse, and crt commands. The key design decision is shared normalization: every residue printed by the CLI lies in its canonical interval, and every failed query has an unambiguous result. A production-shaped case is a service with interchangeable backends: a bounded SPF cache for common values, a fallback factorizer for larger ones, overflow-safe arithmetic in fixed-width deployments, typed error results, and metrics that reveal which threshold actually minimizes latency.

Exercises should climb in scope. Recall: state the existence condition for a modular inverse and a generalized CRT merge. Derivation: prove the reduced CRT equation and its lcm modulus. Implementation: add Euler's phi from the factor list without refactoring the public contracts. Debugging: find a signed-residue bug and a zero-modulus bug with property tests. Extension: add JSON CLI output, benchmark sieve crossover points, and write a verification report covering every exported function.

Deliverable

Apply Number-Theory Toolkit

Case study 1 — trace a small session. Build spf = smallest_prime_factors(100). Record the outputs for factorize(84), divisors(12), inverse_mod(4, 5), and CRT on (2, 6) with (5, 9). Checkpoint: obtain [(2, 2), (3, 1), (7, 1)], [1, 2, 3, 4, 6, 12], 4, and CRTSolution(14, 18). Verify every divisor divides its source and substitute $(14)$ into both congruences.

Case study 2 — production-shaped batch. Generate 10,000 deterministic queries with values at most $(10^6)$, 60 percent factor queries, 25 percent power or inverse queries, and 15 percent two-way CRT queries. Compare one shared SPF table with trial division, recording preprocessing time, median query time, peak memory, and mismatches against brute-force checks on a small sampled prefix. Checkpoint: determine the measured query count at which preprocessing wins; success requires zero oracle mismatches and explicit counts for invalid and unsatisfiable inputs.

Hands-on deliverable. Package the module with a CLI supporting gcd, inverse, power, factor, divisors, and crt. Add help text containing domains and complexity, then run the embedded tests plus randomized properties: Bézout substitution, inverse multiplication, factor-product reconstruction, divisor completeness, and CRT substitution. Success means every command has a documented failure form, all boundary tests pass, and a fresh user can reproduce the trace from one command.

Choose the simplest method that satisfies the domain and workload: direct Euclid for isolated arithmetic, precomputation for repeated bounded queries, and generalized CRT only after compatibility checks. Do not mistake “works on positive toy inputs” for an API contract; use the prerequisite lessons on congruences and sieves to justify invariants, then extend the toolkit through new functions that consume the same canonical representations.

Cookie Consent