crossall.py

"""crossall.py -- the crossings to 10^5 read at the horizon 10^6, lazily.

QUESTION. crosswide.py proves c0 >= 11251/3373 by the crossing test
over its 1922 crossings with T(d) <= 10^5 at the horizon H = 10^5,
failing just below it at the crossing 2 alone. The bound rose past 10^4
through an old crossing's candidates read deeper (11's crossing 11^4
entering the horizon), not through the new decade. At H = 10^6, does
some crossing of the 1922 fail at 11251/3373, and so lift c0?

DESIGN.
(1) The crossings, single.py's World.crossings without F: d = r^e, T(d)
= min(r^(e+1), p(d)) at odd r, min(2^(e+2), p(d)) at 2, p(d) the least
prime = 1 mod d, kept when T(d) <= the crossing bound B; World's order,
2 then the odd primes ascending, e ascending within r.
(2) The crossing test at c and the horizon H over that list, read
lazily: the crossing d passes when some prime a = 1 mod d, a <= X(d) =
floor(cF(max(T(d), 4))), has Y*(a) <= X(d), the a read ascending and
stopping at the first passing one (World.test's order); the test stops
at the first failing crossing. Y*(a) is World.ystar's at the horizon H:
0 for a > H; else the max over the crossings s = a^f <= H of need_f(a),
need_f(a) = 0 if some prime b = 1 mod s has b <= Z = floor(cF(s)) and
Y*(b) <= Z, else the least max(b, Y*(b)), the scan of b = s + 1,
2s + 1, ... stopping once b reaches the least need so far. Read by
memoized recursion from the candidates the test reaches, fresh per c
and shared by every crossing, since Y*(a) does not depend on the
crossing that reads it. Uncapped: World caps Y* at its sieve's W + 1,
and by downward induction its value is min(Y*, W + 1), every Z, X <=
7F(H) < W + 1, so every comparison the test makes is unchanged (C6
checks it at 10^4).
(3) F(s) = CB(T+(s)) by one pstar.Engine at T+(s) per prime power read,
memoized across c and horizons (pstar.py's feeder formula). Primality
by a sieve to 2 * 10^7 and growth.is_prime_mr past it.
(4) SOUND ON A SUBSET (property, single.py's Lemma D read one crossing
at a time): Lemma D dams every c-run at a crossing that fails, whatever
the others do, and Y*_H over-approximates viability at every H (it reads
only the crossings a^f <= H and counts a prime past H viable), so one
failing crossing of the 1922 at H proves c0 above that c, given F by the
feeder formula and the lemmas. The free primes past H inflate only the
pass arm.
(5) MONOTONE IN H (crosswide.py's (3)) on a fixed list: the test at
10^6 failing at c' implies nothing at 10^5, but passing at 10^6 implies
passing at 10^5, so c_low(10^6) >= c_low(10^5) = 11251/3373.
(6) c_low EXACTLY. The test is monotone in c (Y* nonincreasing, X
nondecreasing). Bisection over exact fractions in [lo, hi] to 10^-7;
then every ratio n/F(s) in the window (lo, hi], s over the caps read at
both ends (each read crossing's own cap max(T(d), 4) among them),
binary-searched to an adjacent pair c' failing, c* passing. An
evaluation at c is a deterministic function of the caps it reads, so if
no ratio n/F(s), s read at c', lies in (c', c*), every c in [c', c*)
reads c''s caps and fails, and c_low = c* exactly; else the missing
ratios join the list and the search repeats. The read set depends on
the order of (2); it is kept.
(7) Printed at the horizon: the failing crossings just below c_low (the
test read with every crossing, not stopping), each with its candidates
(a, Y*(a)); the crossing 2's candidates; the primes read and the Engines
built for the horizon; how many crossings pass by a prime past H, and
how many crossings a^f read inside Y* are passed by one (the free arm
of (4), on which a pass at H may rest).
(8) THE FRONTIER: the crossings with T(d) in (10^5, 10^6] are counted,
the distinct caps T+(max(T(d), 4)) their X(d) needs and not yet built
counted, and ten of them, spread evenly over the sorted list, built and
timed. They are not read.

PREDICTIONS (frozen before the code).
P1 (forced) (1) gives 309 crossings at B = 10^4 and 1922 at B = 10^5.
P2 (forced) at H = B = 10^4, c_low = 739/233, failing just below at the
   crossing 41 alone (single.py).
P3 (forced) at H = B = 10^5, c_low = 11251/3373, failing just below at
   the crossing 2 alone (crosswide.py).
P4 (guess, 35%) at H = 10^6, B = 10^5, the test fails at 11251/3373:
   the decade deepens towers like 11's at 10^5, and the 1922 hold
   many old crossings whose single passing candidate can rise.
P5 (guess, 55%) the Engines built for the horizon 10^6, over the
   whole search, number fewer than 5000.
P6 (guess, 60%) between 8,000 and 14,000 crossings have T(d) in (10^5,
   10^6] (309 to 1922 is a factor 6.2 a decade).

THE KILL, frozen. Over the 1922 crossings with T(d) <= 10^5 at H = 10^6,
the script prints "THE CROSSINGS TO 10^5 LIFT c0 TO <c> AT 10^6" or
"KILL: THE CROSSINGS TO 10^5 LEAVE c0's BOUND AT 11251/3373 AT 10^6".
A lift is a proof by (4), on the fail arm, which the free primes cannot
inflate; the kill is of P4, and a pass is read as a range, never a bound.

CONTROLS, read first. C1 (1)'s lists equal World.crossings's (d, T(d))
at B = 10^4 and 10^5, 309 and 1922 long (World.crossings reads no F;
at 10^5 its sieve is sized by a stand-in F(H)). C2 is P2 and C3 is P3,
each by (6)'s search from [1, 7] and (7)'s failing list at the c' it
returns. C4 at H = 10^5 the test fails at c = 1 and passes at c = 7.
C5 (6)'s exactness, broken: at H = B = 10^4 the first ratio list built
from the crossing 2's own cap F(4) = 5 alone, which holds no n/5 in the
window, so the check must catch a missing ratio and c_low still land at
739/233. C6 at H = B = 10^4, F by CB on 4..10^4 (crosswide.f_feeder),
the lazy Y* of (2) equals World.ystar, capped at W + 1, at every prime
the recursion reads, at c = 739/233 and at c = 3.3.

FINDINGS (from the printed run; the controls passed first: C1 the
lists equal, 309 and 1922 long; C6 over 655 (c, a) pairs; C2 739/233
by 2 points in 1 round, failing just below at the crossing 41 alone,
its one candidate (83, 997); C5 caught, 2 rounds, landing at 739/233;
C4; C3 11251/3373 by 2 points in 1 round, failing just below at the
crossing 2 alone, its candidates (3, 31), (5, 11251), (7, 491),
(11, 673487), (13, 131), as crosswide.py printed).
At H = 10^6 the 1922 crossings PASS at 11251/3373, none failing, 3262
primes read and 3392 caps, 1367 Engines built for the horizon; no
crossing of the list is passed by a candidate past H. The crossing 2
passes by 5 alone, Y*(5) = 0, 3, 7, 11 and 13 keeping Y* 31, 491,
673487 and 131, past X(2) = 16; 5's crossing 5^8 = 390,625 is passed by
9,375,001, past H. Inside Y*, 1560 of the 3387 crossings a^f read are
passed by a prime past H, at 1559 of the 3262 primes, the least such
a^f 25,589: the pass at 10^6 is the arm (4) inflates, so it is a range
for the test, and the bound is crosswide.py's. THE PRINT: KILL: THE
CROSSINGS TO 10^5 LEAVE c0's BOUND AT 11251/3373 AT 10^6. P4 FAILS,
P5 HOLDS (1367 < 5000). THE FRONTIER: 13254 crossings have T(d) <=
10^6, 11332 of them past 10^5 (P6 holds); their caps need 10028
Engines not yet built, ten timed from 0.04 s at T 100049 to 0.43 s at
T 999961, mean 0.22 s, about 38 minutes.
RUN RECORD. Rehearsal --H 200000: 5/5, 52.1 s, 142.6 MB, passing at
11251/3373 (337 Engines for the horizon). Full: 5/5, 265.8 s, peak
160.3 MB under the memory guard. Then the crossing 2's script was
folded in: its lazy Y* control as C6, its reading as the free arm's
counts of (7), C1
comparing lists where it had compared counts: rehearsal --H 200000
6/6, 53.0 s, 144.0 MB; full 6/6, 264.8 s, peak 160.3 MB, every earlier
print unchanged.

    python crossall.py [--H N] [--B N]   (horizon, default 10^6; crossing
                                          bound, default 10^5)
"""

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

