#!/usr/bin/env python3
"""
typed_connection_proof.py

Exact executable audit of the following mathematical connections:

  affine contraction
      <-> displacement scaling
      -> grid quantization
      -> integer floor dynamics
      == rational floor dynamics on a certified finite domain
      + ordered Euclidean remainders
      -> exact finite-trajectory reconstruction

and, independently but compatibly,

  positive homothety in an affine real space
      -- affine-linear scalar observable -->
  scalar affine contraction.

The file also checks three separate auxiliary branches:

  * the integer-pair automorphism with positive eigenvalue phi;
  * a finite center-relative Euclidean configuration under homothety;
  * a radial function, its negative derivative, and two complementary
    normalized radial coordinates.

Evidence levels are kept distinct:

  SYMBOLIC IDENTITY
      An exact algebraic identity checked by symbolic simplification.

  ANALYTIC ARGUMENT
      A displayed mathematical proof whose hypotheses are stated here.

  EXACT FINITE CERTIFICATE
      Exhaustive integer/rational verification on the printed domain.

  REGRESSION CHECK
      A finite test of an already proved general identity.

The enclosure of lambda uses only Python integers and fractions.Fraction.
It uses integer-square-root bounds for sqrt(5), alternating-series bounds
for arctangent, and the Machin identity for pi.  No binary floating-point
number enters the floor certificate.

Dependency:
    sympy >= 1.12

This is an executable mathematical audit accompanied by explicit proofs.  It
is not a formal proof-assistant kernel.
"""

from __future__ import annotations

from collections.abc import Callable, Iterable, Sequence
from dataclasses import dataclass
from fractions import Fraction
from itertools import combinations
from math import isqrt
from random import Random
from typing import NoReturn

try:
    import sympy as sp
except ModuleNotFoundError as exc:  # pragma: no cover - environment diagnostic
    raise SystemExit(
        "This exact audit requires SymPy. Install it with: "
        "python3 -m pip install 'sympy>=1.12'"
    ) from exc


# =============================================================================
# 0. Exact constants and declared finite domain
# =============================================================================

phi = (1 + sp.sqrt(5)) / 2
lam = sp.simplify(phi / (phi + sp.pi))
mu = sp.simplify(1 - lam)

S = 10**15
Q = 339_949_771_344_778
q = Fraction(Q, S)

DEFAULT_D_MIN = 0
DEFAULT_D_MAX = 32_767
DEFAULT_STEPS = 23


def fail(message: str) -> NoReturn:
    raise AssertionError(message)


def assert_zero(expr: sp.Expr, label: str) -> None:
    reduced = sp.simplify(expr)
    if reduced != 0:
        fail(f"{label}: expected 0, obtained {reduced}")


def assert_symbolic_equal(lhs: sp.Expr, rhs: sp.Expr, label: str) -> None:
    assert_zero(lhs - rhs, label)


assert_symbolic_equal(lam + mu, sp.Integer(1), "lambda+mu=1")
assert_symbolic_equal(mu, sp.pi / (phi + sp.pi), "mu=pi/(phi+pi)")


# =============================================================================
# 1. Affine contraction, inverse, and displacement conjugacy
# =============================================================================

# The abstract positive symbol Lambda proves the algebra once.  The rational
# enclosure in Section 5 certifies that the concrete lambda satisfies
# 0 < lambda < 1, so it is a valid contracting instance.
Lambda = sp.symbols("Lambda", positive=True, finite=True)
v, w, a, b, z = sp.symbols("v w a b z", real=True, finite=True)


def symbolic_T(x: sp.Expr) -> sp.Expr:
    return 1 + Lambda * (x - 1)


def symbolic_T_inverse(x: sp.Expr) -> sp.Expr:
    return 1 + (x - 1) / Lambda


assert_symbolic_equal(symbolic_T(1), sp.Integer(1), "T(1)=1")
assert_symbolic_equal(sp.diff(symbolic_T(v), v), Lambda, "T'(v)=Lambda")
assert_symbolic_equal(symbolic_T_inverse(symbolic_T(v)), v, "T^{-1} composed with T")
assert_symbolic_equal(symbolic_T(symbolic_T_inverse(w)), w, "T composed with T^{-1}")
assert_symbolic_equal(
    symbolic_T(a) - symbolic_T(b),
    Lambda * (a - b),
    "affine difference scaling",
)

