#!/usr/bin/env python3
"""Exact educational geometry kernels. No third-party packages or floating epsilon.
The randomized bounds use ideal independent uniform draws, not a PRNG security claim.
"""
from fractions import Fraction as F
from itertools import combinations, product, permutations
from random import Random
from math import ceil
import json
from pathlib import Path


def check(test, message):
    if not test:
        raise AssertionError(message)


def point(p):
    if len(p) != 2 or any(not isinstance(x, (int, F)) or isinstance(x, bool) for x in p):
        raise ValueError('coordinates must be integers or Fractions')
    return tuple(F(x) for x in p)


def dot(a, b):
    return a[0]*b[0] + a[1]*b[1]


def sub(a, b):
    return (a[0]-b[0], a[1]-b[1])


def cross(a, b):
    return a[0]*b[1]-a[1]*b[0]


def distance2(a, b):
    v = sub(a, b)
    return dot(v, v)


def boundary_circle(points):
    """Minimum radius circle with ALL <=3 given points on its boundary.
    Three collinear distinct points have no such circle. This is not MEC(points).
    """
    if not points:
        return None
    if len(points) == 1:
        return (points[0], F(0))
    if len(points) == 2:
        a, b = points
        center = ((a[0]+b[0])/2, (a[1]+b[1])/2)
    elif len(points) == 3:
        a, b, c = points
        u, v = sub(b, a), sub(c, a)
        determinant = 2*cross(u, v)
        if determinant == 0:
            raise ValueError('three forced boundary points are collinear')
        du, dv = dot(u, u), dot(v, v)
        center = (a[0]+(du*v[1]-dv*u[1])/determinant,
                  a[1]+(u[0]*dv-v[0]*du)/determinant)
    else:
        raise ValueError('at most three boundary points')
    return center, distance2(center, points[0])


def inside(circle, p):
    return circle is not None and distance2(circle[0], p) <= circle[1]


def circle_certificate(points, circle):
    """Linear scan, fixed-radius boundary directions. Returns IDs and barycentric weights."""
    if circle is None:
        check(not points, 'empty circle needs empty input')
        return ()
    center, radius2 = circle
    boundary = [i for i, p in enumerate(points) if distance2(p, center) == radius2]
    check(boundary, 'nonempty boundary')
    first = boundary[0]
    if radius2 == 0:
        return ((first, F(1)),)
    a = sub(points[first], center)
    upper = lower = None
    for i in boundary[1:]:
        v = sub(points[i], center)
        sign, projection = cross(a, v), dot(a, v)
        if sign == 0 and projection < 0:
            return ((first, F(1, 2)), (i, F(1, 2)))
        if sign > 0 and (upper is None or projection < upper[0]):
            upper = projection, i
        if sign < 0 and (lower is None or projection < lower[0]):
            lower = projection, i
    check(upper is not None and lower is not None, 'center lies in boundary convex hull')
    ids = first, upper[1], lower[1]
    a, b, c = [sub(points[i], center) for i in ids]
    determinant = cross(sub(b, a), sub(c, a))
    weights = (cross(b, c)/determinant, cross(c, a)/determinant, cross(a, b)/determinant)
    check(min(weights) >= 0 and sum(weights) == 1, 'nonnegative barycentric certificate')
    return tuple((i, w) for i, w in zip(ids, weights) if w)


def verify_circle_certificate(points, circle, certificate):
    if circle is None:
        return not points and not certificate
    center, radius2 = circle
    if not (1 <= len(certificate) <= 3 and all(inside(circle, p) for p in points)):
        return False
    if sum(w for i, w in certificate) != 1:
        return False
    if any(w <= 0 or not 0 <= i < len(points) or distance2(center, points[i]) != radius2
           for i, w in certificate):
        return False
    return all(sum(w*points[i][k] for i, w in certificate) == center[k] for k in range(2))


