• Primes
• Factorization
• Arithmetic Functions
• Modular Arithmetic
• Nonlinear Congruences
• Diophantine Equations
• Algebraic Systems
• Lattices
• Appendix

is_prime(n: int) -> bool

[source]

Test if a given integer n is prime.

Uses a combination of trial division, the Miller-Rabin primality test with deterministic bases, or the extra-strong variant of the Baillie-PSW primality test (this variant has no known pseudoprimes in any range, and has been computationally verified to have no counterexamples for all n < 2^64).

See: https://www.techneon.com/download/is.prime.32.base.data (MR hash for n < 2^32) See: https://miller-rabin.appspot.com (other deterministic MR base sets) See: https://ntheory.org/pseudoprimes.html (BPSW verification up to 2^64)

next_prime(n: int) -> int

[source]

Get the smallest prime number greater than n.

random_prime(num_bits: int, *, safe: bool=False) -> int

[source]

Generate a random prime with the given number of bits.

primes(*, low: int=2, high: int | None=None, count: int | None=None) -> Iterator[int]

[source]

Generate at most count primes in increasing order within the range [low, high].

Uses the sieve of Eratosthenes, with a segmented approach for large or unbounded ranges.

count_primes(x: int) -> int

[source]

Prime counting function π(x). Returns the number of primes p ≤ x.

Uses the Lagarias-Miller-Odlyzko (LMO) extension of the Meissel-Lehmer algorithm.

sum_primes(x: int, f: Callable[[int], Number] | None=None, f_prefix_sum: Callable[[int], Number] | None=None) -> Number

[source]

Compute F(x) as the sum of f(p) over all primes p ≤ x, where f is a completely multiplicative function (by default, f(n) = n).

Uses a generalized version of the LMO prime counting algorithm. Ideally f() and f_prefix_sum() can be calculated efficiently in O(1) time via closed-form expression.

perfect_power(n: int) -> tuple[int, int]

[source]

Find integers a, b such that a^b = n.

Returns the solution (a, b) with minimal b > 1 if there are any such solutions, otherwise returns the trivial solution (n, 1).

prime_factors(n: int) -> tuple[int, ...]

[source]

Get all prime factors of n in sorted order (with multiplicity).

Uses a combination of trial division, Fermat's factorization method, Brent's variant of Pollard's rho, Lenstra's elliptic curve method (ECM), and a self-initializing quadratic sieve (SIQS).

prime_factorization(n: int) -> dict[int, int]

[source]

Get the prime factorization of n as a dictionary of {prime: exponent}.

divisors(n: int) -> tuple[int, ...]

[source]

Get all positive divisors of n in sorted order.

omega(n: int) -> int

[source]

Compute the value of ω(n), the number of distinct prime factors of n.

big_omega(n: int) -> int

[source]

Compute the value of Ω(n), the number of prime factors of n (with multiplicity).

divisor_count(n: int) -> int

[source]

Compute the value of σ₀(n), the number of divisors of n.

divisor_sum(n: int) -> int

[source]

Compute the value of σ₁(n), the sum of divisors of n.

divisor_function(n: int, k: int=1) -> int

[source]

Compute the value of the divisor function σₖ(n), where σₖ(n) = ∑_{d|n} dᵏ.

partition(n: int, mod: int | None=None, restrict: Callable[[int], bool] | None=None) -> int

[source]

Return the value of the partition function p(n).

radical(n: int) -> int

[source]

Compute rad(n) as the product of the distinct prime factors of n.

mobius(n: int) -> int

[source]

Compute the Mobius function μ(n) for a positive integer n.

totient(n: int) -> int

[source]

Compute Euler's totient function φ(n) for a positive integer n.

carmichael(n: int) -> int

[source]

Compute Carmichael's lambda function λ(n) for a positive integer n.

valuation(n: int, p: int) -> int

[source]

Compute the p-adic valuation νₚ(n), the exponent of p in the prime factorization of n.

multiplicative_range(f: Callable[..., int], N: int, f0: int=1) -> list[int]

[source]

Find the value of a multiplicative function f(n) for each n = 0, 1, 2, ..., N - 1. Uses a sieving approach across the range to calculate function values efficiently.

