#!/usr/bin/env python3
"""Exact verifier for OMEGA geometry, numerical difference, and recovery.

Python 3.10+; standard library only. Run:
    python3 omega_geometry_recovery_verify.py
    python3 -O omega_geometry_recovery_verify.py

All acceptance checks use explicit exceptions. No floats or assert statements
are used. The enclosure of pi is proved using the finite geometric-series
remainder for arctangent and Machin's identity; see the accompanying proof.

The general mathematical proofs are in the proof document. This program
checks exact symbolic shape identities, supplies rational certificates, and
audits the stated finite integer domain. It does not certify physical data
or execute a proof assistant.
"""

from dataclasses import dataclass
from fractions import Fraction as F
from math import gcd, isqrt
from pathlib import Path
import ast
import hashlib
import sys

Q = 339949771344778
S = 10**15
HORIZON = 23
MAX_SEED = 32767
D_STAR = 42189065
ELL = F(339949771344778929576702809321, 10**30)
UPPER = ELL + F(1, 10**30)


def require(condition, message):
    if not condition:
        raise ValueError(message)


def nonnegative_integer(value, name):
    require(type(value) is int and value >= 0, name + " must be a nonnegative integer")


@dataclass(frozen=True)
class Golden:
    """An exact element u + v*phi with phi**2 = phi + 1."""
    u: F = F(0)
    v: F = F(0)

    def __post_init__(self):
        object.__setattr__(self, "u", F(self.u))
        object.__setattr__(self, "v", F(self.v))

    @staticmethod
    def cast(value):
        return value if isinstance(value, Golden) else Golden(F(value), F(0))

    def __add__(self, other):
        other = self.cast(other)
        return Golden(self.u + other.u, self.v + other.v)

    __radd__ = __add__

    def __neg__(self):
        return Golden(-self.u, -self.v)

    def __sub__(self, other):
        return self + (-self.cast(other))

    def __rsub__(self, other):
        return self.cast(other) - self

    def __mul__(self, other):
        other = self.cast(other)
        return Golden(
            self.u * other.u + self.v * other.v,
            self.u * other.v + self.v * other.u + self.v * other.v,
        )

    __rmul__ = __mul__


@dataclass(frozen=True)
class Poly:
    """Univariate polynomial with exact Golden coefficients."""
    coeffs: tuple

    def __post_init__(self):
        data = [Golden.cast(x) for x in self.coeffs]
        while len(data) > 1 and data[-1] == Golden():
            data.pop()
        if not data:
            data = [Golden()]
        object.__setattr__(self, "coeffs", tuple(data))

    @staticmethod
    def cast(value):
        return value if isinstance(value, Poly) else Poly((Golden.cast(value),))

    def __add__(self, other):
        other = self.cast(other)
        size = max(len(self.coeffs), len(other.coeffs))
        return Poly(tuple(
            (self.coeffs[i] if i < len(self.coeffs) else Golden())
            + (other.coeffs[i] if i < len(other.coeffs) else Golden())
            for i in range(size)
        ))

    __radd__ = __add__

    def __neg__(self):
        return Poly(tuple(-x for x in self.coeffs))

    def __sub__(self, other):
        return self + (-self.cast(other))

    def __rsub__(self, other):
        return self.cast(other) - self

    def __mul__(self, other):
        other = self.cast(other)
        data = [Golden()] * (len(self.coeffs) + len(other.coeffs) - 1)
        for i, x in enumerate(self.coeffs):
            for j, y in enumerate(other.coeffs):
                data[i+j] = data[i+j] + x*y
        return Poly(tuple(data))

    __rmul__ = __mul__

    def __pow__(self, n):
        nonnegative_integer(n, "polynomial exponent")
        answer = Poly((1,))
        factor = self
        while n:
            if n & 1:
                answer = answer * factor
            factor = factor * factor
            n //= 2
        return answer


def shape_certificate():
    phi = Golden(0, 1)
    variable = Poly((0, 1))
    A, B = variable + phi*phi, variable + 1
    z = (-A, -B, -variable, variable, B, A)
    expected = (Poly((phi,)), Poly((1,)), 2*variable, Poly((1,)), Poly((phi,)))
    require(phi*phi == phi+1, "golden defining identity")
    require(tuple(z[i+1]-z[i] for i in range(5)) == expected, "symbolic adjacent gaps")
    c2 = 1 - variable*variable
    require(c2 + variable*variable == Poly((1,)), "seed radius identity")
    gram = tuple(tuple(c2+x*y for y in z) for x in z)
    minor = gram[0][0]*gram[5][5]-gram[0][5]*gram[5][0]
    require(minor == 4*c2*A*A, "symbolic nonzero Gram minor formula")
    print("PASS 1: symbolic six-mark gaps, seed radii, and Gram minor")
    print("        Coefficient identities hold in Q(phi)[s]; positivity uses the stated geometry assumptions.")


