"""THE DRIP'S SHAPE — does Van Eck's zero-gap hump have a limiting FORM?

Companion / sharpening to THE EVEN DRIP and THE ALPHABET THROTTLE.

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 < ...; gaps g_k = z_k - z_{k-1}.

THE EVEN DRIP recorded that the gap std stays "nearly constant ~1.73" while the
mean marches right (mean = 1/zero_share), and left the WISH: is the std truly
BOUNDED (let alone -> sqrt(3) = 1.7320508...)?  Read at 3 decimals the std
looked flat.  This instrument reads it sharper and asks a stronger question
than boundedness: does the whole RECENTERED distribution converge to a fixed
shape?  If so, the drip translates rigidly and std, skewness, and excess
kurtosis all tend to constants.

Method (exact where it matters): accumulate the raw power sums S1=Sum g,
S2=Sum g^2, S3=Sum g^3, S4=Sum g^4 as EXACT Python integers over the gaps
g_k (k>=2, origin gap excluded so the accounting identity stays clean).  At
each geometric checkpoint convert to central moments by the standard
raw->central formulas.  Because the power sums are exact integers, the only
floating error enters at the final division — variance/skew/kurtosis are
reported to full double precision, not the 3-digit read the earlier tool gave.

Reported per checkpoint:
  meanG   = S1/n
  stdG    = sqrt(S2/n - mean^2)          -- watched against sqrt(3)
  skew    = mu3 / var^1.5                 -- >0 means right-skewed
  exkurt  = mu4 / var^2 - 3               -- tail weight vs. normal
  sqrt3-std = 1.7320508... - stdG         -- the approach, signed

If stdG rises monotonically toward sqrt(3) while skew and exkurt settle to
constants, the hump has a genuine limiting FORM (a fixed fluctuation law), and
"bounded std" is a corollary, not the headline.  If stdG marches through
sqrt(3) without slowing, boundedness itself is in doubt.  Either way the read
is exact; rerun to reproduce.

Evidence class of printed rows: exact integer power-sums, one float divide.
"""

from array import array
from argparse import ArgumentParser
from math import sqrt

SQRT3 = 1.7320508075688772


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
    z1 = None
    last_zero = 0
    n = 0                 # count of gaps g_k, k>=2
    S1 = S2 = S3 = S4 = 0  # EXACT integer power sums of the gaps
    ckset = set(checkpoints)
    rows = []

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

        if value == 0:
            g = i - last_zero
            if z1 is None:
                z1 = i
            else:
                n += 1
                g2 = g * g
                S1 += g
                S2 += g2
                S3 += g2 * g
                S4 += g2 * g2
            last_zero = i

        terms = i + 1
        if terms in ckset and n > 1:
            m1 = S1 / n
            m2 = S2 / n
            m3 = S3 / n
            m4 = S4 / n
            var = m2 - m1 * m1
            mu3 = m3 - 3 * m1 * m2 + 2 * m1 ** 3
            mu4 = m4 - 4 * m1 * m3 + 6 * m1 * m1 * m2 - 3 * m1 ** 4
            std = sqrt(max(var, 0.0))
            skew = mu3 / var ** 1.5 if var > 0 else float("nan")
            exk = mu4 / (var * var) - 3 if var > 0 else float("nan")
            rows.append((terms, n, m1, std, skew, exk, SQRT3 - std))
    return rows


def main():
    ap = ArgumentParser()
    ap.add_argument("--terms", type=int, default=100_000_000)
    args = ap.parse_args()
    cks = [c for c in (1_000_000, 2_000_000, 5_000_000, 10_000_000,
                       20_000_000, 50_000_000, 100_000_000, 200_000_000)
           if c <= args.terms]
    if args.terms not in cks:
        cks.append(args.terms)
    rows = census(args.terms, cks)

    print(f"sqrt(3) = {SQRT3:.10f}")
    print(f"{'terms':>13} {'#gaps':>12} {'meanG':>9} {'stdG':>9}"
          f" {'sqrt3-std':>10} {'skew':>9} {'exkurt':>9}")
    for (t, n, mean, std, skew, exk, d3) in rows:
        print(f"{t:>13,} {n:>12,} {mean:>9.4f} {std:>9.6f}"
              f" {d3:>+10.6f} {skew:>9.5f} {exk:>9.5f}")


if __name__ == "__main__":
    main()