#!/usr/bin/env python3
"""Torn-verify for 011 — Measure, Not Memory. Pure Python, deterministic."""

import math
import random

FAILS = []
A, B, C, D = -1.4, 1.6, 1.0, 0.7


def check(name, condition, detail=""):
    print(("PASS" if condition else "FAIL") + "  " + name +
          (("  [" + detail + "]") if detail else ""))
    if not condition:
        FAILS.append(name)


def step(p):
    x, y = p
    return math.sin(A*y) + C*math.cos(A*x), math.sin(B*x) + D*math.cos(B*y)


def jac(p):
    x, y = p
    return ((-A*C*math.sin(A*x), A*math.cos(A*y)),
            (B*math.cos(B*x), -B*D*math.sin(B*y)))


def burn(p, n=1000):
    for _ in range(n):
        p = step(p)
    return p


p = burn((0.1, 0.0))
q0, q1 = (1.0, 0.0), (0.0, 1.0)
s0 = s1 = sum_log_det = 0.0
n_lyap = 200_000
for _ in range(n_lyap):
    j = jac(p)
    a = (j[0][0]*q0[0] + j[0][1]*q0[1],
         j[1][0]*q0[0] + j[1][1]*q0[1])
    b = (j[0][0]*q1[0] + j[0][1]*q1[1],
         j[1][0]*q1[0] + j[1][1]*q1[1])
    r0 = math.hypot(*a)
    q0 = (a[0]/r0, a[1]/r0)
    dot = q0[0]*b[0] + q0[1]*b[1]
    v1 = (b[0] - dot*q0[0], b[1] - dot*q0[1])
    r1 = math.hypot(*v1)
    q1 = (v1[0]/r1, v1[1]/r1)
    s0 += math.log(abs(r0))
    s1 += math.log(abs(r1))
    det = j[0][0]*j[1][1] - j[0][1]*j[1][0]
    sum_log_det += math.log(abs(det))
    p = step(p)

lambda1, lambda2 = s0/n_lyap, s1/n_lyap
log_det_mean = sum_log_det/n_lyap
check("the leading Lyapunov exponent is positive and near the recorded value",
      0.27 < lambda1 < 0.31, f"lambda1={lambda1:.6f}")
check("the QR exponent sum equals the orbit-average log absolute determinant",
      abs((lambda1 + lambda2) - log_det_mean) < 1e-11,
      f"sum={lambda1+lambda2:.12f}, avg log|det J|={log_det_mean:.12f}")
check("the map expands area on average rather than shrinking it",
      0.14 < log_det_mean < 0.19 and 1.15 < math.exp(log_det_mean) < 1.21,
      f"mean={log_det_mean:.6f}, factor={math.exp(log_det_mean):.6f}")

p = burn((0.1, 0.0))
twin = (p[0] + 1e-12, p[1])
fit = []
for n in range(1, 500):
    p, twin = step(p), step(twin)
    distance = math.dist(p, twin)
    if 1e-10 < distance < 1e-3:
        fit.append((n, math.log(distance)))
mx = sum(x for x, _ in fit)/len(fit)
my = sum(y for _, y in fit)/len(fit)
twin_slope = sum((x-mx)*(y-my) for x, y in fit) / sum((x-mx)**2 for x, _ in fit)
check("a true neighbouring orbit independently measures exponential separation",
      len(fit) > 30 and 0.24 < twin_slope < 0.36,
      f"fit points={len(fit)}, slope={twin_slope:.6f}")
check("the QR and twin instruments agree on the scale of the exponent",
      abs(lambda1 - twin_slope) < 0.08,
      f"QR={lambda1:.6f}, twin={twin_slope:.6f}")

grid, radius = 256, 2.2
seeds = (11, 2026, 777)
orbits = []
for seed in seeds:
    rng = random.Random(seed)
    orbits.append(burn((rng.uniform(-1, 1), rng.uniform(-1, 1))))
hist = [[0]*(grid*grid) for _ in seeds]
checkpoints = {}
first_distance = last_distance = 0.0

for n in range(1, 2_000_001):
    for k, p in enumerate(orbits):
        p = step(p)
        orbits[k] = p
        ix = int((p[0] + radius)/(2*radius)*grid)
        iy = int((p[1] + radius)/(2*radius)*grid)
        if 0 <= ix < grid and 0 <= iy < grid:
            hist[k][iy*grid + ix] += 1
    d = math.dist(orbits[0], orbits[1])
    if n <= 20_000:
        first_distance += d
    if n > 1_980_000:
        last_distance += d
    if n in (20_000, 200_000, 2_000_000):
        def l1(i, j):
            return sum(abs(a-b) for a, b in zip(hist[i], hist[j]))/n
        checkpoints[n] = (l1(0, 1), l1(0, 2), l1(1, 2))

first_distance /= 20_000
last_distance /= 20_000
check("two paths never converge pointwise despite sharing the attractor",
      1.35 < first_distance < 1.80 and 1.35 < last_distance < 1.80,
      f"first={first_distance:.3f}, last={last_distance:.3f}")

l20, l200, l2m = (checkpoints[n][0] for n in (20_000, 200_000, 2_000_000))
check("independent occupancy measures converge as the record grows",
      0.72 < l20 < 0.90 and 0.24 < l200 < 0.34 and 0.075 < l2m < 0.11,
      f"{l20:.3f} -> {l200:.3f} -> {l2m:.3f}")
check("each tenfold of patience shrinks histogram disagreement near sqrt(10)",
      2.5 < l20/l200 < 3.8 and 2.5 < l200/l2m < 3.8,
      f"ratios={l20/l200:.3f}, {l200/l2m:.3f}")
check("a third seed converges to the same measured picture",
      all(0.075 < value < 0.11 for value in checkpoints[2_000_000]),
      "pairwise=" + ", ".join(f"{v:.3f}" for v in checkpoints[2_000_000]))

print()
if FAILS:
    print(f"TORN-VERIFY: {len(FAILS)} recorded claim(s) no longer hold: {FAILS}")
    raise SystemExit(1)
print("TORN-VERIFY: all recorded claims for 011 still hold.")