def welzl(raw_points, rng=None, record_events=False, audit=False):
    """Fresh random choice per recursive node, simulated by an explicit stack.
    Prefix swaps are restored. No prefix slices and no Python recursion.
    """
    pts = [point(p) for p in raw_points]
    ids = list(range(len(pts)))
    rng = Random(0) if rng is None else rng
    # frame = [prefix length, forced boundary IDs, phase, chosen swap index]
    stack = [[len(pts), (), 0, None]]
    circle = None
    events = []
    stats = dict(nodes=0, tests=0, violations=0, swaps=0, max_stack=1)
    while stack:
        k, boundary, phase, j = stack[-1]
        if phase == 0:
            stats['nodes'] += 1
            if k == 0 or len(boundary) == 3:
                circle = boundary_circle([pts[i] for i in boundary])
                if audit:
                    check(all(inside(circle, pts[i]) for i in ids[:k]), 'reachable base covers prefix')
                stack.pop()
                continue
            j = rng.randrange(k)
            ids[j], ids[k-1] = ids[k-1], ids[j]
            stats['swaps'] += 1
            stack[-1][2:] = [1, j]
            stack.append([k-1, boundary, 0, None])
        elif phase == 1:
            p = ids[k-1]
            stats['tests'] += 1
            if inside(circle, pts[p]):
                ids[j], ids[k-1] = ids[k-1], ids[j]
                stack.pop()
            else:
                stats['violations'] += 1
                if record_events:
                    events.append(dict(prefix=k, point=p, forced=list(boundary), before=circle))
                stack[-1][2] = 2
                stack.append([k-1, boundary+(p,), 0, None])
        else:
            if audit:
                check(all(inside(circle, pts[i]) for i in ids[:k]), 'return covers prefix')
                check(all(distance2(circle[0], pts[i]) == circle[1] for i in boundary), 'forced boundary')
            ids[j], ids[k-1] = ids[k-1], ids[j]
            stack.pop()
        stats['max_stack'] = max(stats['max_stack'], len(stack))
    check(ids == list(range(len(pts))), 'all random swaps restored')
    certificate = circle_certificate(pts, circle)
    check(verify_circle_certificate(pts, circle, certificate), 'circle certificate')
    return dict(circle=circle, certificate=certificate, stats=stats, events=events)


def circle_enumeration(raw_points):
    pts = [point(p) for p in raw_points]
    if not pts:
        return None
    best = None
    for size in (1, 2, 3):
        for ids in combinations(range(len(pts)), size):
            try:
                candidate = boundary_circle([pts[i] for i in ids])
            except ValueError:
                continue
            if all(inside(candidate, p) for p in pts) and (best is None or candidate[1] < best[1]):
                best = candidate
    check(best is not None, 'finite set has enclosing circle')
    return best


# A = (slope, intercept) represents slope*M + intercept as M tends to +infinity.
# Lexicographic tuple order is eventual numerical order; only affine arithmetic is used.
ZERO = (F(0), F(0))
ONE_M = (F(1), F(0))


def add(a, b):
    return a[0]+b[0], a[1]+b[1]


def scale(k, a):
    return k*a[0], k*a[1]


def constant(k):
    return F(0), F(k)


def affine_dot(a, xy):
    return add(scale(a[0], xy[0]), scale(a[1], xy[1]))


def constraints(raw):
    answer = []
    for row in raw:
        if len(row) != 3 or any(not isinstance(x, (int, F)) or isinstance(x, bool) for x in row):
            raise ValueError('halfplanes need three exact rational coefficients')
        answer.append(tuple(F(x) for x in row))
    return answer


def box_rows():
    return [(F(1), F(0), ONE_M, -1), (F(-1), F(0), ONE_M, -2),
            (F(0), F(1), ONE_M, -3), (F(0), F(-1), ONE_M, -4)]


def satisfies(xy, row):
    return affine_dot(row[:2], xy) <= constant(row[2])


def symbolic_key(xy, objective):
    return affine_dot(objective, xy) + xy[0] + xy[1]


def boundary_lp(row, previous, objective, stats=None):
    a, b, c, identity = row
    if a == 0 and b == 0:
        check(c < 0, 'only contradictory zero row can violate')
        return None, (identity,)
    p = (c/a, F(0)) if a else (F(0), c/b)
    v = (-b, a)
    low = high = None
    low_id = high_id = None
    rows = box_rows()
    # Iterating a prefix iterator avoids copying a growing prefix at every call.
    for aa, bb, cc, idx in (*rows,):
        coefficient = aa*v[0]+bb*v[1]
        rhs = add(cc, constant(-aa*p[0]-bb*p[1]))
        if coefficient == 0:
            if rhs < ZERO:
                return None, tuple(i for i in (identity, idx) if i >= 0)
        else:
            t = scale(1/coefficient, rhs)
            if coefficient > 0 and (high is None or t < high):
                high, high_id = t, idx
            if coefficient < 0 and (low is None or t > low):
                low, low_id = t, idx
    for aa, bb, cc, idx in previous:
        if stats is not None:
            stats['boundary_rows'] += 1
        coefficient = aa*v[0]+bb*v[1]
        rhs = constant(cc-aa*p[0]-bb*p[1])
        if coefficient == 0:
            if rhs < ZERO:
                return None, tuple(i for i in (identity, idx) if i >= 0)
        else:
            t = scale(1/coefficient, rhs)
            if coefficient > 0 and (high is None or t < high):
                high, high_id = t, idx
            if coefficient < 0 and (low is None or t > low):
                low, low_id = t, idx
    check(low is not None and high is not None, 'formal box bounds every nonconstant line')
    if low > high:
        return None, tuple(i for i in (identity, low_id, high_id) if i >= 0)
    # Objective, x, then y, with no perturbation of the feasible constraints.
    direction = (dot(objective, v), v[0], v[1])
    t, endpoint_id = (high, high_id) if direction > (0, 0, 0) else (low, low_id)
    xy = tuple(add(constant(p[k]), scale(v[k], t)) for k in range(2))
    return xy, tuple(i for i in (identity, endpoint_id) if i >= 0)