def arctan_sum(denominator, terms):
    return sum((F((-1)**k, (2*k+1)*denominator**(2*k+1)) for k in range(terms)), F(0))


def lambda_enclosure():
    # Even 40-term sum is below arctan; odd 41-term sum is above it.
    a5_low, a5_high = arctan_sum(5, 40), arctan_sum(5, 41)
    a239_low, a239_high = arctan_sum(239, 40), arctan_sum(239, 41)
    require(a5_low < a5_high and a239_low < a239_high, "arctangent bounds ordering")
    pi_low = 16*a5_low - 4*a239_high
    pi_high = 16*a5_high - 4*a239_low

    # Exact tangent-addition values proving Machin's identity on the proper branch.
    tan2 = F(2, 5) / (1-F(1, 25))
    tan4 = 2*tan2 / (1-tan2*tan2)
    require(tan2 == F(5, 12) and tan4 == F(120, 119), "Machin double-angle values")
    require((tan4-F(1, 239))/(1+tan4*F(1, 239)) == 1, "Machin final tangent")
    require(4*(F(1, 5)-F(1, 375))-F(1, 239) > 0, "Machin angle positive")
    require(F(4, 5) < 1, "Machin angle below 1")

    scale = 10**60
    root = isqrt(5*scale*scale)
    require(root*root < 5*scale*scale < (root+1)*(root+1), "sqrt5 square certificate")
    sqrt_low, sqrt_high = F(root, scale), F(root+1, scale)
    phi_low, phi_high = (1+sqrt_low)/2, (1+sqrt_high)/2
    low = phi_low/(phi_low+pi_high)
    high = phi_high/(phi_high+pi_low)
    require(3 < pi_low < pi_high, "pi bound above 3")
    require(1 < phi_low < phi_high < 2, "phi interval")
    require(ELL < low < high < UPPER, "claimed lambda enclosure")
    q = F(Q, S)
    require(q < ELL and UPPER < F(Q+1, S), "canonical decimal truncation")
    print("PASS 2: exact Machin / integer-square-root certificate for ell < lambda < u")
    print("        sqrt5 certificate integer:", root)
    print("        ell =", ELL)
    print("        u   =", UPPER)
    print("        0 < q < lambda < 1/2; Q < S*lambda < Q+1")
    return low, high


def threshold_certificate():
    q = F(Q, S)
    left, middle, right = F(846547, 2490212), F(14342163, D_STAR), F(13495616, 39698853)
    require(left < q < middle < ELL < UPPER < right, "Farey interval ordering")
    det = right.numerator*left.denominator-left.numerator*right.denominator
    require(det == 1, "Farey determinant")
    require(left.denominator+right.denominator == D_STAR, "Farey denominator obstruction")
    require((q*D_STAR).__floor__() == 14342162, "q threshold floor")
    require(14342163 < ELL*D_STAR < UPPER*D_STAR < 14342164, "lambda threshold floor")
    print("PASS 3: first-disagreement certificate; determinant 1 and denominator sum 42189065")
    print("        q floor at d* = 14342162; lambda floor at d* = 14342163")


def encode(seed, horizon, numerator=Q, denominator=S):
    nonnegative_integer(seed, "seed")
    nonnegative_integer(horizon, "horizon")
    require(type(numerator) is int and type(denominator) is int
            and 0 < numerator < denominator, "0 < Q < S required")
    states, residuals = [seed], []
    for _ in range(horizon):
        nxt, rho = divmod(numerator*states[-1], denominator)
        states.append(nxt)
        residuals.append(rho)
    return states, residuals


