choosy.py

"""choosy.py -- does a run that seats only what the tide demands pay less
than the record run?

QUESTION. From the void every escaping run pays the depth dam's B(t)
pointwise, and the record run is f-greedy (f >= t) exactly when
f >= max(S, B), S(t) the largest
p(r^e) over its seated primes r with r^e <= t (envelope.py); to 10^4
that is S, with E(10^4) = 359747. Does a run that seats only the least
anchor of the tide's unborn power pay less: does its envelope fall below
E/10 at t = 10^4?

DESIGN. One odd prime per move, exponent 1, 2 never seated. The tide
m0 is a price at one prime r, by the door formulas (as envelope.m0_of),
and the move is:
  r seated: s = m0 = r^(u_r+1), seat p(s);
  r unheld with r - 1 | lambda: s = m0 / r, seat p(s);
  r = 2: s = m0 / 4 = 2^(u_2+1), seat p(s);
  r unheld with r - 1 not dividing lambda (m0 = r), two arms:
    ARM A seats r;
    ARM B seats the least prime P > r with P = 1 mod m*, m* the part of
    r - 1 not dividing lambda, so r is skipped from then on.
Each move brings a new prime (a held p has p - 1 | lambda; s and m* do
not divide lambda), so each arm escapes; E_run(t) = max{seat : m0 <= t}.
Prices only rise and a new seat P has P >= m0 (P is unheld with
P - 1 not dividing lambda, so it priced P), so the tide is
nondecreasing and the run is read until m0 > 10^4.
S_run(t) = max{p(r^e) : r held, e > u_r at r's entry, r^e <= t}.
(1) E_run >= max(B, S_run): the run escapes, so the depth dam gives B;
a held r at exponent 1 crosses each e past its entry u_r from a tide at
most r^e, by a seat at least p(r^e).
(2) Arm A only: the seat is p(m0) <= S_run(m0) (u_r + 1 is past r's
entry u_r), p(m0/r) or p(m0/4), B's terms at the point m0, or r = m0;
so E_run <= max(S_run, B, t), equal to max(S_run, B) where that is at
least t. Arm B's seat P at a tide r has no such bound.

PREDICTIONS (frozen before the code).
P1 (forced) the engine under the record policy gives envelope.py's E.
P2 (forced) E_run >= max(B, S_run) at every t <= 10^4, both arms.
P3 (forced) arm A: E_run <= max(S_run, B, t), equal to max(S_run, B)
   wherever that is >= t.
P4 (a guess, 75%) neither arm's E_run(10^4) is below E(10^4)/10.
P5 (a guess, 60%) arm A's E_run(10^4) is below E(10^4).
P6 (a guess, 55%) arm B's E_run(10^4) exceeds arm A's.

THE PRINT. "BELOW E/10: A", "B", "A and B" or "neither", from
E_run(10^4) against 35975. "RUN KILL" if P2 or P3 fails, an engine
fault rather than an answer.

CONTROLS, read first. C1 is P1, at every t <= 10^4. C2 the arms' first
seats against hand lists: arm A 3, 7, 5, 11, 29, 19, 23 (at lambda =
420 the tide is 9, 3 held with u_3 = 1, before 11); arm B 5, 11, 13,
23, 53, 97 (at the void m* = 2 and 5 is the least prime past 3 that is
1 mod 2; at the tide 7, lambda = 20 and m* = 3; at 17, lambda = 8580
and m* = 16). The lists first frozen (A ending 23, 19; B opening 3, 7)
broke the arms' own rules by hand, and were re-derived before any
verdict was printed. C3 P2's test flags E_run lowered by 1 at the
largest point of S_run.

(3), derived after the first run from its readings: every p(s), s a
prime power, is a record seat, since no prime below p(s) is 1 mod s, so
s divides no lcm of p - 1 over primes below p(s). So arm A's anchor
seats are record seats; its open seats r have no such argument, and the
run's "every held prime a record seat" and "held <= 10^4 the record's
seats <= 10^4" are readings, added with the B-above-E line after the
first run.

FINDINGS (from the printed run; the controls passed first). No run
kill: P2 holds on both arms and P3 on arm A at every t <= 10^4. THE
PRINT: BELOW E/10 (35974.7): neither; P4 holds. Arm A seats only record
seats to 10^4: 519 seats to tide 10^4, every one a record seat,
its 302 held primes <= 10^4 the record run's 302 seats <= 10^4, and
E_run = E = S_run at every t <= 10^4 (P5 fails: equal, not below).
Arm B, which skips each open prime by seating a prime = 1 mod its
unborn part, holds 175 primes <= 10^4 (425 seats, not all record
seats) and is cheaper at 10^4, E_run = 324469 = 0.9019 E (P6 fails),
largest E_run/t 44.0005 = 87869/1997 against E's 60.0003; but it is
above E at 484 t, from 13 to 9466, below at 6816, least E_run/E
0.2538 over 3..10^4, its envelope S_run at 9183 t and neither S_run
nor B at 816 (arm B's open seats, which (2) does not bound).
RUN RECORD. First run: C2 failed on the frozen hand lists, which broke
the arms' own rules; re-derived, 6/6, 5.1 s, 55.4 MB. With the
readings of (3) and the B-above-E line: 6/6, 5.1 s, 55.3 MB. An audit:
the bar was E // 10 under a strict <, missing E_run = 35974, now
10 E_run < E: 6/6, 5.4 s, 54.5 MB, the verdict unchanged.
Moved beside envelope.py, its import path dropped: 6/6, 4.4 s,
55.0 MB, the prints unchanged.
"""

