#!/usr/bin/env python3
"""Exact evidence for OMEGA: Relative Representation and Recoverable Information.

Python 3.10+, standard library only. Run normally and with python -O.
Symbolic identities are checked by sparse rational-polynomial normalization.
Finite audits are separately labeled. No floating-point acceptance tests, no
Lean claim, no assertion that a successful run proves arbitrary analytic facts.
The auxiliary six-mark seed is s=1/2, NOT the canonical sine-defined seed.
"""
from __future__ import annotations
import ast
from dataclasses import dataclass
from fractions import Fraction as F
from itertools import combinations, product
from math import gcd, isqrt
from pathlib import Path
import hashlib
import json

Q, S, N = 339949771344778, 10**15, 23
GROUPS: dict[str, dict] = {}


def require(ok: bool, label: str) -> None:
    if not ok:
        raise ValueError(label)


def natural(x: int, label: str) -> None:
    if type(x) is not int or x < 0:
        raise ValueError(label + ': nonnegative integer required')


# Polynomial coefficients are rational; monomials are tuples of exponents.
# The empty dictionary is the unique zero normal form.
class Poly:
    def __init__(self, terms: dict, width: int):
        self.width = width
        self.terms = {e: F(c) for e, c in terms.items() if c != 0}

    @classmethod
    def constant(cls, c, width):
        return cls({(0,) * width: F(c)}, width)

    @classmethod
    def variable(cls, i, width):
        e = [0] * width
        e[i] = 1
        return cls({tuple(e): F(1)}, width)

    def lift(self, x):
        if isinstance(x, Poly):
            require(x.width == self.width, 'polynomial arity')
            return x
        return Poly.constant(x, self.width)

    def __add__(self, x):
        x = self.lift(x)
        out = dict(self.terms)
        for e, c in x.terms.items():
            out[e] = out.get(e, F(0)) + c
        return Poly(out, self.width)
    __radd__ = __add__

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

    def __sub__(self, x):
        return self + -self.lift(x)

    def __rsub__(self, x):
        return self.lift(x) - self

    def __mul__(self, x):
        x = self.lift(x)
        out = {}
        for e, c in self.terms.items():
            for f, d in x.terms.items():
                h = tuple(a + b for a, b in zip(e, f))
                out[h] = out.get(h, F(0)) + c * d
        return Poly(out, self.width)
    __rmul__ = __mul__

    def __pow__(self, n):
        natural(n, 'polynomial exponent')
        out, b = self.lift(1), self
        while n:
            if n & 1:
                out = out * b
            b = b * b
            n //= 2
        return out

    def zero(self, label):
        require(not self.terms, label)


def symbolic_proofs():
    m, a, t, u, v, lam, ph, pi, den = [Poly.variable(i, 9) for i in range(9)]
    (m + t * ((m + a * u) - m) - (m + t * a * u)).zero('similarity')
    ((ph * v + pi) - ((ph + pi) + ph * (v - 1))).zero('normalization')
    # Inductive step, with u denoting lambda**k as an indeterminate.
    (1 + lam * ((1 + u * (v - 1)) - 1)
     - (1 + lam * u * (v - 1))).zero('affine iteration induction')
    ((v * ph + pi - pi) - ph * v).zero('affine inverse numerator')
    checks = 4
    for n in range(1, 24):
        width = n + 3
        q, s = Poly.variable(0, width), Poly.variable(1, width)
        d = [Poly.variable(i + 2, width) for i in range(n + 1)]
        total = Poly.constant(0, width)
        for k in range(n):
            total += (q * d[k] - s * d[k + 1]) * s**k * q**(n-1-k)
        (q**n * d[0] - s**n * d[n] - total).zero('telescoping')
        checks += 1
    x = Poly.variable(0, 1)
    (x**5 + x + 1 - (x**2 + x + 1)*(x**3 - x**2 + 1)).zero('factorization')
    checks += 1
    require((x + 1).terms != x.terms, 'nonzero polynomial rejection control')
    GROUPS['symbolic_polynomials'] = {
        'identities': checks, 'telescoping_horizons': [1, 23],
        'scope': 'all substitutions into each tested polynomial identity',
        'unbounded_horizon': 'proved by induction in the thesis, not this loop'}


