single.py

"""single.py -- the least multiple of the dam's bar under which one run
escapes from the void.

QUESTION. A c-run is a run from the void whose every move from a state
of tide m0 is at most floor(c F(max(m0, 4))), F the dam's bar; a run's
envelope stays within cF at every t exactly when it is a c-run. c0 is
the infimum of the c at which some c-run escapes. Where is c0?

DESIGN. F from fwide.py's widened tau table, fixed on 4..10^4.
(1) Lemma D, the crossing. For a prime r and e >= 1 write d = r^e and its
point T(d) = min(r^(e+1), p(d)) at odd r, min(2^(e+2), p(d)) at 2. At
a state with v_r(lambda) <= e - 1, m0 <= T(d): r's price is at most the
point, and p(d), unheld, prices itself. An escaping run's tide is
unbounded, so it passes every crossing, by a move from a tide <= T(d), at
most X(d) = floor(cF(T(d))), that brings a prime = 1 mod d, every
earlier prime being at most X(d) too. Raising r's own exponent instead
leaves the tide at most r (odd r past depth e with no anchor), at most
2 (2 at b >= e + 2 with u_2 <= b - 2) or at most 4 (the move 4 at
e = 1), until a prime = 1 mod d enters by a move at most X(d).
(2) Lemma H, a held crossing. A held odd q with v_q(lambda) <= f - 1 has
m0 <= q^f, its price q^max(1, v - a + 2) <= q^(f + 1 - a).
(3) Viability. Y*(a), the least Y such that a, entering with every
earlier prime at most Y, has each crossing a^f (a^f <= 10^4) passable: by
a prime b = 1 mod a^f held before a (b <= Y) or after (b <=
Z = floor(cF(a^f)), by (2)), b's earlier primes at most max(Y, Z). So
need_f(a) = 0 if some such b <= Z has Y*(b) <= Z, else the least
max(b, Y*(b)); Y*(a) = max_f need_f(a), and 0 past 10^4. An
over-approximation of viability, so sound for a dam.
(4) THE CROSSING TEST at c: every crossing d with T(d) <= 10^4 has a prime
a = 1 mod d with a <= X(d) and Y*(a) <= X(d). A failing crossing dams
every c-run: a proof at that c, given F's table. It is monotone in c;
c_low, the least c at which it passes, is found by bisection over
exact fractions and named as the critical ratio n/F(x) above the last
failing c.
(5) THE CROSSING-LED RUN, one prime per move at exponent 1 from the void:
at each state the tide's class (anchor.py's modulus and cls), the
candidates its unheld primes P <= floor(cF(m0)) with Y*(P) <=
floor(cF(m0)); it seats the least, and a state with none stops it.
c_up is the least c on the grid 2.30, 2.35, ..., 6.75 at which its
tide passes 10^4.

PREDICTIONS (frozen before the code).
P1 (proved by hand) the crossing test fails at every c < 53/23: c_low >=
   53/23, and at c = 2.30 the 4-crossing fails with Y*(5), Y*(13), Y*(17)
   > 25.
P2 (guess, 50%) c_low < 3.
P3 (guess, 40%) c_up < 6.7315: some crossing-led run beats the chain-led.
P4 (guess, 30%) c_up <= 1.5 c_low.

CONTROLS, read first. C1 the crossing test fails at c = 1 (F is dammed
from the void at 23). C2 it passes at c = 10831/1609 (the chain-led
run's moves stay within that multiple of F to tide 10^4). C3 Lemma D
over cgreedy's brute box (the primes <= 13) and 2000 random states on
the primes < 300: m0 <= T(d) at every state and every d = r^e <= 10^4,
r < 60, e <= 4, with v_r(lambda) <= e - 1. C4 Y* against a brute recursion
(viable(a, Y) by its definition, bisected on Y) on the primes <= 200
at c = 3.

THE PRINT. "c0 >= <c_low>, PROVED BY THE CROSSING TEST; A c-RUN AT <c_up>
REACHES TIDE 10^4" (frozen as "c0 IN [<c_low>, <c_up>]", which read
the run as an upper bound past 10^4); c_up / c_low above 2 prints
"KILL: THE CROSSING TEST DOES NOT LOCATE c0", and no crossing-led run
below 6.75 prints "KILL: NO CROSSING-LED RUN BELOW THE CHAIN-LED".
(Frozen with "door" where every line now says "crossing".)

FINDINGS (from the printed run; the controls passed first: 309 crossings
with T(d) <= 10^4; C1's first failing crossing at c = 1 is d = 2, T 3, X 5;
C3 over 5780 states and 221,864 (state, crossing) pairs; C4 equal on the
45 primes <= 200).
At c = 2.30 the failing crossings are 2 (X 11), 4 (X 25), 16, 3, 5 and 41,
with Y*(5) = 11251, Y*(13) = 131, Y*(17) = 127739: P1 holds.
c_low = 739/233 = 3.171674, the test failing at 3.1716738 at one crossing
only, d = 41: T(41) = min(1681, 83) = 83, X = 738, and its one
candidate 83 has Y* = 997. While 41 does not divide lambda, 83 prices
itself, so the tide is at most 83 and a prime = 1 mod 41 enters by a
move at most 233c; below 739/233 that is 83 alone. 83's crossing, by a
move at most 738, reaches 167 (Y* 2339, past 167's own reach 2039)
and 499 (Y* 10979: its one anchor in reach, 1997, has Y* 87869 =
p(1997), so the anchor 10979 must be held first); else 997 is held
first. P2 FAILS: c_low is 3.1717, not below 3.
The crossing-led run is STUCK at the void below 2.95 and at tide 83, the
41-crossing, from 2.95 to 3.15; at 3.20 it ESCAPES, tide 10,007 after 503
seats: c_up = 3.20, P3 holds (below 6.7315) and P4 holds (ratio
1.0089). READINGS: at c_low = 739/233 itself the crossing-led run ESCAPES,
tide 10,007 after 503 seats; re-tided by cgreedy.m0 on the full factor
map it AGREES at all 504 states, the void included, and its largest
seat / F(m0) is 739/233 exactly, the seat 739 = 18*41 + 1 taken from
the tide 83.
THE PRINT: c0 >= 739/233, PROVED BY THE CROSSING TEST; A c-RUN AT 3.20
REACHES TIDE 10^4. No kill fires. With the reading, c_low is met: no
run from the void escapes within cF for c < 739/233 (a proof, given
F's table and the lemmas), and the crossing-led run at c = 739/233 stays
within cF until its tide passes 10^4.
RUN RECORD. First run: C3 read p(d) past the sieve at d = r^e > 10^4
and was cut to the crossings the test uses, d <= 10^4; then 6/6, 34.5 s,
60.4 MB. With the c_low run and the Y* readings: 6/6, 34.3 s,
60.7 MB. With the re-tiding: 6/6, 34.7 s, 60.9 MB. With the void
re-tided too (it had been set to 3 by hand) and THE PRINT reworded:
6/6, 35.1 s, 61.3 MB. From code/, THE PRINT giving c_low as a
fraction and the largest seat placed: 6/6, 34.9 s, 60.9 MB. With
"door" renamed "crossing": 6/6, 34.7 s, 60.8 MB. With each Y* reading's
reach and its primes = 1 mod a within it printed (2039 and 1997 had
been worked by hand): 6/6, 34.7 s, 60.9 MB.

    python single.py
"""

