"""Independent exact replay of the Oct 30 P95 submission; imports no submitted code.

Usage: python tools/audit-p95-proof-2026-10-30.py ARCHIVE [OUTPUT]
The near-axis analytic lemma is a manually reviewed premise, not formalized here.
"""
import gzip
import copy
import hashlib
import itertools
import json
import sys
import time
import zipfile
from fractions import Fraction as Q
from functools import lru_cache
from pathlib import Path

ROOT = ((Q(-1), Q(1)), (Q(-1), Q(1)), (Q(0), Q(1, 2)))
PAIRS = ((0, 1), (0, 2), (1, 2))


def need(ok, message):
    if not ok:
        raise ValueError(message)


def intersect(a, b):
    return max(a[0], b[0]) <= min(a[1], b[1])


def quarter_turn_crossing(left, right):
    # alpha = beta + pi/2: beta must be in [-pi/2,0].
    # tan((beta+pi/2)/2)=(1+t)/(1-t), increasing on [-1,0].
    low, high = max(right[0], Q(-1)), min(right[1], Q(0))
    return low <= high and intersect(left, ((1+low)/(1-low), (1+high)/(1-high)))


@lru_cache(None)
def unit(t):
    return (1-t*t)/(1+t*t), 2*t/(1+t*t)


@lru_cache(None)
def support(a, b, left, right):
    values = []
    for x, y in ((left[0], right[1]), (left[1], right[0])):
        c, s = unit(x)
        d, t = unit(y)
        values.append(a*abs(c*d+s*t)+b*abs(s*d-c*t))
    if intersect(left, right):
        values.append(a)
    if quarter_turn_crossing(left, right) or quarter_turn_crossing(right, left):
        values.append(b)
    # +/-pi can occur only at domain endpoints, already evaluated above.
    return min(values)


def mul(interval, scalar):
    return tuple(sorted(x*scalar for x in interval))


def geometry(box, n):
    target = Q(3*n-4)
    sizes = tuple(range(n-2, n+1))
    extents = [(support(Q(k), Q(k, 2), t, (Q(0), Q(0))),
                support(Q(k, 2), Q(k), t, (Q(0), Q(0)))) for k, t in zip(sizes, box)]
    radii = [target/2-h for extent in extents for h in extent]
    positions = [(-max(Q(0), r), max(Q(0), r)) for r in radii]
    positions[4] = (positions[4][0], Q(0))
    rows = []
    for coordinate in range(6):
        for sign in (-1, 1):
            row = [(Q(0), Q(0))]*6
            row[coordinate] = (Q(sign), Q(sign))
            rows.append((row, radii[coordinate]))
    row = [(Q(0), Q(0))]*6
    row[4] = (Q(1), Q(1))
    rows.append((row, Q(0)))
    pairs = {}
    for pair, (i, j) in enumerate(PAIRS):
        for label in range(8):
            owner, other = (i, j) if label < 4 else (j, i)
            axis, polarity = (label % 4)//2, 1 if label % 2 == 0 else -1
            low, high = box[owner]
            c = [unit(low)[0], unit(high)[0]]
            if low <= 0 <= high:
                c.append(Q(1))
            cosine, sine = (min(c), max(c)), (unit(low)[1], unit(high)[1])
            normal = (cosine, sine) if axis == 0 else (mul(sine, -1), cosine)
            normal = [mul(v, polarity) for v in normal]
            a, b = (Q(sizes[other]), Q(sizes[other], 2)) if axis == 0 else (Q(sizes[other], 2), Q(sizes[other]))
            required = Q(sizes[owner], 1 if axis == 0 else 2)+support(a, b, box[owner], box[other])
            row = [(Q(0), Q(0))]*6
            for d in range(2):
                row[2*i+d], row[2*j+d] = normal[d], mul(normal[d], -1)
            possible = sum(max(abs(x) for x in normal[d])*max(Q(0), target-extents[i][d]-extents[j][d]) for d in range(2))
            pairs[pair, label] = (row, -required, possible < required)
    return radii, positions, rows, pairs