import os
import time

os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import growth  # noqa: E402
import envelope as ev  # noqa: E402

T = ev.T
CHECKS = []


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


def factor(n):
    if n <= ev.N:
        return ev.fac(n)
    out, m = {}, n
    for q in ev.PRIMES:
        if q * q > m:
            break
        while m % q == 0:
            out[q] = out.get(q, 0) + 1
            m //= q
    if m > 1:
        out[m] = out.get(m, 0) + 1
    return out


def tide(held, vN):
    """(m0, r, kind): the least price, its prime, and 'seated',
    'skipped', 'two' or 'open' (unheld, r - 1 not dividing lambda)."""
    w = vN.get(2, 0)
    best = (4 if w == 0 else 2 ** (w + 3), 2, "two")
    for r in ev.PRIMES[1:]:
        if r >= best[0]:
            break
        u = vN.get(r, 0)
        if r in held:
            cand = (r ** (u + 1), r, "seated")
        elif not ev.divides(ev.fac(r - 1), vN):
            cand = (r, r, "open")
        else:
            cand = (r ** (u + 2), r, "skipped")
        if cand[0] < best[0]:
            best = cand
    return best


def least_past(r, m):
    """The least prime P > r with P = 1 mod m, m | r - 1."""
    n = r + m
    while not (ev.SPF[n] == n if n <= ev.N else growth.is_prime(n)):
        n += m
    return n


