"""gapbar.py -- the dam's bar F under a gap hypothesis for progressions:
does the largest normalized gap past T, over the moduli the empty
seating reads, stay bounded as T grows?

QUESTION. The feeder formula makes P*(T) a maximum of least primes
p_T(s) (the least prime = 1 mod s at or past T) over prime powers
s < T. A gap bound for the progressions 1 mod s just past T then
bounds F(t) = P*(T+(t)), and with it the dam ratio rho. Is the gap
bound the conditional theorem asks for borne out on F's own moduli?

DESIGN.
(1) THE EMPTY SEATING (proved). Seat no odd prime below T. Then N0 =
lcm(p - 1 : p prime in [T, P]), and the striking rule's bars are, at
every odd prime l < T unseated, (l - 1) | N0 and v_l(N0) >= k_l - 2,
and at 2, v_2(N0) >= max(1, k_2 - 3). (l - 1) | N0 exactly when each
prime power s exactly dividing l - 1 divides some p - 1 with p in
[T, P], that is p_T(s) <= P. So the empty seating is feasible at P
exactly when P >= M_D(T), the largest p_T(s) over D_T, the s exactly
dividing some l - 1 (l an odd prime below T), the l^(k_l - 2) with
k_l >= 3, and 2^max(1, k_2 - 3); and P*(T) <= M_D(T) <= M(T), the
largest p_T(s) over every prime power s < T. On D_T, phi(s) < T/2: an
odd s exactly dividing l - 1 has l - 1 = s m, m even, so s <= (T - 2)/2;
phi(2^j) = 2^(j - 1) < T/2; l^(k_l - 2) <= (T - 1)/l. An odd s exactly
dividing some l - 1 and in (T/4, T/2) has m = 2, so 2s + 1 < T is
prime (the pair bound's r when s is prime); powers of 2 and
3^(k_3 - 2) can lie there too.
(2) THE HYPOTHESIS. H(A): at every large prime power T, every s in D_T
has a prime = 1 mod s in [T, T + A phi(s) ln^2 T]. Then M_D(T) <=
T + (A/2) T ln^2 T; T+(t) <= 2t, so F(t) <= 2t + A t ln^2(2t); and
rho(P) = P+/tau(P) <= F(tau(P))/tau(P), tau(P) <= P, so rho(P) <=
2 + A ln^2(2P) at every large P. GRH alone, through psi(x; q, 1) =
x/phi(q) + O(sqrt(x) ln^2 x), gives only p_T(s) << T + s^2 ln^4 s, so
F(t) << t^2 ln^4 t and rho << sqrt(P) ln^2 P, no better than the
floor's route under GRH (the GRH form not read in a source).
(3) THE STATISTIC. G_D(T) = max over s in D_T of (p_T(s) - T)/(phi(s)
ln^2 T), with its argmax s; G(T), the same over every prime power
s < T, a print; M_D(T)/P*(T), how loose (1) is. Read at the 1278 T =
T+(t), 4 <= t <= 10^4, and at pstar.py's 36 draws (d = 4, 5, 6, the
twelve t = round(10^(d + i/24)), i = -6..5, at T = T+(t)). p_T(s) and
P*(T) are pstar.Engine's (its L and pstar()).
(4) THE NULL (a TRANSPLANT: each wait p_T(s) - T exponential at density
1/(phi(s) ln n), n = T + the wait, independent, over D_T's own s at
T = 10^d, prime counts only; null() below prints it). G_D's median
0.803 at T = 10^4 (10..90% 0.632..1.081), 0.834 at 10^5, 0.858 at
10^6 (0.746..1.036); G's median 1.140, 1.127, 1.118; the
share of G_D's tail from s in (T/4, T/2) 0.26..0.29. So R_D, G_D's
d = 6 window median over its d = 4 one, is 1.07 under a bounded law,
about 1.6 under a growth like ln T (gaps phi(s) ln^3 T). One draw's
log spread is about 0.21 at 10^4 and 0.13 at 10^6, so twelve
independent draws carry log R_D to about 0.09; the draws are not
independent (pstar.py (5)), so more.

PREDICTIONS (frozen before the code).
Q1 (proved) P*(T) <= M_D(T) <= M(T) at every T read.
Q2 (the model, 70%) R_D < 1.3.
Q3 (the model, 60%) each window's median G_D within 0.8..1.5 times the
   model's 0.803, 0.834, 0.858.
Q4 (the model, 60%) G_D's argmax s in (T/4, T/2) at 10%..45% of the
   1278 T.
Q5 (a guess, 60%) M_D/P* median in 1.2..4 in every window.

THE KILL, frozen. The script prints R_D and then "R_D < 1.3: G_D
BOUNDED AT THESE T" or "KILL: R_D >= 1.3, G_D GROWS". The bar sits
between the bounded law's 1.07 and ln T growth's 1.6, past the slow
lift the normalizer's ln n/ln T gives a bounded law, the arm a confound
inflates.

CONTROLS, read first. K1 T = 8 by hand: s = 2, 3, 4, 5, 7 give p_8 =
11, 13, 13, 11, 29, so G(8) = 21/(6 ln^2 8) = 0.809 at s = 7; D_8 =
{2, 3, 4}, G_D(8) = 3/ln^2 8 = 0.694 at s = 2, M_D(8) = 13 = P*(8).
K2 p_T(s) equals a brute least prime = 1 mod s at or past T at every
s of T = 1009 and 1024. K3 the kill's reading on planted windows: G_D
= 0.8 at every draw prints bounded, G_D = ln T/11 prints the KILL. K4
Q1's test flags a planted break, M_D taken over s < T/4 only, at some
T.

FINDINGS (from the printed run; the controls passed first). K1 prints
G 0.8094 at 7, G_D 0.6938 at 2, M_D 13 = P*; K2 holds at 1009 and 1024;
K3 prints R_D 1.000 and 1.503; K4 flags its broken bound at 1070 of
1278 T. THE NULL prints (4)'s medians, quantiles and shares as quoted.
Q1 holds at every T read. THE PRINT: R_D = 1.1161 (log 0.110):
R_D < 1.3, G_D BOUNDED AT THESE T, the kill silent; Q2 holds; G's own
ratio 0.8981. Window medians of G_D 0.6468, 0.6480, 0.7219 (model
0.803, 0.834, 0.858: Q3 holds at d = 4 and 6 and fails at d = 5,
ratio 0.777), the largest over the draws 0.8148 at T = 100003
(s = 49811, s/T 0.498); over the 1278 T, G_D's median 0.538, its
largest at T >= 100 0.7488 at T = 101 (s = 41). The waits run below
the model's level, as correlated waits would (one prime p past T
serves every s dividing p - 1), a level, not a trend. G_D's argmax
lies in (T/4, T/2) at 658/1278 T (0.515: Q4 FAILS), 652 of them a
prime r with 2r + 1 prime, the pair bound's r, and 6 a power of 2, and
at 9, 6 and 7 of 12 draws, each such an r: the pair bound's r carry
G_D's peak at about half the T, past the model's 0.26..0.29 share. One s holds
the argmax across a run of draws (3389 at 10 of d = 4's 12, 49811 at 6
of d = 5's, 443273 at 7 of d = 6's), its p_T(s) fixed while T grows,
so a window's draws are far from independent. M_D/P* medians 1.579,
1.677, 1.540 (Q5 holds), range 1.143..2.711 over the draws, M_D = P*
at 161 of 1278 T: the empty seating's threshold sits within a factor
2.72 of P* at every draw. G over every s < T: medians 0.8837, 0.7671,
0.7936, its argmax past T/2 at 7, 8 and 9 of 12 draws, one s held
across a run of draws (5227 at every d = 4 draw).
RUN RECORD. Rehearsal --top 4: 5/5, 6.6 s, peak 13.5 MB. Full run:
6/6, 26.7 s, peak 157.7 MB. With the null and the
argmax classes printed: 6/6, 27.5 s, peak 156.5 MB, every figure above
unchanged; rerun 26.2 s, peak 158.7 MB under a memory guard.

    python gapbar.py [--top D]     (D the largest window, default 6; a
                                    rehearsal runs --top 4)
"""

