Add exact algebra and combinatorial algorithm examples

README.md

@@ -2,6 +2,9 @@
 
 [![Tests](https://github.com/pralav-25/python-algorithm-lab/actions/workflows/ci.yml/badge.svg)](https://github.com/pralav-25/python-algorithm-lab/actions/workflows/ci.yml)
 
+[![Python 3.11+](https://img.shields.io/badge/Python-3.11%2B-blue)](pyproject.toml)
+[![License: MIT](https://img.shields.io/badge/License-MIT-green)](LICENSE)
+
 Readable implementations of core algorithms, with executable examples and
 independent regression tests. Python 3.11+; no runtime dependencies.
 
@@ -30,7 +33,7 @@ python -m doctest algorithm_lab/*.py
 Optional installation: `python -m pip install -e .`. For linting and formatting,
 install `.[dev]` and run `ruff check .` and `ruff format --check .`.
 
-## Catalog (100 algorithms)
+## Catalog (127 algorithms)
 
 | Module | Topic | What it demonstrates |
 | --- | --- | --- |
@@ -134,6 +137,33 @@ install `.[dev]` and run `ruff check .` and `ruff format --check .`.
 | [n_queens](algorithm_lab/n_queens.py) | Backtracking | Lazy nonattacking queen placements using bit masks |
 | [fractional_knapsack](algorithm_lab/fractional_knapsack.py) | Optimization | Exact rational fractional-knapsack value and item fractions |
 | [weighted_interval_scheduling](algorithm_lab/weighted_interval_scheduling.py) | Optimization | Maximum-value compatible intervals with original-index reconstruction |
+| [modular_inverse](algorithm_lab/modular_inverse.py) | Number theory | Extended-Euclidean modular inverses |
+| [euler_totient](algorithm_lab/euler_totient.py) | Number theory | Count coprime residues by prime factorization |
+| [mobius_sieve](algorithm_lab/mobius_sieve.py) | Number theory | Linear sieve of the Mobius function |
+| [divisor_sieve](algorithm_lab/divisor_sieve.py) | Number theory | Compute every positive divisor list through a bound |
+| [segmented_sieve](algorithm_lab/segmented_sieve.py) | Number theory | Find primes in a half-open interval with segmented marking |
+| [integer_nth_root](algorithm_lab/integer_nth_root.py) | Number theory | Exact floor roots using integer binary search |
+| [continued_fraction](algorithm_lab/continued_fraction.py) | Number theory | Canonical finite continued fractions for rational inputs |
+| [continued_fraction_convergents](algorithm_lab/continued_fraction_convergents.py) | Number theory | Build exact rational convergents with a two-term recurrence |
+| [stern_brocot_path](algorithm_lab/stern_brocot_path.py) | Number theory | Run-length paths to positive rationals in the Stern-Brocot tree |
+| [multiplicative_order](algorithm_lab/multiplicative_order.py) | Number theory | Find the order of an invertible residue |
+| [polynomial_evaluate](algorithm_lab/polynomial_evaluate.py) | Algebra | Horner evaluation of integer polynomials |
+| [polynomial_multiply](algorithm_lab/polynomial_multiply.py) | Algebra | Exact dense polynomial convolution |
+| [polynomial_divmod](algorithm_lab/polynomial_divmod.py) | Algebra | Exact polynomial long division over rational coefficients |
+| [lagrange_interpolation](algorithm_lab/lagrange_interpolation.py) | Algebra | Evaluate the interpolating polynomial with exact rationals |
+| [gaussian_elimination](algorithm_lab/gaussian_elimination.py) | Algebra | Solve square rational linear systems by pivoted elimination |
+| [determinant](algorithm_lab/determinant.py) | Algebra | Compute exact determinants by elimination |
+| [matrix_power](algorithm_lab/matrix_power.py) | Algebra | Integer matrix powers by repeated squaring |
+| [linear_recurrence](algorithm_lab/linear_recurrence.py) | Algebra | Evaluate integer recurrences with a companion matrix |
+| [stirling_second](algorithm_lab/stirling_second.py) | Combinatorics | Count partitions into a specified number of nonempty blocks |
+| [combination_rank](algorithm_lab/combination_rank.py) | Combinatorics | Lexicographic ranking of fixed-size combinations |
+| [combination_unrank](algorithm_lab/combination_unrank.py) | Combinatorics | Select a lexicographic combination without enumerating predecessors |
+| [permutation_rank](algorithm_lab/permutation_rank.py) | Combinatorics | Lehmer-code ranking of permutations of range(n) |
+| [permutation_unrank](algorithm_lab/permutation_unrank.py) | Combinatorics | Factoradic selection of a lexicographic permutation |
+| [multiset_permutations](algorithm_lab/multiset_permutations.py) | Combinatorics | Generate distinct integer permutations in lexicographic order |
+| [bounded_compositions](algorithm_lab/bounded_compositions.py) | Combinatorics | Enumerate fixed-sum tuples under per-position upper bounds |
+| [subset_zeta_transform](algorithm_lab/subset_zeta_transform.py) | Combinatorics | Sum values over every subset of each bit mask |
+| [subset_mobius_transform](algorithm_lab/subset_mobius_transform.py) | Combinatorics | Invert subset sums by Mobius inversion on bit masks |
 
 ## Validation
 

algorithm_lab/bounded_compositions.py

@@ -0,0 +1,44 @@
+"""Enumerate fixed-sum tuples under per-position upper bounds.
+
+total and each bound are nonnegative integers excluding bool. Yield tuples
+in lexicographic order, with 0 <= value[i] <= bounds[i] and the given sum.
+No bounds yields () only for total=0. An explicit stack avoids recursion;
+O(k) working space, and O(k) work per yielded tuple plus visited prefixes.
+
+>>> list(bounded_compositions(2, [1, 2]))
+[(0, 2), (1, 1)]
+"""
+
+from algorithm_lab._validation import integer
+
+
+def bounded_compositions(total, bounds):
+    integer(total, "total", minimum=0)
+    bounds = list(bounds)
+    for bound in bounds:
+        integer(bound, "bound", minimum=0)
+    suffix = [0] * (len(bounds) + 1)
+    for i in range(len(bounds) - 1, -1, -1):
+        suffix[i] = suffix[i + 1] + bounds[i]
+    if not bounds:
+        if total == 0:
+            yield ()
+        return
+
+    def choices(i, remaining):
+        return iter(range(max(0, remaining - suffix[i + 1]), min(bounds[i], remaining) + 1))
+
+    prefix = []
+    stack = [(total, choices(0, total))]
+    while stack:
+        remaining, iterator = stack[-1]
+        value = next(iterator, None)
+        if value is None:
+            stack.pop()
+            if prefix:
+                prefix.pop()
+        elif len(stack) == len(bounds):
+            yield tuple(prefix + [value])
+        else:
+            prefix.append(value)
+            stack.append((remaining - value, choices(len(prefix), remaining - value)))

algorithm_lab/combination_rank.py

@@ -0,0 +1,28 @@
+"""Lexicographic ranking of fixed-size combinations.
+
+n is a nonnegative integer, and combination is a strictly increasing sequence
+of integers from range(n), excluding bool. Return its zero-based index among
+itertools.combinations(range(n), k). Empty combinations rank 0.
+O(n) binomial evaluations and O(k) copied storage.
+
+>>> combination_rank(5, [1, 3])
+5
+"""
+
+from math import comb
+
+from algorithm_lab._validation import integer
+
+
+def combination_rank(n, combination):
+    integer(n, "n", minimum=0)
+    values = list(combination)
+    previous, rank, k = -1, 0, len(values)
+    for i, value in enumerate(values):
+        integer(value, "element", minimum=0)
+        if value <= previous or value >= n:
+            raise ValueError("combination must be strictly increasing within range(n)")
+        for skipped in range(previous + 1, value):
+            rank += comb(n - skipped - 1, k - i - 1)
+        previous = value
+    return rank

algorithm_lab/combination_unrank.py

@@ -0,0 +1,30 @@
+"""Select a lexicographic combination without enumerating predecessors.
+
+Integers satisfy 0 <= k <= n and 0 <= rank < C(n,k); bool is rejected.
+Return a list matching itertools.combinations(range(n), k)[rank]. O(n)
+binomial evaluations and O(k) output space. The empty combination has rank 0.
+
+>>> combination_unrank(5, 2, 5)
+[1, 3]
+"""
+
+from math import comb
+
+from algorithm_lab._validation import integer
+
+
+def combination_unrank(n, k, rank):
+    integer(n, "n", minimum=0)
+    integer(k, "k", minimum=0)
+    integer(rank, "rank", minimum=0)
+    if k > n or rank >= comb(n, k):
+        raise ValueError("combination rank is out of range")
+    result, candidate = [], 0
+    while len(result) < k:
+        count = comb(n - candidate - 1, k - len(result) - 1)
+        if rank < count:
+            result.append(candidate)
+        else:
+            rank -= count
+        candidate += 1
+    return result

algorithm_lab/continued_fraction.py

@@ -0,0 +1,27 @@
+"""Canonical finite continued fractions for rational inputs.
+
+Accept integer numerator and nonzero integer denominator, excluding bool.
+Return [a0, a1, ...], with a0 the floor and later terms positive. The final
+term is > 1 unless the result is integral. O(log(max(|n|, |d|)+1)) Euclidean
+steps and O(number of terms) output space.
+
+>>> continued_fraction(43, 19)
+[2, 3, 1, 4]
+"""
+
+from algorithm_lab._validation import integer
+
+
+def continued_fraction(numerator, denominator):
+    integer(numerator, "numerator")
+    integer(denominator, "denominator")
+    if denominator == 0:
+        raise ValueError("denominator must be nonzero")
+    if denominator < 0:
+        numerator, denominator = -numerator, -denominator
+    result = []
+    while denominator:
+        term, remainder = divmod(numerator, denominator)
+        result.append(term)
+        numerator, denominator = denominator, remainder
+    return result

algorithm_lab/continued_fraction_convergents.py

@@ -0,0 +1,24 @@
+"""Build exact rational convergents with a two-term recurrence.
+
+terms is an iterable of integers excluding bool; only the first may be
+nonpositive. Return Fractions for all prefixes, or [] for empty input.
+O(n) arithmetic operations and O(n) output space.
+
+>>> continued_fraction_convergents([1, 2, 2])
+[Fraction(1, 1), Fraction(3, 2), Fraction(7, 5)]
+"""
+
+from fractions import Fraction
+
+from algorithm_lab._validation import integer
+
+
+def continued_fraction_convergents(terms):
+    p0, p1, q0, q1 = 0, 1, 1, 0
+    result = []
+    for i, term in enumerate(terms):
+        integer(term, "term", minimum=1 if i else None)
+        p0, p1 = p1, term * p1 + p0
+        q0, q1 = q1, term * q1 + q0
+        result.append(Fraction(p1, q1))
+    return result

algorithm_lab/determinant.py

@@ -0,0 +1,36 @@
+"""Compute exact determinants by elimination.
+
+Input is a square matrix of int/Fraction values, excluding bool/float.
+Return a Fraction; the empty determinant is 1. Singular matrices return 0.
+First nonzero row pivoting is deterministic. O(n^3) rational operations and
+O(n^2) copied storage; the input is unchanged.
+
+>>> determinant([[1, 2], [3, 4]])
+Fraction(-2, 1)
+"""
+
+from fractions import Fraction
+
+
+def determinant(matrix):
+    a = [list(row) for row in matrix]
+    n = len(a)
+    if any(len(row) != n for row in a):
+        raise ValueError("matrix must be square")
+    if any(isinstance(v, bool) or not isinstance(v, (int, Fraction)) for row in a for v in row):
+        raise ValueError("entries must be integers or Fractions")
+    a = [[Fraction(v) for v in row] for row in a]
+    result = Fraction(1)
+    for j in range(n):
+        pivot = next((i for i in range(j, n) if a[i][j]), None)
+        if pivot is None:
+            return Fraction(0)
+        if pivot != j:
+            a[j], a[pivot] = a[pivot], a[j]
+            result = -result
+        result *= a[j][j]
+        for i in range(j + 1, n):
+            factor = a[i][j] / a[j][j]
+            for k in range(j + 1, n):
+                a[i][k] -= factor * a[j][k]
+    return result

algorithm_lab/divisor_sieve.py

@@ -0,0 +1,20 @@
+"""Compute every positive divisor list through a bound.
+
+Return lists indexed by integer 0..limit, each in ascending order; index 0
+is empty. limit is a nonnegative integer excluding bool. O(n log n) time
+and output storage, where n = limit.
+
+>>> divisor_sieve(4)
+[[], [1], [1, 2], [1, 3], [1, 2, 4]]
+"""
+
+from algorithm_lab._validation import integer
+
+
+def divisor_sieve(limit):
+    integer(limit, minimum=0)
+    result = [[] for _ in range(limit + 1)]
+    for divisor in range(1, limit + 1):
+        for multiple in range(divisor, limit + 1, divisor):
+            result[multiple].append(divisor)
+    return result

algorithm_lab/euler_totient.py

@@ -0,0 +1,24 @@
+"""Count coprime residues by prime factorization.
+
+n is a positive integer, excluding bool. phi(1) is 1. Trial division takes
+O(sqrt(n)) arithmetic steps and O(1) integer storage.
+
+>>> euler_totient(36)
+12
+"""
+
+from algorithm_lab._validation import integer
+
+
+def euler_totient(n):
+    integer(n, minimum=1)
+    result, remaining, factor = n, n, 2
+    while factor * factor <= remaining:
+        if remaining % factor == 0:
+            result -= result // factor
+            while remaining % factor == 0:
+                remaining //= factor
+        factor += 1
+    if remaining > 1:
+        result -= result // remaining
+    return result

algorithm_lab/gaussian_elimination.py

@@ -0,0 +1,38 @@
+"""Solve square rational linear systems by pivoted elimination.
+
+Accept a square matrix and matching right-hand side of int/Fraction values,
+excluding bool/float. Return exact Fractions without changing inputs. Empty
+systems return []; singular systems raise ValueError. Pivots choose the first
+nonzero remaining row. O(n^3) rational operations, O(n^2) space.
+
+>>> gaussian_elimination([[2, 1], [1, -1]], [5, 1])
+[Fraction(2, 1), Fraction(1, 1)]
+"""
+
+from fractions import Fraction
+
+
+def gaussian_elimination(matrix, rhs):
+    rows, rhs = [list(row) for row in matrix], list(rhs)
+    n = len(rows)
+    if len(rhs) != n or any(len(row) != n for row in rows):
+        raise ValueError("a square matrix and matching rhs are required")
+    if any(
+        isinstance(v, bool) or not isinstance(v, (int, Fraction)) for row in rows for v in row
+    ) or any(isinstance(v, bool) or not isinstance(v, (int, Fraction)) for v in rhs):
+        raise ValueError("entries must be integers or Fractions")
+    a = [[Fraction(v) for v in row] + [Fraction(b)] for row, b in zip(rows, rhs, strict=True)]
+    for column in range(n):
+        pivot = next((i for i in range(column, n) if a[i][column]), None)
+        if pivot is None:
+            raise ValueError("matrix is singular")
+        a[column], a[pivot] = a[pivot], a[column]
+        scale = a[column][column]
+        a[column] = [value / scale for value in a[column]]
+        for i in range(column + 1, n):
+            factor = a[i][column]
+            a[i] = [v - factor * w for v, w in zip(a[i], a[column], strict=True)]
+    result = [Fraction(0)] * n
+    for i in range(n - 1, -1, -1):
+        result[i] = a[i][-1] - sum(a[i][j] * result[j] for j in range(i + 1, n))
+    return result

algorithm_lab/integer_nth_root.py

@@ -0,0 +1,28 @@
+"""Exact floor roots using integer binary search.
+
+value >= 0 and degree >= 1 are integers excluding bool. Return r satisfying
+r**degree <= value < (r+1)**degree, without floats. O(bit_length(value))
+comparisons, each using integer exponentiation; O(1) big integers stored.
+
+>>> integer_nth_root(1000, 3)
+10
+"""
+
+from algorithm_lab._validation import integer
+
+
+def integer_nth_root(value, degree):
+    integer(value, minimum=0)
+    integer(degree, "degree", minimum=1)
+    if value < 2 or degree == 1:
+        return value
+    if degree >= value.bit_length():
+        return 1
+    low, high = 0, 1 << ((value.bit_length() + degree - 1) // degree)
+    while low + 1 < high:
+        middle = (low + high) // 2
+        if middle**degree <= value:
+            low = middle
+        else:
+            high = middle
+    return high if high**degree <= value else low

algorithm_lab/lagrange_interpolation.py

@@ -0,0 +1,32 @@
+"""Evaluate the interpolating polynomial with exact rationals.
+
+Points and x contain int or Fraction values, excluding bool and float.
+Require at least one point and distinct abscissas. Return a Fraction; the
+polynomial has degree less than the number of points. O(n^2) rational
+operations and O(n) space. Point order does not affect the answer.
+
+>>> lagrange_interpolation([(0, 1), (1, 2), (2, 5)], 3)
+Fraction(10, 1)
+"""
+
+from fractions import Fraction
+
+
+def lagrange_interpolation(points, x):
+    def rational(value):
+        if isinstance(value, bool) or not isinstance(value, (int, Fraction)):
+            raise ValueError("coordinates must be integers or Fractions")
+        return Fraction(value)
+
+    x = rational(x)
+    points = [(rational(a), rational(b)) for a, b in points]
+    if not points or len({a for a, _ in points}) != len(points):
+        raise ValueError("points must be nonempty with distinct abscissas")
+    result = Fraction(0)
+    for i, (a, b) in enumerate(points):
+        term = b
+        for j, (other, _) in enumerate(points):
+            if i != j:
+                term *= (x - other) / (a - other)
+        result += term
+    return result