"""Rigorous rate evaluation of the two RS-theta antinorm certificates (prop:rst-certs, prop:rst-rate),
in ball arithmetic.

Data: cert_final_0.json, cert_final_1.json (the exact rational certificates; their hypotheses are
checked in exact rational arithmetic by verify_cert.py, verify_cert_indep.py and v_cert_primal.py).
Only y, e1, e2 are read here.  For certificate i the rate is
    omega_i(sigma) = log2(1 + y) - sigma log2 e1 - (1 - sigma) log2 e2 = a0 + a1 sigma,
affine in sigma, so its maximum on an interval is attained at an endpoint.

In arb ball arithmetic (python-flint; each operation returns a ball provably containing the exact
real result, i.e. interval arithmetic with outward rounding), with exact rational inputs, the script
checks that the upper endpoint of the enclosure is below the bound in each of:
  (R1) omega_0(sigma) < omega* = 0.7827016 at sigma = 0.31 and sigma = 0.4172629;
  (R2) omega_1(sigma) < omega*            at sigma = 0.4172629 and sigma = 0.54;
  (R3) all four endpoint values are < 0.78270152 (the value in the coverage table);
  (R4) H(0.2) < 0.7802 (closed-form coverage of sigma <= 0.15 and sigma >= 0.8 via F* <= H);
  (R5) 2^omega* < 1.7204.
It also encloses the crossing point of the two lines and the maximum 0.7827015136... (information).
The bounds F* <= 0.7807 on [0.15, 0.31] and F* <= 0.7802 on [0.54, 0.80] are checked by rs_closed_form_interval.py.
Usage: python3 rst_rate_interval.py   (paths relative to this script; under a second)
Requires python-flint (tested with 0.9.0).
"""
import json
import os
import sys
from flint import arb, fmpq, ctx

ctx.prec = 256
HERE = os.path.dirname(os.path.abspath(__file__))
LN2 = arb(2).log()


def lg(x):
    return arb(x).log() / LN2


def H(p):
    p = arb(p)
    return -(p * p.log() + (1 - p) * (1 - p).log()) / LN2


def main():
    ok = True
    W = fmpq(7827016, 10 ** 7)
    TAB = fmpq(78270152, 10 ** 8)
    S1, S2, S3 = fmpq(31, 100), fmpq(4172629, 10 ** 7), fmpq(54, 100)
    co = []
    worst = None
    for i, (a, b) in enumerate(((S1, S2), (S2, S3))):
        d = json.load(open(os.path.join(HERE, f"cert_final_{i}.json")))
        y, e1, e2 = (fmpq(*d[k]) for k in ("y", "e1", "e2"))
        a0 = lg(1 + y) - lg(e2)
        a1 = lg(e2) - lg(e1)
        co.append((a0, a1))
        print(f"certificate {i}: y = {y}, e1 = {e1}, e2 = {e2}")
        print(f"   omega_{i}(sigma) = a0 + a1 sigma, a0 in {a0.str(15, radius=True)}, "
              f"a1 in {a1.str(15, radius=True)}")
        for s in (a, b):
            r = a0 + a1 * arb(s)
            good = (r < W) and (r < TAB)
            ok &= good
            print(f"   omega_{i}({s}) in {r.str(15, radius=True)}  < omega* and < 0.78270152: {good}")
            if worst is None or r.upper() > worst.upper():
                worst = r
    (a0, a1), (b0, b1) = co
    x = (b0 - a0) / (a1 - b1)
    print(f"crossing point sigma_x in {x.str(15, radius=True)}; value there "
          f"{(a0 + a1 * x).str(15, radius=True)} (information)")
    print(f"largest endpoint value: {worst.str(15, radius=True)}; "
          f"certified margin omega* - max >= {(arb(W) - worst.upper()).str(5)}")
    h = H(fmpq(1, 5))
    r4 = h < arb(fmpq(7802, 10000))
    print(f"(R4) H(0.2) in {h.str(12, radius=True)} < 0.7802: {r4}")
    base = (arb(W) * LN2).exp()
    r5 = base < arb(fmpq(17204, 10000))
    print(f"(R5) 2^omega* in {base.str(12, radius=True)} < 1.7204: {r5}")
    ok &= r4 and r5
    print("CERTIFICATE OK" if ok else "CERTIFICATE FAILED")
    return 0 if ok else 1


if __name__ == "__main__":
    sys.exit(main())