crosswide.py

"""crosswide.py -- the crossing test and the crossing-led run past 10^4,
F read by the feeder formula.

QUESTION. single.py proves c0 >= 739/233 by the crossing test over the
crossings with T(d) <= 10^4 and meets it with the crossing-led run to
tide 10,007; it stops at 10^4 only because F came from fwide.py's tau
table. The feeder formula gives F(t) = CB(T+(t)) at any t (pstar.py
(6)). At c = 739/233 and the horizon H = 10^5: does the crossing test
pass over every crossing with T(d) <= H, and does the crossing-led run
reach a tide past H within cF?

DESIGN.
(1) F(t) = CB(T+(t)) on 4..H, one pstar.Engine per prime power T, one
spf to H + 1000.
(2) single.py's World, test, ystar and crossing_run by import, the
horizon a parameter (single.py's own run at H = 10^4 unchanged). At the
horizon H the test reads the crossings d with T(d) <= H, viability Y*
reads the crossings b^f <= H, and every prime past H counts viable
(Y* = 0). That over-approximation is still sound for a dam, so a
failing crossing at H proves c0 above that c, given F to H and the
lemmas; it is weaker farther out, primes in (10^4, H] now read and
primes past H free, so a pass extends the run's range, not the bound.
(3) MONOTONE IN H (property). Y*_H(a) is nondecreasing in H: by
downward induction on b, a larger H adds crossings to pass and turns
Y* = 0 past the old horizon into a value at least 0, and need_f only
grows with the Y* it reads. Capped at the smaller horizon's W + 1 the
comparison holds for c <= 7 (so every Z and X lies within that sieve),
since a b read only past the smaller sieve has max(b, Y*(b)) past it. So the test at H passing implies it passes at
10^4, and c_low(H) >= 739/233; a pass at 739/233 makes c_low(H) =
739/233 exactly.
(4) A failing crossing at 739/233: c_low(H) by single.py's bisection
over exact fractions above 739/233 to 10^-7, then a binary search over
every ratio n/F(x) in that window (the only c where a cap moves), and
the failing crossings just below it with their candidates (a, Y*).
(5) The crossing-led run at c_low(H) from the void, until its tide
passes H or no candidate is left; stuck, the grid 3.20, 3.25, ... of
single.py from c_low(H) up, within a 30-minute budget, the cut printed.

PREDICTIONS (frozen before the code).
P1 (forced) CB-built F equals fwide.py's F at every t in 4..10^4.
P2 (forced) at H = 10^4 the test reproduces single.py: 309 crossings,
   passing at 739/233, failing at 3.1716738 at the crossing 41 alone;
   the run at 739/233 escapes at tide 10,007 after 503 seats.
P3 (the random-model heuristic, 80%) the crossing test at 739/233
   passes over every crossing with T(d) <= 10^5. The random model
   fails a Sophie Germain crossing near 5*10^4 with probability about
   exp(-3cF(p)/(p ln p)) ~ 10^-7 there, but the 41 failure came through
   viability chains the model does not price, so 80%, not ~100%.
P4 (a guess, 60%) the crossing-led run at 739/233 passes tide 10^5.
P5 (a guess, 50%) between 1,500 and 3,500 crossings have T(d) <= 10^5.

THE KILL, frozen. The script prints, over the crossings with T(d) <=
10^5 at c = 739/233, either "THE CROSSING TEST PASSES AT 739/233 OVER
<n> CROSSINGS WITH T(d) <= 10^5" or "KILL: THE CROSSING TEST FAILS AT
739/233 TO 10^5; c0 >= <c_low(H)>, PROVED BY THE CROSSING TEST", the
latter naming the first failing crossing (d, T(d), X(d)). (Frozen as "A
CROSSING PAST 10^4 FAILS", which presumed where a failure lands; the
rehearsal at H = 20,000 failed at the crossing 2, its candidates'
viability read past 10^4.) The KILL is of P3,
the heuristic's "finitely many failing crossings at 739/233" read on
its first decade past the table. The confound, (2)'s free primes past
H, inflates the pass arm, so the pass is read as a range for the run,
never a bound, and the bar's proof sits on the fail arm, which the
confound cannot inflate. Then "THE CROSSING-LED RUN AT <c> REACHES TIDE
<tide>" or "... STOPS AT TIDE <tide>".

CONTROLS, read first. C1 is P1. C2 is P2 (the World at H = 10^4 with
CB's F). C3 (3) at c = 739/233: min(Y*_H(a), W_4 + 1) >=
min(Y*_4(a), W_4 + 1) at every odd prime a <= 10^4. C4 the test at H
fails at 3.1716738 (it must, by (3)), the crossing 41 among its
failures. C5 the test at H fails at c = 1.

FINDINGS (from the printed run; the controls passed first: C1 0 of
4..10^4 differ; C2 309 crossings, failing at 3.1716738 at (41, 83, 738)
alone, the run at 10,007 after 503 seats; C3 over 1228 primes, 2 rising;
C4 failing at (2, 3, 15) and (41, 83, 738); C5 at (2, 3, 5)).
F on 4..10^5 from 9698 prime powers T, F(10^5) = 2021219; the sieve to
14,148,540; 1922 crossings with T(d) <= 10^5, 1613 of them past 10^4
(P5 holds).
THE PRINT: KILL: THE CROSSING TEST FAILS AT 739/233 TO 10^5; c0 >=
11251/3373, PROVED BY THE CROSSING TEST (first failing crossing d = 2,
T(d) = 3, X(d) = 15), the one failing crossing; c_low(H) = 3.3356063.
P3 FAILS, and P4 with it: a failing
crossing dams every c-run at 739/233. The failure is not in the new
decade: the tightest crossing past 10^4 needs max(a, Y*(a))/F(T(d)) =
0.9999 (d = 30323, a = 788399), far below 739/233. It is (3) at work:
the 2-crossing's candidates within X(2) are 3, 5, 7, 11 and 13, and
of the 1228 odd primes <= 10^4 only two Y* rise between the horizons,
Y*(11) 0 -> 673487 and Y*(2153) 0 -> 51673 (12919 = p(2153), free at
10^4, is read at H). 11's crossing 11^4 = 14641 enters the horizon at 14641, and
the least prime = 1 mod 11^4 is 673487, past Z = 540221 at 739/233: so
11, the one candidate viable at 10^4 (deduced: the test passed there,
and the other four did not rise), needs 673487 held first. Just
below c_low(H) the others need: Y*(3) = 31 (3's own crossing, Z 16,
its feeders 7 and 13 not viable within it), Y*(5) = 11251 (5^4 = 625,
Z 11250, 11251 = 18*625 + 1), Y*(7) = 491 (7^2, Z 463, its feeder 197
not viable within it), Y*(13) = 131 (13, Z 76, its feeder 53 not):
only 5 and 11 meet a crossing with no prime = 1 mod it within reach. c_low(H) = 11251/3373 is
5's crossing 5^4 priced at F(625) = 3373, exact: it passes and the last
failing c, 5276759/1581949, fails, 2 ratios n/F(x) and the bisection's
hi searched (printed as 3 points).
The crossing-led run at c_low(H) = 11251/3373 ESCAPES (its tide passes
H), tide 100,043 after 3386 seats; its largest seat / F(m0) is 11251/3373 exactly, the
seat 11251 from the tide 625; re-tided by cgreedy.m0 it AGREES at the
878 states with tide below 19997, where cgreedy's prime table ends,
the 2508 past it resting on choosy.tide alone (its table to 700,000). So
c0 >= 11251/3373 is met to 10^5: no run escapes within cF for c below
it (a proof, given F by the feeder formula and the lemmas), and the
crossing-led run at c = 11251/3373 stays within cF until its tide passes
10^5.
RUN RECORD. Rehearsal --H 20000: 5/5, 41.2 s, 66.0 MB, the test
failing at the crossing 2 (the kill's line reworded, above). First full
run: 5/5, 167.7 s, 84.7 MB; its c_low(H), 5276759/1581949, was the
least ratio n/F(x) above the bisection's last failing c, single.py's
rule, and failed, so the run sat at the void: replaced by a binary
search over every ratio in the bisection's window, checked. Then 6/6,
167.2 s, 83.8 MB. With the run's largest seat and its re-tiding (cut to
the tides below growth's prime table after the first try re-tided the
final tide 20123 as 22201 = 149^2, a prime past the table unseen): 6/6,
167.9 s, 84.4 MB, every figure above unchanged. With the searched count
labelled points: 6/6, 169.5 s, 83.4 MB, unchanged. From code/: 6/6,
166.5 s, 84.4 MB, unchanged.

    python crosswide.py [--H N]     (N the horizon, default 10^5; a
                                     rehearsal runs --H 20000)
"""

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

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

