"""THE EVEN DRIP — the shape of Van Eck's zero-gap distribution.

Van Eck: a(0)=0; a(n+1)=0 if a(n) is a first occurrence, else n - (last prior
index of a(n)).  Zeros stand at 0 = z_0 < z_1 < ...; the gaps are
g_k = z_k - z_{k-1}.

This instrument produces exact integer computation of the gap DISTRIBUTION —
not just its density (THE COUNTING SHADOW) or its inexhaustibility (THE
SELF-FEEDING SPRING), but the shape of the hump the gaps make.  Three things
are measured and one identity is verified live at every zero:

  (1) SELF-RECORDING EQUATION, checked directly: a(z_k + 1) == g_k, i.e. the
      term written immediately after a zero equals the gap that just closed.
      Asserted at every zero; a violation raises.  (PROVED in THE SELF-FEEDING
      SPRING; this is its running witness.)

  (2) THE ACCOUNTING IDENTITY (exact, hence PROVED): the mean gap is forced.
      Sum of the gaps g_2..g_K is exactly z_K - z_1; there are K-1 of them; so
      mean gap = (z_last - z_1)/(#gaps) = N/Z(N) * (1 + o(1)) = 1/zero_share.
      The center of the hump is nailed to THE COUNTING SHADOW.

  (3) THE DRIP STEADIES (COMPUTED): the standard deviation of the gaps stays
      nearly constant (~1.73) while the mean marches right, so the coefficient
      of variation CV = std/mean = std * zero_share falls monotonically,
      ~ std / (2 sqrt(ln N)) -> 0.  The zeros become RELATIVELY evenly spaced
      even as the spring slows.

  (4) THE LEFT EDGE CLOSES (COMPUTED): small gaps freeze.  g=2 occurs a fixed
      handful of times and stops; g=3, g=4 saturate; while the right tail
      (max gap) still creeps outward.  Mechanism (PROVED biconditional):
      g_{k+1} = 2  <=>  a(z_k+1) = g_k is a FIRST occurrence at z_k+1.  A post-
      zero term is a small integer (a gap value); a first occurrence in Van Eck
      almost always carries a LARGE value; the two rarely coincide, so a zero
      two steps after a zero is structurally rare and confined to the early
      sequence.

Evidence class of the printed rows: exact integer computation.  Deterministic;
rerun to reproduce every carved figure.
"""

from array import array
from argparse import ArgumentParser
from math import log, sqrt
from collections import Counter


def census(size, checkpoints):
    if size < 3:
        raise ValueError("size must be >= 3")
    last = array("l", [-1]) * (size + 1)   # last-seen index per value (< size)
    value = 0
    nz = 1                 # zeros seen so far (z_0 counts)
    z1 = None              # position of the first zero after the origin
    last_zero = 0
    sum_gap = 0
    sumsq_gap = 0
    max_gap = 0
    min_gap = 1 << 60
    hist = Counter()
    expect_postzero = None  # g_k, to be matched against a(z_k + 1)
    ckset = set(checkpoints)
    rows = []

    for i in range(1, size):
        prev = last[value]
        last[value] = i - 1
        value_prev = value
        value = 0 if prev == -1 else (i - 1) - prev

        # a(i) has just been computed as `value`.  If a(i-1) was a zero (and not
        # the origin), then a(i) is the post-zero term a(z_k + 1); check it.
        if value_prev == 0 and (i - 1) != 0:
            assert expect_postzero is not None and value == expect_postzero, (
                "SELF-RECORDING EQUATION violated at post-zero position", i)
            expect_postzero = None

        if value == 0:
            g = i - last_zero
            nz += 1
            if z1 is None:
                z1 = i
            else:
                # gaps g_k for k >= 2 (g_1 involves the origin; excluded so the
                # accounting identity is clean).
                sum_gap += g
                sumsq_gap += g * g
                if g > max_gap:
                    max_gap = g
                if g < min_gap:
                    min_gap = g
                hist[g] += 1
            last_zero = i
            expect_postzero = g   # this g_k should reappear as a(z_k + 1)

        terms = i + 1
        if terms in ckset:
            ng = nz - 2
            mean = sum_gap / ng
            std = sqrt(max(sumsq_gap / ng - mean * mean, 0.0))
            zshare = nz / terms
            rows.append((terms, nz, zshare, mean, std, std / mean,
                         min_gap, max_gap, sum_gap, ng, z1,
                         hist.get(2, 0), hist.get(3, 0), hist.get(4, 0)))
    return rows, hist


def main():
    ap = ArgumentParser()
    ap.add_argument("--terms", type=int, default=50_000_000)
    ap.add_argument("--histogram", action="store_true",
                    help="print the full gap histogram at the final N")
    args = ap.parse_args()
    cks = [c for c in (1_000_000, 5_000_000, 10_000_000, 20_000_000,
                       35_000_000, 50_000_000) if c <= args.terms]
    if args.terms not in cks:
        cks.append(args.terms)
    rows, hist = census(args.terms, cks)

    print("SELF-RECORDING EQUATION a(z_k+1)=g_k asserted at every zero: PASSED")
    print(f"{'terms':>12} {'zshare':>9} {'meanG':>7} {'stdG':>7} {'CV':>7}"
          f" {'min':>4} {'max':>4} {'g=2':>5} {'g=3':>5} {'g=4':>7}")
    for (t, nz, zs, mean, std, cv, mn, mx, sg, ng, z1, c2, c3, c4) in rows:
        print(f"{t:>12,} {zs:>9.5f} {mean:>7.3f} {std:>7.3f} {cv:>7.4f}"
              f" {mn:>4} {mx:>4} {c2:>5} {c3:>5} {c4:>7,}")

    print("\nACCOUNTING IDENTITY (exact): mean gap = (z_last - z_1)/#gaps"
          " = 1/zero_share")
    for (t, nz, zs, mean, std, cv, mn, mx, sg, ng, z1, c2, c3, c4) in rows:
        # sg/ng is mean by construction; show it equals 1/zshare to rounding.
        assert abs(sg / ng - mean) < 1e-9
        print(f"  N={t:>12,}  (z_last-z_1)/#gaps={sg/ng:9.5f}"
              f"  1/zero_share={1/zs:9.5f}  diff={sg/ng - 1/zs:+.5f}")

    if args.histogram:
        tot = sum(hist.values())
        print(f"\nGap histogram at N={args.terms:,}:")
        for v in range(min(hist), max(hist) + 1):
            c = hist.get(v, 0)
            bar = "#" * int(60 * c / tot)
            print(f"  {v:>3}: {c:>10,} {c / tot:8.4%}  {bar}")


if __name__ == "__main__":
    main()