import os
import sys
import time
from math import exp, log

os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import pstar  # noqa: E402

CHECKS = []
MODEL = {4: 0.803, 5: 0.834, 6: 0.858}   # the null's G_D medians, (4)
BAR = 1.3


def check(name, ok, detail=""):
    CHECKS.append((name, bool(ok)))
    print(f"  [{'PASS' if ok else 'FAIL'}] {name}"
          + (f" -- {detail}" if detail else ""))


def dset(E):
    """D_T as (r, j) keys of E.L, (1)."""
    k = E.k
    D = {(2, max(1, k[2] - 3))}
    for l in E.odd:
        D.update(E.fac[l])
        if k[l] >= 3:
            D.add((l, k[l] - 2))
    return D


def phi(r, j):
    return r ** j - r ** (j - 1)


def gmax(E, keys):
    """(the largest (p_T(s) - T)/(phi(s) ln^2 T) over keys, its s)."""
    L2 = log(E.T) ** 2
    best = (-1.0, None)
    for r, j in keys:
        g = (E.L[(r, j)] - E.T) / (phi(r, j) * L2)
        if g > best[0]:
            best = (g, r ** j)
    return best


def read(T, spf, primes):
    t0 = time.time()
    E = pstar.Engine(T, spf, primes)
    D = dset(E)
    gd, sd = gmax(E, D)
    kd = kind(sd, spf)
    ga, sa = gmax(E, E.L)
    md = max(E.L[s] for s in D)
    m = max(E.L.values())
    quarter = [s for s in D if s[0] ** s[1] < T / 4]
    mq = max(E.L[s] for s in quarter) if quarter else 0
    P = E.pstar()[0]
    return {"T": T, "P": P, "GD": gd, "sD": sd, "kD": kd, "G": ga,
            "sG": sa,
            "MD": md, "M": m, "MQ": mq, "nD": len(D),
            "s": time.time() - t0}


