#!/usr/bin/env python3
"""Paper VIII exact finite certificates. Python 3.10+, standard library only.

Run: python3 verify_paper8.py --out fresh_certificate_directory
The destination must not exist. No result directory is published until all
checks and rejection controls succeed. This is an executable arithmetic
checker, not a proof-assistant verification of the analytic theorems.
"""
from __future__ import annotations

import argparse
from collections import Counter
from copy import deepcopy
from dataclasses import dataclass
from fractions import Fraction as F
from hashlib import sha256
import json
from math import gcd, isqrt
from pathlib import Path
import tempfile

Q = 339949771344778
S = 10**15
A = 3 * Q - S
N = 32767
DSTAR = 42189065
PERIOD = Q // gcd(Q, S)
ENC_DEN = 10**30
ENC_LO = 339949771344778929576702809321
EXPECTED_DEPTHS = [1, 2, 6, 18, 53, 156, 459, 1350, 3971, 11681, 15071]
EXPECTED_B = [0, 2, 8, 26, 79, 235, 694, 2044, 6015, 17696,
              52057, 153134, 450463, 1325089, 3897899, 11466105,
              33728823, 99217081]


def require(ok, message):
    if not ok:
        raise ValueError(message)


def integer(x):
    return type(x) is int


def canonical(x):
    return (json.dumps(x, sort_keys=True, separators=(",", ":"),
                       ensure_ascii=True) + "\n").encode("ascii")


def frac(x):
    x = F(x)
    return [x.numerator, x.denominator]


def step(d):
    return Q * d // S


def left(n):
    return (S * n + Q - 1) // Q


def fiber_size(n):
    return left(n + 1) - left(n)


def make_row(seed):
    d, residuals = seed, []
    while d:
        d, r = divmod(Q * d, S)
        residuals.append(r)
    return {"seed": seed, "terminal": 0, "remainders": residuals}