def replay(records, n):
    target = 3*n-4
    header = records[0]
    need(header['schema'] == 'continuous-angle-strict-target-v1', 'schema')
    need(Q(header['T']) == target and header['pieces'] == list(range(n-2, n+1)), 'target/pieces')
    need(tuple(tuple(map(Q, t)) for t in header['root']) == ROOT, 'root')
    stack = [(ROOT, set(itertools.product(range(8), repeat=3)))]
    counts = dict(nodes=0, duals=0, analytic=0, walls=0, pair_prunes=0)
    for record in records[1:-1]:
        need(bool(stack), 'extra node')
        box, branches = stack.pop()
        counts['nodes'] += 1
        need(record['node'] == counts['nodes'] and record['incoming'] == len(branches), 'node/inheritance')
        analytic = 'analytic_pair35' if n == 5 else 'analytic_consecutive_pair'
        if analytic in record:
            need(record[analytic] is True and set(record) == {'node', 'incoming', analytic}, 'analytic schema')
            need(box[2][1] <= Q(1, 200) and (box[0][0] >= Q(199, 201) or box[0][1] <= -Q(199, 201)), 'analytic domain')
            counts['analytic'] += 1
            continue
        radii, positions, walls, pairs = geometry(box, n)
        if 'wall' in record:
            k = record['wall']
            need(type(k) is int and 0 <= k < 3 and min(radii[2*k:2*k+2]) < 0, 'wall exclusion')
            counts['walls'] += 1
            continue
        for pair, label in record['pair_prunes']:
            need((pair, label) in pairs and pairs[pair, label][2], 'pair exclusion')
            remaining = {branch for branch in branches if branch[pair] != label}
            counts['pair_prunes'] += len(branches)-len(remaining)
            branches = remaining
        for certificate in record['duals']:
            branch = tuple(certificate['branch'])
            need(branch in branches, 'uncovered/duplicate branch')
            residual = [[Q(0), Q(0)] for _ in range(6)]
            rhs = mass = Q(0)
            for index, raw_weight in certificate['dual']:
                weight = Q(raw_weight)
                need(type(index) is int and 0 <= index < 16 and weight >= 0, 'dual index/sign')
                row, bound = walls[index] if index < 13 else pairs[index-13, branch[index-13]][:2]
                mass += weight
                rhs += weight*bound
                for d in range(6):
                    for end in range(2):
                        residual[d][end] += weight*row[d][end]
            lower = sum(min(a*b for a in r for b in z) for r, z in zip(residual, positions))
            need(mass > 0 and rhs-lower <= Q(certificate['upper']) < 0, 'dual contradiction')
            branches.remove(branch)
            counts['duals'] += 1
        if branches:
            d, mid = record['split']
            mid = Q(mid)
            need(type(d) is int and 0 <= d < 3 and box[d][0] < mid < box[d][1], 'split')
            left, right = list(box), list(box)
            left[d], right[d] = (box[d][0], mid), (mid, box[d][1])
            stack.extend([(tuple(right), branches.copy()), (tuple(left), branches.copy())])
        else:
            need('split' not in record, 'unnecessary split')
    end = records[-1]
    need(end.get('end') is True and end.get('complete') is True and end.get('nodes') == counts['nodes'] and not stack, 'incomplete tree')
    return dict(n=n, target=target, **counts)


def main():
    start = time.monotonic()
    archive = Path(sys.argv[1])
    result = {'archive_sha256': hashlib.sha256(archive.read_bytes()).hexdigest(), 'trees': []}
    with zipfile.ZipFile(archive) as bundle:
        for n in (5, 6, 7):
            raw = bundle.read(f'P95_n1_n7_real_optima/certificates/n{n}_strict_below_{3*n-4}.jsonl.gz')
            records = [json.loads(line) for line in gzip.decompress(raw).splitlines()]
            result['trees'].append(replay(records, n))
            print(f'n={n}: independent exact replay PASS', flush=True)
            if n == 5:
                controls = []
                broken = copy.deepcopy(records)
                broken.pop(-2)
                controls.append(('missing_final_node', broken))
                broken = copy.deepcopy(records)
                broken.insert(2, copy.deepcopy(broken[1]))
                controls.append(('duplicate_node', broken))
                broken = copy.deepcopy(records)
                first = next(r for r in broken[1:-1] if r.get('duals'))
                for term in first['duals'][0]['dual']:
                    term[1] = '0'
                controls.append(('zero_dual', broken))
                broken = copy.deepcopy(records)
                broken[1] = {'node': 1, 'incoming': 512, 'analytic_pair35': True}
                controls.append(('analytic_leaf_outside_cone', broken))
                result['negative_controls'] = []
                for name, broken in controls:
                    try:
                        replay(broken, n)
                    except (ValueError, KeyError, IndexError) as exc:
                        result['negative_controls'].append({'case': name, 'rejected': True, 'reason': str(exc)})
                    else:
                        raise ValueError('accepted corrupt certificate: '+name)
    result.update(passed=True, seconds=time.monotonic()-start,
                  trust_boundary='Manually reviewed near-axis lemma and SAT; exact certificate replay, not proof-assistant formalization.')
    if len(sys.argv) > 2:
        Path(sys.argv[2]).write_text(json.dumps(result, indent=2)+'\n', encoding='utf-8')
    print(json.dumps(result, indent=2))


if __name__ == '__main__':
    main()