def kind(s, spf):
    """An argmax's class: 'sg' a prime s with 2s + 1 prime, 'two' a
    power of 2, else 'other'."""
    if s & (s - 1) == 0:
        return "two"
    if spf[s] == 0 and spf[2 * s + 1] == 0:
        return "sg"
    return "other"


def null(T, spf, primes):
    """(4): the model's G_D quantiles (10%, 50%, 90%), the share of its
    tail at the median from s in (T/4, T/2), and G's median, at T."""
    low = [p for p in primes if p < T]
    k = {}
    for r in low:
        j, x = 1, r
        while x < T:
            j, x = j + 1, x * r
        k[r] = j
    D = {2 ** max(1, k[2] - 3)}
    for l in low[1:]:
        D.update(r ** e for r, e in pstar.factor(l - 1, spf).items())
        if k[l] >= 3:
            D.add(l ** (k[l] - 2))
    A = [r ** j for r in low for j in range(1, k[r])]
    L2 = log(T) ** 2

    def phis(s):
        r = spf[s] or s
        return s - s // r

    def tail(g, ss):
        out = 0.0
        for s in ss:
            y = g * phis(s) * L2
            lam = y / (phis(s) * log(T + y))
            out += -log(1 - exp(-lam)) if lam < 30 else exp(-lam)
        return out

    def solve(target, ss):
        a, b = 0.02, 5.0
        for _ in range(40):
            m = (a + b) / 2
            if tail(m, ss) > target:
                a = m
            else:
                b = m
        return m

    D = sorted(D)
    med = solve(log(2), D)
    band = [s for s in D if T / 4 < s < T / 2]
    return (solve(log(10), D), med, solve(log(1 / 0.9), D),
            tail(med, band) / tail(med, D), solve(log(2), A))