H = int(sys.argv[sys.argv.index("--H") + 1]) if "--H" in sys.argv else 10 ** 5
T4 = single.T
C = Fr(739, 233)
LO4 = Fr(31716738, 10 ** 7)       # single.py's last failing c at 10^4
BUDGET = 30 * 60                  # (5)'s grid, seconds
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 f_feeder(top):
    """F on 4..top by (1): (F, the number of prime powers T read)."""
    spf = pstar.spf_to(top + 1000)
    primes = [p for p in range(2, top + 1001) if spf[p] == 0]
    F, memo = [None] * (top + 1), {}
    for t in range(4, top + 1):
        T = pstar.next_pp(t, spf)
        if T not in memo:
            memo[T] = pstar.Engine(T, spf, primes).feeder_bound()[0]
        F[t] = memo[T]
    return F, len(memo)


def world(F, top):
    w = single.World(F, top)
    w.CROSSINGS = w.crossings()
    return w


def c_low(w, lo, hi):
    """single.py's bisection to 10^-7, then a binary search over every
    ratio n/F(x) in (lo, hi], the only c where a cap moves: (c_low, the
    largest c read below it, which fails, and the number of ratios
    searched)."""
    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 = set(w.F[t] for t in range(4, w.T + 1))
    br = sorted({Fr(n, f) for f in fv
                 for n in range(int(lo * f) + 1, int(hi * f) + 1)} | {hi})
    i, j = -1, len(br) - 1             # fails at lo, passes at br[j]
    while j - i > 1:
        m = (i + j) // 2
        if w.test(br[m])[0]:
            j = m
        else:
            i = m
    return br[j], (br[i] if i >= 0 else lo), len(br)


