"""Rigorous two-sided bound on the exponential rate of Prop. rst-caterpillar(c), in exact
rational arithmetic (Python fractions). Floats are used only to *propose* vectors; every
inequality that the conclusion rests on is checked exactly.

Notation: u_y and M_y as in the proof of Prop. rst-caterpillar(c), so W(y) = fin(M_y^{L-1} u_y),
fin(P) = P1 + P2 + P3/2, n = 7L.  All entries of u_y and M_y are polynomials in y with
nonnegative coefficients, hence nondecreasing in y >= 0.

LOWER BOUND (what the algorithm lower bound needs).  For every L >= 1 and 0 <= t <= n:
      mu_t(S*) / C(n,t) <= K (n+1) B_lo^{-n}.
  Proof: mu_t = 0 for t > 3L.  For t <= 3L put y_t = t/(n-t) in [0, 3/4].  Then
  mu_t y_t^t <= W(y_t) and C(n,t) y_t^t >= (1+y_t)^n/(n+1) (proof of lem:fw-gf), so
  mu_t / C(n,t) <= (n+1) W(y_t)/(1+y_t)^n.  Cover [0,3/4] by cells [a,b].  On a cell, by
  monotonicity, W(y) <= fin(M_b^{L-1} u_b).  With a positive rational vector w and rational
  lam such that w^T M_b <= lam w^T (componentwise, checked exactly),
      fin(M_b^{L-1} u_b) <= c (w^T M_b^{L-1} u_b) <= c lam^{L-1} (w^T u_b),  c = max_i fin_i / w_i,
  and (1+y)^n >= (1+a)^{7L}.  So the ratio is <= (n+1) [c (w^T u_b)/lam] (lam/(1+a)^7)^L, and
  the cell condition  B_lo^7 lam <= (1+a)^7  (checked exactly) gives the claim with
  K = max over cells of c (w^T u_b)/lam.

UPPER BOUND (the value is not larger).  At one rational y*, with v > 0 and lam' such that
  M v >= lam' v (exactly), and c' = min_i u_i / v_i:  W(y*) >= c' fin(v) lam'^{L-1}.  By the
  first inequality of lem:fw-gf (min_t C/mu_t <= (1+y)^n / W(y) for every y),
      min_t C(n,t)/mu_t <= (1+y*)^{7L} / (c' fin(v) lam'^{L-1}),
  and  (1+y*)^7 <= B_hi^7 lam'  (checked exactly) bounds its n-th root by B_hi times a factor
  that tends to 1.
"""
from fractions import Fraction as Fr
import sys

F = (Fr(1), Fr(1), Fr(1, 2))


def uM(y):
    u1 = y / 2 * (y + Fr(1, 4)) ** 2
    u2 = y * y / 4 + y / 8
    u3 = Fr(1, 64)
    M = ((2 * u1, u1, u1), (u2 + u3, u2, u2 / 2), (0, u3, u3 / 2))
    return (u1, u2, u3), M


def perron(Mf, left, it=400):
    v = [1.0, 1.0, 1.0]
    for _ in range(it):
        if left:
            nv = [sum(v[i] * Mf[i][j] for i in range(3)) for j in range(3)]
        else:
            nv = [sum(Mf[i][j] * v[j] for j in range(3)) for i in range(3)]
        s = max(nv); v = [x / s for x in nv]
    return v


def rat(x, den=10 ** 15):
    return Fr(round(x * den), den)


def cell_check(a, b, B7):
    u, M = uM(b)
    Mf = [[float(x) for x in r] for r in M]
    w = [rat(x) for x in perron(Mf, left=True)]
    if min(w) <= 0:
        return None
    wM = [sum(w[i] * M[i][j] for i in range(3)) for j in range(3)]
    lam = max(wM[j] / w[j] for j in range(3))
    if B7 * lam > (1 + a) ** 7:
        return None
    c = max(F[i] / w[i] for i in range(3))
    K = c * sum(w[i] * u[i] for i in range(3)) / lam
    return K


def lower(B, ymax=Fr(3, 4), init=64, maxdepth=40):
    B7 = B ** 7
    cells = [(ymax * i / init, ymax * (i + 1) / init, 0) for i in range(init)]
    done, Kmax, widths = 0, Fr(0), []
    while cells:
        a, b, d = cells.pop()
        K = cell_check(a, b, B7)
        if K is None:
            if d >= maxdepth:
                return None
            m = (a + b) / 2
            cells += [(a, m, d + 1), (m, b, d + 1)]
            continue
        done += 1; Kmax = max(Kmax, K); widths.append(b - a)
    return done, Kmax, min(widths)


def upper(B, ystar):
    u, M = uM(ystar)
    Mf = [[float(x) for x in r] for r in M]
    v = [rat(x) for x in perron(Mf, left=False)]
    Mv = [sum(M[i][j] * v[j] for j in range(3)) for i in range(3)]
    lam = min(Mv[i] / v[i] for i in range(3))
    ok = (1 + ystar) ** 7 <= B ** 7 * lam
    cprime = min(u[i] / v[i] for i in range(3))
    return ok, lam, cprime


if __name__ == '__main__':
    B_lo = Fr(sys.argv[1]) if len(sys.argv) > 1 else Fr('1.71927')
    B_hi = Fr(sys.argv[2]) if len(sys.argv) > 2 else Fr('1.71928')
    ystar = Fr(sys.argv[3]) if len(sys.argv) > 3 else Fr('0.3453')
    r = lower(B_lo)
    ok = r is not None and r[1] <= 1
    if r is None:
        print(f"LOWER: could not certify B_lo = {B_lo}")
    else:
        done, K, wmin = r
        print(f"LOWER: for all L>=1 and t, mu_t/C(n,t) <= K (n+1) B^-n with B = {B_lo} (exact rational checks on"
              f" {done} cells of [0,3/4], narrowest {float(wmin):.2e}); K = {float(K):.6g}, K <= 1 exactly: {K <= 1}")
    up, lam, cp = upper(B_hi, ystar)
    ok = ok and up
    print(f"UPPER at y* = {ystar}: (1+y*)^7 <= B^7 lam' with B = {B_hi}: {up}"
          f"  (lam' = {float(lam):.12g}, c' = {float(cp):.6g})")
    print("CERTIFICATE OK" if ok else "CERTIFICATE FAILED")
    sys.exit(0 if ok else 1)