# A(v)=v-1 and A^{-1}(z)=1+z.
assert_symbolic_equal(
    symbolic_T(1 + z) - 1,
    Lambda * z,
    "A composed with T composed with A^{-1}",
)


# =============================================================================
# 2. Continuous one-parameter group
# =============================================================================

time_t, time_s = sp.symbols("time_t time_s", real=True, finite=True)


def flow_factor(time: sp.Expr) -> sp.Expr:
    # For Lambda>0 this equals Lambda**time without branch ambiguity.
    return sp.exp(time * sp.log(Lambda))


def symbolic_T_time(time: sp.Expr, x: sp.Expr) -> sp.Expr:
    return 1 + flow_factor(time) * (x - 1)


assert_symbolic_equal(
    symbolic_T_time(time_t, symbolic_T_time(time_s, v)),
    symbolic_T_time(time_t + time_s, v),
    "T_t composed with T_s equals T_{t+s}",
)
assert_symbolic_equal(symbolic_T_time(0, v), v, "T_0 is the identity")
assert_symbolic_equal(symbolic_T_time(1, v), symbolic_T(v), "T_1=T")
assert_symbolic_equal(
    symbolic_T_time(-time_t, symbolic_T_time(time_t, v)),
    v,
    "T_{-t} is the inverse of T_t",
)

kappa = -sp.log(Lambda)
assert_symbolic_equal(
    sp.diff(symbolic_T_time(time_t, v), time_t),
    -kappa * (symbolic_T_time(time_t, v) - 1),
    "flow differential equation",
)


# =============================================================================
# 3. Exact floor arithmetic and grid quantization
# =============================================================================


def floor_fraction(x: Fraction) -> int:
    """Mathematical floor for a rational number, including negative inputs."""
    return x.numerator // x.denominator


def quantize_rational(y: Fraction, h: Fraction) -> Fraction:
    """Q_h(y)=h*floor(y/h), for exact rational y and positive rational h."""
    if h <= 0:
        raise ValueError("the quantization scale h must be positive")
    return h * floor_fraction(y / h)


def F_q(d: int) -> int:
    """F_q(d)=floor(Q*d/S), with mathematical floor for every integer d."""
    return (Q * d) // S


def rho_q(d: int) -> int:
    """Canonical Euclidean remainder Q*d=S*F_q(d)+rho_q(d)."""
    remainder = Q * d - S * F_q(d)
    if not 0 <= remainder < S:
        fail("Euclidean remainder condition failed")
    return remainder


def verify_rational_grid_quantization(
    d_min: int, d_max: int, h: Fraction = Fraction(7, 13)
) -> int:
    """Exact regression check of Q_h(q*h*d)=h*F_q(d)."""
    if d_min > d_max:
        raise ValueError("empty integer domain")
    if h <= 0:
        raise ValueError("the quantization scale h must be positive")

    for d in range(d_min, d_max + 1):
        lhs = quantize_rational(q * h * d, h)
        rhs = h * F_q(d)
        if lhs != rhs:
            fail(f"rational grid quantization failed at d={d}")
    return d_max - d_min + 1


# =============================================================================
# 4. Self-contained exact rational enclosures of sqrt(5), pi, and lambda
# =============================================================================


def sqrt5_bounds(decimal_digits: int = 90) -> tuple[Fraction, Fraction]:
    """Return strict exact rational bounds lo < sqrt(5) < hi."""
    if decimal_digits < 1:
        raise ValueError("decimal_digits must be positive")

    scale10 = 10**decimal_digits
    lower_integer = isqrt(5 * scale10 * scale10)
    lo = Fraction(lower_integer, scale10)
    hi = Fraction(lower_integer + 1, scale10)

    # Equality lo**2=5 would make sqrt(5) rational, which is impossible.
    if not lo * lo < 5 < hi * hi:
        fail("sqrt(5) enclosure failed")
    return lo, hi


