#!/usr/bin/env python3
"""OMEGA NAUGHT: exact symbolic certificates and finite audits.

Python 3.10+, standard library only. No floating-point arithmetic or assert
statements. General theorems are proved in the accompanying Markdown;
finite audits are reported with their actual domains.
"""
from __future__ import annotations

from dataclasses import dataclass
from fractions import Fraction as F
from math import isqrt
from typing import Callable, Iterable
import sys

Q = 339_949_771_344_778
S = 1_000_000_000_000_000
MAX_SEED = 32_767
HORIZON = 23
MAX_FIBER = 50_000


def require(condition: bool, message: str) -> None:
    """Always active, including under python -O."""
    if not condition:
        raise RuntimeError(message)


def natural(value: int, name: str) -> None:
    if type(value) is not int or value < 0:
        raise ValueError(f"{name} must be a nonnegative integer")


def ceil_div(n: int, d: int) -> int:
    natural(n, "numerator")
    if type(d) is not int or d <= 0:
        raise ValueError("denominator must be a positive integer")
    return (n + d - 1) // d


def boundary(y: int) -> int:
    natural(y, "output")
    return ceil_div(S * y, Q)


@dataclass(frozen=True)
class Step:
    output: int
    remainder: int
    rank: int


def forward(d: int) -> Step:
    natural(d, "state")
    y, r = divmod(Q * d, S)
    w = d - boundary(y)
    require(0 <= r < S and 0 <= w <= 2, "forward invariant")
    return Step(y, r, w)


def inverse_remainder(y: int, r: int) -> int:
    natural(y, "output")
    natural(r, "remainder")
    if r >= S:
        raise ValueError("remainder outside [0,S)")
    d, residue = divmod(S * y + r, Q)
    if residue != 0:
        raise ValueError("record is not divisible by Q")
    if divmod(Q * d, S) != (y, r):
        raise ValueError("inconsistent remainder record")
    return d


def inverse_rank(y: int, w: int) -> int:
    natural(y, "output")
    natural(w, "rank witness")
    if w > 2:
        raise ValueError("rank witness must be 0, 1, or 2")
    d = boundary(y) + w
    if Q * d // S != y:
        raise ValueError("rank is not present in this output fiber")
    return d


def recover_history(terminal: int, witnesses: Iterable[int],
                    use_rank: bool = False) -> list[int]:
    natural(terminal, "terminal state")
    inverse = inverse_rank if use_rank else inverse_remainder
    history = [terminal]
    for witness in reversed(list(witnesses)):
        history.append(inverse(history[-1], witness))
    return list(reversed(history))


def distances(a: F, x: F, y: F) -> tuple[F, F, F]:
    return ((x + a)**2 + y*y, (x - a)**2 + y*y,
            x*x + (y - a)**2)


def recover_point(a: F, delta: tuple[F, F, F]) -> tuple[F, F]:
    if not isinstance(a, F) or any(not isinstance(z, F) for z in delta):
        raise ValueError("exact Fraction inputs required")
    if a <= 0 or len(delta) != 3 or any(z < 0 for z in delta):
        raise ValueError("positive scale and three squared distances required")
    d0, d1, d2 = delta
    x = (d0 - d1) / (4*a)
    y = (d0 + d1 - 2*d2) / (4*a)
    if distances(a, x, y) != tuple(delta):
        raise ValueError("inconsistent planar distance record")
    return x, y


def encode_message(message: bytes) -> int:
    if not isinstance(message, bytes):
        raise ValueError("message must be bytes")
    return int.from_bytes(b"\x01" + message, "big")


