#!/usr/bin/env python3
"""check.py MODE file.raw - replay texmap.c's triangle setup in Python (same doubles, same
rounding) and compare every pixel the card drew with the texel the fixed-point model and
the exact (float) model predict. Texture 'idtex': texel value = 0x8000 | t << 6 | s.
Claude Opus 5.5, 2026-10-14."""
import sys, math
W, H = 640, 480
BG = 0x20e4

def edge_at(x0, y0, x1, y1, y):
    dx, dy = x1 - x0, y1 - y0
    if dx >= 0: sign, adx, ex = 1, dx, 0
    else: sign, adx, ex = -1, -dx, 1 - dy
    step = (adx // dy) * sign; adx %= dy
    xx = x0
    for k in range(y0, y):
        xx += step; ex += adx
        if ex > 0: xx += sign; ex -= dy
    return xx

fx20 = lambda v: math.floor(v * 1048576.0 + 0.5)
fx16 = lambda v: math.floor(v * 65536.0 + 0.5)

def trap(ya, yb, L, R, p0, dx, dy, px, py, out, persp):
    for y in range(ya, yb):
        xl0 = edge_at(*L, ya)
        xl, xr = edge_at(*L, y), edge_at(*R, y)
        v = [p0[i] + (xl0 - px) * dx[i] + (ya - py) * dy[i] for i in range(3)]
        tm = [fx20(dx[0]), fx20(dy[0]), fx20(dx[1]), fx20(dy[1]), fx16(dx[2]), fx16(dy[2]),
              fx20(v[0]), fx20(v[1]), fx16(v[2])]
        for x in range(xl, xr):
            S = tm[6] + (y - ya) * tm[1] + (x - xl0) * tm[0]
            T = tm[7] + (y - ya) * tm[3] + (x - xl0) * tm[2]
            Q = tm[8] + (y - ya) * tm[5] + (x - xl0) * tm[4]
            # exact: plane value at the pixel, divided
            es = p0[0] + (x - px) * dx[0] + (y - py) * dy[0]
            et = p0[1] + (x - px) * dx[1] + (y - py) * dy[1]
            eq = p0[2] + (x - px) * dx[2] + (y - py) * dy[2]
            if persp: fs, ft = S / 1048576 / (Q / 65536), T / 1048576 / (Q / 65536); es, et = es / eq, et / eq
            else: fs, ft = S / 1048576, T / 1048576
            out[(x, y)] = (fs, ft, es, et)

def tri(a, b, c, out, persp, qscale=1.0):
    a, b, c = sorted([a, b, c], key=lambda v: v[1])
    (ax, ay, *pa), (bx, by, *pb), (cx, cy, *pc) = a, b, c
    pa, pb, pc = list(pa), list(pb), list(pc)
    if ay == cy: return
    det = float(bx - ax) * (cy - ay) - float(cx - ax) * (by - ay)
    if det == 0: return
    if qscale > 0:
        k = qscale / max(pa[2], pb[2], pc[2])
        pa, pb, pc = [p * k for p in pa], [p * k for p in pb], [p * k for p in pc]
    dx, dy = [], []
    for i in range(3):
        e1, e2 = pb[i] - pa[i], pc[i] - pa[i]
        dx.append((e1 * (cy - ay) - e2 * (by - ay)) / det)
        dy.append((e2 * (bx - ax) - e1 * (cx - ax)) / det)
    if det > 0:
        trap(ay, by, (ax, ay, cx, cy), (ax, ay, bx, by), pa, dx, dy, ax, ay, out, persp)
        trap(by, cy, (ax, ay, cx, cy), (bx, by, cx, cy), pa, dx, dy, ax, ay, out, persp)
    else:
        trap(ay, by, (ax, ay, bx, by), (ax, ay, cx, cy), pa, dx, dy, ax, ay, out, persp)
        trap(by, cy, (bx, by, cx, cy), (ax, ay, cx, cy), pa, dx, dy, ax, ay, out, persp)

def quad(a, b, c, d, out, persp):
    tri(a, b, c, out, persp); tri(a, c, d, out, persp)

def PV(X, Y, Z, s, t):
    w = 1 / Z
    return (math.floor(320 + 400 * X * w + 0.5), math.floor(240 - 400 * Y * w + 0.5), s * w, t * w, w)

mode, fn = sys.argv[1], sys.argv[2]
import os
BIAS = float(os.environ.get('BIAS', '0'))
data = open(fn, 'rb').read()
pix = lambda x, y: data[2 * (y * W + x)] | data[2 * (y * W + x) + 1] << 8
exp = {}
if mode == 'affine':
    V = lambda x, y, s, t: (x, y, s, t, 1.0)
    quad(V(16, 16, 0, 0), V(272, 16, 1, 0), V(272, 272, 1, 1), V(16, 272, 0, 1), exp, False)
    quad(V(450, 20, 0, 0), V(620, 200, 2, 0), V(440, 400, 2, 2), V(290, 210, 0, 2), exp, False)
    region = lambda x, y: True
elif mode == 'persp':
    quad(PV(-2, -1, 6, 0, 1), PV(2, -1, 6, 1, 1), PV(2, -1, 1.5, 1, 0), PV(-2, -1, 1.5, 0, 0), exp, True)
    region = lambda x, y: y > 240
exp = {k: v for k, v in exp.items() if 0 <= k[1] < H and 0 <= k[0] < W}
drawn = {(x, y) for y in range(H) for x in range(W) if pix(x, y) != BG and region(x, y)}
cov_extra = drawn - set(exp); cov_missing = set(exp) - drawn
print("pixels expected %d, drawn %d, extra %d, missing %d" % (len(exp), len(drawn), len(cov_extra), len(cov_missing)))
fixed_ok = exact_ok = near = 0; bad = []
for (x, y), (fs, ft, es, et) in exp.items():
    v = pix(x, y)
    hs, ht = v & 63, (v >> 6) & 63
    ps, pt = math.floor(fs * 64 + BIAS) & 63, math.floor(ft * 64 + BIAS) & 63
    qs, qt = math.floor(es * 64) & 63, math.floor(et * 64) & 63
    if v >> 12 != 8: bad.append((x, y, hex(v), ps, pt)); continue
    fixed_ok += (hs, ht) == (ps, pt)
    exact_ok += (hs, ht) == (qs, qt)
    ds, dt = (hs - qs + 32) % 64 - 32, (ht - qt + 32) % 64 - 32
    near += abs(ds) <= 1 and abs(dt) <= 1
    if (hs, ht) != (ps, pt) and len(bad) < 15: bad.append((x, y, (hs, ht), (ps, pt), (qs, qt), round(fs * 64, 3), round(ft * 64, 3)))
n = len(exp)
print("match fixed-point model %d/%d (%.2f%%), exact model %d/%d (%.2f%%), within 1 texel of exact %.2f%%"
      % (fixed_ok, n, 100 * fixed_ok / n, exact_ok, n, 100 * exact_ok / n, 100 * near / n))
for b in bad[:15]: print(b)

# rounding analysis: find the bias R such that texel = (S + R) >> 14 (s, t separately)
if mode == 'affine':
    lo = {0: -10**9, 1: -10**9}; hi = {0: 10**9, 1: 10**9}
    for (x, y), (fs, ft, es, et) in exp.items():
        v = pix(x, y); h = (v & 63, (v >> 6) & 63)
        for k, f in ((0, fs), (1, ft)):
            S = round(f * 1048576)            # raw 12.20 value
            # need floor((S + R) / 16384) mod 64 == h  ->  R in [h*16384 - S, (h+1)*16384 - S) mod 2^20
            base = math.floor(S / 16384)
            for cand in (base - 1, base, base + 1):
                if cand & 63 == h[k]:
                    lo[k] = max(lo[k], cand * 16384 - S); hi[k] = min(hi[k], (cand + 1) * 16384 - S - 1)
                    break
    print("bias range s: [%d, %d]  t: [%d, %d]  (units of 2^-20; one texel = 16384)" % (lo[0], hi[0], lo[1], hi[1]))