"""danger.py -- at the takeover theorem's own e(r), does the exact centred test put every single-odd-part shape above F_e's circle mean? QUESTION. J* = (4e - 1) l - 1, l = r - 2, e = floor(sqrt(l/2)) at odd r, floor(sqrt(l)/2) at even r, e = 3..11 (odd r = 21..289, even r = 38..576). The shapes k=0 2^(r-g) g, g odd >= 5: is the exact centred test, best over p, above F_e's mean at every one? BY HAND (re-derived before the code). al = (r - g + 1 + q)/2, be = (J + p)/2, q = (r - g) mod 2, a = l/2, b = (J + 1)/2: X/mean = B1 |E cos gu| / A, B1 = 2^(p+q+3-g)/|w|_1 B(al, be)/B(a, b), u drawn as cos^(2be-1) sin^(2al-1) on [0, pi/2], A = E|1 - 4eU|, U ~ Beta(a, b). The hand form H replaces E cos gu by cos(g ubar) exp(-g^2 Var u/2), ubar and Var u read off the shape's own Beta(al, be) with Jensen's correction (partfive.py's C5h at general g). The class read off u0 = arcsin(1/(2 sqrt e)) dies past the first crossing: the shape's own mean moves g u by O(1) once g ~ sqrt(l) (e = 5, r = 53, g = 21: g u0 = 4.736, g ubar = 3.834). Scale: X/mean ~ 0.63 l |cos g ubar| near the first crossing, the per-rank least O(1) with no e-trend; p rescues at most ~0.33 g/sqrt(8e). PREDICTIONS (frozen before the code). P1 H/exact in [0.95, 1.05] at every shape with g <= r/2 and |cos g ubar| >= 0.2, at its exact best p. P2 per rank the least exact X/mean sits at a part of the first or second crossing, in [1, 4] at most ranks, with no rising trend in e. P3 at e = 5 the binding shape is the part-7 shape, least in [1.2, 2.2]. P4 somewhere in e = 3..11 a shape falls below the mean. KILL. Printed per rank: the least exact X/mean over the family, its g and p, and H at the same (g, p). H/exact outside [0.9, 1.1] at a rank's binding shape fires: the hand form cannot carry a uniform margin. CONTROLS, read first. C1 at e = 3, odd r = 21..33: the part-5 shape's exact X/mean equals partfive's printed 2.6738 .. 2.7100 (4 digits). C2 at e = 4, odd r = 39..51: the part-7 shape reads exacttest's 9.8599 .. 14.0320. PLANT: H with cos((g - 2) ubar) leaves [0.9, 1.1] of exact at the part-5 shape at every r of C1. C3 psi1 equals exacttest.psi at r = 53, every g; at r = 53, 54, 101, 102, 289 the pipeline's least and below-list equal the full read's (best_x at every shape; the below-lists are empty at all five). The pipeline (changed after the first run, which read exact only the shapes whose hand best fell below 6 and each rank's three least by hand): every shape is read exact at its hand-best p; exacttest.best_x, a scan over p that stops 40 steps past its best (not the maximum: at r = 258, g = 251 the hand-best p reads 6.3 times its value), is read where that read is below the running least or the mean. Every shape then has some p whose exact X/mean is at least the printed least, and none skipped is below the mean; that the least equals the full read's is C3's, at its five ranks. "Least" is the least best_x read, a lower bound on the least maximum over p. FINDINGS (copied from the printed output). (1) The controls hold. C1: the part-5 shape reads 2.6738 .. 2.7100 at r = 21..33, H/exact 0.9981 .. 0.9991, the cos 3u plant 3.4751 .. 4.1417 (caught at every r). C2: 9.8599 (r = 39), 14.0320 (r = 51). C3: psi1 differs at 0 of 25 shapes; the pipeline's least equals the full read at r = 53 (11.1125, g = 7), 54 (2.3361, g = 5), 101 (38.5160, g = 5), 102 (5.5378, g = 7), 289 (124.0961, g = 5). (2) The KILL misses: at all 405 ranks H/exact at the binding shape lies in 0.9981 .. 1.0002. P1 at its frozen p, best_x's, is read only at the binding shapes (the KILL's H/exact); read at each shape's hand-best p it fails at 33 shapes, all at even r = 192..572 with g >= 0.46 r, H/exact 0.942 .. 1.067. (3) All 51165 shapes are read exact at their hand-best p, 405 (one a rank) by best_x: no shape falls below the mean at any rank (P4 fails). The least per (e, parity), at r and g: e odd r even r 3 2.6617 (25, 5) 2.1747 (38, 5) 4 7.8844 (35, 5) 7.8243 (66, 5) 5 8.5955 (73, 7) 4.1169 (144, 7) 6 25.8052 (75, 5) 27.3201 (146, 5) 7 38.5160 (101, 5) 41.1192 (198, 5) 8 32.1444 (163, 9) 57.5718 (258, 9) 9 71.4459 (165, 5) 76.8094 (326, 5) 10 91.6668 (203, 5) 98.6990 (402, 5) 11 114.3926 (245, 5) 123.2842 (486, 5) The binding part is 5 at 358 ranks, 7 at 33 (all of e = 5), 9 at 14 (of e = 8: odd r = 139..163 and r = 258). P2 fails: the least rises with e, with dips at e = 5 and 8. P3 fails: e = 5 binds at the part-7 shape as predicted, but at 4.1169 (even) and 8.5955 (odd), above 2.2. (4) Within e = 5 the part-7 shape falls in r: 11.1125 (r = 53) to 8.5955 (r = 73), 5.5378 (r = 102) to 4.1169 (r = 144). READS. The dips sit where the first-crossing part's g u0 is just above pi/2: 7u0 = 1.579 at e = 5 (0.008 above), 9u0 = 1.599 at e = 8. The shape's own downward shift, about g(g - 3)/(4 sqrt(e) l) (its mean alone, odd r, p = 0; Jensen adds a further g/(8 sqrt(e) l) downward, and p moves it further at even r), holds g ubar below the zero, and it shrinks as l grows, so the margin falls in r at fixed e (at e = 5 the mean alone would close it near l ~ 400, with Jensen near l ~ 450, argued, not printed; e = 5's window ends at l = 142). Away from a dip the binding shape is the part 5, where 5 ubar falls toward 0 as e grows: X/mean over l reads 0.47 (r = 245) to 0.43 (r = 289) at e = 11's odd ranks, 0.25 to 0.23 at its even ones. BY HAND's scale 0.63 l |cos g ubar| held the part 5's prefactor; B1 grows with g, and the part 7 stands above that scale at e = 5. A dip is an accident of e: one window of l moves the shift by O(e^-2.5), so the margin at e is set by how near g u0, read mod pi, sits to pi/2 + shift: a Diophantine condition on pi/(2 arcsin(1/(2 sqrt e))) against the odd integers, which the centred test alone cannot bound uniformly. The design variable is e itself: J_hi >= c r^(3/2) needs only some e of order sqrt(l) per rank, not the theorem's. RUN RECORD. 15.7 s, peak commit 29.0 MB. The full read (best_x at every shape) costs 20.4 s at r = 576 alone, about 25 min over the range (a separate timing run). """ import math import os import sys import time os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import circlemean # noqa: E402 import compbound # noqa: E402 import exacttest # noqa: E402 import takeph # noqa: E402 LN2 = math.log(2) def setup(r): e, J = takeph.theorem_e(r), takeph.theorem_J(r) l = r - 2 lm = circlemean.log_mean(r, J, e) logA = float(circlemean.mp.log(circlemean.mad( circlemean.mp.mpf(l) / 2, circlemean.mp.mpf(J + 1) / 2, 4 * e))) return e, l, J, lm, logA def lbeta(x, y): return math.lgamma(x) + math.lgamma(y) - math.lgamma(x + y) def hand(r, J, logA, g, p, shift=0): """log H and cos(g ubar) at (g, p); shift: the plant's g - 2.""" l = r - 2 q = (r - g) % 2 a, b = l / 2, (J + 1) / 2 al, be = (r - g + 1 + q) / 2, (J + p) / 2 lB1 = ((p + q + 3 - g) * LN2 - compbound.log_w1(p, q) + lbeta(al, be) - lbeta(a, b)) U = al / (al + be) VU = al * be / ((al + be) ** 2 * (al + be + 1)) upp = -(1 - 2 * U) / (4 * (U * (1 - U)) ** 1.5) ub = math.asin(math.sqrt(U)) + upp * VU / 2 Vu = VU / (4 * U * (1 - U)) c = math.cos((g - shift) * ub) if c == 0: return -math.inf, c return lB1 + math.log(abs(c)) - g * g * Vu / 2 - logA, c def hand_best(r, J, logA, g, patience=40): p = J % 2 best, since = (-math.inf, p), 0 while since < patience: v = hand(r, J, logA, g, p)[0] if v > best[0]: best, since = (v, p), 0 else: since += 1 p += 2 return best def exact(r, J, lm, g): v, p = exacttest.best_x(J, 0, (g,) + (2,) * (r - g)) return v - lm, p def psi1(r, g): """exacttest.psi of the shape k=0 2^(r-g) g: 2^(r-g) T_g, built without r - g list products (C3 reads it against psi).""" return [v << (r - g) for v in exacttest._cheb(exacttest._T, g)] def controls(): ok = True want1 = {21: 2.6738, 23: 2.6630, 25: 2.6617, 27: 2.6673, 29: 2.6780, 31: 2.6925, 33: 2.7100} for r, w in want1.items(): e, l, J, lm, logA = setup(r) v, p = exact(r, J, lm, 5) x = math.exp(v) h = math.exp(hand(r, J, logA, 5, p)[0] - v) pl = math.exp(hand(r, J, logA, 5, p, shift=2)[0] - v) good = e == 3 and abs(x - w) < 5e-5 and not 0.9 <= pl <= 1.1 ok &= good print("C1 r = %d e = %d: part 5 X/mean %.4f (want %.4f) p = %d," " H/exact %.4f, plant %.4f %s %s" % (r, e, x, w, p, h, pl, "caught" if not 0.9 <= pl <= 1.1 else "MISSED", "holds" if good else "FAILS"), flush=True) want2 = {39: 9.8599, 51: 14.0320} for r, w in want2.items(): e, l, J, lm, logA = setup(r) x = math.exp(exact(r, J, lm, 7)[0]) good = e == 4 and abs(x - w) < 5e-5 ok &= good print("C2 r = %d e = %d: part 7 X/mean %.4f (want %.4f) %s" % (r, e, x, w, "holds" if good else "FAILS"), flush=True) return ok def rank(r, full=False): """(e, (least exact X/mean, g, p, H/exact at it), shapes read at every p, shapes below the mean, P1 misses, shapes). Every shape is read exact at its hand-best p; best_x is read where that read is below the running least or the mean (every shape when full).""" e, l, J, lm, logA = setup(r) rows = [] for g in range(5, r + 1, 2): hp = hand_best(r, J, logA, g)[1] lo = exacttest.log_x(J, 0, r - g, 1, psi1(r, g), hp) - lm rows.append((lo, g, hp)) rows.sort() least, below, p1, n = (math.inf, 0, 0, 0.0), [], [], 0 for lo, g, hp in rows: hv, c = hand(r, J, logA, g, hp) ratio = math.exp(hv - lo) if g <= r / 2 and abs(c) >= 0.2 and not 0.95 <= ratio <= 1.05: p1.append((g, round(ratio, 3))) if not (full or lo < max(least[0], 0.0)): continue v, p = exact(r, J, lm, g) n += 1 if v < 0: below.append((g, p, round(math.exp(v), 4))) if v < least[0]: least = (v, g, p, math.exp(hand(r, J, logA, g, p)[0] - v)) return e, least, n, below, p1, len(rows) def main(): t0 = time.time() ok = controls() bad = sum(psi1(53, g) != exacttest.psi((g,) + (2,) * (53 - g)) for g in range(5, 54, 2)) ok &= bad == 0 print("C3 r = 53: psi1 against exacttest.psi, %d of 25 shapes differ %s" % (bad, "holds" if bad == 0 else "FAILS"), flush=True) for r in (53, 54, 101, 102, 289): a, b = rank(r), rank(r, full=True) good = (a[1][1] == b[1][1] and abs(a[1][0] - b[1][0]) < 1e-12 and a[3] == b[3]) ok &= good print("C3 r = %d: pipeline least %.4f at g = %d (%d of %d read at" " by best_x), full %.4f at g = %d %s" % (r, math.exp(a[1][0]), a[1][1], a[2], a[5], math.exp(b[1][0]), b[1][1], "holds" if good else "FAILS"), flush=True) print("controls %s [%.1f s]" % ("hold" if ok else "FAIL", time.time() - t0)) if not ok: return False print() print(" r e least g p H/exact read of below / P1 misses") kill, per_e, tot = [], {}, [0, 0] rs = [r for r in range(21, 290, 2)] + [r for r in range(38, 577, 2)] for r in rs: e, (v, g, p, ratio), n, below, p1, m = rank(r) tot[0] += n tot[1] += m x = math.exp(v) if not 0.9 <= ratio <= 1.1: kill.append(r) cur = per_e.setdefault((e, r % 2), [math.inf, 0, 0, 0, 0]) if x < cur[0]: cur[:3] = [x, r, g] cur[3] += len(below) cur[4] += len(p1) print("%3d %2d %.4f %3d %4d %.4f %4d %3d %s %s" % (r, e, x, g, p, ratio, n, m, below or "", p1 or ""), flush=True) print() print(" e parity least at r g shapes below P1 misses") for (e, par), (x, r, g, nb, np1) in sorted(per_e.items()): print("%2d %s %.4f %4d %3d %d %d" % (e, "odd " if par else "even", x, r, g, nb, np1)) print("%d shapes over %d ranks, each read exact at its hand-best p;" " %d read by best_x" % (tot[1], len(rs), tot[0])) print("KILL %s" % ("fires at r = %s" % kill if kill else "misses")) print("%.1f s" % (time.time() - t0)) return ok if __name__ == "__main__": sys.exit(0 if main() else 1)