def arctan_recip_bounds(m: int, terms: int) -> tuple[Fraction, Fraction]:
    r"""Exact adjacent-partial-sum bounds for arctan(1/m).

    For m>=1,

        arctan(1/m) = sum_{j>=0} (-1)^j / ((2j+1)m^(2j+1)).

    The term magnitudes decrease strictly to zero, so the alternating-series
    theorem places the sum between any partial sum and the next partial sum.
    """
    if m < 1 or terms < 1:
        raise ValueError("m and terms must be positive integers")

    x = Fraction(1, m)
    partial = Fraction(0)
    for j in range(terms):
        magnitude = x ** (2 * j + 1) / (2 * j + 1)
        partial += magnitude if j % 2 == 0 else -magnitude

    j = terms
    next_magnitude = x ** (2 * j + 1) / (2 * j + 1)
    next_partial = partial + (next_magnitude if j % 2 == 0 else -next_magnitude)
    return min(partial, next_partial), max(partial, next_partial)


def verify_machin_identity_algebra() -> None:
    r"""Verify the exact algebraic core of Machin's identity.

    Put A=arctan(1/5), B=arctan(1/239).  Exact tangent arithmetic gives

        tan(4A-B)=1.

    Also

        A > 1/5 - 1/(3*5^3),   A < 1/5,   0 < B < 1/239,

    so 0 < 4A-B < 4/5 < 3/2 < pi/2, using the classical bound pi>3.
    Tangent is injective on (-pi/2,pi/2); hence 4A-B=pi/4 and
    pi=16A-4B.
    """

    def tan_double(x: Fraction) -> Fraction:
        denominator = 1 - x * x
        if denominator == 0:
            fail("undefined tangent double-angle denominator")
        return 2 * x / denominator

    tan_A = Fraction(1, 5)
    tan_2A = tan_double(tan_A)
    tan_4A = tan_double(tan_2A)
    tan_B = Fraction(1, 239)
    denominator = 1 + tan_4A * tan_B
    tan_difference = (tan_4A - tan_B) / denominator

    if tan_2A != Fraction(5, 12):
        fail("first tangent double-angle identity failed")
    if tan_4A != Fraction(120, 119):
        fail("second tangent double-angle identity failed")
    if tan_difference != 1:
        fail("Machin tangent identity failed")

    lower_A = Fraction(1, 5) - Fraction(1, 3 * 5**3)
    upper_A = Fraction(1, 5)
    upper_B = Fraction(1, 239)
    if not 4 * lower_A - upper_B > 0:
        fail("Machin angle lower-range certificate failed")
    if not 4 * upper_A < Fraction(3, 2):
        fail("Machin angle upper-range certificate failed")


def pi_bounds(terms_5: int = 75, terms_239: int = 24) -> tuple[Fraction, Fraction]:
    r"""Return exact rational bounds for pi using Machin's identity."""
    verify_machin_identity_algebra()

    a5_lo, a5_hi = arctan_recip_bounds(5, terms_5)
    a239_lo, a239_hi = arctan_recip_bounds(239, terms_239)

    lo = 16 * a5_lo - 4 * a239_hi
    hi = 16 * a5_hi - 4 * a239_lo
    if not 3 < lo < hi < 4:
        fail("pi enclosure failed")
    return lo, hi


def lambda_bounds() -> tuple[Fraction, Fraction]:
    r"""Return exact rational bounds for lambda=phi/(phi+pi).

    The map (phi,pi) -> phi/(phi+pi) is increasing in phi and decreasing
    in pi on the positive quadrant.  Therefore

        lambda_lo = phi_lo/(phi_lo+pi_hi),
        lambda_hi = phi_hi/(phi_hi+pi_lo).
    """
    root5_lo, root5_hi = sqrt5_bounds()
    pi_lo, pi_hi = pi_bounds()

    phi_lo = (1 + root5_lo) / 2
    phi_hi = (1 + root5_hi) / 2
    lo = phi_lo / (phi_lo + pi_hi)
    hi = phi_hi / (phi_hi + pi_lo)

    if not 0 < lo < hi < 1:
        fail("lambda enclosure or contraction condition failed")
    return lo, hi


