pstar.py
Checks: IMAGE: The c-greedy walk is dammed theorem.
Cited on: IMAGE.
To run it: save pstar.py, with bridges.py, cgreedy.py, choosy.py, dambar.py, envelope.py, fwide.py, growth.py, slack.py beside it; standard library only.
"""pstar.py -- the dam's bar F as one round of striking, the feeder
formula (6), and how fast F grows past 10^4: like t ln t, or t ln^2 t?
QUESTION. F(t) = max{P+ : tau(P) <= t} is the least nondecreasing slack
every escape from the void pays; fwide.py fixes it on 4..10^4, where
F(10^4)/10^4 = 15.3757, one point. Whether F(t)/t grows like ln t or
like ln^2 t decides, heuristically, whether the single run's crossings
alone make c0 infinite. Read F at T from 10^4 to past 10^6.
DESIGN.
(1) F(t) = P*(T+(t)). T+(t) is the least prime power past t, and P*(T)
the least prime P at which cgreedy.feasible(T, P) holds. Feasibility is
monotone in P (a prime of [T, P] only enlarges N, and every bar is
monotone in N) and downward in T (a state with every price at least T
has every price at least T' < T). So tau(P) > t, tau a prime power,
exactly when T+(t) is feasible at P, and the largest P with
tau(P) <= t is the prime before P*(T+(t)).
(2) THE ENGINE. In cgreedy.feasible the primes of [T, P] enter only
through vhi[r] = max v_r(p - 1) over them, r < T, and every bar reads
a valuation only up to k_r - 1 (k_r least with r^(k_r) >= T): the
seated bar asks v_l(N) >= k_l - 1, the unseated v_l(N) >= k_l - 2 and
(l - 1) | N, where v_r(l - 1) <= k_r - 1 since r^v <= l - 1 < T, and
2's bar asks v_2(N) >= max(1, k_2 - 3). So capped at k_r - 1, vhi[r] =
max{j : L_T(r^j) <= P}, L_T(s) the least prime = 1 mod s at or past T.
P*(T) is then one of the L_T values (feasibility changes only where P
crosses one; P* > T at a prime power T that is not prime, since
tau(P) <= P from P = 5 on), found by a binary search over their sorted
distinct values, each step one striking pass. The striking pass keeps,
per prime r, the count of members l of X with v_r(l - 1) = e for each
e, so N's valuation is read off the counts plus vhi, and striking l
re-reads only the members whose own valuation dropped (a worklist).
L_T(s) is found over the candidates m*s + 1 (m*s even) from T on, by
a sieve of [T, 30 T] and Miller-Rabin (growth.is_prime_mr) past it.
(3) THE PAIR BOUND (proved). Let l = 2r + 1 with l and r prime, r odd,
and (T + 1)/2 < l < T. Every feasible seated set at (T, P) has
L_T(l) <= P or L_T(r) <= P. Seated, l's bar asks v_l(N) >= k_l - 1 >= 1
(k_l >= 2 as l < T); a seated prime p with l | p - 1 is = 1 mod 2l, so
p >= 2l + 1 > T, outside the seated primes below T; so some prime of
[T, P] is = 1 mod l. Unseated, l's bar asks 2r | N; a seated prime p
with r | p - 1, p < T and p != l, is 2rm + 1 with m >= 2, so
p >= 4r + 1 = 2l - 1 > T, impossible; so some prime of [T, P] is
= 1 mod r. Hence P*(T) >= PB(T) = max over such pairs of
min(L_T(l), L_T(r)).
(4) THE HAND COUNT (a TRANSPLANT: least primes in residue classes as
independent exponential waits at Cramer's density, re-derived here).
The candidates = 1 mod 2q past T, q an odd prime, are spaced 2q and each
is prime with probability about 2q/(phi(2q) ln x) ~ 2/ln x, so
L_T(q) - T ~ q ln x * Exp(1); for a pair the two waits have rates
1/(l ln x) and 2/(l ln x), so min(L_T(l), L_T(r)) - T ~ (l ln x/3)
Exp(1); the pairs, r in (T/4, T/2), number about 1.3203 * integral of
dr/(ln r ln 2r) (Hardy-Littlewood's twin-type count), and the median of
their maximum, solved numerically with x = P* itself, gives P*/T =
13.8, 24.5, 38.7 at T = 10^4, 10^5, 10^6 (45, 280, 1908 pairs): F/(t ln t)
1.49, 2.12, 2.80, F/(t ln^2 t) 0.162, 0.185, 0.203, rising slowly
toward 1/3 as ln(pairs) catches ln T. A truth like C t ln t holds
F/(t ln t) flat (falling, by the 1/ln t of P* > T). The count says R,
the ratio of F/(t ln t)'s median at 10^6 to that at 10^4, is 1.87 (2.80
over 1.49 before rounding; P3 below used 1.88); t ln t says 1.0, a pure
C t ln^2 t 1.50. One draw's spread: the maximum of the waits is Gumbel,
its relative spread about 1.28/(ln(pairs) + 0.58), 30% at 10^4 and 16%
at 10^6, so twelve independent draws would carry log R to about 0.12,
against log 1.87 = 0.63.
(5) THE WINDOWS. For d = 4, 5, 6, the twelve t = round(10^(d + i/24)),
i = -6..5, each read at T+(t), F(t)/t taken at the drawn t. Draws a
factor 10^(1/24) apart share most of their least primes, each wait
spanning about T ln x/3, far past the spacing, and P* is nondecreasing
in T: the twelve are not independent, and the window's log se, computed
as if they were, understates the spread.
(6) THE FEEDER FORMULA (proved after the first run). For s = r^j < T let
C(s), the feeders of s, be the primes below T that are = 1 mod s, each
larger than s, and call an odd p < T UNFED when C(p^(k_p - 1)) is
empty. Put nu(s) = min(L_T(s), min over p in C(s) of nu(p^(k_p - 1))),
well founded as p^(k_p - 1) >= p > s. Striking at P keeps exactly the p with
P >= nu(p^(k_p - 1)): downward from the largest, p's seated bar is met
by a prime of [T, P] or by a kept member of C(p^(k_p - 1)), all larger
than p, and a struck p fails its bar at the final N, N only shrinking.
So v_r(N) >= j exactly when P >= nu(r^j). An odd l, kept or not, then
meets its unseated bar exactly when P >= nu(s) at each prime power s
exactly dividing l - 1 and at s = l^(k_l - 2) (k_l >= 3), since each
such nu(s) is at most nu(l^(k_l - 1)) (l in C(s), or C(l^(k_l - 1)) inside
C(l^(k_l - 2))); 2's bar asks nu(2^max(1, k_2 - 3)) <= P. nu is monotone
along divisibility, so P*(T) is the largest nu(s) over the s with C(s)
nonempty and the forced s (2's, and l^(k_l - 2)) with C(s) empty. And
nu(s) <= nu(p^(k_p - 1)) at every p in C(s), and C(p^(k_p - 1)) is
nonempty when p is not UNFED: following such members, each larger, ends
at a nonempty C(s) of UNFED members, where nu(s) = min(L_T(s), min over p in
C(s) of L_T(p^(k_p - 1))) outright. So P*(T) = CB(T), the largest of
L_T(s) over the forced s with C(s) empty and of that minimum over the
s whose C(s) is nonempty and UNFED throughout: F(t) = CB(T+(t)), one
round of striking. The pair bound is the case C(r) = {2r + 1}.
PREDICTIONS (frozen before the code).
P1 (forced) P*(T+(t)) equals fwide.py's F at every t in 4..10^4.
P2 (proved) PB(T) <= P*(T) at every T read.
P3 (the count, 70%) R >= 1.35, the geometric midpoint of 1.0 and 1.88.
P4 (the count, 50%) each window's median F/t within 0.7..1.5 times the
count's 13.8, 24.5, 38.7.
P5 (a guess, 55%) PB/P* has median >= 0.8 in every window.
P6 (a guess, 50%) at the L_T value below P*, a failing unseated bar
belongs to a pair's l at half the draws of the 10^6 window or more.
THE KILL, frozen. The script prints R, the median F/(t ln t) over the
d = 6 window over that over the d = 4 window, with its spread, and then
"R >= 1.35: F/t RISES PAST ln t" or "KILL: R < 1.35, F/t NO FASTER THAN
ln t". KILL of "F grows like t ln^2 t" (the count's prediction) at
R < 1.35. The bar sits on R, which the 1/ln t of P* > T deflates under
a t ln t truth, so a confound pushes toward the kill, not past it.
CONTROLS, read first. C1 is P1. C2 L_T(s) equals a brute least prime
= 1 mod s at or past T, by trial division, at every s the engine reads
for T = 1009 and T = 1024, and at T = 1009 with the sieve cut to
[T, 2T], so that Miller-Rabin decides some (added after the first run).
C3 the kill's reading on planted windows:
F = 2 t ln t at every draw prints the KILL; F = t ln^2 t/5 prints
R >= 1.35. C4 a planted break of the pair bound, the larger of each
pair's least primes in place of the smaller, is flagged by P2's test at
some draw. C5 (added after the first run) is (6): CB(T) = P*(T) at
every T of C1. C6 CB over one-member C(s) only falls below P* at
some T of C1, so C5 can fail.
FINDINGS (from the printed run; the controls passed first). C1 holds at
1278 T (0 differ), C2 at its three readings (108 values by
Miller-Rabin), C3 prints R 1.000 and 1.503, C4 flags its broken bound
at all 36 draws, C6 at 54 T, 2243..2663. P2 holds at every T read, and
(6) CB = P* at every T of
C1 and every draw. THE PRINT: R = 1.8329, log R 0.606 against the
count's log 1.87 = 0.626, log se 0.107: R >= 1.35, F/t RISES PAST ln t,
the kill silent; P3 holds. Window medians of F/t 12.74, 19.47, 34.61
(count 13.8, 24.5, 38.7: P4 holds, ratios 0.923, 0.795, 0.894); of
F/(t ln t) 1.3526, 1.6773, 2.4791; of F/(t ln^2 t) 0.1441, 0.1457,
0.1764. The draws are not twelve independent ones: P* takes 8, 8 and 7
distinct values in the three windows (19424281 and 58168417 three times
each at d = 6), so the printed log se understates R's spread, and R is
read as a value, its bar at 1.35, not as a significance. PB = P* at 12,
11 and 9 of 12 draws (P5 holds, PB/P* least
0.789); where PB < P*, CB's witness is (157, 2) at t = 10^5 (98597 =
4 * 157^2 + 1, the one prime below T = 1 mod 157^2) and (161869, 1) at
the three draws from 681292 (647477 = 4 * 161869 + 1); a pair's l fails
below P* at 9 of 12 at d = 6 (P6 holds). Below 10^4, P* = PB at 1139 of
1278 T, the last differing T 2879; one-member C(s) fall short at
54 T, 2243..2663 (C6), and at T = 2243 P* = 15667 is L_T(373), C(373) =
{1493, 2239}, both UNFED. The largest F/t over the draws is
43.6202 at t = 1333521, so rho >= 43.62 at the prime before 58168417
(tau there at most t), against 17.998, the table's record to 153757;
PB = P* there, so the floor rests on the pair bound alone. So F(t)/t
rises faster than ln t: in log, R sits 0.20 from a pure C ln^2 t law's
1.50 (C3) and 0.02 from the pair count's 1.87, but the printed log se,
0.107, understates the draws' dependence, so three windows do not
separate the two; C = F/(t ln^2 t), 0.1441, 0.1457 and 0.1764, rises
by that same excess (0.1764/0.1441 = 1.224 = 1.8329/1.503).
RUN RECORD. Rehearsal --top 4: 4/4, 40.5 s, 65.6 MB. First full run:
5/5, 61.3 s, 194.6 MB. With the unfed-prime reading: 7/7, 63.0 s,
193.7 MB. With the feeder formula as checks: 7/7, 63.4 s, 204.2 MB.
With the later controls (C2's Miller-Rabin reading, C4 on the
computation, C6) and r odd in the pairs: 8/8, 63.8 s, 201.7 MB, every
figure above unchanged.
python pstar.py [--top D] (D the largest window, default 6; a
rehearsal runs --top 4)
"""
import os
import sys
import time
from array import array
from math import log, sqrt
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import growth # noqa: E402
CHECKS = []
SEG = 30 # the sieve reads [T, SEG * T]
COUNT = {4: 13.8, 5: 24.5, 6: 38.7} # the hand count's P*/T, (4)
def check(name, ok, detail=""):
CHECKS.append((name, bool(ok)))
print(f" [{'PASS' if ok else 'FAIL'}] {name}"
+ (f" -- {detail}" if detail else ""))
def spf_to(n):
"""Smallest prime factor of every m <= n, 0 at the primes."""
spf = array("i", [0]) * (n + 1)
for p in range(int(n ** 0.5), 1, -1):
cnt = len(range(p * p, n + 1, p))
spf[p * p::p] = array("i", [p]) * cnt
return spf
def factor(m, spf):
out = {}
while m > 1:
q = spf[m] or m
out[q] = out.get(q, 0) + 1
m //= q
return out
def next_pp(t, spf):
"""T+(t), the least prime power past t."""
n = t + 1
while True:
q = spf[n] or n
m = n
while m % q == 0:
m //= q
if m == 1:
return n
n += 1
class Engine:
"""Everything P*(T) reads at one T."""
def __init__(self, T, spf, primes, seg=SEG):
self.T = T
self.low = [l for l in primes if l < T]
self.k = {}
for r in self.low:
k, x = 1, r
while x < T:
k, x = k + 1, x * r
self.k[r] = k
hi = seg * T
seg = bytearray([1]) * (hi - T + 1)
for p in primes:
if p * p > hi:
break
s = max(p * p, -(-T // p) * p)
seg[s - T::p] = bytes(len(range(s, hi + 1, p)))
self.seg, self.hi = seg, hi
self.L = {} # (r, j) -> L_T(r^j)
for r in self.low:
for j in range(1, self.k[r]):
self.L[(r, j)] = self.least(r ** j)
self.odd = [l for l in self.low if l > 2]
self.fac = {l: [(r, min(e, self.k[r] - 1))
for r, e in factor(l - 1, spf).items()]
for l in self.odd}
self.vals = sorted(set(self.L.values()))
def isp(self, n):
if n <= self.hi:
return self.seg[n - self.T] == 1
return growth.is_prime_mr(n)
def least(self, s):
"""L_T(s): the least prime = 1 mod s at or past T."""
step = s if s % 2 == 0 else 2 * s
n = -(-(self.T - 1) // step) * step + 1
while not self.isp(n):
n += step
return n
def vhi(self, P):
out = {}
for (r, j), v in self.L.items():
if v <= P and j > out.get(r, 0):
out[r] = j
return out
def strike(self, P):
"""The striking pass at P: (X, the failing bars), each bar read
as cgreedy.feasible reads it, valuations capped at k_r - 1."""
k, fac = self.k, self.fac
vh = self.vhi(P)
cnt = {}
for l in self.odd:
for r, e in fac[l]:
cnt.setdefault(r, [0] * k[r])[e] += 1
def vN(r):
top = vh.get(r, 0)
c = cnt.get(r)
if c:
for e in range(len(c) - 1, top, -1):
if c[e]:
return e
return top
X = set(self.odd)
work = list(self.odd)
while work:
l = work.pop()
if l not in X or vN(l) >= k[l] - 1:
continue
X.discard(l)
for r, e in fac[l]:
c = cnt[r]
c[e] -= 1
if c[e] == 0 and r in X:
work.append(r)
fail = [l for l in self.odd if l not in X and
(any(vN(r) < e for r, e in fac[l]) or vN(l) < k[l] - 2)]
if vN(2) < max(1, k[2] - 3):
fail.append(2)
return X, fail
def feasible(self, P):
return not self.strike(P)[1]
def pstar(self):
"""(P*(T), the L_T value below it, None at the least)."""
v = self.vals
assert self.feasible(v[-1])
lo, hi = -1, len(v) - 1 # infeasible at v[lo], feasible at v[hi]
while hi - lo > 1:
mid = (lo + hi) // 2
if self.feasible(v[mid]):
hi = mid
else:
lo = mid
return v[hi], (v[lo] if lo >= 0 else None)
def pairs(self):
ps = set(self.low)
return [(l, (l - 1) // 2) for l in self.low
if l > 5 and 2 * l > self.T + 1 and (l - 1) // 2 in ps]
def pair_bound(self, pick=min):
"""(PB(T), its pair), (0, None) with no pair; pick=max is C4's
broken bound."""
best = (0, None)
for l, r in self.pairs():
m = pick(self.L[(l, 1)], self.L[(r, 1)])
if m > best[0]:
best = (m, (l, r))
return best
def feeders(self):
"""{(r, j): C(r^j)}, the primes below T = 1 mod r^j, nonempty."""
C = {}
for p in self.odd:
for r, e in self.fac[p]:
for j in range(1, e + 1):
C.setdefault((r, j), []).append(p)
return C
def feeder_bound(self, single=False):
"""(CB(T), its witness), added after the first run. C(s) is the
primes below T = 1 mod s, s = r^j < T; p is UNFED when C(p^(k_p
- 1)) is empty, so only a prime past T meets its seated bar.
Proved as (3) is: every bar needing s | N is met only by a
member of C(s) seated or by L_T(s) <= P. So P*(T) is at least
L_T(s) when C(s) is empty and s | N is forced: s = 2^(max(1,
k_2 - 3)) (2's bar), or s = l^(k_l - 2) at odd l, k_l >= 3
(l's bar, seated or not);
min(L_T(s), min over p in C(s) of L_T(p^(k_p - 1))) when C(s)
is nonempty and every member is UNFED (an unseated member
asks s | N, a seated one its own bar)."""
k, L = self.k, self.L
C = self.feeders()
best = (0, None)
forced = [(2, max(1, k[2] - 3))]
forced += [(l, k[l] - 2) for l in self.odd if k[l] >= 3]
for s in forced:
if s not in C and L[s] > best[0]:
best = (L[s], s)
for s, ps in C.items():
if ((len(ps) == 1 or not single) and
all((p, k[p] - 1) not in C for p in ps)):
m = min([L[s]] + [L[(p, k[p] - 1)] for p in ps])
if m > best[0]:
best = (m, s)
return best
def brute_least(s, T):
n = T
while not ((n - 1) % s == 0 and growth.is_prime(n)):
n += 1
return n
def read(T, spf, primes):
t0 = time.time()
E = Engine(T, spf, primes)
P, below = E.pstar()
PB, pair = E.pair_bound()
fail = E.strike(below)[1] if below is not None else []
pl = {l for l, _ in E.pairs()}
CB, wit = E.feeder_bound()
return {"T": T, "P": P, "PB": PB, "PBx": E.pair_bound(max)[0],
"pair": pair, "nfail": len(fail),
"pairfail": any(l in pl for l in fail), "CB": CB,
"wit": wit, "s": time.time() - t0}
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 lnse(xs):
"""The standard error of a median's log: 1.2533 sd / sqrt(n)."""
ys = [log(x) for x in xs]
m = sum(ys) / len(ys)
sd = sqrt(sum((y - m) ** 2 for y in ys) / (len(ys) - 1))
return 1.2533 * sd / sqrt(len(ys))
def verdict(win):
"""R, its log se and the line, from {d: [F/(t ln t)]}: the top d's
median over d = 4's."""
top = max(win)
R = median(win[top]) / median(win[4])
se = sqrt(lnse(win[4]) ** 2 + lnse(win[top]) ** 2)
line = ("R >= 1.35: F/t RISES PAST ln t" if R >= 1.35 else
"KILL: R < 1.35, F/t NO FASTER THAN ln t")
return R, se, line
def p2_bad(draws):
return [x["t"] for x in draws if x["PB"] > x["P"]]
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 = 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")
bad = []
past = 0
for T, sg in ((1009, SEG), (1024, SEG), (1009, 2)):
E = Engine(T, spf, primes, sg)
bad += [(T, r, j) for (r, j), v in E.L.items()
if v != brute_least(r ** j, T)]
past += sg == 2 and sum(v > E.hi for v in E.L.values())
check("C2 L_T(s) equals the brute least prime = 1 mod s at or past T "
"at every s read, T = 1009 and 1024, and T = 1009 with the "
"sieve cut to [T, 2T]", not bad and past > 0,
f"{bad[:3]}; {past} values by Miller-Rabin")
import cgreedy # noqa: E402
import dambar # noqa: E402
import fwide # noqa: E402
cgreedy.tau_table()
fwide.widen(fwide.PMAX)
tab = fwide.tau_wide(fwide.PMAX, fwide.T)
F = dambar.f_table(tab)
t1 = time.time()
memo, pb, lb, lb1, off = {}, {}, {}, {}, []
for t in range(4, fwide.T + 1):
T = next_pp(t, spf)
if T not in memo:
E = Engine(T, spf, primes)
memo[T], pb[T] = E.pstar()[0], E.pair_bound()[0]
lb[T] = E.feeder_bound()[0]
lb1[T] = E.feeder_bound(single=True)
if memo[T] != F[t]:
off.append((t, memo[T], F[t]))
check(f"C1 = P1 P*(T+(t)) equals fwide.py's F at every t in "
f"4..{fwide.T}", not off,
f"{len(memo)} T, {len(off)} differ {off[:3]} "
f"({time.time() - t1:.1f} s)")
# added after the rehearsal: where the pair bound is P* itself
eq = [T for T in memo if pb[T] == memo[T]]
ne = sorted(T for T in memo if pb[T] != memo[T])
print(f" reading: P* = PB at {len(eq)} of the {len(memo)} T; the "
f"largest T where they differ {ne[-3:]}, the least where they "
f"agree {min(eq) if eq else None}")
# added after the first run: the feeder formula (6), proved
ne = sorted(T for T in memo if lb[T] != memo[T])
check("C5 = (6) CB(T) = P*(T) at every T of C1", not ne, f"{ne[:3]}")
ne = sorted(T for T in memo if lb1[T][0] != memo[T])
check("C6 CB over one-member C(s) only falls below P* at some T of "
"C1 (C5 can fail)", ne and all(lb1[T][0] < memo[T] for T in ne),
f"{len(ne)} T, {ne[0] if ne else None}..{ne[-1] if ne else None}")
if ne:
T = ne[0]
E = Engine(T, spf, primes)
C = E.feeders()
wit = [(s, ps) for s, ps in C.items()
if len(ps) > 1 and E.L[s] == memo[T] and
all((p, E.k[p] - 1) not in C for p in ps)]
print(f" reading: at T = {T}, P* = {memo[T]} is L_T(s) at (s, its "
f"feeders) {wit}, all UNFED")
if top > 4:
a = verdict({d: [2.0] * len(ts[d]) for d in ts})
b = verdict({d: [log(t) / 5 for t in ts[d]] for d in ts})
check("C3 the kill's reading: F = 2 t ln t prints the KILL, "
"F = t ln^2 t/5 prints R >= 1.35",
a[2].startswith("KILL") and b[2].startswith("R >= 1.35"),
f"R {a[0]:.3f}, {b[0]:.3f}")
if not all(ok for _, ok in CHECKS):
print("CONTROL FAILED: nothing below is read.")
raise SystemExit(1)
print("THE DRAWS")
draws = []
for d in sorted(ts):
for t in ts[d]:
x = read(next_pp(t, spf), spf, primes)
x["t"], x["d"] = t, d
draws.append(x)
ft = x["P"] / t
print(f" d {d} t {t} T {x['T']}: P* {x['P']}, F/t {ft:.3f}, "
f"F/(t ln t) {ft / log(t):.4f}, F/(t ln^2 t) "
f"{ft / log(t) ** 2:.4f}, PB/P* {x['PB'] / x['P']:.3f} at "
f"{x['pair']}; below P* {x['nfail']} bars fail, a pair's l "
f"among them {x['pairfail']} ({x['s']:.1f} s)", flush=True)
plant = [dict(x, PB=x["PBx"]) for x in draws]
check("C4 P2's test flags the broken pair bound (the larger least "
"prime of each pair) at some draw", p2_bad(plant) != [],
f"{len(p2_bad(plant))} draws")
bad = p2_bad(draws)
check("P2 PB(T) <= P*(T) at every T read", not bad, f"{bad[:3]}")
# added after the first run: the feeder formula at the draws
ne = [(x["t"], x["wit"]) for x in draws if x["CB"] != x["P"]]
check("(6) CB(T) = P*(T) at every draw", not ne, f"{ne[:3]}")
print(f" reading: CB's witness where PB < P*: "
f"{[(x['t'], x['wit']) for x in draws if x['PB'] < x['P']]}")
tm = max(draws, key=lambda x: x["P"] / x["t"])
print(f" reading: the largest F/t over the draws, {tm['P'] / tm['t']:.4f}"
f" at t = {tm['t']}, is a floor for rho at the prime before "
f"{tm['P']}")
print("THE WINDOWS")
win = {}
for d in sorted(ts):
ws = [x for x in draws if x["d"] == d]
ft = [x["P"] / x["t"] for x in ws]
f1 = [f / log(x["t"]) for f, x in zip(ft, ws)]
f2 = [f / log(x["t"]) ** 2 for f, x in zip(ft, ws)]
sh = [x["PB"] / x["P"] for x in ws]
win[d] = f1
print(f" d {d}: {len({x['P'] for x in ws})} distinct P* in "
f"{len(ws)} draws")
print(f" d {d}: F/t median {median(ft):.2f} (range {min(ft):.2f}"
f"..{max(ft):.2f}; count {COUNT.get(d)}), F/(t ln t) "
f"{median(f1):.4f} (log se {lnse(f1):.3f}), F/(t ln^2 t) "
f"{median(f2):.4f}; PB/P* median {median(sh):.3f}, least "
f"{min(sh):.3f}; P* = PB at "
f"{sum(x['P'] == x['PB'] for x in ws)}/{len(ws)}; a pair's l "
f"failing below P* at {sum(x['pairfail'] for x in ws)}"
f"/{len(ws)}")
if d in COUNT:
r = median(ft) / COUNT[d]
print(f" P4 at d {d}: median/count {r:.3f}, "
f"{'HOLDS' if 0.7 <= r <= 1.5 else 'FAILS'}")
print(f" P5 at d {d}: {'HOLDS' if median(sh) >= 0.8 else 'FAILS'}")
if top >= 6:
n = sum(x["pairfail"] for x in draws if x["d"] == 6)
print(f" P6: {'HOLDS' if 2 * n >= len(ts[6]) else 'FAILS'} "
f"({n}/{len(ts[6])})")
print("THE PRINT")
if top > 4:
R, se, line = verdict(win)
print(f" R = {R:.4f} (log R {log(R):.3f}, log se {se:.3f}; the "
f"count's log 1.87 = {log(1.87):.3f})")
print(f" P3: {'HOLDS' if R >= 1.35 else 'FAILS'}")
print(" " + line)
else:
print(" rehearsal: one window, no R")
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()