@dataclass(frozen=True)
class Golden:
    """a+b*phi in Q(phi), phi**2=phi+1; exact, canonical coefficients."""
    a: F = F(0)
    b: F = F(0)

    def __post_init__(self):
        for x in (self.a, self.b):
            if type(x) is not int and not isinstance(x, F):
                raise TypeError('integer or Fraction required')
        object.__setattr__(self, 'a', F(self.a))
        object.__setattr__(self, 'b', F(self.b))

    @staticmethod
    def lift(x):
        return x if isinstance(x, Golden) else Golden(x)

    def __add__(self, x):
        x = self.lift(x)
        return Golden(self.a+x.a, self.b+x.b)
    __radd__ = __add__

    def __neg__(self):
        return Golden(-self.a, -self.b)

    def __sub__(self, x):
        return self + -self.lift(x)

    def __rsub__(self, x):
        return self.lift(x) - self

    def __mul__(self, x):
        x = self.lift(x)
        return Golden(self.a*x.a+self.b*x.b,
                      self.a*x.b+self.b*x.a+self.b*x.b)
    __rmul__ = __mul__

    def __truediv__(self, x):
        x = self.lift(x)
        norm = x.a*x.a+x.a*x.b-x.b*x.b
        if norm == 0:
            raise ZeroDivisionError('zero in Q(phi)')
        return self * Golden((x.a+x.b)/norm, -x.b/norm)

    def __pow__(self, n):
        natural(n, 'algebraic exponent')
        out, b = Golden(1), self
        while n:
            if n & 1:
                out = out*b
            b = b*b
            n //= 2
        return out

    def record(self):
        return [str(self.a), str(self.b)]


ZERO, ONE, PHI = Golden(), Golden(1), Golden(0, 1)


def inverse(matrix):
    n = len(matrix)
    require(n > 0 and all(len(r) == n for r in matrix), 'square matrix')
    a = [[Golden.lift(x) for x in row] +
         [ONE if i == j else ZERO for j in range(n)]
         for i, row in enumerate(matrix)]
    for col in range(n):
        pivot = next((r for r in range(col, n) if a[r][col] != ZERO), None)
        if pivot is None:
            raise ValueError('singular matrix')
        a[col], a[pivot] = a[pivot], a[col]
        p = a[col][col]
        a[col] = [x/p for x in a[col]]
        for r in range(n):
            if r != col:
                t = a[r][col]
                a[r] = [x-t*y for x, y in zip(a[r], a[col])]
    return [row[n:] for row in a]


def mv(a, v):
    require(all(len(row) == len(v) for row in a), 'matrix-vector dimensions')
    return [sum((x*y for x, y in zip(row, v)), ZERO) for row in a]