@dataclass(frozen=True)
class FloorCertificate:
    d_min: int
    d_max: int
    count: int
    lambda_lo: Fraction
    lambda_hi: Fraction
    q_is_strict_lower_bound: bool


def certified_floor_lambda(d: int, lambda_lo: Fraction, lambda_hi: Fraction) -> int:
    """Return floor(lambda*d) when the supplied enclosure proves it."""
    if not lambda_lo < lambda_hi:
        raise ValueError("invalid lambda enclosure")

    if d >= 0:
        lower = lambda_lo * d
        upper = lambda_hi * d
    else:
        lower = lambda_hi * d
        upper = lambda_lo * d

    floor_lower = floor_fraction(lower)
    floor_upper = floor_fraction(upper)
    if floor_lower != floor_upper:
        fail(f"lambda*d enclosure crosses an integer boundary at d={d}")
    return floor_lower


def certify_floor_equality(d_min: int, d_max: int) -> FloorCertificate:
    r"""Certify floor(lambda*d)=floor(Q*d/S) on every declared state."""
    if d_min > d_max:
        raise ValueError("empty integer domain")

    lambda_lo, lambda_hi = lambda_bounds()
    for d in range(d_min, d_max + 1):
        certified = certified_floor_lambda(d, lambda_lo, lambda_hi)
        rational_value = F_q(d)
        if certified != rational_value:
            fail(
                f"floor mismatch at d={d}: "
                f"certified={certified}, rational={rational_value}"
            )

    return FloorCertificate(
        d_min=d_min,
        d_max=d_max,
        count=d_max - d_min + 1,
        lambda_lo=lambda_lo,
        lambda_hi=lambda_hi,
        q_is_strict_lower_bound=q < lambda_lo,
    )


# =============================================================================
# 5. Ordered Euclidean remainders and exact finite reconstruction
# =============================================================================


@dataclass(frozen=True)
class TrajectoryRecord:
    terminal: int
    remainders: tuple[int, ...]
    states: tuple[int, ...]


def forward_with_ledger(d0: int, steps: int) -> TrajectoryRecord:
    if steps < 0:
        raise ValueError("steps must be nonnegative")

    current = d0
    remainders: list[int] = []
    states = [current]

    for _ in range(steps):
        next_state = F_q(current)
        remainder = Q * current - S * next_state
        if not 0 <= remainder < S:
            fail("forward step produced a noncanonical remainder")
        remainders.append(remainder)
        current = next_state
        states.append(current)

    return TrajectoryRecord(
        terminal=current,
        remainders=tuple(remainders),
        states=tuple(states),
    )


def reconstruct_from_terminal(
    terminal: int, remainders: Sequence[int]
) -> tuple[int, ...]:
    current = terminal
    reverse_states = [current]

    for reverse_index, remainder in enumerate(reversed(remainders)):
        if not isinstance(remainder, int) or isinstance(remainder, bool):
            raise TypeError("every remainder must be an integer")
        if not 0 <= remainder < S:
            raise ValueError(
                f"remainder outside [0,S) at reverse index {reverse_index}"
            )

        numerator = S * current + remainder
        if numerator % Q != 0:
            raise ValueError(
                f"remainder has no integer predecessor at reverse index {reverse_index}"
            )

        predecessor = numerator // Q
        if F_q(predecessor) != current or rho_q(predecessor) != remainder:
            fail("backward step failed the forward-map consistency check")

        current = predecessor
        reverse_states.append(current)

    return tuple(reversed(reverse_states))


def verify_reconstruction(seeds: Iterable[int], steps: int = DEFAULT_STEPS) -> int:
    count = 0
    for seed in seeds:
        record = forward_with_ledger(seed, steps)
        rebuilt = reconstruct_from_terminal(record.terminal, record.remainders)
        if rebuilt != record.states:
            fail(f"reconstruction failed for seed {seed}")
        count += 1
    return count


# =============================================================================
# 6. General quotient-residual encoding
# =============================================================================


@dataclass(frozen=True)
class Encoded:
    d: tuple[int, ...]
    rho: tuple[Fraction, ...]

    def __post_init__(self) -> None:
        if len(self.d) != len(self.rho):
            raise ValueError("encoded quotient and residual dimensions differ")