import os
import random
import time
from fractions import Fraction as Fr

os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import cgreedy  # noqa: E402
import envelope as ev  # noqa: E402
import choosy  # noqa: E402
import dambar  # noqa: E402
import fwide  # noqa: E402
import anchor  # 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 sieve(n):
    s = bytearray([1]) * (n + 1)
    s[0] = s[1] = 0
    for i in range(2, int(n ** 0.5) + 1):
        if s[i]:
            s[i * i::i] = bytearray(len(range(i * i, n + 1, i)))
    return s


class World:
    def __init__(self, F):
        self.F = F
        self.W = 7 * F[T] + 7
        self.S = sieve(self.W)
        self.small = [q for q in range(3, T + 1) if self.S[q]]

    def Fx(self, t):
        return self.F[max(t, 4)]

    def cap(self, c, t):
        return (c.numerator * self.Fx(t)) // c.denominator

    def p(self, s):
        n = s + 1
        while not self.S[n]:
            n += s
        return n

    def ystar(self, c):
        """{a: Y*(a)} for the odd primes a <= T; W + 1 stands for no
        Y within the sieve."""
        Y = {}
        INF = self.W + 1
        for a in reversed(self.small):
            worst = 0
            s = a
            while s <= T:
                Z = self.cap(c, s)
                need, b = INF, s + 1
                while b <= self.W:
                    if self.S[b]:
                        yb = Y.get(b, 0)
                        if b <= Z and yb <= Z:
                            need = 0
                            break
                        need = min(need, max(b, yb))
                    if b >= need:
                        break
                    b += s
                worst = max(worst, need)
                if worst >= INF:
                    break
                s *= a
            Y[a] = worst
        return Y

    def crossings(self):
        out = []
        for r in [2] + self.small:
            e = 1
            while True:
                d = r ** e
                if d > T:
                    break
                pt = 2 ** (e + 2) if r == 2 else r ** (e + 1)
                Td = min(pt, self.p(d))
                if Td <= T:
                    out.append((d, r, e, Td))
                e += 1
        return out

    def test(self, c, Y=None, why=False):
        Y = self.ystar(c) if Y is None else Y
        bad = []
        for d, r, e, Td in self.CROSSINGS:
            X = self.cap(c, Td)
            ok = False
            a = d + 1
            while a <= X:
                if self.S[a] and Y.get(a, 0) <= X:
                    ok = True
                    break
                a += d
            if not ok:
                bad.append((d, Td, X))
                if not why:
                    return False, bad
        return not bad, bad


