# verify_dist.py
# Author: David Niedbala Giraudin, ORCID 0009-0009-1526-1178
# Technical assistance: Claude (Anthropic)
#
# Checks, from the certificate DIST_trefoil_*.csv (same folder), the statement
#     there is a closed polygon of knot type trefoil whose distortion is at most D,
# where distortion = sup over pairs of distinct points of (shorter arc length)/(distance).
#   1. P0 = an explicit polygon on the boundary of a square solid torus (R = 8, r = 5,
#      16 cut segments, each split in 6), of slope (2,3): a trefoil (recomputed exactly).
#   2. A list of integer polygons X_0 = P0, ..., X_K (unit 1/U); every linear step is subdivided
#      until each sub-step is certified to pass only through embedded polygons.  Then every edge
#      of X_K is split at its midpoint (exact in the unit 1/(2U)) and a second list of polygons
#      continues the path in that unit.  Hence the last polygon F is a trefoil.
#   3. The distortion of F is at most D, checked exactly (rational arithmetic) for every pair
#      of edges (cert_dist.py).
# No argument is needed.  Plain Python 3 with numpy.
import os, sys, time
import numpy as np
from fractions import Fraction as Fr
HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, HERE)
from certlib import step_ok
from p0 import p0_polygon
import cert_dist

D_CLAIM = Fr(6920, 1000)


def load_cert(prefix):
    """Read the certificate stored as CSV files <prefix>_<array>.csv (exact integers);
    the first line of each file gives the shape of the array."""
    import glob as _glob
    out = {}
    for fn in sorted(_glob.glob(prefix + "_*.csv")):
        key = fn[len(prefix) + 1:-4]
        with open(fn) as fh:
            head = fh.readline()
        assert head.startswith("# shape "), fn
        shape = tuple(int(x) for x in head[len("# shape "):].split(","))
        arr = np.loadtxt(fn, dtype=np.int64, delimiter=",", comments="#", ndmin=2)
        out[key] = arr.reshape(shape)
    return out

def say(*a): print(*a, flush=True)

def check_path(states, U, label):
    t0 = time.time(); nsub = 0
    for k in range(len(states) - 1):
        A = states[k]; B = states[k + 1]
        disp = np.sqrt(((B - A) ** 2).sum(1).astype(np.float64)).max()
        m = max(1, int(np.ceil(disp / (0.08 * U))))
        while True:
            C = [A + ((B - A) * j) // m for j in range(m)] + [B]
            ok = True
            for j in range(m):
                good, info = step_ok(C[j], C[j + 1])
                if not good: ok = False; break
            if ok: break
            m *= 2
            if m > 4096:
                say("REJECTED: step", k, "of", label, "cannot be certified", info); return False
        nsub += m
    say("  %s: %d steps, %d sub-steps, all embedded (%.0f s)" % (label, len(states) - 1, nsub, time.time() - t0))
    return True

def main():
    t_start = time.time()
    z = load_cert(os.path.join(HERE, "DIST_trefoil"))
    p, q, R, r, k, U = [int(v) for v in z["meta"]]
    if (p, q, R, r, k) != (2, 3, 8, 5, 6):
        say("REJECTED: torus parameters differ from those of the note"); say("VERDICT: REJECTED"); return
    S = z["states"].astype(np.int64); SR = z["states_ref"].astype(np.int64)
    P0, _ = p0_polygon(p, q, R, r, k)
    w = [[c * U for c in v] for v in P0]
    if not all(x.denominator == 1 for v in w for x in v):
        say("REJECTED: unit"); say("VERDICT: REJECTED"); return
    if not np.array_equal(np.array([[int(x) for x in v] for v in w], dtype=np.int64), S[0]):
        say("REJECTED: first state is not P0"); say("VERDICT: REJECTED"); return
    say("1. first state equals the polyhedral (2,3) torus-knot polygon (trefoil), %d vertices" % len(S[0]))
    if not check_path(list(S), U, "2a. path"):
        say("VERDICT: REJECTED"); return
    X = S[-1]; D0 = np.empty((2 * len(X), 3), dtype=np.int64)
    D0[0::2] = 2 * X; D0[1::2] = X + np.roll(X, -1, 0)
    if not check_path([D0] + list(SR), 2 * U, "2b. path after midpoint refinement"):
        say("VERDICT: REJECTED"); return
    F = SR[-1]
    say("   final polygon: %d vertices, a trefoil" % len(F))
    say("3. distortion bound, all pairs of edges, exact rational arithmetic ...")
    ok, info = cert_dist.certify([tuple(int(c) for c in v) for v in F], D_CLAIM)
    if not ok:
        say("REJECTED:", info); say("VERDICT: REJECTED"); return
    say("   %d boxes checked" % info["boxes"])
    say("VERDICT: CERTIFIED.  The trefoil has a polygonal realization with distortion <= %s = %.3f"
        % (D_CLAIM, float(D_CLAIM)))
    say("time %.0f s" % (time.time() - t_start))

if __name__ == "__main__":
    main()