def encode(y: Sequence[Fraction], scales: Sequence[Fraction]) -> Encoded:
    if len(y) != len(scales):
        raise ValueError("dimension mismatch")

    quotients: list[int] = []
    residuals: list[Fraction] = []
    for coordinate, scale_i in zip(y, scales):
        if scale_i <= 0:
            raise ValueError("every scale must be positive")

        quotient = floor_fraction(coordinate / scale_i)
        residual = coordinate - scale_i * quotient
        if not Fraction(0) <= residual < scale_i:
            fail("residual outside its canonical half-open cell")

        quotients.append(quotient)
        residuals.append(residual)

    return Encoded(tuple(quotients), tuple(residuals))


def decode(e: Encoded, scales: Sequence[Fraction]) -> tuple[Fraction, ...]:
    if len(scales) != len(e.d) or len(e.d) != len(e.rho):
        raise ValueError("dimension mismatch")

    reconstructed: list[Fraction] = []
    for scale_i, quotient, residual in zip(scales, e.d, e.rho):
        if scale_i <= 0:
            raise ValueError("every scale must be positive")
        if not Fraction(0) <= residual < scale_i:
            raise ValueError("residual outside its canonical half-open cell")
        reconstructed.append(scale_i * quotient + residual)
    return tuple(reconstructed)


def verify_general_encoding(case_count: int = 1_000) -> int:
    """Deterministic exact rational regression checks, including negatives."""
    if case_count < 0:
        raise ValueError("case_count must be nonnegative")

    rng = Random(20_260_909)
    for _ in range(case_count):
        dimension = rng.randint(1, 7)
        y = tuple(
            Fraction(rng.randint(-100_000, 100_000), rng.randint(1, 997))
            for _ in range(dimension)
        )
        scales = tuple(
            Fraction(rng.randint(1, 10_000), rng.randint(1, 997))
            for _ in range(dimension)
        )
        encoded = encode(y, scales)
        if decode(encoded, scales) != y:
            fail("quotient-residual reconstruction failed")
    return case_count


# =============================================================================
# 7. Integer-pair automorphism and its phi eigenline
# =============================================================================

M = sp.Matrix([[0, 1], [1, 1]])
M_inverse = sp.Matrix([[-1, 1], [1, 0]])

if M.det() != -1:
    fail("det(M) must equal -1")
if M * M_inverse != sp.eye(2) or M_inverse * M != sp.eye(2):
    fail("integer matrix inverse failed")

characteristic_variable = sp.symbols("characteristic_variable")
assert_symbolic_equal(
    M.charpoly(characteristic_variable).as_expr(),
    characteristic_variable**2 - characteristic_variable - 1,
    "characteristic polynomial of M",
)
assert_symbolic_equal(phi**2, phi + 1, "golden quadratic identity")

phi_eigenvector = sp.Matrix([1, phi])
for component in M * phi_eigenvector - phi * phi_eigenvector:
    assert_zero(component, "phi eigenpair")


# =============================================================================
# 8. Arbitrary-dimensional homothety and scalar semiconjugacy
# =============================================================================

# Universal analytic proof:
# Let X be an affine real space with translation space V, alpha in X, and
# L:V->R linear.  For
#
#     H_t(x)=alpha+lambda^t(x-alpha),
#     v_L(x)=1+L(x-alpha),
#     T_t(v)=1+lambda^t(v-1),
#
# linearity gives
#
#     v_L(H_t(x))
#       =1+L(lambda^t(x-alpha))
#       =1+lambda^t L(x-alpha)
#       =T_t(v_L(x)).
#
# The following symbolic component calculation is an executable check of the
# same identity in seven arbitrary coordinates.  The written calculation above
# supplies the dimension-independent proof.

homothety_dimension = 7
alpha_symbols = sp.symbols(f"alpha_0:{homothety_dimension}", real=True, finite=True)
x_symbols = sp.symbols(f"x_0:{homothety_dimension}", real=True, finite=True)
functional_symbols = sp.symbols(
    f"functional_0:{homothety_dimension}", real=True, finite=True
)