def realize_symbolic(rows, objective, xy, witness=()):
    if xy is None:
        return dict(status='infeasible', witness=tuple(dict.fromkeys(witness)))
    slope = tuple(q[0] for q in xy)
    intercept = tuple(q[1] for q in xy)
    threshold = F(1)
    for a, b, c in rows:
        s, residual = a*slope[0]+b*slope[1], a*intercept[0]+b*intercept[1]-c
        check(s <= 0 and (s != 0 or residual <= 0), 'eventual feasibility')
        if s < 0:
            threshold = max(threshold, residual/(-s))
    M = ceil(threshold)
    p = tuple(intercept[k]+M*slope[k] for k in range(2))
    check(all(dot(row[:2], p) <= row[2] for row in rows), 'finite feasible realization')
    improvement = dot(objective, slope)
    check(improvement >= 0, 'nested boxes cannot decrease objective')
    status = 'unbounded' if improvement > 0 else 'optimal'
    result = dict(status=status, point=p, M=M, symbolic=xy, witness=tuple(dict.fromkeys(witness)))
    if status == 'unbounded':
        result['ray'] = slope
    else:
        result['value'] = dot(objective, p)
    return result


def seidel(raw_rows, raw_objective, rng=None, record_events=False, audit=False):
    rows, objective = constraints(raw_rows), point(raw_objective)
    rng = Random(0) if rng is None else rng
    order = list(range(len(rows)))
    rng.shuffle(order)
    indexed = [rows[i]+(i,) for i in order]
    xy = tuple((F(1 if coefficient >= 0 else -1), F(0)) for coefficient in objective)
    witness = ()
    events = []
    stats = dict(tests=0, repairs=0, boundary_rows=0)
    for j, row in enumerate(indexed):
        stats['tests'] += 1
        if not satisfies(xy, row):
            stats['repairs'] += 1
            before = xy
            xy, witness = boundary_lp(row, (indexed[k] for k in range(j)), objective, stats)
            if record_events:
                events.append(dict(position=j, constraint=row[3], before=before, after=xy,
                                   active_or_conflict=witness))
            if xy is None:
                break
        if audit:
            check(all(satisfies(xy, indexed[k]) for k in range(j+1)), 'processed constraints')
    return dict(result=realize_symbolic(rows, objective, xy, witness), order=order,
                stats=stats, events=events)


def enumerate_lp_symbolic(raw_rows, raw_objective):
    """Independent small-kernel: enumerate formal box polygon vertices, O((m+4)^3)."""
    rows, objective = constraints(raw_rows), point(raw_objective)
    symbolic_rows = [(a, b, constant(c), i) for i, (a, b, c) in enumerate(rows)] + box_rows()
    best = None
    for first, second in combinations(symbolic_rows, 2):
        a, b, rhs, _ = first
        d, e, other, _ = second
        determinant = a*e-b*d
        if not determinant:
            continue
        xy = (scale(1/determinant, add(scale(e, rhs), scale(-b, other))),
              scale(1/determinant, add(scale(a, other), scale(-d, rhs))))
        if all(affine_dot((aa, bb), xy) <= cc for aa, bb, cc, _ in symbolic_rows):
            if best is None or symbolic_key(xy, objective) > symbolic_key(best, objective):
                best = xy
    return best