def run(policy):
    """[(m0, seat)], the held primes with their u_r at entry, the
    tides, until m0 > T."""
    held, vN, entry = set(), {}, {}
    pairs, nxt = [], 1                       # nxt: the record's read index
    while True:
        m0, r, kind = tide(held, vN)
        if m0 > T:
            return pairs, entry
        if policy == "record":
            while ev.divides(ev.fac(ev.PRIMES[nxt] - 1), vN):
                nxt += 1
            P = ev.PRIMES[nxt]
            nxt += 1
        elif kind == "seated":
            P = ev.p_of(m0)
        elif kind == "skipped":
            P = ev.p_of(m0 // r)
        elif kind == "two":
            P = ev.p_of(m0 // 4)
        elif policy == "A":
            P = r
        else:
            fr = ev.fac(r - 1)
            ms = 1
            for q, e in fr.items():
                if e > vN.get(q, 0):
                    ms *= q ** e
            P = least_past(r, ms)
        assert P not in held and P > 2
        pairs.append((m0, P))
        entry[P] = vN.get(P, 0)
        held.add(P)
        for q, e in factor(P - 1).items():
            if e > vN.get(q, 0):
                vN[q] = e


def bars():
    """B and the prime-power table, as envelope.py builds them."""
    pps = ev.prime_powers(T)
    pv = {s: ev.p_of(s) for s, _, _ in pps}
    bpts = [((4 * s if r == 2 else r * s), pv[s]) for s, r, _ in pps]
    return ev.hull([b for b in bpts if b[0] <= T]), pv


def s_run(entry, pv):
    pts = []
    for r, u in entry.items():
        q, e = r, 1
        while q <= T:
            if e > u:
                pts.append((q, pv[q]))
            q, e = q * r, e + 1
    return ev.hull(pts)


def main():
    t0 = time.time()
    B, pv = bars()
    epairs, eseats, _, _ = ev.record_run()
    E = ev.hull([(m, R) for m, R, _ in epairs])
    print("CONTROLS")
    rp, _ = run("record")
    Er = ev.hull(rp)
    check("C1 the record policy gives envelope.py's E at every t <= 10^4",
          Er[1:] == E[1:], f"E(10^4) = {Er[T]}, {len(rp)} seats")
    runs = {arm: run(arm) for arm in ("A", "B")}
    hand = {"A": [3, 7, 5, 11, 29, 19, 23], "B": [5, 11, 13, 23, 53, 97]}
    got = {arm: [P for _, P in runs[arm][0][:len(hand[arm])]]
           for arm in hand}
    check("C2 the arms' first seats equal the hand lists",
          got == hand, f"A {got['A']}, B {got['B']}")
    env, srun = {}, {}
    for arm, (pairs, entry) in runs.items():
        env[arm] = ev.hull(pairs)
        srun[arm] = s_run(entry, pv)

    def p2(Ex, Sx):
        return [t for t in range(1, T + 1) if Ex[t] < max(B[t], Sx[t])]
    sa = srun["A"]
    tm = min(t for t in range(1, T + 1) if sa[t] == sa[T])
    plant = list(env["A"])
    plant[tm] = sa[tm] - 1
    check("C3 P2's test flags E_run lowered below S_run at its largest"
          " point", tm in p2(plant, sa), f"t = {tm}, S_run {sa[tm]}")
    if not all(ok for _, ok in CHECKS):
        print("CONTROL FAILED: nothing below is read.")
        raise SystemExit(1)

    print("THE ARMS")
    bad = False
    for arm in ("A", "B"):
        f2 = p2(env[arm], srun[arm])
        check(f"P2 arm {arm}: E_run >= max(B, S_run) at every t <= 10^4",
              not f2, f"{len(f2)} t off {f2[:3]}")
        bad = bad or bool(f2)
    Ea, Sa = env["A"], srun["A"]
    over = [t for t in range(2, T + 1) if Ea[t] > max(Sa[t], B[t], t)]
    off = [t for t in range(2, T + 1) if max(Sa[t], B[t]) >= t
           and Ea[t] != max(Sa[t], B[t])]
    check("P3 arm A: E_run <= max(S_run, B, t), equal to max(S_run, B)"
          " where that is >= t", not over and not off,
          f"{len(over)} and {len(off)} t off")
    bad = bad or bool(over or off)

    print("READINGS")
    small = [p for p in eseats if p <= T]
    print(f"  record: E(10^4) = {E[T]}, {len(small)} seats <= 10^4, "
          f"B(10^4) = {B[T]}")
    for arm, (pairs, entry) in runs.items():
        Ex, Sx = env[arm], srun[arm]
        held = sorted(entry)
        low = [p for p in held if p <= T]
        tm = max(range(2, T + 1), key=lambda t: Ex[t] / t)
        kinds = {}
        for t in range(2, T + 1):
            k = "S_run" if Ex[t] == Sx[t] else (
                "B" if Ex[t] == B[t] else "neither")
            kinds[k] = kinds.get(k, 0) + 1
        print(f"  arm {arm}: {len(pairs)} seats, {len(low)} held <= 10^4, "
              f"largest held {held[-1]}; E_run(10^4) = {Ex[T]}, "
              f"S_run(10^4) = {Sx[T]}, E_run/E = {Ex[T] / E[T]:.4f}; "
              f"largest E_run/t {Ex[tm] / tm:.4f} = {Ex[tm]}/{tm}; "
              f"E_run equals {kinds}")
        lt = sum(Ex[t] < E[t] for t in range(2, T + 1))
        gt = sum(Ex[t] > E[t] for t in range(2, T + 1))
        print(f"    against E over 2..10^4: below at {lt} t, above at {gt}")
        up = [t for t in range(2, T + 1) if Ex[t] > E[t]]
        if up:
            print(f"    above E from t = {up[0]} to {up[-1]}; least E_run/E "
                  f"{min(Ex[t] / E[t] for t in range(3, T + 1)):.4f} over 3..10^4")
        sub = set(entry) <= set(eseats)
        same = low == [p for p in eseats if p <= T]
        print(f"    every held prime a record seat: {sub}; held <= 10^4 "
              f"the record's seats <= 10^4: {same}")
    bar = E[T] / 10
    below = [arm for arm in ("A", "B") if 10 * env[arm][T] < E[T]]
    print("THE PRINT")
    print(f"  BELOW E/10 ({bar:.1f}): "
          + (" and ".join(below) if below else "neither"))
    print("  " + ("RUN KILL" if bad else "no run kill"))
    print(f"P5 arm A below E at 10^4: {env['A'][T] < E[T]}; "
          f"P6 arm B above arm A at 10^4: {env['B'][T] > env['A'][T]}")
    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()