def check_row(row, seed):
    require(type(row) is dict and set(row) == {"seed", "terminal", "remainders"},
            "row schema")
    require(integer(row["seed"]) and row["seed"] == seed, "row seed")
    require(integer(row["terminal"]) and row["terminal"] == 0, "terminal")
    rs = row["remainders"]
    require(type(rs) is list and len(rs) <= 10, "ordered remainder list")
    require(all(integer(r) and 0 <= r < S for r in rs), "remainder bounds")
    d = 0
    for r in reversed(rs):
        previous, remainder = divmod(S * d + r, Q)
        require(remainder == 0 and previous > d, "reverse divisibility or descent")
        require(Q * previous // S == d, "forward quotient")
        d = previous
    require(d == seed, "reconstructed seed")
    return len(rs)


def check_block(rows, job):
    require(integer(job) and 0 <= job < 32, "job index")
    require(type(rows) is list and len(rows) == 1024, "block coverage")
    return [check_row(row, 1024 * job + i) for i, row in enumerate(rows)]


def atan_bounds(b):
    lower = sum((F((-1)**j, (2*j + 1) * b**(2*j + 1))
                 for j in range(64)), F(0))
    return lower, lower + F(1, 129 * b**129)


def enclosure():
    scale = 10**80
    root = isqrt(5 * scale**2)
    require(root**2 < 5 * scale**2 < (root + 1)**2, "sqrt enclosure")
    phi_lo, phi_hi = (1 + F(root, scale))/2, (1 + F(root + 1, scale))/2
    a, b = atan_bounds(5)
    c, d = atan_bounds(239)
    pi_lo, pi_hi = 16*a - 4*d, 16*b - 4*c
    lam_lo, lam_hi = phi_lo/(phi_lo + pi_hi), phi_hi/(phi_hi + pi_lo)
    lo, hi = F(ENC_LO, ENC_DEN), F(ENC_LO + 1, ENC_DEN)
    require(F(Q, S) < lo < lam_lo < lam_hi < hi, "lambda enclosure")
    return lo, hi


def floor_sum(n, m, a, b=0):
    value = 0
    while True:
        u, a = divmod(a, m)
        v, b = divmod(b, m)
        value += u*n*(n-1)//2 + v*n
        k, r = divmod(a*n + b, m)
        if k == 0:
            return value
        n, m, a, b = k, a, m, r


def arithmetic():
    lo, hi = enclosure()
    lower = F(846547, 2490212)
    upper = F(13495616, 39698853)
    crossing = F(14342163, DSTAR)
    require(13495616*2490212 - 846547*39698853 == 1, "Farey determinant")
    require(crossing == F(846547 + 13495616, 2490212 + 39698853), "mediant")
    require(lower < F(Q, S) < crossing < lo < hi < upper, "Farey ordering")
    sums = [floor_sum(DSTAR, S, Q), floor_sum(DSTAR, ENC_DEN, ENC_LO+1)]
    require(sums == [302541195258184]*2, "floor sum agreement")
    require(step(DSTAR) == 14342162, "rational first crossing")
    require(lo*DSTAR // 1 == hi*DSTAR // 1 == 14342163, "real first crossing")
    require(lo*(DSTAR+1) // 1 == hi*(DSTAR+1) // 1 == step(DSTAR+1),
            "agreement resumes at the next state")
    require(gcd(Q, S) == 2 and A == 19849314034334, "frozen constants")
    require(((PERIOD-1)*S) % Q == A and (PERIOD*S) % Q == 0, "remote residues")
    require(fiber_size(PERIOD-1) == 2 and fiber_size(PERIOD) == 3, "remote endpoints")
    require(((A//2)*Q + A-1)//A - 1 == PERIOD-1, "Beatty endpoint")
    for n in range(50001):
        r = n*S % Q
        predicate = 2 if 0 < r <= A else 3
        mechanical = 3 - ((n+1)*A//Q - n*A//Q)
        require(fiber_size(n) == predicate == mechanical, "fiber regression")
    bs = [0]
    for _ in range(17):
        bs.append((S*(bs[-1]+1)-1)//Q)
    require(bs == EXPECTED_B, "depth endpoints")
    low_margin, high_margin = F(0), F(1)
    for d in range(1, N+1):
        y = step(d)
        low_margin = max(low_margin, F(y, d))
        high_margin = min(high_margin, F(y+1, d))
    require((low_margin, high_margin) == (F(4196, 12343), F(7039, 20706)),
            "maximal finite-domain parameter interval")
    require(F(step(12343),12343) == F(step(24686),24686) == low_margin,
            "lower interval witnesses")
    require(F(step(20706)+1,20706) == high_margin, "upper interval witness")
    counts = Counter(step(d) for d in range(N+1))
    fiber_hist = Counter(counts.values())
    require(len(counts) == 11140 and fiber_hist == {1: 1, 2: 650, 3: 10489},
            "restricted fiber histogram")
    require([left(11139), left(11140)-1] == [32767, 32769], "clipped fiber")
    overflow = (2**63-1)//Q + 1
    require(overflow == 27132, "signed 64-bit overflow boundary")
    return {"lambda_enclosure": [frac(lo), frac(hi)], "depth_endpoints": bs,
            "first_disagreement": DSTAR, "floor_sums": sums,
            "parameter_interval": [frac(low_margin), frac(high_margin)],
            "image_size": len(counts), "restricted_fiber_histogram": dict(fiber_hist),
            "rational_period": PERIOD, "period_size_two_count": A//2,
            "signed_64_first_overflow": overflow,
            "fiber_regression_bases": 50001,
            "remote_boundary_bases": [PERIOD-1, PERIOD]}


@dataclass(frozen=True)
class K5:
    a: F = F(0)
    b: F = F(0)

    def __add__(self, other):
        other = other if isinstance(other, K5) else K5(F(other))
        return K5(self.a+other.a, self.b+other.b)

    __radd__ = __add__

    def __neg__(self):
        return K5(-self.a, -self.b)

    def __sub__(self, other):
        return self + (-other if isinstance(other, K5) else -F(other))

    def __mul__(self, other):
        other = other if isinstance(other, K5) else K5(F(other))
        return K5(self.a*other.a+5*self.b*other.b,
                  self.a*other.b+self.b*other.a)

    __rmul__ = __mul__

    def __truediv__(self, other):
        other = other if isinstance(other, K5) else K5(F(other))
        den = other.a**2 - 5*other.b**2
        require(den != 0, "field divisor")
        return self * K5(other.a/den, -other.b/den)

    def __pow__(self, n):
        require(integer(n) and n >= 0, "field exponent")
        result = K5(F(1))
        for _ in range(n):
            result = result*self
        return result

    def record(self):
        return [frac(self.a), frac(self.b)]


def odd_factorial(n):
    value = 1
    for k in range(1, n+1, 2):
        value *= k
    return value


def cubature():
    phi = K5(F(1, 2), F(1, 2))
    zero, one = K5(), K5(F(1))
    vertices = []
    for i in [-1, 1]:
        for j in [-1, 1]:
            vertices.extend([(zero, i*one, j*phi),
                             (i*one, j*phi, zero), (j*phi, zero, i*one)])
    radius2 = phi + 2
    require(len(set(vertices)) == 12, "vertex inventory")
    require(all(sum((x*x for x in v), zero) == radius2 for v in vertices), "node norms")
    moments = []
    for total in range(6):
        for i in range(total+1):
            for j in range(total-i+1):
                k = total-i-j
                raw = sum((v[0]**i * v[1]**j * v[2]**k for v in vertices), zero)/12
                expected = F(0) if any(x % 2 for x in (i,j,k)) else F(
                    odd_factorial(i-1)*odd_factorial(j-1)*odd_factorial(k-1),
                    odd_factorial(total+1))
                if total % 2:
                    require(raw == zero and expected == 0, "odd moment")
                    value = zero
                else:
                    value = raw / radius2**(total//2)
                    require(value == K5(expected), "even moment")
                moments.append({"powers": [i,j,k], "value": value.record()})
    sixth = sum((v[0]**6 for v in vertices), zero)/(12*radius2**3)
    require(len(moments) == 56 and sixth == K5(F(2,15)), "degree inventory and obstruction")
    return {"nodes": 12, "weight_per_node": [1,12], "degree": 5,
            "monomials_checked": 56, "moments": moments,
            "sixth_moment": sixth.record(), "sixth_spherical_moment": [1,7]}


def rejection_controls():
    cases = []
    base = make_row(399)
    for key, value in [("seed", True), ("seed", 398), ("terminal", 1),
                       ("terminal", False), ("remainders", []),
                       ("remainders", tuple(base["remainders"]))]:
        row = deepcopy(base)
        row[key] = value
        cases.append(("invalid_"+key+"_"+str(len(cases)), lambda row=row: check_row(row,399)))
    row = deepcopy(base)
    del row["remainders"]
    cases.append(("missing_field", lambda: check_row(row,399)))
    for label, transform in [
        ("changed_remainder", lambda r: r.__setitem__(0,r[0]+1)),
        ("reordered_remainders", lambda r: r.reverse()),
        ("upper_remainder_endpoint", lambda r: r.__setitem__(0,S)),
        ("negative_remainder", lambda r: r.__setitem__(0,-1)),
        ("boolean_remainder", lambda r: r.__setitem__(0,True)),
        ("trailing_zero_step", lambda r: r.append(0))]:
        bad = deepcopy(base)
        transform(bad["remainders"])
        cases.append((label, lambda bad=bad: check_row(bad,399)))
    block = [make_row(d) for d in range(1024)]
    cases.extend([
        ("empty_block", lambda: check_block([],0)),
        ("missing_seed", lambda: check_block(block[:-1],0)),
        ("duplicate_seed", lambda: check_block([block[0]]+block[:-1],0)),
        ("wrong_block_index", lambda: check_block(block,1)),
    ])
    for label, action in cases:
        try:
            action()
        except ValueError:
            continue
        raise ValueError("negative control accepted: "+label)
    return [label for label, _ in cases]


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--out", type=Path, required=True)
    args = parser.parse_args()
    destination = args.out.resolve()
    require(not destination.exists(), "destination must not exist")
    destination.parent.mkdir(parents=True, exist_ok=True)
    arithmetic_result = arithmetic()
    cubature_result = cubature()
    controls = rejection_controls()
    with tempfile.TemporaryDirectory(prefix="paper8_staging_", dir=destination.parent) as temp:
        staging = Path(temp)/"certificate"
        staging.mkdir()
        jobs, histogram, unique_images = [], Counter(), Counter()
        for job in range(32):
            rows = [make_row(d) for d in range(1024*job,1024*(job+1))]
            depths = check_block(rows,job)
            data = b"".join(canonical(row) for row in rows)
            require(check_block([json.loads(line) for line in data.splitlines()],job) == depths,
                    "serialized block replay")
            file_name = f"block_{job:02d}.jsonl"
            (staging/file_name).write_bytes(data)
            images = len({step(row["seed"]) for row in rows})
            jobs.append({"index":job,"start":1024*job,"end":1024*(job+1)-1,
                         "rows":1024,"images":images,"file":file_name,
                         "sha256":sha256(data).hexdigest()})
            histogram.update(depths)
            unique_images.update([images])
        require([histogram[i] for i in range(11)] == EXPECTED_DEPTHS, "depth histogram")
        require(unique_images == {348:9,349:23}, "block image counts")
        report = {"schema":"OMEGA_PAPER_VIII_CERTIFICATE/1.0", "result":"PASS",
                  "map":{"Q":Q,"S":S},"domain":{"first":0,"last":N},
                  "checker_sha256":sha256(Path(__file__).read_bytes()).hexdigest(),
                  "arithmetic":arithmetic_result,"cubature":cubature_result,
                  "depth_histogram":EXPECTED_DEPTHS,
                  "ledger_rows":N+1,"ledger_steps":sum(k*v for k,v in histogram.items()),
                  "rejected_controls":controls,"blocks":jobs}
        payload = canonical(report)
        (staging/"CERTIFICATE.json").write_bytes(payload)
        staging.rename(destination)
    print("RESULT=PASS")
    print("CERTIFICATE_SHA256="+sha256(payload).hexdigest())
    print("ROWS=32768; MONOMIALS=56; REJECTED_CONTROLS="+str(len(controls)))


if __name__ == "__main__":
    main()