def median(xs):
    xs = sorted(xs)
    n = len(xs)
    return xs[n // 2] if n % 2 else (xs[n // 2 - 1] + xs[n // 2]) / 2


def verdict(win):
    top = max(win)
    R = median(win[top]) / median(win[4])
    line = (f"R_D < {BAR}: G_D BOUNDED AT THESE T" if R < BAR else
            f"KILL: R_D >= {BAR}, G_D GROWS")
    return R, line


def main():
    top = 6
    if "--top" in sys.argv:
        top = int(sys.argv[sys.argv.index("--top") + 1])
    t0 = time.time()
    ts = {d: [round(10 ** (d + i / 24)) for i in range(-6, 6)]
          for d in range(4, top + 1)}
    nmax = max(max(v) for v in ts.values()) + 1000
    spf = pstar.spf_to(nmax)
    primes = [p for p in range(2, nmax + 1) if spf[p] == 0]
    print(f"sieve to {nmax}: {len(primes)} primes ({time.time() - t0:.1f} s)")

    print("CONTROLS")
    x = read(8, spf, primes)
    check("K1 T = 8 by hand: G 0.809 at 7, G_D 0.694 at 2, M_D 13 = P*",
          abs(x["G"] - 21 / (6 * log(8) ** 2)) < 1e-12 and x["sG"] == 7
          and abs(x["GD"] - 3 / log(8) ** 2) < 1e-12 and x["sD"] == 2
          and x["MD"] == 13 == x["P"] and round(x["G"], 3) == 0.809
          and round(x["GD"], 3) == 0.694,
          f"G {x['G']:.4f} at {x['sG']}, G_D {x['GD']:.4f} at {x['sD']}, "
          f"M_D {x['MD']}, P* {x['P']}")
    bad = []
    for T in (1009, 1024):
        E = pstar.Engine(T, spf, primes)
        bad += [(T, r, j) for (r, j), v in E.L.items()
                if v != pstar.brute_least(r ** j, T)]
    check("K2 p_T(s) equals the brute least prime = 1 mod s at or past T "
          "at every s, T = 1009 and 1024", not bad, f"{bad[:3]}")
    if top > 4:
        a = verdict({d: [0.8] * len(ts[d]) for d in ts})
        b = verdict({d: [log(pstar.next_pp(t, spf)) / 11 for t in ts[d]]
                     for d in ts})
        check("K3 the kill's reading: G_D = 0.8 prints bounded, G_D = "
              "ln T/11 prints the KILL",
              not a[1].startswith("KILL") and b[1].startswith("KILL"),
              f"R_D {a[0]:.3f}, {b[0]:.3f}")

    print("THE NULL (4)")
    for d in sorted(ts):
        q10, q50, q90, sh, ga = null(10 ** d, spf, primes)
        print(f"  T = 10^{d}: G_D median {q50:.3f} (10..90% {q10:.3f}.."
              f"{q90:.3f}), share from (T/4, T/2) {sh:.2f}; G median "
              f"{ga:.3f}", flush=True)

    t1 = time.time()
    small = {}
    for t in range(4, 10 ** 4 + 1):
        T = pstar.next_pp(t, spf)
        if T not in small:
            small[T] = read(T, spf, primes)
    rows = list(small.values())
    plant = [r["T"] for r in rows if r["P"] > r["MQ"]]
    check("K4 Q1's test flags the planted break (M_D over s < T/4 only) "
          "at some T", plant != [], f"{len(plant)} of {len(rows)} T")
    if not all(ok for _, ok in CHECKS):
        print("CONTROL FAILED: nothing below is read.")
        raise SystemExit(1)

    print(f"THE 1278 T ({len(rows)} T to {max(small)}, "
          f"{time.time() - t1:.1f} s)")
    bad = [r["T"] for r in rows if not r["P"] <= r["MD"] <= r["M"]]
    check("Q1 P*(T) <= M_D(T) <= M(T) at every small T", not bad,
          f"{bad[:3]}")
    big = [r for r in rows if r["T"] >= 100]
    top_d = max(big, key=lambda r: r["GD"])
    top_a = max(big, key=lambda r: r["G"])
    band = sum(r["T"] / 4 < r["sD"] < r["T"] / 2 for r in rows)
    half = sum(r["sD"] > r["T"] / 2 for r in rows)
    inb = [r["kD"] for r in rows if r["T"] / 4 < r["sD"] < r["T"] / 2]
    print(f"  G_D's argmax in (T/4, T/2) by class: "
          f"{ {c: inb.count(c) for c in ('sg', 'two', 'other')} }")
    print(f"  G_D: median {median([r['GD'] for r in rows]):.3f}, at T >= 100 "
          f"largest {top_d['GD']:.4f} at T = {top_d['T']} (s = "
          f"{top_d['sD']}, s/T {top_d['sD'] / top_d['T']:.3f})")
    print(f"  G: median {median([r['G'] for r in rows]):.3f}, at T >= 100 "
          f"largest {top_a['G']:.4f} at T = {top_a['T']} (s = "
          f"{top_a['sG']}, s/T {top_a['sG'] / top_a['T']:.3f})")
    print(f"  G_D's argmax s in (T/4, T/2) at {band}/{len(rows)} T, past "
          f"T/2 (powers of 2) at {half}")
    print(f"    Q4: {'HOLDS' if 0.10 <= band / len(rows) <= 0.45 else 'FAILS'}"
          f" ({band / len(rows):.3f})")
    print(f"  M_D/P*: median {median([r['MD'] / r['P'] for r in rows]):.3f}, "
          f"M_D = P* at {sum(r['MD'] == r['P'] for r in rows)}/{len(rows)}")

    print("THE DRAWS")
    draws = []
    for d in sorted(ts):
        for t in ts[d]:
            x = read(pstar.next_pp(t, spf), spf, primes)
            x["t"], x["d"] = t, d
            draws.append(x)
            print(f"  d {d} t {t} T {x['T']}: G_D {x['GD']:.4f} at s "
                  f"{x['sD']} (s/T {x['sD'] / x['T']:.3f}), G {x['G']:.4f} "
                  f"at s {x['sG']} (s/T {x['sG'] / x['T']:.3f}), |D| "
                  f"{x['nD']}, P* {x['P']}, M_D/P* {x['MD'] / x['P']:.3f}, "
                  f"M/P* {x['M'] / x['P']:.3f} ({x['s']:.1f} s)", flush=True)
    bad = [x["t"] for x in draws if not x["P"] <= x["MD"] <= x["M"]]
    check("Q1 P*(T) <= M_D(T) <= M(T) at every draw", not bad, f"{bad[:3]}")

    print("THE WINDOWS")
    win = {}
    for d in sorted(ts):
        ws = [x for x in draws if x["d"] == d]
        gd = [x["GD"] for x in ws]
        win[d] = gd
        md = median([x["MD"] / x["P"] for x in ws])
        inb = [x for x in ws if x["T"] / 4 < x["sD"] < x["T"] / 2]
        print(f"  d {d}: G_D median {median(gd):.4f} (range {min(gd):.4f}"
              f"..{max(gd):.4f}; model {MODEL[d]}), G median "
              f"{median([x['G'] for x in ws]):.4f}, M_D/P* median {md:.3f}, "
              f"argmax s in (T/4, T/2) at "
              f"{len(inb)}/{len(ws)}, its classes "
              f"{sorted(x['kD'] for x in inb)}; G's argmax past T/2 at "
              f"{sum(x['sG'] > x['T'] / 2 for x in ws)}/{len(ws)}")
        r = median(gd) / MODEL[d]
        print(f"    Q3 at d {d}: median/model {r:.3f}, "
              f"{'HOLDS' if 0.8 <= r <= 1.5 else 'FAILS'}")
        print(f"    Q5 at d {d}: {'HOLDS' if 1.2 <= md <= 4 else 'FAILS'}")

    print("THE PRINT")
    if top > 4:
        R, line = verdict(win)
        Ra = (median([x["G"] for x in draws if x["d"] == top]) /
              median([x["G"] for x in draws if x["d"] == 4]))
        print(f"  R_D = {R:.4f} (log {log(R):.3f}; the null's 1.07, ln T "
              f"growth's 1.6); G's own ratio, a print, {Ra:.4f}")
        print(f"  Q2: {'HOLDS' if R < BAR else 'FAILS'}")
        print("  " + line)
    else:
        print("  rehearsal: one window, no R_D")
    print(f"{sum(ok for _, ok in CHECKS)}/{len(CHECKS)} checks passed, "
          f"{time.time() - t0:.1f} s")
    raise SystemExit(0 if all(ok for _, ok in CHECKS) else 1)


if __name__ == "__main__":
    main()