fix(math): preflight exact result growth

Author: morluto

src/jacobian/math/finite_stochastic_processes/_models.py

@@ -53,6 +53,12 @@ def require_same_space(self) -> Self:
             raise ValueError(
                 "random variable and sigma algebra must share the same probability space"
             )
+        for value in self.rv.values:
+            require_bounded_rational(
+                value,
+                max_digits=256,
+                label="random-variable value",
+            )
         return self
 
 

src/jacobian/math/finite_stochastic_processes/values.py

@@ -62,12 +62,6 @@ class FiniteRandomVariable(StrictModel):
     def require_valid_rv(self) -> Self:
         if len(self.values) != len(self.space.samples):
             raise ValueError("values must have one entry per sample")
-        for value in self.values:
-            require_bounded_rational(
-                value,
-                max_digits=256,
-                label="random-variable value",
-            )
         return self
 
 

src/jacobian/math/galois_theory/_models.py

@@ -162,6 +162,8 @@ def require_reconstruction_certificate(self) -> Self:
 
 
 def _require_factor_residues(result: GaloisFactorResult) -> None:
+    from sympy import GF, Poly, Symbol
+
     prime = result.field_order
     if result.unit >= prime:
         raise ValueError("factorization unit must be a canonical nonzero residue")
@@ -172,6 +174,13 @@ def _require_factor_residues(result: GaloisFactorResult) -> None:
             raise ValueError("factor coefficients must be canonical field residues")
         if factor.coefficients[-1] != 1:
             raise ValueError("finite-field factors must be monic")
+        polynomial = Poly(
+            list(reversed(factor.coefficients)),
+            Symbol("x"),
+            domain=GF(prime),
+        )
+        if not polynomial.is_irreducible:
+            raise ValueError("every finite-field factor must be irreducible")
 
 
 def _reconstruct_factorization(result: GaloisFactorResult) -> tuple[int, ...]:

src/jacobian/math/matrices/symbolic/_models.py

@@ -89,6 +89,35 @@ def require_square(self) -> Self:
         return self
 
 
+class SymbolicDeterminantRequest(SquareSymbolicMatrixRequest):
+    """A square matrix whose exact determinant fits the public result type."""
+
+    @model_validator(mode="after")
+    def require_representable_determinant(self) -> Self:
+        # The input sparsity budget bounds backend work, but rational-function
+        # expansion can still exceed the canonical result representation.  Run
+        # the bounded exact kernel at admission so an accepted request cannot
+        # fail later while constructing its typed result.
+        from jacobian.math.matrices.symbolic import symbolic_determinant
+
+        symbolic_determinant(self.matrix.entries, self.matrix.variables)
+        return self
+
+
+class SymbolicCharacteristicPolynomialRequest(SquareSymbolicMatrixRequest):
+    """A square matrix whose characteristic polynomial fits the result type."""
+
+    @model_validator(mode="after")
+    def require_representable_characteristic_polynomial(self) -> Self:
+        from jacobian.math.matrices.symbolic import symbolic_characteristic_polynomial
+
+        symbolic_characteristic_polynomial(
+            self.matrix.entries,
+            self.matrix.variables,
+        )
+        return self
+
+
 class SymbolicDeterminantResult(StrictModel):
     """The exact determinant in the matrix's rational-function field."""
 

src/jacobian/math/matrices/symbolic/_operations.py

