"""Exact verification of the Directed Steiner Tree certificates (prop:fa-dst-certs).

Reads dst_certificates.json: certificates (y, e1, e2, Z) with rational entries, and rational weights
lambda_j >= 0 with sum 1.  With R = conv(Z) + R_{>=0}^2 and
    (a, b) o (a', b') = (a (a' + b'), b (a'/2 + b')),
it checks in exact rational arithmetic (fractions.Fraction):
    (C0) every z in Z is a nonnegative point with positive first coordinate, y >= 0, e1, e2 > 0;
    (C1) z o z' in R for all z, z' in Z;
    (C2) (1, 0)/e2 in R;
    (C3) ((y, 1) o z)/e1 in R for all z in Z.
Membership P in R is certified by a witness: a point z of Z, or two points z, z' of Z and a rational
mu in [0, 1], with mu z + (1 - mu) z' <= P coordinatewise.  (In the plane this search is complete.)
It also reports c_Z = min_z z_1 exactly.

Rates.  For certificate j put L_j = log2(1 + y_j) and
    a_j = (L_j - log2 e1_j,  -log2 e2_j,  L_j)     (coefficients of sigma_s, sigma_p, sigma_r).
The logarithms and the weighted sums sum_j lambda_j a_j[i] (exact rational lambda_j) are enclosed in
arb ball arithmetic (python-flint, 256-bit; every operation returns a ball provably containing the
exact real result, i.e. interval arithmetic with outward rounding).  The script checks that the upper
endpoint of the enclosure of sum_j lambda_j a_j[i] is below 0.457 for i = 0, 1, 2, and that
0.457 < log2(1.3727) (lower endpoint of its enclosure).  Then for every point sigma of the simplex,
min_j a_j . sigma <= sum_j lambda_j a_j . sigma < 0.457 < log2(1.3727).
(Version of 2026-10-17; before, the rates were evaluated in 50-digit decimal arithmetic without
outward rounding.)  Requires python-flint (tested with 0.9.0).

Usage: python3 dst_certificate_verify.py [dst_certificates.json]
"""
import json
import sys
import os
from fractions import Fraction as F
from flint import arb, fmpq, ctx

ctx.prec = 256
LN2 = arb(2).log()


def comp(P, Q):
    return (P[0] * (Q[0] + Q[1]), P[1] * (Q[0] / 2 + Q[1]))


def mu_interval(z, w, P):
    """Rationals mu in [0,1] with mu z + (1-mu) w <= P coordinatewise: return one, or None."""
    lo, hi = F(0), F(1)
    for k in range(2):
        # mu (z_k - w_k) <= P_k - w_k
        d, r = z[k] - w[k], P[k] - w[k]
        if d > 0:
            hi = min(hi, r / d)
        elif d < 0:
            lo = max(lo, r / d)
        elif r < 0:
            return None
    return lo if lo <= hi else None


def member(P, Z):
    for z in Z:
        if z[0] <= P[0] and z[1] <= P[1]:
            return True
    Zs = sorted(Z)
    for z, w in zip(Zs, Zs[1:]):
        if mu_interval(z, w, P) is not None:
            return True
    for i in range(len(Z)):
        for j in range(i + 1, len(Z)):
            if mu_interval(Z[i], Z[j], P) is not None:
                return True
    return False


def log2(x):
    """Ball enclosure of log2 of the positive rational x."""
    return arb(fmpq(x.numerator, x.denominator)).log() / LN2


def main(path):
    data = json.load(open(path))
    certs = data["certificates"]
    lam = [F(s) for s in data["lambda"]]
    ok = len(lam) == len(certs) and all(l >= 0 for l in lam) and sum(lam) == 1
    print(f"weights: {len(lam)}, nonnegative and summing to 1: {ok}")
    coeffs = []
    for j, c in enumerate(certs):
        y, e1, e2 = F(c["y"]), F(c["e1"]), F(c["e2"])
        Z = [(F(a), F(b)) for a, b in c["Z"]]
        c0 = y >= 0 and e1 > 0 and e2 > 0 and all(a > 0 and b >= 0 for a, b in Z)
        c1 = all(member(comp(z, w), Z) for z in Z for w in Z)
        c2 = member((1 / e2, F(0)), Z)
        c3 = all(member(tuple(x / e1 for x in comp((y, F(1)), z)), Z) for z in Z)
        cz = min(a for a, b in Z)
        L = log2(1 + y)
        a = (L - log2(e1), -log2(e2), L)
        coeffs.append(a)
        print(f"certificate {j}: y = {float(y):.6f}, e1 = {float(e1):.12f}, e2 = {float(e2):.12f}, "
              f"|Z| = {len(Z)}, c_Z = {float(cz):.6f}; C0 {c0} C1 {c1} C2 {c2} C3 {c3}")
        print(f"   plane coefficients (sigma_s, sigma_p, sigma_r): "
              f"{a[0].str(12)} {a[1].str(12)} {a[2].str(12)}")
        ok = ok and c0 and c1 and c2 and c3
    omega = log2(F(13727, 10000))
    w0 = arb(fmpq(457, 1000))
    sums = [sum((arb(fmpq(l.numerator, l.denominator)) * a[i] for l, a in zip(lam, coeffs)), arb(0))
            for i in range(3)]
    worst = max(sums, key=lambda x: x.upper())
    for i, x in enumerate(sums):
        print(f"sum_j lambda_j a_j[{i}] in {x.str(20, radius=True)}  < 0.457: {x < w0}")
    print(f"log2(1.3727) in {omega.str(20, radius=True)}  > 0.457: {omega > w0}")
    print(f"certified margins: 0.457 - max >= {(w0 - worst.upper()).str(5)}, "
          f"log2(1.3727) - max >= {(omega.lower() - worst.upper()).str(5)}")
    ok = ok and all(x < w0 for x in sums) and omega > w0
    print("CERTIFICATE OK" if ok else "CERTIFICATE FAILED")


if __name__ == "__main__":
    main(sys.argv[1] if len(sys.argv) > 1 else
         os.path.join(os.path.dirname(os.path.abspath(__file__)), "dst_certificates.json"))