def lp_certificate(raw_rows, raw_objective, result):
    """Additional deterministic certificate extraction, not part of Seidel's linear core.
    Uses the core's <=3 infeasibility witness or at most two active normal directions.
    Full active-pair search has O(n^2) arithmetic cost.
    """
    rows, objective = constraints(raw_rows), point(raw_objective)
    if result['status'] == 'infeasible':
        ids = result['witness']
        check(len(ids) <= 3 and enumerate_lp_symbolic([rows[i] for i in ids], objective) is None,
              'small infeasible witness')
        return dict(infeasible_rows=ids)
    p = result['point']
    if result['status'] == 'unbounded':
        ray = result['ray']
        check(all(dot(row[:2], ray) <= 0 for row in rows) and dot(objective, ray) > 0,
              'strict improving feasible ray')
        return dict(point=p, ray=ray)
    if objective == (0, 0):
        return dict(point=p, multipliers=())
    active = [i for i, (a, b, c) in enumerate(rows) if a*p[0]+b*p[1] == c]
    for i in active:
        normal = rows[i][:2]
        if cross(normal, objective) == 0 and dot(normal, objective) > 0:
            k = 0 if normal[0] else 1
            return dict(point=p, multipliers=((i, objective[k]/normal[k]),))
    for i, j in combinations(active, 2):
        u, v = rows[i][:2], rows[j][:2]
        determinant = cross(u, v)
        if determinant:
            alpha, beta = cross(objective, v)/determinant, cross(u, objective)/determinant
            if alpha >= 0 and beta >= 0:
                return dict(point=p, multipliers=tuple((k, w) for k, w in ((i, alpha), (j, beta)) if w))
    raise AssertionError('bounded optimum has nonnegative normal certificate')


def clarkson(raw_rows, raw_objective, rng=None, record_events=False, max_trials=None):
    """Clarkson-style weighted inner algorithm, iid WITH replacement each trial.
    LP-type dimension bound 3 includes infeasible planar systems. A trial budget
    reports UNKNOWN, never a guessed answer. Integer tickets avoid float rounding.
    """
    rows, objective = constraints(raw_rows), point(raw_objective)
    rng = Random(0) if rng is None else rng
    if max_trials is not None and (not isinstance(max_trials, int) or isinstance(max_trials, bool) or max_trials < 0):
        raise ValueError('trial budget must be a nonnegative integer')
    n, delta = len(rows), 3
    sample_size = 6*delta*delta
    if not rows:
        xy = enumerate_lp_symbolic([], objective)
        return dict(status='complete', symbolic=xy, result=realize_symbolic(rows, objective, xy),
                    weights=[], events=[], stats=dict(trials=0, accepted=0, rejected=0, tests=0))
    weights = [1]*n
    events = []
    stats = dict(trials=0, accepted=0, rejected=0, tests=0)
    while max_trials is None or stats['trials'] < max_trials:
        stats['trials'] += 1
        prefix, total = [], 0
        for weight in weights:
            total += weight
            prefix.append(total)
        # bisect_right picks the first cumulative weight strictly beyond ticket.
        from bisect import bisect_right
        draws = [bisect_right(prefix, rng.randrange(total)) for _ in range(sample_size)]
        # Boolean marker makes deduplication deterministic O(n+r), no hash assumption.
        marked = [False]*n
        for i in draws:
            marked[i] = True
        selected = [i for i in range(n) if marked[i]]
        xy = enumerate_lp_symbolic([rows[i] for i in selected], objective)
        violators = [] if xy is None else [i for i, row in enumerate(rows) if not satisfies(xy, row)]
        stats['tests'] += 0 if xy is None else n
        violating_weight = sum(weights[i] for i in violators)
        accepted = 3*delta*violating_weight <= total
        event = dict(draws=draws, selected=selected, total_weight=total,
                     violating_weight=violating_weight, violators=violators, accepted=accepted)
        if record_events:
            events.append(event)
        if not violators:
            witness = []
            if xy is None:
                witness = selected[:]
                # Deletion-only reduction; the subset kernel remains deterministic.
                for identity in selected:
                    candidate = [i for i in witness if i != identity]
                    if enumerate_lp_symbolic([rows[i] for i in candidate], objective) is None:
                        witness = candidate
                check(len(witness) <= delta, 'Helly-sized infeasibility certificate')
            return dict(status='complete', symbolic=xy,
                        result=realize_symbolic(rows, objective, xy, witness), weights=weights,
                        events=events, stats=stats)
        if accepted:
            stats['accepted'] += 1
            for i in violators:
                weights[i] *= 2
        else:
            stats['rejected'] += 1
    return dict(status='unknown', weights=weights, stats=stats, events=events)