os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import crosswide  # noqa: E402
import growth  # 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 ** 6
B = int(sys.argv[sys.argv.index("--B") + 1]) if "--B" in sys.argv else 10 ** 5
LO5 = Fr(11251, 3373)              # c_low(10^5) by crosswide.py
LO4 = Fr(739, 233)                 # c_low(10^4) by single.py
SIEVE = 2 * 10 ** 7
CHECKS = []
sys.setrecursionlimit(100000)


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


class Lazy:
    """F and primality by (3); the test of (2) over self.CR."""

    CR = ()

    def __init__(self, top):
        self.spf = pstar.spf_to(top + 1000)
        self.primes = [p for p in range(2, top + 1001) if self.spf[p] == 0]
        self.S = single.sieve(SIEVE)
        self.FT, self.Ts = {}, {}

    def isp(self, n):
        return self.S[n] == 1 if n <= SIEVE else growth.is_prime_mr(n)

    def F(self, t):
        T = self.Ts.get(t)
        if T is None:
            T = self.Ts[t] = pstar.next_pp(t, self.spf)
        if T not in self.FT:
            self.FT[T] = pstar.Engine(T, self.spf,
                                      self.primes).feeder_bound()[0]
        return self.FT[T]

    def crossings(self, bound):
        """(1): [(d, T(d))] with T(d) <= bound, in World's order."""
        out = []
        for r in self.primes:
            if r > bound:
                break
            e, d = 1, r
            while d <= bound:
                pt = 2 ** (e + 2) if r == 2 else r ** (e + 1)
                b = d + 1
                while b < pt and b <= bound and not self.isp(b):
                    b += d
                Td = min(pt, b) if self.isp(b) else pt
                if Td <= bound:
                    out.append((d, Td))
                e, d = e + 1, d * r
        return out

    def evaluate(self, c, h, every=False):
        """(2) at c and the horizon h: (passes, the t whose cap F(t) was
        read, Y*, the dams {a: [(s, Z, need, b, Y*(b))]}, the failing
        crossings [(d, T(d), X)], the crossings passed by a prime past
        h). every: read every crossing, not stop at the first failing."""
        num, den = c.numerator, c.denominator
        Y, dams, read = {}, {}, set()

        def ys(a):
            if a > h:
                return 0
            v = Y.get(a)
            if v is not None:
                return v
            worst, s, rec = 0, a, []
            while s <= h:
                t = max(s, 4)
                read.add(t)
                Z = num * self.F(t) // den
                need, arg, b = None, None, s + 1
                while True:
                    if self.isp(b):
                        yb = ys(b)
                        if b <= Z and yb <= Z:
                            need, arg = 0, b
                            break
                        if need is None or max(b, yb) < need:
                            need, arg = max(b, yb), b
                    if need is not None and b >= need:
                        break
                    b += s
                rec.append((s, Z, need, arg, ys(arg)))
                worst = max(worst, need)
                s *= a
            Y[a], dams[a] = worst, rec
            return worst

        bad, free = [], 0
        for d, Td in self.CR:
            t = max(Td, 4)
            read.add(t)
            X = num * self.F(t) // den
            ok, a = False, d + 1
            while a <= X:
                if self.isp(a) and ys(a) <= X:
                    ok, free = True, free + (a > h)
                    break
                a += d
            if not ok:
                bad.append((d, Td, X))
                if not every:
                    break
        return not bad, read, Y, dams, bad, free


