#!/usr/bin/env python3
"""Torn-verify for 014 — The Smell Is the Only Thing That Gets Out. Pure Python.

Claim under test: a black hole is supposed to let nothing out, yet it radiates.
Page's theorem says what that radiation carries. For a Haar-random pure state on
m (x) n, with m <= n, the average entanglement entropy of the m-dimensional part is

    <S> = sum_{j=n+1}^{mn} 1/j  -  (m-1)/(2n)      [Page 1993; proved Foong-Kanno 1994]

Emit k qubits from an N-qubit hole: m = 2^k (radiation), n = 2^(N-k) (hole).
Two independent computations of the same quantity — exact summation and an
asymptotic expansion of the harmonic numbers — must agree.
"""

from fractions import Fraction
from math import log, fsum

LN2 = log(2)
GAMMA = 0.57721566490153286060651209008240243104215933593992

FAILS = []


def check(name, condition, detail=""):
    print(("PASS" if condition else "FAIL") + "  " + name +
          (("  [" + detail + "]") if detail else ""))
    if not condition:
        FAILS.append(name)


# --- method A: exact / direct summation (feasible while mn is small) ----------
def page_exact(m, n):
    """Exact Fraction value of Page's formula. Only for tiny m*n."""
    return sum(Fraction(1, j) for j in range(n + 1, m * n + 1)) - Fraction(m - 1, 2 * n)


def page_direct(m, n):
    """Float value by direct summation of the harmonic difference."""
    return fsum(1.0 / j for j in range(n + 1, m * n + 1)) - (m - 1) / (2.0 * n)


# --- method B: asymptotic expansion of H(x) ----------------------------------
def H(x):
    """Harmonic number via Euler-Maclaurin; x is large."""
    return (log(x) + GAMMA + 1.0 / (2 * x) - 1.0 / (12 * x * x)
            + 1.0 / (120 * x ** 4) - 1.0 / (252 * x ** 6))


def page_asymptotic(m, n):
    return H(m * n) - H(n) - (m - 1) / (2.0 * n)


def S_of_k(N, k, method):
    """Average entropy of the radiation after k of N qubits have been emitted."""
    a, b = k, N - k
    m, n = (2 ** a, 2 ** b) if a <= b else (2 ** b, 2 ** a)   # Page needs m <= n
    return method(m, n)


# --- 1. the textbook anchor: two qubits, one each side ------------------------
two_qubits = page_exact(2, 2)
check("m=n=2 reproduces the known exact average, <S> = 1/3 nats",
      two_qubits == Fraction(1, 3), f"got {two_qubits}")

# --- 2. the two methods must agree -------------------------------------------
N = 20
worst = max(abs(S_of_k(N, k, page_direct) - S_of_k(N, k, page_asymptotic))
            for k in range(1, N))
check("direct summation and asymptotic expansion agree across the whole curve",
      worst < 1e-12, f"N={N}, max |A-B| = {worst:.2e} nats")