scalar_displacement = sum(
    coefficient * (coordinate - center)
    for coefficient, coordinate, center in zip(
        functional_symbols, x_symbols, alpha_symbols
    )
)
homothetic_scalar_displacement = flow_factor(time_t) * scalar_displacement

assert_symbolic_equal(
    1 + homothetic_scalar_displacement,
    symbolic_T_time(time_t, 1 + scalar_displacement),
    "homothety/scalar semiconjugacy",
)


# =============================================================================
# 9. Finite center-relative Euclidean configuration under homothety
# =============================================================================

epsilon = sp.pi / 648_000
delta = epsilon / 2
sin_delta = sp.sin(delta)
cos_delta = sp.cos(delta)

x0 = sin_delta
x1 = 1 + sin_delta
x2 = phi**2 + sin_delta

# Coordinates in an orthonormal two-plane, relative to the distinguished
# center.  The six endpoints themselves lie on the affine line x=cos(delta),
# while the six center-relative vectors span the two-plane.
points = tuple(
    sp.Matrix([cos_delta, sign * offset]) for offset in (x0, x1, x2) for sign in (-1, 1)
)

affine_difference_matrix = sp.Matrix.hstack(
    *(point - points[0] for point in points[1:])
)
center_relative_matrix = sp.Matrix.hstack(*points)
gram_matrix = sp.simplify(center_relative_matrix.T * center_relative_matrix)

if affine_difference_matrix.rank() != 1:
    fail("the six endpoints must have affine dimension one")
if center_relative_matrix.rank() != 2:
    fail("the center-relative vectors must span dimension two")
if gram_matrix.rank() != 2:
    fail("the center-relative Gram matrix must have rank two")

ordinate_vector = sp.Matrix([point[1] for point in points])
expected_gram = (
    cos_delta**2 * sp.ones(len(points), len(points))
    + ordinate_vector * ordinate_vector.T
)
for entry in gram_matrix - expected_gram:
    assert_zero(entry, "center-relative Gram formula")

positive_scale = sp.symbols("positive_scale", positive=True, finite=True)
scaled_points = tuple(positive_scale * point for point in points)

pair_count = 0
for i, j in combinations(range(len(points)), 2):
    original_difference = points[i] - points[j]
    scaled_difference = scaled_points[i] - scaled_points[j]
    original_squared_distance = original_difference.dot(original_difference)
    scaled_squared_distance = scaled_difference.dot(scaled_difference)
    assert_symbolic_equal(
        scaled_squared_distance,
        positive_scale**2 * original_squared_distance,
        f"squared-distance scaling for pair ({i},{j})",
    )
    pair_count += 1

if pair_count != 15:
    fail("all fifteen unordered point pairs were not checked")

scaled_center_relative_matrix = positive_scale * center_relative_matrix
scaled_gram = sp.simplify(
    scaled_center_relative_matrix.T * scaled_center_relative_matrix
)
for entry in scaled_gram - positive_scale**2 * gram_matrix:
    assert_zero(entry, "Gram scaling")


# =============================================================================
# 10. Radial function and complementary normalized coordinates
# =============================================================================

r, ell, lambda_star, coupling = sp.symbols(
    "r ell lambda_star coupling", positive=True, finite=True
)

V = sp.exp(-lambda_star / r) / r
expected_V_derivative = (lambda_star - r) * sp.exp(-lambda_star / r) / r**3
assert_symbolic_equal(sp.diff(V, r), expected_V_derivative, "radial derivative")

U = coupling * V
expected_negative_U_derivative = (
    coupling / r**2 * (1 - lambda_star / r) * sp.exp(-lambda_star / r)
)
assert_symbolic_equal(
    -sp.diff(U, r),
    expected_negative_U_derivative,
    "negative radial derivative",
)

