"""JAM 5, line 5 self-audit (founder): the raw conditional-entropy
probe showed an H-drop growing sharply with M (0.057 -> 0.522 ->
1.506 bits) exactly as avg-samples-per-cell shrank (686 -> 172 ->
76). That correlation is the signature of the classic plug-in
entropy estimator's negative bias, not necessarily real structure:
with few samples spread over up to M possible outcomes, the empirical
distribution looks artificially peaked. The standard first-order
(Miller-Madow) correction estimates this bias per cell as
(K-1)/(2*n*ln2) bits, K = alphabet size (M), n = samples in that
cell. Applying it before trusting the line-5 result."""
import math
from collections import defaultdict

N_MAX = 2_000_000
TEST_M = [54, 108, 162]


def build_q(n_max: int):
    Q = [0] * (n_max + 1)
    Q[1] = 1
    if n_max >= 2:
        Q[2] = 1
    for n in range(3, n_max + 1):
        Q[n] = Q[n - Q[n - 1]] + Q[n - Q[n - 2]]
    return Q


def entropy_bits(counts):
    total = sum(counts.values())
    h = 0.0
    for c in counts.values():
        p = c / total
        h -= p * math.log2(p)
    return h


def main():
    print("building Q(1..%d) ..." % N_MAX)
    Q = build_q(N_MAX)
    print("\n  M | raw H-drop | Miller-Madow bias est | corrected H-drop")
    for M in TEST_M:
        marg = defaultdict(int)
        joint = defaultdict(lambda: defaultdict(int))
        for n in range(3, N_MAX + 1):
            a, b, c = Q[n - 2] % M, Q[n - 1] % M, Q[n] % M
            marg[c] += 1
            joint[(a, b)][c] += 1
        H_marg = entropy_bits(marg)
        total_n = N_MAX - 2
        H_cond = 0.0
        bias_total = 0.0
        for (a, b), dist in joint.items():
            n_ab = sum(dist.values())
            H_cond += (n_ab / total_n) * entropy_bits(dist)
            k_observed = len(dist)  # observed distinct outcomes in this cell
            bias_total += (n_ab / total_n) * ((k_observed - 1) / (2.0 * n_ab * math.log(2)))
        raw_drop = H_marg - H_cond
        corrected_drop = raw_drop - bias_total
        print("  %3d | %10.4f | %22.4f | %17.4f"
              % (M, raw_drop, bias_total, corrected_drop))
    print("\nIf corrected H-drop collapses toward ~0 (or goes slightly negative,")
    print("which is a sign the correction itself is noisy at low counts), the raw")
    print("signal was sampling bias, not structure - the honest verdict this jam")
    print("needs before any claim reaches Grok.")


if __name__ == "__main__":
    main()