# --- 3. before the Page time the radiation is thermal to exponential accuracy -
N = 40
deficits = {}
ok = True
for k in range(1, N // 2):
    S = S_of_k(N, k, page_asymptotic)
    deficit = k * LN2 - S            # how far below maximally mixed, in nats
    deficits[k] = deficit
    if not (0 <= deficit <= 2.0 ** (2 * k - N)):
        ok = False
check("every emission before the Page time is thermal within 2^(2k-N) nats",
      ok, f"N={N}, k=1..{N//2-1}")

first = deficits[1] / LN2
check("the first emitted qubit is thermal to twelve decimal places",
      first < 1e-11, f"information in whiff #1 = {first:.3e} bits")

# --- 4. the curve turns over exactly at half ---------------------------------
curve = [S_of_k(N, k, page_asymptotic) for k in range(N + 1)]
check("the entropy of the radiation peaks at the Page time k = N/2",
      max(range(N + 1), key=lambda k: curve[k]) == N // 2, f"peak at k={N//2}")

sym = max(abs(curve[k] - curve[N - k]) for k in range(N + 1))
check("the curve is its own mirror: S(k) = S(N-k), the global state stays pure",
      sym < 1e-12, f"max asymmetry = {sym:.2e} nats")

peak_deficit = N // 2 * LN2 - curve[N // 2]
check("at the peak the shortfall is exactly half a nat (Page's 1/2)",
      abs(peak_deficit - 0.5) < 1e-9, f"shortfall = {peak_deficit:.12f} nats")

# --- 5. after the Page time information leaves at twice the rate --------------
# FAILED EXPECTATION, KEPT: the first draft asserted the 2-bit slope from k=N/2+2
# onward and it failed at 1.9662 bits. The turnover is not a kink -- it is rounded,
# and the asymptotic rate only sets in a few qubits past the Page time. The width
# of that rounding is now a measured quantity instead of a hidden one.
info = [(k * LN2 - curve[k]) / LN2 for k in range(N + 1)]   # bits recovered
slopes = [info[k + 1] - info[k] for k in range(N)]

settle = next(k for k in range(N // 2, N) if abs(slopes[k] - 2.0) < 1e-9)
check("far past the Page time each further qubit returns 2 bits, not 1",
      all(abs(slopes[k] - 2.0) < 1e-9 for k in range(settle, N)),
      f"exactly 2 bits/qubit from k={settle} on")

# The width of that rounding is not a property of the curve. It is a property of
# how strictly you interrogate it: the shortfall decays like 2^(2k-N), so demanding
# the rate to within eps costs log2(1/eps)/2 further qubits. Falsifiable, so tested.
def width(eps, N=100):
    """Qubits past the Page time before the rate is 2 to within eps. Needs a hole
    big enough that the answer fits inside it -- N=40 is already too small at 1e-12."""
    c = [S_of_k(N, k, page_asymptotic) for k in range(N + 1)]
    s = [((k + 1) * LN2 - c[k + 1] - (k * LN2 - c[k])) / LN2 for k in range(N)]
    return next(k for k in range(N // 2, N) if abs(s[k] - 2.0) < eps) - N // 2


measured = {eps: width(eps) for eps in (1e-3, 1e-6, 1e-9, 1e-12)}
predicted = {eps: log(1 / eps, 2) / 2 for eps in measured}
check("the turnover is rounded, and the width is set by the tolerance: log2(1/eps)/2",
      all(abs(measured[e] - predicted[e]) <= 1.0 for e in measured),
      "  ".join(f"eps={e:.0e}: {measured[e]} vs {predicted[e]:.1f}" for e in measured))

check("so the kink is exponentially sharp: one more qubit buys two more bits of proof",
      all(measured[1e-3] < measured[1e-6] < measured[1e-9] < measured[1e-12] for _ in [0]),
      f"slope at N/2+2 is {slopes[N//2+2]:.6f}, not 2 -- the failure that found this")

check("and the ledger closes: everything thrown in comes back out",
      abs(info[N] - N) < 1e-9, f"recovered {info[N]:.9f} of {N} bits")

# --- 6. third leg: do not believe Page, sample him ---------------------------
# The two computations above are both analytic -- one theorem wearing two faces.
# This one draws actual Haar-random pure states and measures their entanglement.
# Optional: the file stays runnable on a machine without numpy.
try:
    import numpy as np
except ImportError:
    print("SKIP  Haar-random cross-check (no numpy on this machine)")
else:
    rng = np.random.default_rng(20260726)
    Nq, trials = 12, 400
    worst_z = 0.0
    detail = []
    for k in range(1, Nq):
        a, b = min(k, Nq - k), max(k, Nq - k)
        m, n = 2 ** a, 2 ** b
        ent = np.empty(trials)
        for t in range(trials):
            psi = rng.normal(size=(m, n)) + 1j * rng.normal(size=(m, n))
            psi /= np.linalg.norm(psi)               # Haar-random pure state
            p = np.linalg.svd(psi, compute_uv=False) ** 2
            p = p[p > 1e-15]
            ent[t] = -np.sum(p * np.log(p))
        mean, sem = ent.mean(), ent.std(ddof=1) / trials ** 0.5
        pred = page_asymptotic(m, n)
        z = abs(mean - pred) / sem
        worst_z = max(worst_z, z)
        detail.append(f"k={k}: {mean:.4f}+-{sem:.4f} vs {pred:.4f} (z={z:.1f})")
    check("sampled Haar-random states reproduce Page's formula at every k",
          worst_z < 4.0,
          f"N={Nq}, {trials} states per point, worst z-score = {worst_z:.2f}")
    for line in detail:
        print("      " + line)

print()
print("the curve (N=40), bits:")
for k in range(0, N + 1, 4):
    print(f"  emitted {k:2d}   S = {curve[k]/LN2:7.4f}   recovered = {info[k]:7.4f}")

print()
print("FAILS:", FAILS if FAILS else "none")
raise SystemExit(1 if FAILS else 0)