assert_symbolic_equal(
    expected_V_derivative.subs(r, lambda_star),
    sp.Integer(0),
    "critical radius",
)
assert_symbolic_equal(
    V.subs(r, lambda_star),
    1 / (sp.E * lambda_star),
    "radial maximum value",
)
assert_symbolic_equal(
    sp.diff(V, r, 2).subs(r, lambda_star),
    -1 / (sp.E * lambda_star**3),
    "negative second derivative at the critical radius",
)
assert_symbolic_equal(sp.limit(V, r, 0, dir="+"), sp.Integer(0), "radial limit at zero")
assert_symbolic_equal(sp.limit(V, r, sp.oo), sp.Integer(0), "radial limit at infinity")

# Since exp(-lambda_star/r)/r**3 is positive for r>0, the derivative has
# the sign of lambda_star-r.  Hence V increases on (0,lambda_star), decreases
# on (lambda_star,infinity), and the critical point is the unique global
# maximum.  This is the analytic sign proof; the identities and endpoint
# limits used by it are checked above.

v_minus = ell / (r + ell)
v_plus = r / (r + ell)
assert_symbolic_equal(v_minus + v_plus, sp.Integer(1), "v_-+v_+=1")

Phi_minus = sp.simplify(mu + lam * v_minus)
Phi_plus = sp.simplify(mu + lam * v_plus)

assert_symbolic_equal(
    Phi_minus,
    1 + lam * (v_minus - 1),
    "Phi_- is the affine contraction evaluated at v_-",
)
assert_symbolic_equal(
    Phi_plus,
    1 + lam * (v_plus - 1),
    "Phi_+ is the affine contraction evaluated at v_+",
)
assert_symbolic_equal(Phi_minus + Phi_plus, 1 + mu, "Phi_-+Phi_+=1+mu")
assert_symbolic_equal(v_minus.subs(r, ell), sp.Rational(1, 2), "v_-(ell)=1/2")
assert_symbolic_equal(v_plus.subs(r, ell), sp.Rational(1, 2), "v_+(ell)=1/2")
assert_symbolic_equal(sp.limit(Phi_minus, r, 0, dir="+"), sp.Integer(1), "Phi_-(0+)")
assert_symbolic_equal(sp.limit(Phi_minus, r, sp.oo), mu, "Phi_-(infinity)")
assert_symbolic_equal(sp.limit(Phi_plus, r, 0, dir="+"), mu, "Phi_+(0+)")
assert_symbolic_equal(sp.limit(Phi_plus, r, sp.oo), sp.Integer(1), "Phi_+(infinity)")

# A separately declared multiplier changes the derivative by the product rule.
radial_multiplier = sp.Function("radial_multiplier")
display_function = V * radial_multiplier(r)
expected_display_negative_derivative = expected_negative_U_derivative.subs(
    coupling, 1
) * radial_multiplier(r) - V * sp.diff(radial_multiplier(r), r)
assert_symbolic_equal(
    -sp.diff(display_function, r),
    expected_display_negative_derivative,
    "radial product-rule separation",
)


# =============================================================================
# 11. Negative validation tests
# =============================================================================


def expect_exception(
    expected_type: type[Exception], operation: Callable[[], object], label: str
) -> None:
    try:
        operation()
    except expected_type:
        return
    fail(f"{label}: expected {expected_type.__name__}, but no exception occurred")


def verify_negative_validation_paths() -> int:
    tests: tuple[tuple[type[Exception], Callable[[], object], str], ...] = (
        (
            ValueError,
            lambda: quantize_rational(Fraction(1), Fraction(0)),
            "zero quantization scale",
        ),
        (
            ValueError,
            lambda: encode((Fraction(1),), (Fraction(1), Fraction(2))),
            "encoder dimension mismatch",
        ),
        (
            ValueError,
            lambda: decode(Encoded((1, 2), (Fraction(0), Fraction(0))), (Fraction(1),)),
            "decoder dimension mismatch",
        ),
        (
            ValueError,
            lambda: decode(Encoded((0,), (Fraction(1),)), (Fraction(1),)),
            "decoder noncanonical residual",
        ),
        (
            ValueError,
            lambda: reconstruct_from_terminal(0, (S,)),
            "out-of-range ledger remainder",
        ),
        (
            ValueError,
            lambda: reconstruct_from_terminal(0, (1,)),
            "nonintegral ledger predecessor",
        ),
        (
            ValueError,
            lambda: forward_with_ledger(0, -1),
            "negative trajectory length",
        ),
    )

    for expected_type, operation, label in tests:
        expect_exception(expected_type, operation, label)
    return len(tests)