@@ -11,8 +11,9 @@
     symbolic_rank,
 )
 from jacobian.math.matrices.symbolic._models import (
-    SquareSymbolicMatrixRequest,
+    SymbolicCharacteristicPolynomialRequest,
     SymbolicCharacteristicPolynomialResult,
+    SymbolicDeterminantRequest,
     SymbolicDeterminantResult,
     SymbolicEigenvaluesResult,
     SymbolicMatrixRequest,
@@ -21,7 +22,7 @@
 
 
 def compute_symbolic_determinant(
-    request: SquareSymbolicMatrixRequest,
+    request: SymbolicDeterminantRequest,
 ) -> SymbolicDeterminantResult:
     determinant = symbolic_determinant(
         request.matrix.entries,
@@ -41,7 +42,7 @@ def compute_symbolic_rank(
 
 
 def compute_symbolic_characteristic_polynomial(
-    request: SquareSymbolicMatrixRequest,
+    request: SymbolicCharacteristicPolynomialRequest,
 ) -> SymbolicCharacteristicPolynomialResult:
     degree, coeffs = symbolic_characteristic_polynomial(
         request.matrix.entries,
@@ -54,7 +55,7 @@ def compute_symbolic_characteristic_polynomial(
 
 
 def compute_symbolic_eigenvalues(
-    request: SquareSymbolicMatrixRequest,
+    request: SymbolicCharacteristicPolynomialRequest,
 ) -> SymbolicEigenvaluesResult:
     entries = request.matrix.entries
     variables = request.matrix.variables

src/jacobian/math/matrices/symbolic/_tools.py

@@ -7,8 +7,9 @@
 from jacobian.catalog._examples import example
 from jacobian.catalog.models import MathTool, OperationExample
 from jacobian.math.matrices.symbolic._models import (
-    SquareSymbolicMatrixRequest,
+    SymbolicCharacteristicPolynomialRequest,
     SymbolicCharacteristicPolynomialResult,
+    SymbolicDeterminantRequest,
     SymbolicDeterminantResult,
     SymbolicEigenvaluesResult,
     SymbolicMatrixRequest,
@@ -97,7 +98,7 @@ def _generic_two_by_two() -> dict[str, Any]:
         "matrix.symbolic.determinant.compute",
         "Compute an exact symbolic matrix determinant (det) over QQ(t_1, ..., t_n)",
         "Compute the determinant of a square matrix whose entries are rational functions in declared algebraically independent variables, using SymPy's exact fraction-free Bareiss algorithm.",
-        SquareSymbolicMatrixRequest,
+        SymbolicDeterminantRequest,
         SymbolicDeterminantResult,
         compute_symbolic_determinant,
         "matrix",
@@ -142,7 +143,7 @@ def _generic_two_by_two() -> dict[str, Any]:
         "matrix.symbolic.characteristic_polynomial.compute",
         "Compute an exact symbolic characteristic polynomial",
         "Compute the dense monic coefficients of det(lambda I - A) for a square symbolic matrix whose entries are rational functions in declared algebraically independent variables.",
-        SquareSymbolicMatrixRequest,
+        SymbolicCharacteristicPolynomialRequest,
         SymbolicCharacteristicPolynomialResult,
         compute_symbolic_characteristic_polynomial,
         "matrix",
@@ -165,7 +166,7 @@ def _generic_two_by_two() -> dict[str, Any]:
         "matrix.symbolic.eigenvalues.compute",
         "Compute exact symbolic eigenvalues",
         "Compute the exact eigenvalues with algebraic multiplicities of a square symbolic matrix using SymPy's eigenvals. Entries may be rational functions in declared algebraically independent variables; eigenvalues are returned as canonical SymPy expression strings.",
-        SquareSymbolicMatrixRequest,
+        SymbolicCharacteristicPolynomialRequest,
         SymbolicEigenvaluesResult,
         compute_symbolic_eigenvalues,
         "matrix",

src/jacobian/math/polynomial_interpolation_ops/_kernel.py

@@ -0,0 +1,45 @@
+"""Shared exact kernels for polynomial interpolation contracts and operations."""
+
+from __future__ import annotations
+
+from fractions import Fraction
+
+from jacobian._exact import CanonicalRational
+
+
+def divided_difference_coefficients(
+    nodes: tuple[CanonicalRational, ...],
+    values: tuple[CanonicalRational, ...],
+) -> tuple[Fraction, ...]:
+    """Return the exact Newton coefficients for pairwise-distinct nodes."""
+
+    node_values = tuple(node.as_fraction() for node in nodes)
+    row = [value.as_fraction() for value in values]
+    coefficients = [row[0]]
+    for width in range(1, len(node_values)):
+        row = [
+            (row[index + 1] - row[index])
+            / (node_values[index + width] - node_values[index])
+            for index in range(len(node_values) - width)
+        ]
+        coefficients.append(row[0])
+    return tuple(coefficients)
+
+
+def evaluate_newton_form(
+    nodes: tuple[CanonicalRational, ...],
+    coefficients: tuple[CanonicalRational, ...],
+    point: CanonicalRational,
+) -> Fraction:
+    """Evaluate one Newton form exactly with nested multiplication."""
+
+    node_values = tuple(node.as_fraction() for node in nodes)
+    coefficient_values = tuple(value.as_fraction() for value in coefficients)
+    point_value = point.as_fraction()
+    result = coefficient_values[-1]
+    for index in range(len(coefficient_values) - 2, -1, -1):
+        result = coefficient_values[index] + (point_value - node_values[index]) * result
+    return result
+
+
+__all__ = ["divided_difference_coefficients", "evaluate_newton_form"]

src/jacobian/math/polynomial_interpolation_ops/_models.py

@@ -6,8 +6,16 @@
 
 from pydantic import Field, model_validator
 
-from jacobian._exact import CanonicalRational, require_bounded_rational
+from jacobian._exact import (
+    MAX_CANONICAL_RATIONAL_DIGITS,
+    CanonicalRational,
+    require_bounded_rational,
+)
 from jacobian._models import StrictModel
+from jacobian.math.polynomial_interpolation_ops._kernel import (
+    divided_difference_coefficients,
+    evaluate_newton_form,
+)
 
 MAX_POINTS = 32
 _MAX_RATIONAL_DIGITS = 256
@@ -46,6 +54,12 @@ def require_well_formed_samples(self) -> Self:
         _require_distinct(self.nodes)
         _require_bounded(self.nodes, "interpolation node")
         _require_bounded(self.values, "interpolation value")
+        for coefficient in divided_difference_coefficients(self.nodes, self.values):
+            require_bounded_rational(
+                CanonicalRational.from_fraction(coefficient),
+                max_digits=MAX_CANONICAL_RATIONAL_DIGITS,
+                label="derived Newton coefficient",
+            )
         return self
 
 
@@ -75,7 +89,6 @@ def require_basis_shape(self) -> Self:
             raise ValueError("Newton nodes and coefficients must have the same length")
         _require_distinct(self.nodes)
         _require_bounded(self.nodes, "Newton node")
-        _require_bounded(self.coefficients, "Newton coefficient")
         return self
 
 
@@ -90,6 +103,17 @@ def require_bounded_point(self) -> Self:
             max_digits=_MAX_RATIONAL_DIGITS,
             label="evaluation point",
         )
+        require_bounded_rational(
+            CanonicalRational.from_fraction(
+                evaluate_newton_form(
+                    self.newton_form.nodes,
+                    self.newton_form.coefficients,
+                    self.evaluation_point,
+                )
+            ),
+            max_digits=MAX_CANONICAL_RATIONAL_DIGITS,
+            label="derived Newton evaluation",
+        )
         return self
 
 

src/jacobian/math/polynomial_interpolation_ops/_operations.py

@@ -5,6 +5,10 @@
 from fractions import Fraction
 
 from jacobian._exact import CanonicalRational
+from jacobian.math.polynomial_interpolation_ops._kernel import (
+    divided_difference_coefficients,
+    evaluate_newton_form,
+)
 from jacobian.math.polynomial_interpolation_ops._models import (
     DividedDifferencesRequest,
     DividedDifferencesResult,
@@ -15,39 +19,22 @@
 )
 
 
-def _divided_difference_coefficients(
-    nodes: tuple[CanonicalRational, ...],
-    values: tuple[CanonicalRational, ...],
-) -> tuple[Fraction, ...]:
-    node_values = tuple(node.as_fraction() for node in nodes)
-    row = [value.as_fraction() for value in values]
-    coefficients = [row[0]]
-    for width in range(1, len(node_values)):
-        row = [
-            (row[index + 1] - row[index])
-            / (node_values[index + width] - node_values[index])
-            for index in range(len(node_values) - width)
-        ]
-        coefficients.append(row[0])
-    return tuple(coefficients)
-
-
 def _canonical(values: tuple[Fraction, ...]) -> tuple[CanonicalRational, ...]:
     return tuple(CanonicalRational.from_fraction(value) for value in values)
 
 
 def compute_divided_differences(
     request: DividedDifferencesRequest,
 ) -> DividedDifferencesResult:
-    coefficients = _divided_difference_coefficients(
+    coefficients = divided_difference_coefficients(
         request.samples.nodes,
         request.samples.values,
     )
     return DividedDifferencesResult(coefficients=_canonical(coefficients))
 
 
 def compute_newton_form(request: NewtonFormRequest) -> NewtonForm:
-    coefficients = _divided_difference_coefficients(
+    coefficients = divided_difference_coefficients(
         request.samples.nodes,
         request.samples.values,
     )
@@ -58,15 +45,15 @@ def compute_newton_form(request: NewtonFormRequest) -> NewtonForm:
 
 
 def compute_newton_evaluate(request: NewtonEvaluateRequest) -> NewtonEvaluateResult:
-    nodes = tuple(node.as_fraction() for node in request.newton_form.nodes)
-    coefficients = tuple(
-        coefficient.as_fraction() for coefficient in request.newton_form.coefficients
+    return NewtonEvaluateResult(
+        result=CanonicalRational.from_fraction(
+            evaluate_newton_form(
+                request.newton_form.nodes,
+                request.newton_form.coefficients,
+                request.evaluation_point,
+            )
+        )
     )
-    point = request.evaluation_point.as_fraction()
-    result = coefficients[-1]
-    for index in range(len(coefficients) - 2, -1, -1):
-        result = coefficients[index] + (point - nodes[index]) * result
-    return NewtonEvaluateResult(result=CanonicalRational.from_fraction(result))
 
 
 __all__ = [

tests/math/finite_stochastic_processes/test_finite_stochastic_processes.py

@@ -123,6 +123,20 @@ def test_discrete_sigma(self) -> None:
         # E[X | {H}] = 1, E[X | {T}] = 0
         assert result.values == (_q(1), _q(0))
 
+    def test_exact_result_may_grow_beyond_input_digit_bound(self) -> None:
+        left = 10**255 + 19
+        right = 10**255 + 21
+        rv = FiniteRandomVariable(
+            space=_coin_space(), values=(_q(1, left), _q(1, right))
+        )
+        sigma = FiniteSigmaAlgebra(space=_coin_space(), blocks=(("H", "T"),))
+
+        result = compute_conditional_expectation(
+            ConditionalExpectationRequest(rv=rv, sigma=sigma)
+        )
+
+        assert len(result.values[0].den) > 256
+
 
 # ---------------------------------------------------------------------------
 # Filtration