def decode_message(n: int) -> bytes:
    natural(n, "message integer")
    if n == 0:
        raise ValueError("zero has no leading marker")
    raw = n.to_bytes((n.bit_length() + 7) // 8, "big")
    if raw[0] != 1:
        raise ValueError("invalid message marker")
    return raw[1:]


def encode_bounded(message: bytes) -> F:
    return F(1) + F(1, encode_message(message) + 1)


def decode_bounded(a: F) -> bytes:
    if not isinstance(a, F):
        raise ValueError("exact Fraction scale required")
    if not F(1) < a <= F(3, 2):
        raise ValueError("bounded scale must lie in (1,3/2]")
    n = 1/(a - 1) - 1
    if n.denominator != 1:
        raise ValueError("scale is not a message code")
    return decode_message(n.numerator)


class Poly:
    """Sparse Z-polynomials; a monomial is a sorted string of variables.

    Example: 'aax' denotes a^2*x. This intentionally small normalizer
    checks algebraic identities, not just numerical sample evaluations.
    """
    def __init__(self, terms: dict[str, int] | None = None):
        self.terms = {m: c for m, c in (terms or {}).items() if c != 0}

    @staticmethod
    def number(n: int) -> Poly:
        return Poly({"": n})

    @staticmethod
    def variable(name: str) -> Poly:
        if len(name) != 1 or not name.isalpha():
            raise ValueError("one-letter variable required")
        return Poly({name: 1})

    @staticmethod
    def cast(value: Poly | int) -> Poly:
        return value if isinstance(value, Poly) else Poly.number(value)

    def __add__(self, other: Poly | int) -> Poly:
        other = Poly.cast(other)
        result = self.terms.copy()
        for monomial, coefficient in other.terms.items():
            result[monomial] = result.get(monomial, 0) + coefficient
        return Poly(result)

    __radd__ = __add__

    def __neg__(self) -> Poly:
        return Poly({m: -c for m, c in self.terms.items()})

    def __sub__(self, other: Poly | int) -> Poly:
        return self + (-Poly.cast(other))

    def __rsub__(self, other: Poly | int) -> Poly:
        return Poly.cast(other) + (-self)

    def __mul__(self, other: Poly | int) -> Poly:
        other = Poly.cast(other)
        result: dict[str, int] = {}
        for left, c in self.terms.items():
            for right, d in other.terms.items():
                monomial = "".join(sorted(left + right))
                result[monomial] = result.get(monomial, 0) + c*d
        return Poly(result)

    __rmul__ = __mul__

    def __pow__(self, n: int) -> Poly:
        natural(n, "polynomial exponent")
        result = Poly.number(1)
        for _ in range(n):
            result = result * self
        return result


def symbolic_certificates() -> int:
    a, x, y, u, v, p, q, d, e, h, f, c, r, l, z, j = (
        Poly.variable(name) for name in "axyuvpqdehfcrlzj")
    t = q*d - e
    nxt = q*t - j
    lam = q + h
    ds = lambda X, Y: ((X+a)**2 + Y**2, (X-a)**2 + Y**2,
                      X**2 + (Y-a)**2)
    dx = ds(x, y)
    du = ds(u, v)
    first = dx[0] - dx[1] - du[0] + du[1]
    second = dx[0] + dx[1] - 2*dx[2] - du[0] - du[1] + 2*du[2]
    certificates = [
        a*f*(f+c) - f*(a*f+a*c),
        (a*x-a*u)**2 + (a*y-a*v)**2 - a**2*((x-u)**2+(y-v)**2),
        (r+l*a) - (r+l*((r+a)-r)),
        q*d*p - (t*p+e*p),
        lam*d*p - (t*p+(e+h*d)*p),
        dx[0]-dx[1]-4*a*x,
        dx[0]+dx[1]-2*dx[2]-4*a*y,
        ((x+a)**2+y**2+z**2)-dx[0]-z**2,
        16*a**2*((x-u)**2+(y-v)**2)-first**2-second**2,
        q**2*d-nxt-q*e-j,
        lam**2*d-nxt-lam*(e+h*d)-(j+h*t),
    ]
    for index, expression in enumerate(certificates, 1):
        require(not expression.terms, f"symbolic certificate {index}")
    # A nonidentity must not normalize to zero.
    require(bool((x*x + 1 - x*x).terms), "normalizer negative control")
    return len(certificates)


def atan_sum(inverse: int, terms: int) -> F:
    return sum((F((-1)**k, (2*k+1)*inverse**(2*k+1))
                for k in range(terms)), F(0))


def lambda_bracket() -> tuple[F, F]:
    b = 10**30
    n = 5*b*b
    root = isqrt(n)
    require(root*root < n < (root+1)**2, "strict sqrt(5) bracket")
    phi_low = F(b+root, 2*b)
    phi_high = F(b+root+1, 2*b)
    pi_low = 16*atan_sum(5, 18)-4*atan_sum(239, 9)
    pi_high = 16*atan_sum(5, 19)-4*atan_sum(239, 8)
    require(F(3) < pi_low < pi_high < F(22, 7), "pi bracket")
    low = phi_low/(phi_low+pi_high)
    high = phi_high/(phi_high+pi_low)
    require(F(Q, S) < low < high < F(Q, S)+F(1, 10**14),
            "lambda bracket ordering")
    require(F(9295, 10**19) < low-F(Q, S)
            and high-F(Q, S) < F(9297, 10**19), "delta bracket")
    return low, high


def audit_history(seed: int, horizon: int) -> tuple[list[int], list[int], list[int]]:
    natural(seed, "seed")
    natural(horizon, "horizon")
    states, remainders, ranks = [seed], [], []
    q_power = s_power = 1
    weighted_sum = 0
    for _ in range(horizon):
        step = forward(states[-1])
        require(Q*states[-1] == S*step.output+step.remainder, "step identity")
        require(inverse_remainder(step.output, step.remainder) == states[-1],
                "one-step remainder inverse")
        require(inverse_rank(step.output, step.rank) == states[-1],
                "one-step rank inverse")
        weighted_sum = Q*weighted_sum + step.remainder*s_power
        q_power *= Q
        s_power *= S
        require(q_power*seed == s_power*step.output+weighted_sum,
                "full telescoping prefix identity")
        states.append(step.output)
        remainders.append(step.remainder)
        ranks.append(step.rank)
    require(recover_history(states[-1], remainders) == states,
            "complete remainder-history inverse")
    require(recover_history(states[-1], ranks, True) == states,
            "complete rank-history inverse")
    return states, remainders, ranks


def expect_rejection(action: Callable[[], object]) -> None:
    try:
        action()
    except ValueError:
        return
    raise RuntimeError("malformed input was accepted")


def run() -> dict[str, int]:
    result: dict[str, int] = {}
    result["symbolic_polynomial_identities"] = symbolic_certificates()
    require(3*Q > S and 2*Q < S, "canonical ratio range")
    low, high = lambda_bracket()
    for d in range(MAX_SEED+1):
        t = Q*d//S
        require(t*low.denominator <= low.numerator*d, "lambda lower floor")
        require(high.numerator*d < (t+1)*high.denominator, "lambda upper floor")
    result["certified_lambda_floor_states"] = MAX_SEED+1
    # A disagreement outside the certified domain; no minimality claim.
    counterexample = 42_189_065
    t = Q*counterexample//S
    require(t == 14_342_162, "counterexample rational floor")
    require((t+1)*low.denominator <= low.numerator*counterexample,
            "counterexample lower bound")
    require(high.numerator*counterexample < (t+2)*high.denominator,
            "counterexample upper bound")
    result["out_of_domain_disagreement_examples"] = 1

    predecessor_count = 0
    first_two = None
    for y in range(MAX_FIBER+1):
        lo, hi = boundary(y), boundary(y+1)
        require(hi-lo in (2, 3), "floor fiber cardinality")
        if hi-lo == 2 and first_two is None:
            first_two = y
        for d in range(lo, hi):
            step = forward(d)
            require(step.output == y, "fiber membership")
            require(inverse_rank(y, step.rank) == d, "fiber rank inverse")
            require(inverse_remainder(y, step.remainder) == d, "fiber remainder inverse")
            predecessor_count += 1
    require([forward(d).output for d in range(3)] == [0, 0, 0],
            "two-bit lower-bound example")
    result["consecutive_output_fibers"] = MAX_FIBER+1
    result["fiber_predecessors"] = predecessor_count
    result["maximum_fiber_cardinality_proved"] = 3
    result["optimal_fixed_width_witness_bits_proved"] = 2

    max_extinction = 0
    points = [(-1, 0), (1, 0), (0, 1), (3, 5)]
    for seed in range(MAX_SEED+1):
        step = forward(seed)
        for px, py in points:
            require(Q*seed*px == S*step.output*px+step.remainder*px,
                    "geometric x discrepancy")
            require(Q*seed*py == S*step.output*py+step.remainder*py,
                    "geometric y discrepancy")
        if seed > 0:
            a, x, y = F(seed), F(2*seed, 3), F(-5*seed, 7)
            require(recover_point(a, distances(a, x, y)) == (x, y),
                    "positive-scale planar recovery")
        states, _, _ = audit_history(seed, HORIZON)
        require(states[-1] == 0, "canonical extinction by horizon 23")
        max_extinction = max(max_extinction, states.index(0))
    result["geometric_scale_cases_including_zero"] = MAX_SEED+1
    result["marked_point_vector_identities"] = (MAX_SEED+1)*len(points)
    result["positive_scale_planar_recoveries"] = MAX_SEED
    result["canonical_transition_checks"] = (MAX_SEED+1)*HORIZON
    result["canonical_telescoping_prefix_checks"] = (MAX_SEED+1)*HORIZON
    result["canonical_full_history_roundtrips_per_witness_type"] = MAX_SEED+1
    result["maximum_extinction_horizon_observed"] = max_extinction

    for seed in (0, 1, MAX_SEED, 2**256+123):
        audit_history(seed, 0)
    result["zero_horizon_roundtrips"] = 4

    message_count = 0
    for length in range(3):
        for value in range(256**length):
            message = value.to_bytes(length, "big")
            require(decode_message(encode_message(message)) == message, "integer message code")
            require(decode_bounded(encode_bounded(message)) == message, "bounded-scale code")
            message_count += 1
    result["messages_of_lengths_zero_through_two_per_codec"] = message_count
    samples = [b"", b"\0", b"\0\0", b"\0\1\xff", b"OMEGA",
               bytes(range(256)), bytes(64)]
    for message in samples:
        n = encode_message(message)
        states, remainders, ranks = audit_history(n, HORIZON)
        require(decode_message(recover_history(states[-1], remainders)[0]) == message,
                "message remainder pipeline")
        require(decode_message(recover_history(states[-1], ranks, True)[0]) == message,
                "message rank pipeline")
    result["composed_message_history_roundtrips_per_witness_type"] = len(samples)

    require(first_two is not None, "two-element fiber encountered")
    rejects = [lambda: forward(-1), lambda: inverse_remainder(-1, 0),
               lambda: inverse_remainder(0, -1), lambda: inverse_remainder(0, S),
               lambda: inverse_remainder(0, 1), lambda: inverse_rank(-1, 0),
               lambda: inverse_rank(0, -1), lambda: inverse_rank(0, 3),
               lambda: inverse_rank(first_two, 2),
               lambda: recover_point(F(0), (F(1), F(1), F(1))),
               lambda: recover_point(F(1), (F(0), F(0), F(0))),
               lambda: decode_message(0), lambda: decode_message(2),
               lambda: decode_bounded(F(1)), lambda: decode_bounded(F(2))]
    for action in rejects:
        expect_rejection(action)
    result["malformed_input_cases_rejected"] = len(rejects)
    return result


if __name__ == "__main__":
    if len(sys.argv) != 1:
        raise SystemExit("Usage: python omega_unified_proof.py")
    report = run()
    print("OMEGA_UNIFIED_EXACT_CHECK_V2")
    for key, value in report.items():
        print(f"{key}={value}")
    print("status=PASS")