# verify_C22a.py
# Author: David Niedbala Giraudin, ORCID 0009-0009-1526-1178
# Technical assistance: Claude (Anthropic)
#
# Checks, with exact integer / rational arithmetic (numpy is used only on exact
# integers; floating point serves only to choose how finely a step is subdivided,
# which does not affect the validity of the checks), the statement
#     L(T(p,q)) <= L_up,   hence   C_22a <= L_up / Cr(T(p,q))**(3/4),
# for every certificate C22a_T*_meta.csv (+ _states, _final, ...) in the same folder:
#   1. P0 = an explicit polygon on the boundary of a square solid torus, of
#      slope (p,q); it is a (p,q) torus knot (recomputed here exactly).
#   2. A list of integer polygons X_0 = P0, X_1, ..., X_K; every linear step
#      X_k -> X_{k+1} is subdivided until each sub-step is certified to pass
#      only through embedded polygons.  Hence X_K is a (p,q) torus knot.
#   3. The scaled final polygon F (also joined to X_K by a certified step),
#      smoothed by unit-radius circular arcs at its vertices, has curvature <= 1
#      and every pair of points at arc distance >= pi is at distance >= 2;
#      hence its thickness (unit radius) is >= 1, and its length is <= L_up.
# No argument is needed.  Plain Python 3 with numpy.

import os, sys, time
import numpy as np
from fractions import Fraction as Fr
from math import isqrt, gcd

HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, HERE)
from certlib import step_ok, thickness_cert, grid_pairs
from p0 import p0_polygon

U_GLOBAL = 1

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

def check_path(states, 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_GLOBAL))))
        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, "cannot be certified", info); return False
        nsub += m
        if k % 50 == 0:
            say("  %s step %d/%d certified (%d sub-steps so far, %.0f s)" % (label, k, len(states) - 1, nsub, time.time() - t0))
    say("  %s: %d steps, %d sub-steps, all embedded" % (label, len(states) - 1, nsub))
    return True

def cr34_lo(c):
    # rational lower bound for c**(3/4): largest x = a/2**20 with x**4 <= c**3
    lo, hi = 0, c << 20
    while hi - lo > 1:
        mid = (lo + hi) // 2
        if mid ** 4 <= (c ** 3) << 80: lo = mid
        else: hi = mid
    return Fr(lo, 1 << 20)


def load_cert(prefix):
    """Read a certificate stored as CSV files <prefix>_<array>.csv (exact integers);
    the first line of each file gives the shape of the array.  For a sequence of polygons
    marked "delta", the file holds the first polygon and then the differences between
    consecutive polygons, and the polygons are recovered by cumulative sums."""
    import glob as _glob
    import re as _re
    out = {}
    groups = {}                        # a large array may be split into parts <key>_p01, _p02, ...
    for fn in sorted(_glob.glob(prefix + "_*.csv")):
        key = fn[len(prefix) + 1:-4]
        mm = _re.fullmatch(r"(.+)_p(\d\d)", key)
        groups.setdefault(mm.group(1) if mm else key, []).append(fn)
    for key, fns in groups.items():
        with open(fns[0]) as fh:
            head = fh.readline()
        delta = head.startswith("# delta shape ")
        assert delta or head.startswith("# shape "), fns[0]
        shape = tuple(int(x) for x in head.split("shape ")[1].split(","))
        arr = np.concatenate([np.loadtxt(f, dtype=np.int64, delimiter=",", comments="#", ndmin=2)
                              for f in fns]).reshape(shape)
        if delta:                    # first polygon, then differences between consecutive polygons
            arr = np.cumsum(arr, axis=0)
        out[key] = arr
    return out

def main(CERT):
    z = load_cert(CERT)
    global U_GLOBAL
    p, q, U, Rr, rr, kk = [int(v) for v in z["meta"]]
    U_GLOBAL = U
    states = z["states"].astype(np.int64)
    F = z["final"].astype(np.int64)
    Cr = min(p * (q - 1), q * (p - 1))
    if not (2 <= p < q and gcd(p, q) == 1 and Rr == 2 * q + 6 and rr == q + 2 and kk == 16):
        say("REJECTED: parameters (p, q, R, r, k) differ from those of the note"); say("VERDICT: REJECTED"); return
    say("certificate: T(%d,%d), crossing number %d, unit 1/%d, %d states of %d vertices"
        % (p, q, Cr, U, len(states), states.shape[1]))
    # 1. exact P0
    P0, _ = p0_polygon(p, q, Rr, rr, kk)
    P0i = []
    for v in P0:
        w = [c * U for c in v]
        assert all(x.denominator == 1 for x in w), "unit does not clear denominators"
        P0i.append([int(x) for x in w])
    P0i = np.array(P0i, dtype=np.int64)
    if not np.array_equal(P0i, states[0]):
        say("REJECTED: first state is not P0"); say("VERDICT: REJECTED"); return
    say("1. first state equals the polyhedral (%d,%d) torus-knot polygon" % (p, q))
    # 2. path (optionally followed by an exact refinement: every edge is split at its
    #    midpoint, coordinates expressed in the unit 1/(2U); the polygon is unchanged)
    if "states_ref" in z:
        if not check_path(list(states), "2a. path"):
            say("VERDICT: REJECTED"); return
        Xl = states[-1]
        D0 = np.empty((2 * len(Xl), 3), dtype=np.int64)
        D0[0::2] = 2 * Xl; D0[1::2] = Xl + np.roll(Xl, -1, 0)
        say("   refinement: %d -> %d vertices (same polygon, unit 1/%d)" % (len(Xl), len(D0), 2 * U))
        U = 2 * U
        globals()["U_GLOBAL"] = U
        if not check_path([D0] + list(z["states_ref"].astype(np.int64)) + [F], "2b. path"):
            say("VERDICT: REJECTED"); return
    else:
        if not check_path(list(states) + [F], "2. path"):
            say("VERDICT: REJECTED"); return
    # 3. thickness and length of the smoothed final polygon
    ok, L = thickness_cert(F, U)
    if not ok:
        say("REJECTED:", L); say("VERDICT: REJECTED"); return
    c34 = cr34_lo(Cr)
    alpha = L / c34
    say("3. smoothed final curve: thickness >= 1, length <= %.6f" % float(L))
    say("VERDICT: CERTIFIED.  L(T(%d,%d)) <= %.6f and C_22a <= %.6f" % (p, q, float(L), float(alpha) + 1e-6))

if __name__ == "__main__":
    import glob
    files = sorted(f[:-len("_meta.csv")] for f in glob.glob(os.path.join(HERE, "C22a_T*_meta.csv")))
    if not files:
        say("no certificate C22a_T*_meta.csv found next to this script")
    for f in files:
        say("=" * 70); say(os.path.basename(f)); t_start = time.time()
        main(f)
        say("time %.0f s" % (time.time() - t_start))