def decode(terminal, residuals, horizon=None, numerator=Q, denominator=S):
    nonnegative_integer(terminal, "terminal")
    require(type(numerator) is int and type(denominator) is int
            and 0 < numerator < denominator, "0 < Q < S required")
    if horizon is not None:
        nonnegative_integer(horizon, "horizon")
        require(len(residuals) == horizon, "record length disagrees with horizon")
    for rho in residuals:
        require(type(rho) is int and 0 <= rho < denominator, "remainder outside [0,S)")
    reverse_states = [terminal]
    for rho in reversed(residuals):
        previous, failure = divmod(denominator*reverse_states[-1]+rho, numerator)
        require(failure == 0, "backward numerator is not divisible by Q")
        require(previous >= 0, "negative reconstructed state")
        require(divmod(numerator*previous, denominator) == (reverse_states[-1], rho),
                "forward recheck failed")
        reverse_states.append(previous)
    return list(reversed(reverse_states))


def master_numerator(terminal, residuals, numerator=Q, denominator=S):
    n = len(residuals)
    return denominator**n*terminal + sum(
        rho*denominator**k*numerator**(n-1-k) for k, rho in enumerate(residuals)
    )


def must_reject(terminal, residuals, horizon=None):
    try:
        decode(terminal, residuals, horizon=horizon)
    except ValueError:
        return
    raise ValueError("invalid record was accepted")