def exact_measurements():
    require(PHI*PHI == PHI+1, 'golden relation')
    seed = Golden(F(1, 2))  # Explicit auxiliary geometry, not canonical s.
    nodes = [-seed-PHI**2, -seed-1, -seed, seed, seed+1, seed+PHI**2]
    require(len(set(nodes)) == 6, 'distinct nodes')
    point_minors = interval_minors = 0
    for degree in range(6):
        for ids in combinations(range(6), degree+1):
            a = [[nodes[i]**k for k in range(degree+1)] for i in ids]
            coefficients = [Golden(k-2, F(k, 7)) for k in range(degree+1)]
            require(mv(inverse(a), mv(a, coefficients)) == coefficients,
                    'point interpolation')
            point_minors += 1
    for degree in range(5):
        for ids in combinations(range(5), degree+1):
            a = [[(nodes[i+1]**(k+1)-nodes[i]**(k+1))/(k+1)
                  for k in range(degree+1)] for i in ids]
            coefficients = [Golden(k-1, F(k+1, 5)) for k in range(degree+1)]
            require(mv(inverse(a), mv(a, coefficients)) == coefficients,
                    'interval interpolation')
            interval_minors += 1
    # One explicit affine word, all pairs of corrupted locations.
    values = [3-2*z for z in nodes]
    error_pairs = 0
    for bad in combinations(range(6), 2):
        observed = list(values)
        for i in bad:
            observed[i] += 7+i
        candidates = set()
        for i, j in combinations(range(6), 2):
            slope = (observed[j]-observed[i])/(nodes[j]-nodes[i])
            offset = observed[i]-slope*nodes[i]
            if sum(offset+slope*z != y for z, y in zip(nodes, observed)) <= 2:
                candidates.add((offset, slope))
        require(candidates == {(Golden(3), Golden(-2))}, 'two-error correction')
        error_pairs += 1
    GROUPS['six_mark_measurements'] = {
        'auxiliary_seed': '1/2', 'point_interpolation_subsets': point_minors,
        'interval_interpolation_subsets': interval_minors,
        'two_error_position_pairs_one_word': error_pairs,
        'scope': 'exact finite instances; general distance proof is in thesis'}
    return nodes, values


def forward(d, horizon, q=Q, s=S):
    natural(d, 'seed'); natural(horizon, 'horizon')
    if type(q) is not int or type(s) is not int or not 0 < q < s:
        raise ValueError('require integers 0<Q<S')
    residues = []
    for _ in range(horizon):
        d, r = divmod(q*d, s)
        residues.append(r)
    return d, tuple(residues)


def backward(terminal, residues, horizon, q=Q, s=S):
    natural(terminal, 'terminal'); natural(horizon, 'horizon')
    if type(q) is not int or type(s) is not int or not 0 < q < s:
        raise ValueError('require integers 0<Q<S')
    if not isinstance(residues, (tuple, list)) or len(residues) != horizon:
        raise ValueError('wrong residue count')
    d = terminal
    for r in reversed(residues):
        natural(r, 'residue')
        require(r < s, 'residue bound')
        previous, rem = divmod(s*d+r, q)
        require(rem == 0, 'reverse divisibility')
        require(divmod(q*previous, s) == (d, r), 'reverse consistency')
        d = previous
    return d