def ratios(lz, ts, lo, hi, closed=True):
    """Every n/F(t), t in ts, in (lo, hi] (closed) or (lo, hi)."""
    out = set()
    for t in ts:
        f = lz.F(t)
        top = (hi.numerator * f) // hi.denominator
        if not closed and Fr(top, f) == hi:
            top -= 1
        out.update(Fr(n, f) for n in range((lo.numerator * f) //
                                           lo.denominator + 1, top + 1))
    return out


def exact(lz, h, lo, hi, use="both"):
    """(6): (c_low, the c' below it, ratios searched, rounds, evals).
    use = "four" builds the first list from the crossing 2's own cap
    alone (C5)."""
    ev = 0
    while hi - lo > Fr(1, 10 ** 7):
        mid = Fr(round((lo + hi) / 2 * 10 ** 8), 10 ** 8)
        ev += 1
        if lz.evaluate(mid, h)[0]:
            hi = mid
        else:
            lo = mid
    if use == "four":
        ts = {4}
    else:
        ts = lz.evaluate(lo, h)[1] | lz.evaluate(hi, h)[1]
        ev += 2
    rounds, nbr = 0, 0
    while True:
        rounds += 1
        br = sorted(ratios(lz, ts, lo, hi) | {hi})
        nbr += len(br)
        i, j = -1, len(br) - 1
        while j - i > 1:
            m = (i + j) // 2
            ev += 1
            if lz.evaluate(br[m], h)[0]:
                j = m
            else:
                i = m
        cs = br[j]
        cp = br[i] if i >= 0 else lo
        readp = lz.evaluate(cp, h)[1]
        ev += 1
        miss = ratios(lz, readp, cp, cs, closed=False)
        if not miss:
            return cs, cp, nbr, rounds, ev
        ts |= readp
        lo, hi = cp, cs