def crossing_run(w, c, Y, log=None):
    """The crossing-led run: (escaped, final tide, seats); log gathers
    (m0, seat) per move."""
    held, vN, seats = set(), {}, 0
    while True:
        m0, r, kind = choosy.tide(held, vN)
        if m0 > T:
            return True, m0, seats
        s = anchor.modulus(m0, r, kind, vN)
        X = w.cap(c, m0)
        P = None
        n = s + 1
        while n <= X:
            if w.S[n] and n not in held and Y.get(n, 0) <= X:
                P = n
                break
            n += s
        if P is None:
            return False, m0, seats
        held.add(P)
        seats += 1
        if log is not None:
            log.append((m0, P))
        for q, e in choosy.factor(P - 1).items():
            if e > vN.get(q, 0):
                vN[q] = e


def brute_y(w, c, top):
    """Y* by its definition on the primes <= top, recursing upward."""
    memo = {}

    def viable(a, Yb):
        if a > T:
            return True
        key = (a, Yb)
        if key in memo:
            return memo[key]
        memo[key] = False
        s, ok = a, True
        while s <= T and ok:
            Z = w.cap(c, s)
            M = max(Yb, Z)
            ok = any(w.S[b] and viable(b, M) for b in range(s + 1, M + 1, s))
            s *= a
        memo[key] = ok
        return ok

    out = {}
    for a in [q for q in w.small if q <= top]:
        lo, hi = 0, w.W + 1
        if viable(a, 0):
            out[a] = 0
            continue
        if not viable(a, w.W):
            out[a] = w.W + 1
            continue
        lo, hi = 0, w.W
        while hi - lo > 1:
            mid = (lo + hi) // 2
            if viable(a, mid):
                hi = mid
            else:
                lo = mid
        out[a] = hi
    return out