def finite_ledger_audit():
    powers_q = [Q**k for k in range(N+1)]
    powers_s = [S**k for k in range(N+1)]
    g = gcd(Q, S); modulus = S//g
    inv = pow(Q//g, -1, modulus)
    max_extinction = transitions = 0
    for seed in range(32768):
        d, witness, residues = seed, 0, []
        extinction = 0 if seed == 0 else None
        for k in range(N):
            d, r = divmod(Q*d, S)
            residues.append(r)
            require(0 <= r < S, 'Euclidean range')
            witness = Q*witness+r*powers_s[k]
            require(powers_q[k+1]*seed == powers_s[k+1]*d+witness,
                    'prefix telescoping invariant')
            if d == 0 and extinction is None:
                extinction = k+1
            transitions += 1
        require(extinction is not None, 'finite extinction')
        max_extinction = max(max_extinction, extinction)
        require(backward(d, residues, N) == seed, 'full reverse')
        require(((residues[0]//g)*inv) % modulus == seed, 'first-residue inverse')
    terminal, good = forward(399, N)
    bad = list(good); bad[0] += 1
    invalid = [(terminal, bad, N), (terminal, good[:-1], N),
               (terminal, [S]+list(good[1:]), N),
               (terminal, [True]+list(good[1:]), N)]
    rejected = 0
    for t, r, n in invalid:
        try:
            backward(t, r, n)
        except ValueError:
            rejected += 1
    require(rejected == len(invalid), 'invalid records rejected')
    require(backward(399, [], 0) == 399, 'zero horizon')
    GROUPS['canonical_ledger'] = {
        'seeds': 32768, 'horizons': [1, N], 'prefix_identities': transitions,
        'full_reversals': 32768, 'first_residue_inversions': 32768,
        'maximum_extinction': max_extinction, 'invalid_records_rejected': rejected}


def floor_sum(n, m, a, b=0):
    """sum_{i=0}^{n-1} floor((a*i+b)/m), all arguments nonnegative, m>0."""
    for value in (n, a, b):
        natural(value, 'floor sum input')
    require(type(m) is int and m > 0, 'positive modulus')
    total = 0
    while True:
        qa, a = divmod(a, m); qb, b = divmod(b, m)
        total += qa*n*(n-1)//2 + qb*n
        y = a*n+b
        if y < m:
            return total
        n, b = divmod(y, m)
        m, a = a, m


def atan_bounds(q, terms):
    lo_or_hi = sum((F((-1)**k, (2*k+1)*q**(2*k+1))
                    for k in range(terms)), F(0))
    adjacent = lo_or_hi+F((-1)**terms, (2*terms+1)*q**(2*terms+1))
    return min(lo_or_hi, adjacent), max(lo_or_hi, adjacent)


def certified_floor_prefix():
    small = 0
    for n, m, a, b in product(range(9), range(1, 7), range(9), range(5)):
        require(floor_sum(n,m,a,b) == sum((a*i+b)//m for i in range(n)),
                'independent small floor sum')
        small += 1
    al, au = atan_bounds(5, 50); bl, bu = atan_bounds(239, 18)
    pi_l, pi_u = 16*al-4*bu, 16*au-4*bl
    m = 10**80; r = isqrt(5*m*m)
    require(r*r < 5*m*m < (r+1)*(r+1), 'strict sqrt bracket')
    ph_l, ph_u = (1+F(r,m))/2, (1+F(r+1,m))/2
    low, high = ph_l/(ph_l+pi_u), ph_u/(ph_u+pi_l)
    require(F(Q,S) < low < high, 'lambda brackets and rational slope')
    dstar = 42189065
    sums = [floor_sum(dstar, r.denominator, r.numerator)
            for r in (F(Q,S), low, high)]
    require(sums == [302541195258184]*3, 'prefix sum equality')
    real_ends = [(r*dstar).__floor__() for r in (low, high)]
    rational_end = Q*dstar//S
    require(real_ends == [14342163]*2 and rational_end == 14342162,
            'first-disagreement endpoint')
    GROUPS['floor_prefix_certificate'] = {
        'small_direct_crosschecks': small, 'first_disagreement': dstar,
        'common_floor_sum': sums[0], 'real_floor_at_boundary': real_ends[0],
        'rational_floor_at_boundary': rational_end,
        'method': 'rational interval bounds plus exact lattice floor sums'}


def icosahedral_certificate():
    vertices = []
    for e, f in product((-1, 1), repeat=2):
        vertices += [(ZERO, Golden(e), f*PHI),
                     (Golden(e), f*PHI, ZERO), (f*PHI, ZERO, Golden(e))]
    radius2 = PHI+2
    require(all(sum((x*x for x in p), ZERO) == radius2 for p in vertices),
            'one sphere')
    moments = 0
    for total in range(6):
        for i in range(total+1):
            for j in range(total-i+1):
                k = total-i-j
                got = sum((x**i*y**j*z**k for x,y,z in vertices), ZERO)/12
                if any(e % 2 for e in (i,j,k)):
                    expected = ZERO
                elif total == 0:
                    expected = ONE
                elif total == 2:
                    expected = radius2/3
                else:
                    expected = radius2**2/(5 if 4 in (i,j,k) else 15)
                require(got == expected, 'exact spherical monomial')
                moments += 1
    features = [[ONE,x,y,z,x*y,x*z,y*z,x*x-y*y,2*z*z-x*x-y*y]
                for x,y,z in vertices]
    diag = [ONE]+[radius2/3]*3+[radius2**2/15]*3+[4*radius2**2/15,4*radius2**2/5]
    for i,j in product(range(9), repeat=2):
        entry = sum((r[i]*r[j] for r in features), ZERO)/12
        require(entry == (diag[i] if i == j else ZERO), 'harmonic Gram entry')
    coeff = [Golden(i-3, F(i,11)) for i in range(9)]
    samples = mv(features, coeff)
    recovered = [sum((r[j]*y for r,y in zip(features,samples)), ZERO)/(12*diag[j])
                 for j in range(9)]
    require(recovered == coeff, 'nine-mode recovery')
    defect = sum((x**6 for x,y,z in vertices), ZERO)/12-radius2**3/7
    require(defect == -(4*PHI+3)/21 and defect != ZERO, 'degree-six obstruction')
    GROUPS['icosahedral_cubature'] = {
        'nodes': 12, 'exact_monomials_degree_at_most_5': moments,
        'gram_entries': 81, 'recovered_harmonic_modes': 9,
        'sixth_moment_defect_coefficients_in_1_phi': defect.record()}


def relational_examples(nodes, values):
    # Kuratowski representation of a declared finite line metric.
    points = [F(-3), F(-1,2), F(0), F(2), F(7)]
    base = points[2]
    enc = [[abs(x-z)-abs(base-z) for z in points] for x in points]
    for i,j in product(range(len(points)), repeat=2):
        require(max(abs(a-b) for a,b in zip(enc[i],enc[j])) == abs(points[i]-points[j]),
                'metric representation')
    adjunctions = 0
    for n,d in product(range(30), range(100)):
        left = (S*n+Q-1)//Q
        require((left <= d) == (n <= Q*d//S), 'order adjunction')
        adjunctions += 1
    # One observer, two times: x=(position,velocity), step (p,v)->(p+v,v).
    for p,v in product(range(-5,6), repeat=2):
        y0,y1 = p,p+v
        require((y0,y1-y0) == (p,v), 'two-time observer')
    payload = json.dumps([x.record() for x in values], separators=(',',':')).encode('ascii')
    packets = [forward(b,N) for b in payload]
    recovered = bytes(backward(t,r,N) for t,r in packets)
    require(recovered == payload, 'exact byte reconstruction')
    decoded = [Golden(F(a),F(b)) for a,b in json.loads(recovered.decode('ascii'))]
    slope = (decoded[-1]-decoded[0])/(nodes[-1]-nodes[0])
    offset = decoded[0]-slope*nodes[0]
    require((offset,slope) == (Golden(3),Golden(-2)), 'end-to-end field reconstruction')
    GROUPS['relational_and_end_to_end'] = {
        'metric_pairs': 25, 'order_adjunction_instances': adjunctions,
        'two_time_observer_states': 121, 'payload_bytes': len(payload),
        'byte_ledgers': len(packets), 'field': 'p(z)=3-2z',
        'recovered_offset': offset.record(), 'recovered_slope': slope.record()}


def main():
    source = Path(__file__).read_bytes()
    tree = ast.parse(source)
    require(not any(isinstance(x, ast.Assert) for x in ast.walk(tree)), 'no disabled asserts')
    require(not any(isinstance(x, ast.Constant) and isinstance(x.value,float)
                    for x in ast.walk(tree)), 'no floating constants')
    symbolic_proofs()
    nodes, values = exact_measurements()
    finite_ledger_audit()
    certified_floor_prefix()
    icosahedral_certificate()
    relational_examples(nodes, values)
    report = {'status': 'PASS', 'source_sha256': hashlib.sha256(source).hexdigest(),
              'groups': GROUPS,
              'limits': ['Not a proof-assistant kernel certificate.',
                         'Finite examples do not prove universal analytic statements.',
                         'No physical or neural-performance experiment is claimed.']}
    print(json.dumps(report, indent=2, sort_keys=True))


if __name__ == '__main__':
    main()