egcd(a: int, b: int) -> tuple[int, int, int]

[source]

Extended Euclidean algorithm.

Returns

Complexity: O(log min(a, b)) time

crt(congruences: Iterable[tuple[int, int]]) -> int | None

[source]

Solve a system of linear congruences x ≡ aᵢ (mod nᵢ) via the Chinese Remainder Theorem.

Returns a solution to the system of congruences, mod the LCM of the moduli, or None if no solution exists.

Supports non-coprime moduli.

coprimes(n: int) -> Iterator[int]

[source]

Generate all integers k in the range [0, n) that are coprime to n.

Returns the reduced residue system modulo n, i.e., the unit group (Z/nZ)×. The size of this set is φ(n) (Euler's totient function).

For small n, uses an O(n) space sieve for speed. For large n, uses an O(1) space generator that checks gcd(k, n) = 1 for each k.

Complexity: O(n * ω(n)) time and O(n) space for n ≤ 10⁷ (sieve approach). O(n log n) time and O(1) space for n > 10⁷ (gcd approach).

multiplicative_order(a: int, mod: int) -> int

[source]

Compute ordₘ(a), the smallest positive integer such that a^k ≡ 1 (mod m).

primitive_root(n: int) -> int | None

[source]

Find a primitive root modulo n.

Use Bach's primitive root finding algorithm to search for candidates.

legendre(a: int, p: int) -> int

[source]

Compute the Legendre symbol (a | p), where p is an odd prime.

jacobi(a: int, n: int) -> int

[source]

Compute the Jacobi symbol (a | n), where n is an odd positive integer.

Complexity: O(log a log n) time

kronecker(a: int, n: int) -> int

[source]

Compute the Kronecker symbol (a | n).

dirichlet_character(m: int, k: int) -> Callable[[int], Number]

[source]

Return the Dirichlet character χₘ‚ₖ : ℤ → ℂ under Conrey labeling, where m is the modulus and k is an index such that gcd(m, k) = 1.

See: https://www.lmfdb.org/knowledge/show/character.dirichlet.conrey

Returns

hensel(coefficients: Sequence[int], p: int, k: int, initial: Iterable[int] | None=None) -> tuple[int, ...]

[source]

Find all solutions to the polynomial congruence f(x) ≡ 0 (mod pᵏ).

Assumes f(x) = a₀ + a₁x + a₂x² + a₃x³ ... is a polynomial. Uses Hensel lifting to find solutions.

Complexity: O(ksd) arithmetic operations, where s is total number of solutions and d = deg(f). O(pd) to find initial solutions if not provided.

polynomial_roots(coefficients: Sequence[int], mod: int) -> tuple[int, ...]

[source]

Find all roots x of a univariate polynomial f(x) ≡ 0 (mod m).

Factors m into prime powers, finds roots modulo each p^e via Cantor-Zassenhaus + Hensel lifting, then combines solutions with CRT.

nth_roots(a: int, n: int, mod: int) -> tuple[int, ...]

[source]

Find all solutions x to x^n ≡ a (mod m).

Uses the Tonelli-Shanks / Adleman-Manders-Miller to find roots modulo primes, Hensel lifting to roots modulo prime powers, and the Chinese Remainder Theorem to combine solutions.

discrete_log(target: int, base: int, mod: int) -> int | None

[source]

Find the smallest non-negative integer x such that target ≡ base^x (mod m). Returns None if no such integer exists.

Uses the Pohlig-Hellman algorithm, with either baby-step giant-step or Pollard's rho for discrete logarithms on the prime-order sub-problems.

bezout(a: int, b: int, c: int) -> Iterator[tuple[int, int]]

[source]

Generate all integer solutions to the linear Diophantine equation ax + by = c.

Uses the extended Euclidean algorithm to find a pair of Bézout coefficients, and then generate an infinite family of solutions.

Yields

Complexity: O(log(min(a, b))) time to find initial solution, O(1) per additional solution.

cornacchia(d: int, m: int) -> Iterator[tuple[int, int]]

[source]

Generate all unique positive integer solutions to the equation x² + dy² = m where 0 < d < m and gcd(d, m) = 1.

Yields

Complexity: O(f(m) + τ(m) g(m) + τ(m) log m) time, where f, τ, g are the cost of factorization, divisor count, and cost of modular roots respectively.

pell(D: int, N: int=1) -> Iterator[tuple[int, int]]

[source]

Generate all unique positive integer solutions to the generalized Pell equation x² - Dy² = N, where D is not a perfect square.

Yields infinite positive integer solutions x, y > 0 in order of increasing x. Uses the Lagrange-Matthews-Mollin (LMM) algorithm.

See: https://cjhb.site/Files.php/Books/math/B3.4/pell.pdf See: http://www.numbertheory.org/PDFS/patz_improved.pdf

Yields

Complexity: O(L(D) + f(|N|) + τ(|N|) * (g(|N|) + L(D))) time, where L, f, τ, g are the continued-fraction period length, cost of factoring, divisor count, and cost of modular roots respectively.

conic(a: int, b: int, c: int, d: int, e: int, f: int) -> Iterator[tuple[int, int]]

[source]

Generate all unique integer solutions (x, y) to the binary quadratic Diophantine conic equation ax² + bxy + cy² + dx + ey + f = 0.

Uses the theory of binary quadratic forms, classifying by discriminant Δ = b² - 4ac:

Δ < 0 (ellipse): Lagrange reduction, finite solutions
Δ = 0 (parabola): parametric families via modular square roots
Δ > 0 (hyperbola): reduction to Pell equation, infinite solutions
Degenerate cases: factorization into linear forms

Yields

pythagorean_triples(max_c: float | None=None, max_sum: float | None=None) -> Iterator[tuple[int, int, int]]

[source]

Generate positive integer solutions to the equation a² + b² = c².

Uses Euclid's formula to generate unique Pythagorean triples (a, b, c) where a ≤ b ≤ c.

If no bounds are specified, infinitely generates triples in order of increasing c. When bounds are specified, no order is guaranteed.

pillai(a: int, b: int, c: int) -> Iterator[tuple[int, int]]

[source]

Generate all positive integer solutions (x, y) to the exponential Diophantine Pillai equation aˣ - bʸ = c, where a, b >= 2 and x, y > 0.

Yields

solve_linear_system(A: Matrix[int], b: Vector[int] | None=None, *, nullspace: bool=False) -> tuple[Vector[int] | None, list[Vector[int]] | None]

[source]

Find integer solutions to the system of linear equations given by Ax = b.

Returns

solve_polynomial_system(polynomials: list[Polynomial[int]], bounds: tuple[int, ...]) -> tuple[tuple[int, ...], ...]

[source]

Find integer solutions to a system of multivariate polynomial equations f₁(x₁, x₂, ...) = f₂(x₁, x₂, ...) = ... = fₖ(x₁, x₂, ...) = 0, where |xᵢ| < bounds[i] for each variable.

Polynomials are represented as dictionaries mapping monomial tuples to coefficients, e.g. {(2, 0): 3, (0, 1): -5, (0, 0): 7} represents 3x² - 5y + 7.

Uses the Gao-Volny-Wang (GVW) algorithm to find a Grobner grevlex basis over Fₚ, the Faugère-Gianni-Lazard-Mora (FGLM) algorithm to convert to a lex basis, and either extracts univariate polynomials or solves the system via backtracking. Prunes by directly solving the linear subsystem, if one exists.

lll_reduce(B: Matrix[int]) -> Matrix[int]

[source]

Lenstra-Lenstra-Lovász (LLL) lattice basis reduction.

Returns a reduced basis with shorter, more orthogonal vectors, satisfying:

Size-reduction: |μ_{i,j}| ≤ 0.5 for all i > j
Lovász condition: δ‖b*_k‖² ≤ ‖b*_{k+1}‖² + μ_{k+1,k}² ‖b*_k‖²

Uses floating-point arithmetic for speed, with automatic escalation to exact rational arithmetic if precision issues are detected.

See: https://www.cs.cmu.edu/~avrim/451f11/lectures/lect1129_LLL.pdf

Complexity: O(n⁵d log³B) time for n × d matrix with max entry size B, O(n² + nd) space

bkz_reduce(B: Matrix[int], block_size: int=20) -> Matrix[int]

[source]

BKZ (Block Korkine-Zolotarev) lattice basis reduction.

BKZ generalizes LLL by applying an SVP (Shortest Vector Problem) oracle to sliding blocks of consecutive basis vectors.

Uses Schnorr-Euchner enumeration for the SVP oracle.

See: https://www.sciencedirect.com/science/article/pii/0304397587900648

Complexity: O(2^(0.25β²)) per block, O(β^β) worst case

closest_vector(B: Matrix[int], target: Vector[int]) -> Vector[int]

[source]

Find the (approximate) closest vector to the target in the lattice with basis given by rows of matrix B.

Uses Babai nearest-plane algorithm for approximate closest vector.

Complexity: O(n²d) time, O(n² + nd) space for n × d matrix (n vectors of dimension d)

small_roots(f: Polynomial[int], mod: int, bounds: tuple[int, ...] | None=None, *, m: int | None=None) -> list[tuple[int, ...]]

[source]

Find small integer roots of a multivariate polynomial f(x₁, x₂, ...) ≡ 0 (mod M).

Uses the Jochemsz-May multivariate generalization of Coppersmith's method.

See: https://www.iacr.org/archive/asiacrypt2006/42840270/42840270.pdf See: https://cr.yp.to/bib/2001/howgrave-graham.pdf See: https://link.springer.com/chapter/10.1007/3-540-68339-9_14

Complexity: Lattice path is dominated by LLL on an H × W matrix, about O(H⁵W log³A) time and O(H² + HW) space, where A is the max lattice

integers() -> Iterator[int]

[source]

Generate all integers (0, 1, -1, 2, -2, ...) in an infinite generator.

integer_pairs() -> Iterator[tuple[int, int]]

[source]

Generate all integer pairs (x, y) via diagonal enumeration.

alternating(*iterables: Iterable) -> Iterator

[source]

Visit input iterables in a cycle until each is exhausted.

below(f: Callable[[int], int], upper_bound: int, start: int=0) -> Iterable[int]

[source]

Yield consecutive values of n >= start as long f(n) < upper_bound.

lower_bound(f: Callable[[int], int], f_min: int, low: int=0, high: int | None=None) -> int

[source]

Given a monotonically increasing function f, find where it first reaches f_min. Returns the smallest integer n in [low, high] such that f(n) >= f_min.

permutation(n: int, master_key: bytes | None=None) -> Iterator[int]

[source]

Generate a pseudorandom permutation of the integers 0, 1, ..., n - 1.

is_square(n: int) -> bool

[source]

Check if an integer n is a square.

iroot(x: int, n: int) -> int

[source]

Find the integer n-th root of x. Returns the largest integer a such that a^n ≤ x. Uses Newton's method.

ilog(a: int, b: int=2) -> int

[source]

Find the integer logarithm of a with base b. Returns the largest integer n such that b^n ≤ a. Uses repeated squaring and binary search.

fibonacci(n: int, mod: int | None=None) -> int

[source]

Return the n-th Fibonacci number.

fibonacci_index(n: int) -> int

[source]

Find the index of n in the Fibonacci sequence. Returns the largest integer i such that F(i) <= n.

Complexity: O(log² n) time for logarithmic search with Fibonacci evaluations

polygonal(s: int, i: int) -> int

[source]

Return the i-th s-gonal number.

polygonal_index(s: int, n: int) -> int

[source]

Find the index of n in the s-gonal numbers. Returns the largest integer i such that P(s, i) ≤ n.

periodic_continued_fraction(D: int, P: int=0, Q: int=1) -> tuple[Iterator[int], int, int]

[source]

Compute coefficients for the periodic continued fraction (P + sqrt(D)) / Q = a₀ + 1 / (a₁ + 1 / (a₂ + ...)).

Returns

convergents(coefficients: Iterable[int], num: int | None=None) -> Iterator[tuple[int, int]]

[source]

Return convergents of the continued fraction with the given coefficients.

Yields

polynomial(coefficients: Sequence[Number], mod: int | None=None) -> Callable[[Number], Number]

[source]

Create a univariate polynomial function with the given coefficients (a₀, ..., aₙ). Uses Horner's method for polynomial evaluation.