"""Rate certificate for Paired Domination (Section "Further applications", lem:fa-pd-rate).

For an integer m >= 2 put beta_m = 2^(-m/2) and
    h_m(y) = log2(1 + y) - log2(y^2 + 4 beta_m y + beta_m^2) / m .
Let y* = (3 sqrt(3) - 1)/16, the unique positive critical point of h_6.

Arithmetic (version of 2026-10-17; before, 50-digit decimal arithmetic without outward rounding was
used): arb ball arithmetic from python-flint at 256 bits.  Every operation (sqrt, exp, log, division,
multiplication) returns a ball that provably contains the exact real result, i.e. interval arithmetic
with outward rounding; y* itself is enclosed by a ball.  A strict inequality between two enclosures
is accepted only if it holds for every pair of points of the two balls.  The script checks
  (1) h_m(y*) < h_6(y*) for m = 2,...,5 and m = 7,...,M0 (M0 = 40), and reports a certified lower
      bound for the smallest gap h_6(y*) - h_m(y*);
  (2) the tail bound log2(1 + y*) + (2/(M0+1)) log2(1/y*) < 0.431 < h_6(y*), which covers m > M0,
      because y^2 + 4 beta_m y + beta_m^2 >= y^2 and 0 < y* < 1 give
      h_m(y*) <= log2(1 + y*) + (2/m) log2(1/y*), a quantity decreasing in m;
  (3) y* is a critical point of h_6: 4 y*^2 + y*/2 - 13/32 = 0 holds exactly, since
      y*^2 = (28 - 6 sqrt 3)/256 (also checked: the enclosure of the residual contains 0);
      and h_6(y*) < 0.70502 and 2^0.70502 < 1.6302 (so 2^(h_6(y*)) < 1.6302).
Usage: python3 paired_rate_certificate.py   (well under a second)
Requires python-flint (tested with 0.9.0).
"""
from flint import arb, fmpq, ctx

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


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


def beta(m):
    return (-(arb(m) / 2) * LN2).exp()


def h(m, y):
    b = beta(m)
    return log2(1 + y) - log2(y * y + 4 * b * y + b * b) / m


ys = (3 * arb(3).sqrt() - 1) / 16
h6 = h(6, ys)
ok = True
print(f"y* in {ys.str(25, radius=True)}")
print(f"h_6(y*) in {h6.str(25, radius=True)}")
print(f"2^h_6(y*) in {(h6 * LN2).exp().str(15, radius=True)}")
crit = 4 * ys * ys + ys / 2 - arb(fmpq(13, 32))
print(f"critical-point residual 4y^2 + y/2 - 13/32 in {crit.str(5, radius=True)} (contains 0: {crit.contains(0)})")
ok &= crit.contains(0)
worst = None
for m in list(range(2, 6)) + list(range(7, M0 + 1)):
    v = h(m, ys)
    gap = h6 - v
    good = gap > 0
    ok &= good
    if worst is None or gap.lower() < worst[0].lower():
        worst = (gap, m)
    if m <= 12:
        print(f"m = {m:2d}: h_m(y*) in {v.str(20, radius=True)}  gap to h_6 in {gap.str(12, radius=True)}")
print(f"smallest gap over 2 <= m <= {M0}, m != 6: >= {worst[0].lower().str(12)} at m = {worst[1]}")
tail = log2(1 + ys) + arb(2) / (M0 + 1) * log2(1 / ys)
t_ok = (tail < arb(fmpq(431, 1000))) and (arb(fmpq(431, 1000)) < h6)
print(f"tail bound for m > {M0}: {tail.str(15, radius=True)} < 0.431 < h_6(y*): {t_ok}")
ok &= t_ok
w = fmpq(70502, 10 ** 5)
w_ok = (h6 < arb(w)) and ((arb(w) * LN2).exp() < arb(fmpq(16302, 10 ** 4)))
print(f"h_6(y*) < 0.70502 and 2^0.70502 in {(arb(w) * LN2).exp().str(15, radius=True)} < 1.6302: {w_ok}")
ok &= w_ok
print(f"certified margin log2(1.6302) - h_6(y*) >= {(log2(fmpq(16302, 10 ** 4)).lower() - h6.upper()).str(5)}")
print("CERTIFICATE OK" if ok else "CERTIFICATE FAILED")