def failing(lz, h, c, show=6):
    """(7) at c: the failing crossings [(d, T(d), X)], each printed with
    its candidates (a, Y*(a)), the first `show` of them."""
    ok, read, Y, dams, bad, free = lz.evaluate(c, h, every=True)
    print(f"    at {c} = {float(c):.7f}: {len(bad)} failing of {len(lz.CR)}"
          f", {len(Y)} primes read, {len(read)} caps, {free} crossings "
          f"passed by a prime past H")
    for d, Td, X in bad:
        cand = [(a, Y[a]) for a in range(d + 1, X + 1, d) if lz.isp(a)]
        print(f"      crossing {d}: T(d) {Td}, X {X}, {len(cand)} candidates"
              f" {cand[:show]}{' ...' if len(cand) > show else ''}")
        for a, _ in cand[:show]:
            top = max(dams[a], key=lambda x: x[2])
            print(f"        {a}: its largest need at s = {top[0]}, Z "
                  f"{top[1]}, need {top[2]} by b = {top[3]}, Y*(b) {top[4]}")
    return [d for d, _, _ in bad]


def search(lz, h, lo, hi, label, use="both"):
    t1 = time.time()
    e0 = len(lz.FT)
    cs, cp, nbr, rounds, ev = exact(lz, h, lo, hi, use=use)
    print(f"  {label}: c_low = {cs} = {float(cs):.7f}, c' = {cp}, {nbr} "
          f"points, {rounds} rounds, {ev} evaluations, {len(lz.FT) - e0} "
          f"Engines built ({time.time() - t1:.1f} s)")
    return cs, cp, rounds


def free_arm(lz, h, c):
    """(7) at c: the crossing 2's candidates, and the crossings a^f read
    inside Y* that are passed by a prime past h."""
    _, _, Y, dams, _, _ = lz.evaluate(c, h, every=True)
    X = (c.numerator * lz.F(4)) // c.denominator
    cand = [(a, Y[a]) for a in range(3, X + 1, 2) if lz.isp(a)]
    print(f"  the crossing 2 at {c}: X(2) = {X}, candidates {cand}")
    rows = [(a, s, b) for a, rec in dams.items() for s, _, need, b, _ in rec
            if need == 0 and b > h]
    reads = sum(len(rec) for rec in dams.values())
    print(f"  inside Y*: {len(rows)} of {reads} crossings a^f read passed by "
          f"a prime past H, at {len({a for a, _, _ in rows})} of {len(Y)} "
          f"primes; the least such a^f {min((s for _, s, _ in rows), default=None)}")
    five = [(s, b) for a, s, b in rows if a == 5]
    print(f"  5's own: {five}")