def margins(w, c, Y):
    """Per crossing: the least max(a, Y*(a)) / F(T(d)) over its primes
    a = 1 mod d within X(d) with Y*(a) <= X(d), Y* read at c; None when
    none passes."""
    out = []
    for d, r, e, Td in w.CROSSINGS:
        X = w.cap(c, Td)
        best, by = None, None
        for a in range(d + 1, X + 1, d):
            if w.S[a] and Y.get(a, 0) <= X:
                m = Fr(max(a, Y.get(a, 0)), w.Fx(Td))
                if best is None or m < best:
                    best, by = m, a
        out.append((d, Td, X, best, by))
    return out


def why_y(w, c, Y, a):
    """Y*(a)'s crossings at c, as ystar reads them: [(s = a^f, Z,
    need_f(a), the b giving it, Y*(b))]."""
    out, s = [], a
    while s <= w.T:
        Z = w.cap(c, s)
        need, arg, b = w.W + 1, None, s + 1
        while b <= w.W:
            if w.S[b]:
                yb = Y.get(b, 0)
                if b <= Z and yb <= Z:
                    need, arg = 0, b
                    break
                if max(b, yb) < need:
                    need, arg = max(b, yb), b
            if b >= need:
                break
            b += s
        out.append((s, Z, need, arg, Y.get(arg, 0) if arg else None))
        s *= a
    return out