# =============================================================================
# 12. Exact theorem statement
# =============================================================================

END_TO_END_THEOREM = r"""
THEOREM (typed contraction-to-ledger connection)

Let X be an affine real space with translation space V, let alpha belong to X,
and let L:V->R be linear.  Let 0<lambda<1 and define, for real t,

    H_t(x) = alpha + lambda^t (x-alpha),
    v_L(x) = 1 + L(x-alpha),
    T_t(v) = 1 + lambda^t(v-1).

Then

    v_L composed with H_t = T_t composed with v_L.

At t=1, suppose L(x-alpha)=h*d, where h>0 and d is an integer.  Then

    Q_h(L(H_1(x)-alpha))
      = h floor(lambda*d)
      = h F_lambda(d).

On any declared integer domain D on which

    F_lambda(d)=F_q(d),

the same quantized state equals h*F_q(d).

For a finite F_q-trajectory d_0,...,d_n, define the ordered canonical
remainders

    rho_j = Q*d_j - S*d_{j+1},       0 <= rho_j < S.

Every predecessor is then uniquely recovered by

    d_j = (S*d_{j+1}+rho_j)/Q.

Consequently, the terminal integer d_n together with the ordered remainder
sequence (rho_0,...,rho_{n-1}) reconstructs the entire finite trajectory
exactly.  The statement requires no geometric or application-specific
interpretation.
"""


# =============================================================================
# 13. Full audit
# =============================================================================


def main() -> None:
    print("EXACT SYMBOLIC IDENTITIES: PASS")
    print("  affine fixed point, derivative, inverse, and difference scaling")
    print("  displacement conjugacy and real one-parameter group law")
    print("  flow differential equation")
    print("  integer-pair automorphism, characteristic polynomial, and eigenpair")
    print("  affine-homothety/scalar semiconjugacy")
    print("  all 15 marked-point distance identities and Gram scaling")
    print("  affine dimension 1; center-relative and Gram ranks 2")
    print("  radial derivative, limits, unique-maximum ingredients, and product rule")
    print("  complementary radial normalization identities and limits")

    lambda_lo, lambda_hi = lambda_bounds()
    if not 0 < lambda_lo < lambda_hi < 1:
        fail("the concrete lambda was not certified to lie in (0,1)")
    print("EXACT CONSTANT ENCLOSURES: PASS")
    print("  sqrt(5), Machin identity, pi, and 0<lambda<1")

    encoding_count = verify_general_encoding(case_count=1_000)
    print(
        "QUOTIENT-RESIDUAL REGRESSION: PASS "
        f"({encoding_count} deterministic exact rational cases)"
    )

    grid_count = verify_rational_grid_quantization(DEFAULT_D_MIN, DEFAULT_D_MAX)
    print(f"RATIONAL GRID-QUANTIZATION REGRESSION: PASS ({grid_count} states)")

    floor_certificate = certify_floor_equality(DEFAULT_D_MIN, DEFAULT_D_MAX)
    if not floor_certificate.q_is_strict_lower_bound:
        fail("q<lambda was not certified")
    print(
        "FINITE FLOOR CERTIFICATE: PASS for every integer "
        f"d in [{floor_certificate.d_min},{floor_certificate.d_max}] "
        f"({floor_certificate.count} states; q<lambda certified)"
    )

    trajectory_count = verify_reconstruction(
        range(DEFAULT_D_MIN, DEFAULT_D_MAX + 1),
        steps=DEFAULT_STEPS,
    )
    print(
        "EXHAUSTIVE ORDERED-REMAINDER RECONSTRUCTION: PASS "
        f"({trajectory_count} seeds x {DEFAULT_STEPS} steps)"
    )

    negative_count = verify_negative_validation_paths()
    print(
        "NEGATIVE VALIDATION PATHS: PASS "
        f"({negative_count} malformed-input classes rejected)"
    )

    print()
    print(END_TO_END_THEOREM)


if __name__ == "__main__":
    main()