def jsonify(value):
    if isinstance(value, F):
        return str(value)
    if isinstance(value, dict):
        return {str(k): jsonify(v) for k, v in value.items()}
    if isinstance(value, (list, tuple)):
        return [jsonify(v) for v in value]
    return value


def main():
    rng = Random(20261009)
    counts = dict(circles=0, linear_programs=0, certificates=0, clarkson=0)
    grid = list(product(range(-1, 2), repeat=2))
    for size in range(6):
        for ids in combinations(range(9), size):
            pts = [grid[i] for i in ids]
            expected = circle_enumeration(pts)
            for seed in range(3):
                actual = welzl(pts, Random(seed), audit=True)
                check(actual['circle'] == expected, 'grid MEC oracle')
                counts['circles'] += 1
    for _ in range(350):
        pts = [(rng.randrange(-6, 7), rng.randrange(-6, 7)) for _ in range(rng.randrange(12))]
        check(welzl(pts, rng, audit=True)['circle'] == circle_enumeration(pts), 'duplicate/random MEC')
        counts['circles'] += 1
    special = [([(1, 0, 0), (0, -1, -100)], (1, 0)), ([], (0, 0)), ([], (1, 0)), ([(1, 0, 0)], (1, 0)),
               ([(1, 0, 0), (-1, 0, 0)], (0, 1)), ([(0, 0, -1)], (0, 0)),
               ([(1, 0, 0), (-1, 0, -1)], (1, 1)),
               ([(1, 0, 0), (0, 1, 0), (-1, -1, -1)], (1, 2))]
    programs = special + [([(rng.randrange(-3, 4), rng.randrange(-3, 4), rng.randrange(-5, 6))
                            for _ in range(rng.randrange(9))],
                           (rng.randrange(-2, 3), rng.randrange(-2, 3))) for _ in range(600)]
    for rows, objective in programs:
        expected = enumerate_lp_symbolic(rows, objective)
        for seed in range(3):
            run = seidel(rows, objective, Random(seed), audit=True)
            result = run['result']
            check((None if result['status'] == 'infeasible' else result['symbolic']) == expected,
                  'symbolic LP independent vertex oracle')
            certificate = lp_certificate(rows, objective, result)
            if 'multipliers' in certificate:
                weights = certificate['multipliers']
                check(all(w >= 0 for i, w in weights), 'nonnegative multipliers')
                check(tuple(sum(w*F(rows[i][k]) for i, w in weights) for k in range(2)) == objective,
                      'dual normal sum')
                check(sum(w*F(rows[i][2]) for i, w in weights) == result['value'], 'dual bound attained')
            counts['linear_programs'] += 1
            counts['certificates'] += 1
    for rows, objective in programs[:30]:
        actual = clarkson(rows, objective, rng)
        check(actual['symbolic'] == enumerate_lp_symbolic(rows, objective), 'Clarkson full optimum')
        counts['clarkson'] += 1
    demo_points = [(0, 0), (4, 0), (1, 3), (1, 1), (4, 0), (2, 0)]
    demo_rows = [(1, 0, 4), (0, 1, 3), (-1, 0, 0), (0, -1, 0), (1, 1, 5), (2, 1, 8)]
    weighted_rows = [(1, 0, i) for i in range(128)]
    rejected_seed = None
    for seed in range(100000):
        rr = Random(seed)
        if min(rr.randrange(128) for _ in range(54)) >= 15:
            rejected_seed = seed
            break
    check(rejected_seed is not None, 'deterministic rejection example exists')
    weighted_demo = clarkson(weighted_rows, (1, 0), Random(rejected_seed), record_events=True)
    check(weighted_demo['stats']['rejected'] >= 1 and weighted_demo['stats']['accepted'] >= 1,
          'both weighted branches exercised')
    counts['clarkson'] += 1
    check(clarkson(weighted_rows, (1, 0), Random(0), max_trials=0)['status'] == 'unknown',
          'zero budget is not an answer')
    result = dict(status='PASS', counts=counts, weighted_seed=rejected_seed, weighted=weighted_demo,
                  circle=welzl(demo_points, Random(17), record_events=True, audit=True),
                  lp=seidel(demo_rows, (3, 2), Random(7), record_events=True, audit=True),
                  boundaries=[seidel(rows, c, Random(0)) for rows, c in special])
    target = Path(__file__).with_name('algorithms-small-basis-optimization-results.json')
    target.write_text(json.dumps(jsonify(result), ensure_ascii=False, indent=2)+'\n')
    print(json.dumps(jsonify(result), ensure_ascii=False, indent=2))


if __name__ == '__main__':
    main()