def symbolic_difference_checks():
    q, X = F(Q, S), Poly((0, 1))
    seeds = (0, 1, 2, 3, 17, 26, MAX_SEED, D_STAR, D_STAR+1, S//2, 10**60)
    count = 0
    for seed in seeds:
        for horizon in (0, 1, 2, 7, HORIZON):
            states, residues = encode(seed, horizon)
            eta = [F(rho, S) for rho in residues]
            for k in range(horizon):
                require(X*states[k]-states[k+1] == (X-q)*states[k]+eta[k],
                        "symbolic one-step difference")
            rational_loss = sum((q**(horizon-1-k)*eta[k] for k in range(horizon)), F(0))
            require(q**horizon*seed-states[-1] == rational_loss, "rational accumulated loss")
            lhs = X**horizon*seed-states[-1]
            separated = (X**horizon-q**horizon)*seed+rational_loss
            local = sum((X**(horizon-1-k)*((X-q)*states[k]+eta[k])
                         for k in range(horizon)), Poly((0,)))
            require(lhs == separated == local, "symbolic accumulated difference")
            power_diff = (X-q)*sum((X**(horizon-1-j)*q**j for j in range(horizon)), Poly((0,)))
            require(X**horizon-q**horizon == power_diff, "power-difference identity")
            if horizon:
                require(0 <= rational_loss < (1-q**horizon)/(1-q), "quantization bound")
            require(decode(states[-1], residues, horizon) == states, "example round trip")
            require(master_numerator(states[-1], residues) == Q**horizon*seed, "example master identity")
            # Recovered ideal-scale coefficient is exactly seed*X**j for every tested j.
            recovered = F(master_numerator(states[-1], residues), Q**horizon)
            require(recovered == seed, "seed from global identity for valid record")
            for j in (0, 1, horizon, horizon+1):
                require(X**j*recovered == X**j*seed, "ideal geometry reconstruction coefficient")
            count += 1
    print("PASS 4: symbolic local and accumulated differences, bounds, and ideal-scale recovery")
    print("       ", count, "seed/horizon cases; X is an indeterminate, never a floating approximation to lambda")


def admissibility_checks():
    rho0, rho1 = 2*Q-S//2, Q//2
    require(0 <= rho0 < S and 0 <= rho1 < S, "counterexample remainder ranges")
    require(master_numerator(0, [rho0, rho1]) == 2*Q**2, "counterexample global identity")
    require(F(rho1, Q) == F(1, 2), "counterexample nonintegral intermediate")
    must_reject(0, [rho0, rho1], 2)
    must_reject(-1, [])
    must_reject(0, [-1])
    must_reject(0, [S])
    must_reject(0, [True])
    must_reject(0, [F(Q)])
    must_reject(0, [], 1)
    states, residues = encode(17, 3)
    residues[0] += 1
    must_reject(states[-1], residues, 3)
    require(decode(0, [Q], 1) == [1, 0], "valid seed 1")
    require(decode(0, [2*Q], 1) == [2, 0], "distinct valid seed 2")
    print("PASS 5: fail-closed record validation; global identity alone is insufficient")
    print("        false record: terminal 0; residuals", rho0, rho1, "; last inverse is 1/2")


def information_and_extinction_checks():
    period = S//gcd(Q, S)
    inverse = pow(Q//2, -1, period)
    require(gcd(Q, S) == 2 and period == 5*10**14, "canonical residue period")
    require(MAX_SEED < period, "canonical injectivity interval")
    require(Q**10*MAX_SEED < S**10, "canonical extinction bound")
    require(3 < ELL**2*26 < UPPER**2*26 < 4, "single-round ideal floor 3")
    require(encode(26, 2)[0] == [26, 8, 2], "repeated quantization example")
    for endpoint in range(1000):
        lo = (S*endpoint+Q-1)//Q
        hi = (S*(endpoint+1)+Q-1)//Q-1
        require(hi-lo+1 in (2, 3), "canonical fiber size")
        for d in range(lo, hi+1):
            require((Q*d)//S == endpoint, "fiber membership")
        if lo:
            require((Q*(lo-1))//S != endpoint, "fiber lower endpoint")
        require((Q*(hi+1))//S != endpoint, "fiber upper endpoint")
    # An example proving that ordering can distinguish two valid histories.
    require(decode(0, [1, 0], 2, 1, 2) == [1, 0, 0], "ordered example 1")
    require(decode(0, [0, 1], 2, 1, 2) == [2, 1, 0], "ordered example 2")
    print("PASS 6: fibers, residue period, ordered-history example, and canonical extinction bound")
    print("        first-residue inverse =", inverse, "mod", period)
    print("        repeated quantization 26 -> 8 -> 2; floor(lambda^2 * 26) = 3")
    return inverse, period


def canonical_audit(inverse, period):
    transition_count = 0
    max_extinction = 0
    for seed in range(MAX_SEED+1):
        product = F(Q*seed, S)
        require(product.__floor__() == (ELL*seed).__floor__() == (UPPER*seed).__floor__(),
                "canonical floor agreement")
        first_remainder = (Q*seed) % S
        require(first_remainder % 2 == 0, "residual parity")
        require((inverse*(first_remainder//2)) % period == seed, "first-residue seed inverse")
        states, residues = [seed], []
        ledger_sum, q_power, s_power = 0, 1, 1
        extinction = 0 if seed == 0 else None
        for k in range(HORIZON):
            previous = states[-1]
            nxt, rho = divmod(Q*previous, S)
            require(0 <= rho < S and Q*previous == S*nxt+rho, "canonical local ledger identity")
            require(nxt <= previous and (previous == 0 or nxt < previous), "state contraction")
            require(F(S*nxt+rho, Q) == previous, "canonical local inverse")
            ledger_sum = Q*ledger_sum+rho*s_power
            q_power *= Q
            s_power *= S
            require(q_power*seed == s_power*nxt+ledger_sum, "canonical prefix invariant")
            states.append(nxt)
            residues.append(rho)
            if nxt == 0 and extinction is None:
                extinction = k+1
            transition_count += 1
        require(extinction is not None and extinction <= 10, "canonical extinction")
        max_extinction = max(max_extinction, extinction)
        require(decode(states[-1], residues, HORIZON) == states, "canonical full history round trip")
        # Every prefix obeys the master invariant; every transition has an exact inverse.
        require(states[-1] == 0, "canonical terminal")
    require(transition_count == 753664 and max_extinction == 10, "canonical audit totals")
    print("PASS 7: exhaustive canonical audit")
    print("        states:", MAX_SEED+1)
    print("        horizon:", HORIZON)
    print("        transitions / prefix invariants:", transition_count)
    print("        complete ledger round trips:", MAX_SEED+1)
    print("        first-residue seed inversions:", MAX_SEED+1)
    print("        maximum extinction time:", max_extinction)


def main():
    print("OMEGA geometry, numerical difference, and exact recovery verifier v1.0")
    print("Q =", Q, "; S =", S)
    shape_certificate()
    lambda_enclosure()
    threshold_certificate()
    symbolic_difference_checks()
    admissibility_checks()
    inverse, period = information_and_extinction_checks()
    canonical_audit(inverse, period)
    source = Path(__file__).read_bytes()
    require(not any(isinstance(node, ast.Assert) for node in ast.walk(ast.parse(source))),
            "source contains assert statements")
    print("PASS 8: verifier contains no assert statements")
    print("SOURCE SHA-256:", hashlib.sha256(source).hexdigest())
    print("ALL 8 GROUPS PASS")


if __name__ == "__main__":
    try:
        main()
    except Exception as error:
        print("FAIL:", str(error), file=sys.stderr)
        raise