"""Lower bound for RS-theta: on the caterpillar C_L (spine of L vertices without outside
neighbours, two teeth per spine vertex, each tooth with two pendant leaves; n = 7L,
gamma_c = s = 3L, the minimum CDS is unique) |T(X)| = Phi(X) for X subset S*, and RS-theta
succeeds only through X subset S*.  The expected number of successes with budget 2^{Rn} is
at most poly(n) 2^{Rn} max_t mu_t / C(n,t), so RS-theta needs
   R >= min_t [log2 C(n,t) - log2 mu_t] / n,    mu_t = sum_{X in C(S*,t)} 2^{-|T(X)|}.
mu_t is computed exactly (float polynomials) by the tree DP with polynomial messages."""
import math
import numpy as np

def padd(a, b):
    L = max(len(a), len(b)); r = np.zeros(L); r[:len(a)] += a; r[:len(b)] += b; return r

def caterpillar_mu(L):
    # message = (A, B, C) each a polynomial in y (numpy arrays), with a common log2 scale
    def init(m):
        return [np.array([0.0, 0.5]), np.array([0.0]), np.array([2.0 ** (-m)])]
    def merge(P, Q):
        A, B, C = P; a, b, c = Q
        A2 = np.convolve(A, padd(padd(2 * a, b), c))
        B2 = padd(np.convolve(B, padd(padd(a, b), c / 2)), np.convolve(C, a))
        C2 = np.convolve(C, padd(b, c / 2))
        return [A2, B2, C2]
    def norm(P):
        s = max(x.max() for x in P if len(x))
        return [x / s for x in P], math.log2(s)
    tooth = init(2)
    msg = None; lg = 0.0
    for i in range(L):
        P = init(0)
        P = merge(P, tooth); P = merge(P, tooth)
        if msg is not None:
            P = merge(P, msg)
        msg, l = norm(P); lg += l
    A, B, C = msg
    W = padd(padd(A, B), C / 2)
    return W, lg

if __name__ == '__main__':
    for L in (20, 50, 100, 200, 400):
        W, lg = caterpillar_mu(L)
        n, s = 7 * L, 3 * L
        best = None
        for t in range(1, s + 1):
            if W[t] <= 0: continue
            lc = (math.lgamma(n + 1) - math.lgamma(t + 1) - math.lgamma(n - t + 1)) / math.log(2)
            v = (lc - (math.log2(W[t]) + lg)) / n
            if best is None or v < best[0]: best = (v, t)
        print(f"L={L} n={n}: min_t [log2 C(n,t) - log2 mu_t]/n = {best[0]:.6f}  -> base {2 ** best[0]:.5f} (t/n = {best[1]/n:.4f})")