def main():
    t0 = time.time()
    F, nT = f_feeder(H)
    print(f"F by the feeder formula on 4..{H}: {nT} prime powers T, "
          f"F({H}) = {F[H]} ({time.time() - t0:.1f} s)")

    print("CONTROLS")
    fwide.widen(fwide.PMAX)
    F4 = dambar.f_table(fwide.tau_wide(fwide.PMAX, T4))
    off = [t for t in range(4, T4 + 1) if F4[t] != F[t]]
    check(f"C1 = P1 CB-built F equals fwide.py's F on 4..{T4}", not off,
          f"{len(off)} differ {off[:3]}")
    t1 = time.time()
    w4 = world(F, T4)
    ok4, _ = w4.test(C, why=True)
    okl, badl = w4.test(LO4, why=True)
    Y4 = w4.ystar(C)
    esc, tide, seats = single.crossing_run(w4, C, Y4)
    check(f"C2 = P2 at H = {T4}: 309 crossings, the test passing at "
          f"739/233 and failing at {float(LO4)} at 41 alone; the run "
          f"escaping at tide 10007 after 503 seats",
          len(w4.CROSSINGS) == 309 and ok4 and not okl
          and [b[0] for b in badl] == [41] and esc and tide == 10007
          and seats == 503,
          f"{len(w4.CROSSINGS)} crossings, pass {ok4}, failing {badl}, "
          f"run {esc} {tide} {seats} ({time.time() - t1:.1f} s)")

    t1 = time.time()
    w = world(F, H)
    print(f"  the world at H = {H}: sieve to {w.W}, {len(w.CROSSINGS)} "
          f"crossings with T(d) <= {H} ({time.time() - t1:.1f} s)")
    t1 = time.time()
    Y = w.ystar(C)
    print(f"  Y* at 739/233 to H: {time.time() - t1:.1f} s")
    cap4 = w4.W + 1
    low = [a for a in w4.small if min(Y[a], cap4) < min(Y4[a], cap4)]
    up = sum(min(Y[a], cap4) > min(Y4[a], cap4) for a in w4.small)
    check("C3 (3) Y*_H >= Y*_4 at every odd prime <= 10^4, at 739/233, "
          "capped at W_4 + 1", not low,
          f"{len(w4.small)} primes, {up} rise, lower at {low[:5]}")
    rise = [a for a in w4.small if min(Y[a], cap4) > min(Y4[a], cap4)]
    for a in rise[:8]:
        top = max(why_y(w, C, Y, a), key=lambda x: x[2])
        print(f"    rises: Y*({a}) {Y4[a]} -> {Y[a]}, from its crossing "
              f"s = {top[0]}: Z {top[1]}, need {top[2]} by b = {top[3]}, "
              f"Y*(b) {top[4]}")
    t1 = time.time()
    okH, badH = w.test(LO4, why=True)
    check(f"C4 the test at H fails at {float(LO4)}, 41 among its failures",
          not okH and 41 in [b[0] for b in badH],
          f"{len(badH)} failing, first {badH[:4]} ({time.time() - t1:.1f} s)")
    ok1, bad1 = w.test(Fr(1))
    check("C5 the test at H fails at c = 1", not ok1, f"{bad1[:1]}")
    if not all(ok for _, ok in CHECKS):
        print("CONTROL FAILED: nothing below is read.")
        raise SystemExit(1)

    print("THE CROSSING TEST AT 739/233 TO 10^5")
    t1 = time.time()
    ok, bad = w.test(C, Y, why=True)
    past = [x for x in w.CROSSINGS if x[3] > T4]
    print(f"  {len(w.CROSSINGS)} crossings, {len(past)} with T(d) in "
          f"({T4}, {H}]; {len(bad)} failing ({time.time() - t1:.1f} s)")
    print(f"  P5 (guess) 1500..3500 crossings: "
          f"{'HOLDS' if 1500 <= len(w.CROSSINGS) <= 3500 else 'FAILS'}")
    mg = margins(w, C, Y)
    tight = sorted((m for m in mg if m[3] is not None and m[1] > T4),
                   key=lambda m: -m[3])
    print("  the tightest crossings past 10^4, (d, T(d), X(d), "
          "max(a, Y*(a))/F(T(d)), a):")
    for d, Td, X, m, a in tight[:6]:
        print(f"    {d}, {Td}, {X}, {float(m):.4f}, {a}")
    free = [m for m in mg if m[3] is not None and m[4] > H]
    print(f"  crossings whose least-need prime lies past H (free by (2)): "
          f"{len(free)}, the least d "
          f"{min((m[0] for m in free), default=None)}")
    if ok:
        cl = C
    else:
        t1 = time.time()
        hi = Fr(7)
        if not w.test(hi)[0]:
            print("  the test fails at c = 7 too; (4) stops")
            cl = None
        else:
            cl, lo, nbr = c_low(w, C, hi)
            _, badc = w.test(lo, why=True)
            Yl = w.ystar(lo)
            check("c_low(H) passes and the largest c read below it fails",
                  w.test(cl)[0] and bool(badc),
                  f"c_low(H) = {cl} = {float(cl):.7f}, last failing "
                  f"{lo} = {float(lo):.7f}, {nbr} points searched "
                  f"({time.time() - t1:.1f} s)")
            for d, Td, X in badc[: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 "
                      f"(a, Y*) {cands[:10]}")
                for a, ya in cands[:5]:
                    top = max(why_y(w, lo, Yl, a), key=lambda x: x[2])
                    print(f"      Y*({a}) = {ya} from its crossing s = "
                          f"{top[0]}: Z {top[1]}, need {top[2]} by b = "
                          f"{top[3]}, Y*(b) {top[4]}")
    print(f"  P3 (the heuristic) the test passes at 739/233 to 10^5: "
          f"{'HOLDS' if ok else 'FAILS'}")

    print("THE CROSSING-LED RUN")
    runs = []
    if cl is not None:
        t1 = time.time()
        Yc = Y if cl == C else w.ystar(cl)
        log = []
        runs.append((cl,) + single.crossing_run(w, cl, Yc, log))
        print(f"  c = {cl}: {'ESCAPES' if runs[-1][1] else 'STUCK'} at "
              f"tide {runs[-1][2]} after {runs[-1][3]} seats "
              f"({time.time() - t1:.1f} s)")
        if log:
            t2 = time.time()
            worst = max(log, key=lambda x: Fr(x[1], w.Fx(x[0])))
            # cgreedy.m0 reads growth.PRIMES, a table to 20,000: it
            # re-tides the states whose tide lies below that table only
            top = cgreedy.growth.PRIMES[-1]
            held, agree, n = [], True, 0
            for m0, P in log:
                if m0 < top:
                    agree = agree and cgreedy.m0({p: 1 for p in held}) == m0
                    n += 1
                held.append(P)
            print(f"    largest seat / F(m0) {Fr(worst[1], w.Fx(worst[0]))}"
                  f" at (tide, seat) {worst}; re-tided by cgreedy.m0 at "
                  f"the {n} states with tide below {top}: "
                  f"{'AGREES' if agree else 'DIFFERS'} "
                  f"({time.time() - t2:.1f} s)")
        grid = [Fr(320 + 5 * i, 100) for i in range(72)]
        for c in [g for g in grid if g > cl]:
            if runs[-1][1] or time.time() - t1 > BUDGET:
                break
            runs.append((c,) + single.crossing_run(w, c, w.ystar(c)))
            print(f"  c = {float(c):.2f}: "
                  f"{'ESCAPES' if runs[-1][1] else 'STUCK'} at tide "
                  f"{runs[-1][2]} after {runs[-1][3]} seats")
        if not runs[-1][1] and time.time() - t1 > BUDGET:
            print(f"  GRID CUT at c = {float(runs[-1][0]):.2f}, budget "
                  f"spent")
        if cl == C:
            print(f"  P4 (guess) the run at 739/233 passes 10^5: "
                  f"{'HOLDS' if runs[0][1] else 'FAILS'}")

    print("THE PRINT")
    if ok:
        print(f"  THE CROSSING TEST PASSES AT 739/233 OVER "
              f"{len(w.CROSSINGS)} CROSSINGS WITH T(d) <= 10^5")
    else:
        d, Td, X = bad[0]
        print(f"  KILL: THE CROSSING TEST FAILS AT 739/233 TO 10^5; "
              f"c0 >= {cl}, "
              f"PROVED BY THE CROSSING TEST (first failing crossing d = "
              f"{d}, T(d) = {Td}, X(d) = {X})")
    for c, esc, tide, seats in runs[-1:]:
        print(f"  THE CROSSING-LED RUN AT {c} "
              f"{'REACHES' if esc else 'STOPS AT'} TIDE {tide}")
    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()