def main():
    t0 = time.time()
    fwide.widen(fwide.PMAX)
    F = dambar.f_table(fwide.tau_wide(fwide.PMAX, T))
    assert all(F[t] is not None for t in range(4, T + 1))
    w = World(F)
    w.CROSSINGS = w.crossings()
    print(f"F fixed on 4..{T}; {len(w.CROSSINGS)} crossings with T(d) <= {T}; "
          f"sieve to {w.W} ({time.time() - t0:.1f} s)")

    print("CONTROLS")
    ok1, bad1 = w.test(Fr(1), why=True)
    check("C1 the crossing test fails at c = 1", not ok1,
          f"first failing crossing {bad1[0] if bad1 else None}")
    cch = Fr(10831, 1609)
    ok2, bad2 = w.test(cch, why=True)
    check("C2 the crossing test passes at c = 10831/1609", ok2,
          f"failing {bad2[:3]}")
    rng = random.Random(1)
    states = cgreedy.brute_box(13)
    pr = [p for p in cgreedy.PR if p < 300]
    for _ in range(2000):
        mf = {2: rng.randint(0, 8)} if rng.random() < 0.7 else {}
        for q in rng.sample(pr[1:], rng.randint(1, 12)):
            mf[q] = rng.randint(1, 3)
        states.append({q: a for q, a in mf.items() if a})
    viol, tested = [], 0
    rs = [r for r in cgreedy.PR if r < 60]
    for mf in states:
        lf = cgreedy.lam(mf)
        m0 = cgreedy.m0(mf)
        for r in rs:
            for e in range(1, 5):
                d = r ** e
                if d <= T and lf.get(r, 0) <= e - 1:
                    pt = 2 ** (e + 2) if r == 2 else r ** (e + 1)
                    tested += 1
                    if m0 > min(pt, w.p(d)):
                        viol.append((mf, d))
    check("C3 Lemma D: m0 <= T(d) wherever v_r(lambda) <= e - 1",
          not viol, f"{len(states)} states, {tested} (state, crossing) pairs, "
          f"first violation {viol[:1]}")
    Y3 = w.ystar(Fr(3))
    B3 = brute_y(w, Fr(3), 200)
    diff = [a for a in B3 if min(B3[a], w.W + 1) != min(Y3[a], w.W + 1)]
    check("C4 Y* against the brute recursion, primes <= 200 at c = 3",
          not diff, f"{len(B3)} primes, differing {diff[:5]}")
    if not all(ok for _, ok in CHECKS):
        print("CONTROL FAILED: nothing below is read.")
        raise SystemExit(1)

    print("THE CROSSING TEST")
    c23 = Fr(23, 10)
    Y = w.ystar(c23)
    ok, bad = w.test(c23, Y, why=True)
    four = [b for b in bad if b[0] == 4]
    ys = {a: Y.get(a) for a in (5, 13, 17)}
    p1a = (not ok and bool(four) and four[0][2] == 25
           and all(v > 25 for v in ys.values()))
    print(f"  c = 2.30: failing crossings {bad[:6]}; Y*(5, 13, 17) = {ys}")
    lo, hi = Fr(1), cch
    while hi - lo > Fr(1, 10 ** 7):
        mid = (lo + hi) / 2
        mid = Fr(round(mid * 10 ** 8), 10 ** 8)
        if w.test(mid)[0]:
            hi = mid
        else:
            lo = mid
    fv = sorted(set(F[t] for t in range(4, T + 1)))
    crit = min(Fr(int(lo * f) + 1, f) for f in fv)
    okc = crit <= hi and w.test(crit)[0] and not w.test(lo)[0]
    check("the critical ratio above the last failing c passes", okc,
          f"c_low = {crit} = {float(crit):.6f}, last failing {float(lo):.7f}")
    Yl = w.ystar(lo)
    _, badl = w.test(lo, Yl, why=True)
    check("P1 c_low >= 53/23, the 4-crossing failing at 2.30 with "
          "Y*(5), Y*(13), Y*(17) > 25", p1a and crit >= Fr(53, 23))
    print(f"  crossings failing just below c_low: {badl[:8]}")
    for d, Td, X in badl[:4]:
        cands = [(a, Yl.get(a, 0)) for a in range(d + 1, X + 1, d)
                 if w.S[a]]
        print(f"    crossing {d}: T {Td}, X {X}, candidates (a, Y*) "
              f"{cands[:10]}")
    print(f"  P2 (guess) c_low < 3: "
          f"{'HOLDS' if crit < 3 else 'FAILS'}, c_low {float(crit):.4f}")

    print("THE CROSSING-LED RUN")
    grid = [Fr(230 + 5 * i, 100) for i in range(int((675 - 230) / 5) + 1)]
    c_up, rows = None, []
    for c in grid:
        Yc = w.ystar(c)
        esc, tide, seats = crossing_run(w, c, Yc)
        rows.append((c, esc, tide, seats))
        if esc:
            c_up = c
            break
    for c, esc, tide, seats in rows[::max(1, len(rows) // 12)] + rows[-1:]:
        print(f"  c = {float(c):.2f}: {'ESCAPES' if esc else 'STUCK'} at "
              f"tide {tide} after {seats} seats")
    if c_up is not None:
        print(f"  P3 (guess) c_up < 6.7315: "
              f"{'HOLDS' if c_up < cch else 'FAILS'}, c_up {float(c_up):.2f}")
        print(f"  P4 (guess) c_up <= 1.5 c_low: "
              f"{'HOLDS' if c_up <= Fr(3, 2) * crit else 'FAILS'}, "
              f"ratio {float(c_up / crit):.4f}")

    print("READINGS")
    log = []
    esc, tide, seats = crossing_run(w, crit, w.ystar(crit), log)
    print(f"  the crossing-led run at c_low = {crit} itself: "
          f"{'ESCAPES' if esc else 'STUCK'} at tide {tide} after {seats} "
          f"seats")
    held, worst, at, agree = set(), Fr(0), None, True
    for m0, P in log:
        m_full = cgreedy.m0({p: 1 for p in held})
        agree = agree and m_full == m0
        if Fr(P, w.Fx(m_full)) > worst:
            worst, at = Fr(P, w.Fx(m_full)), (m_full, P)
        held.add(P)
    agree = agree and cgreedy.m0({p: 1 for p in held}) == tide
    print(f"  re-tided by cgreedy.m0 on the full factor map: "
          f"{'AGREES' if agree else 'DIFFERS'} at all {len(log) + 1} states; "
          f"largest seat / F(m0) {worst} = {float(worst):.6f}, within "
          f"c_low: {worst <= crit}, at (tide, seat) {at}")
    for a in (83, 167, 499, 1997):
        reach = w.cap(lo, a)
        anchors = [b for b in range(2 * a + 1, reach + 1, 2 * a) if w.S[b]]
        print(f"  Y*({a}) just below c_low: {Yl.get(a)}; its reach "
              f"{reach}, the primes = 1 mod {a} within it {anchors}")

    print("THE PRINT")
    if c_up is None:
        print("  KILL: NO CROSSING-LED RUN BELOW THE CHAIN-LED")
    else:
        print(f"  c0 >= {crit}, PROVED BY THE CROSSING TEST; A "
              f"c-RUN AT {float(c_up):.2f} REACHES TIDE 10^4")
        if c_up / crit > 2:
            print("  KILL: THE CROSSING TEST DOES NOT LOCATE c0")
    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()