def main():
    t0 = time.time()
    top = max(H, 10 ** 6)
    lz = Lazy(top)
    print(f"spf to {top + 1000}, sieve to {SIEVE} ({time.time() - t0:.1f} s)")

    print("CONTROLS")
    t1 = time.time()
    F4, _ = crosswide.f_feeder(10 ** 4)
    w4 = single.World(F4, 10 ** 4)
    w5 = single.World({10 ** 5: SIEVE // 7}, 10 ** 5)
    cr4, cr5 = lz.crossings(10 ** 4), lz.crossings(10 ** 5)
    same = all([(d, Td) for d, _, _, Td in w.crossings()] == cr
               for w, cr in ((w4, cr4), (w5, cr5)))
    del w5
    check("C1 (1)'s lists equal World.crossings at B = 10^4 and 10^5, 309 "
          "and 1922 long", same and (len(cr4), len(cr5)) == (309, 1922),
          f"{len(cr4)}, {len(cr5)} ({time.time() - t1:.1f} s)")
    lz.CR = cr4
    t1 = time.time()
    agree, n = True, 0
    for c in (LO4, Fr(33, 10)):
        Yw = w4.ystar(c)
        Yl = lz.evaluate(c, 10 ** 4, every=True)[2]
        agree = agree and all(min(Yl[a], w4.W + 1) == Yw[a] for a in Yl)
        n += len(Yl)
    del w4
    check("C6 at H = 10^4 the lazy Y* equals World.ystar (capped) at every "
          "prime read, c = 739/233 and 3.3", agree,
          f"{n} (c, a) pairs ({time.time() - t1:.1f} s)")
    cs, cp, _ = search(lz, 10 ** 4, Fr(1), Fr(7), "H = B = 10^4")
    bad = failing(lz, 10 ** 4, cp)
    check("C2 = P2 at H = B = 10^4 c_low = 739/233, failing just below at "
          "the crossing 41 alone", cs == LO4 and bad == [41], f"{bad}")
    cs, _, rounds = search(lz, 10 ** 4, Fr(1), Fr(7), "C5, from F(4) alone",
                           use="four")
    check("C5 the list from the crossing 2's cap alone: the check catches a "
          "missing ratio and c_low lands at 739/233",
          rounds > 1 and cs == LO4, f"{rounds} rounds")
    lz.CR = cr5
    f1 = lz.evaluate(Fr(1), 10 ** 5)[0]
    p7 = lz.evaluate(Fr(7), 10 ** 5)[0]
    check("C4 at H = 10^5 the test fails at c = 1 and passes at 7",
          not f1 and p7)
    cs, cp, _ = search(lz, 10 ** 5, Fr(1), Fr(7), "H = B = 10^5")
    bad = failing(lz, 10 ** 5, cp)
    check("C3 = P3 at H = B = 10^5 c_low = 11251/3373, failing just below "
          "at the crossing 2 alone", cs == LO5 and bad == [2], f"{bad}")
    if not all(ok for _, ok in CHECKS):
        print("CONTROL FAILED: nothing below is read.")
        raise SystemExit(1)

    lz.CR = lz.crossings(B)
    print(f"THE {len(lz.CR)} CROSSINGS WITH T(d) <= {B} AT H = {H}")
    t1 = time.time()
    e0 = len(lz.FT)
    at = lz.evaluate(LO5, H)
    print(f"  at 11251/3373: {'PASSES' if at[0] else 'FAILS'}, "
          f"{len(at[2])} primes read, {len(lz.FT) - e0} Engines built "
          f"({time.time() - t1:.1f} s)")
    if at[0]:
        cs = LO5
    else:
        print(f"  first failing crossing {at[4][0]}")
        if not lz.evaluate(Fr(7), H)[0]:
            print("  the test fails at c = 7 too; (6) stops")
            raise SystemExit(1)
        cs, cp, _ = search(lz, H, LO5, Fr(7), f"H = {H}")
        failing(lz, H, cp)
    failing(lz, H, cs)
    free_arm(lz, H, cs)
    built = len(lz.FT) - e0
    print(f"  Engines built for this horizon {built}, {len(lz.FT)} in all "
          f"({time.time() - t1:.1f} s)")
    print(f"  P4 (guess) the test fails at 11251/3373: "
          f"{'HOLDS' if cs > LO5 else 'FAILS'}")
    print(f"  P5 (guess) fewer than 5000 Engines for the horizon: "
          f"{'HOLDS' if built < 5000 else 'FAILS'}")
    hs = f"10^{len(str(H)) - 1}" if str(H).strip("0") == "1" else str(H)
    print("THE PRINT")
    if cs > LO5:
        print(f"  THE CROSSINGS TO 10^5 LIFT c0 TO {cs} AT {hs}")
    else:
        print(f"  KILL: THE CROSSINGS TO 10^5 LEAVE c0's BOUND AT 11251/3373 "
              f"AT {hs}")

    print("THE FRONTIER (8)")
    t1 = time.time()
    cr6 = lz.crossings(10 ** 6)
    new = [(d, Td) for d, Td in cr6 if Td > 10 ** 5]
    caps = sorted({pstar.next_pp(max(Td, 4), lz.spf) for _, Td in new}
                  - set(lz.FT))
    print(f"  {len(cr6)} crossings with T(d) <= 10^6, {len(new)} of them "
          f"past 10^5; {len(caps)} caps T+ not yet built "
          f"({time.time() - t1:.1f} s)")
    print(f"  P6 (guess) 8,000 to 14,000 past 10^5: "
          f"{'HOLDS' if 8000 <= len(new) <= 14000 else 'FAILS'}")
    if caps:
        pick = [caps[i * (len(caps) - 1) // 9] for i in range(10)]
        times = []
        for T in pick:
            t1 = time.time()
            pstar.Engine(T, lz.spf, lz.primes).feeder_bound()
            times.append(time.time() - t1)
        print("  ten Engines timed: " + ", ".join(
            f"T {T} {s:.2f} s" for T, s in zip(pick, times)))
        print(f"  mean {sum(times) / len(times):.2f} s, so about "
              f"{sum(times) / len(times) * len(caps) / 60:.0f} minutes for "
              f"the {len(caps)}")
    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()