"""takeph.py -- at the takeover theorem's depth, is the least pure product the champion, and does F_e beat it? QUESTION. At J* = (4e - 1) l - 1, l = r - 2, e = floor(sqrt(l/2)) at odd r and floor(sqrt(l)/2) at even r, the takeover theorem gives height(F_e) < H_ch. The cell (J* + r, J*) fails when height(F_e) < ph. Is ph = H_ch there, which shape attains ph, and is F_e below it? r = 5..25. BY HAND (re-derived before the code). An admissible pure product is (x - 1)^(J+k) prod [g_i], [g] = 1 + x + ... + x^(g-1), parts g_i >= 2, k + sum(g_i - 1) <= r - 1, k >= 0, at most J + k parts >= 2. Its order at -1 is s = #even parts <= r - 1. s = r - 3 has the five cofactors 1, 1 + x, (1 + x)^2, 1 + x^2, 1 - x + x^2 (after x -> -x); the last is F_1, so at e = 1 the theorem's witness is itself pure and the cell need not fail. PREDICTIONS (frozen before the code). P1 at the e = 1 ranks (r = 5, 7, 9 and even 6..16) ph < H_ch (proved: F_1 is pure), ph attained by the part-3 shape at most of them. P2 at the e >= 2 ranks (odd 11..25, even 18..24) height(F_e) < ph, ph/height(F_e) >= 1.05 at each. P3 there ph = H_ch at most ranks, else the part-3 shape, ph/H_ch in [0.9, 1]. P4 the least order-(r - 2) shape over H_ch within 30% of sqrt(J*/(4(r - 1))) (TRANSPLANT: Gaussian, M_(s-1)/M_s ~ s^(-1/2)). KILL. K1 an e >= 2 rank where height(F_e) >= ph is printed. CONTROLS, read first. C1 the restricted ph against flatten.ph_of( pure_table(23)) at every cell M <= 24; its PLANT: k pinned to 0 (frozen), and the parts capped at 3 (added after the run, the k plant read nothing). C2 the champion in the scan against circlemean's exact_Hch. C3 the lattice route at (J* + r, J*): the witness clears J, h <= height(F_e), the enumeration complete by prunerr's bound (3E + 2E^2 + 2^-51 < 1e-9; its coordinate and overflow conditions are not printed). A failure h < ph stands on its cleared witness and the scan's ph (C1); completeness carries only h = ph. FINDINGS (copied from the printed output). (1) C1 0 of 276 cells differ. The parts-at-3 plant is caught at 80 of 276; the k plant at none: no cell of width <= 24 needs k > 0, so C1 reads nothing for k (k > 0 attains no ph and no least non-champion below either). C2 holds. C3 holds at every rank, E <= 9.48e-14. (2) K1 misses: ph/height(F_e) >= 1.8895 (r = 21) over the e >= 2 arm (P2 holds). The lattice route reads h < ph at the cell (J* + r, J*) at r = 11, 13..25: J_hi(r) >= J* there, computed. At r = 5..10 and 12 (e = 1) the cell does not fail, h = ph; at r = 14, 16 (e = 1) it fails, h/ph = 0.7407 and 0.6655. (3) P1: ph/H_ch = 0.1103 .. 0.4000 at the e = 1 ranks; ph is the part-3 shape F_1 at r = 5, 7, 9 only; one part 3 with an excess unspent (order r - 4) at r = 6, and two or three parts 3 at r = 8, 10, 12, 14, 16: P1's shape half fails. (4) P3 fails: ph = H_ch at r = 11, 13 (e = 2, the part 3 tying it, e - 1 = 1) and 21, 23, 25 (e = 3); at r = 15 and 17..24 (e = 2) ph is the part-5 shape k = 0 2^(r-5) 5, of order r - 5, ph/H_ch = 0.9724, 0.8890, 0.8621, 0.8393, 0.8403, 0.8236, 0.8109 at r = 15, 17, 18, 19, 20, 22, 24. (5) P4 holds: the order-(r - 2) shape over H_ch reads 1.2169 .. 1.3423 at e = 2 and 1.6594 .. 1.6701 at e = 3, within 5.1% of the Gaussian 1.2450 .. 1.6202; at e = 1 0.7857 .. 0.9667 against 0.7071 .. 0.8266. (6) At odd r the five order-(r - 3) shapes stand at exactly e - 1 (the part 3), e (cofactor 1), 2e - 1 (the part 4), 2e - 1/2 (k = 1) and 4e - 1 (k = 2) times H_ch, exact Fractions, at every odd rank read. By hand, Q = (1 + x)^J (1 - x)^(r-3): J* makes the Beta mean exactly 1/(4e), so the central coefficient of x Q is -e times the champion's, F_t's is (t - e) times it, and F_e's vanishes; the k = 1 shape's middle pair is 2e(1 - mu) = 2e - 1/2. That each height sits at its centre is read off the run. READS (not verdicts). The middle coefficient bounds no class: it is 0 at some shape at every rank read. The test family (p <= 16, q = 0, 1, k within one of centring) puts every non-champion pure product above F_e's circle mean at e = 3 (least test/H_ch 1.3174, 1.3189, 1.3201 at r = 21, 23, 25 against means 0.7370, 0.7016, 0.6709, the least at the order-(r - 2) shape) and fails at every e = 2 rank read (0.1350 .. 0.4928 against 0.4581 .. 0.6905, the least at part-5 shapes). RUN RECORD. 27.3 s, peak commit 23.5 MB (the reads 21 s of it). """ import math import os import sys import time from fractions import Fraction os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import mpmath as mp # noqa: E402 import flatten # noqa: E402 import circlemean # noqa: E402 import prunerr # noqa: E402 def mul_block(p, g): """p * (1 + x + ... + x^(g-1)), by a running window.""" n = len(p) + g - 1 out = [0] * n acc = 0 for i in range(n): if i < len(p): acc += p[i] if i - g >= 0: acc -= p[i - g] out[i] = acc return out def pure_scan(J, r, kmax=None, xmax=None): """Every admissible pure product at (J + r, J): yields (height, k, parts) with parts the multiset of g_i >= 2, nonincreasing (kmax, xmax: the PLANTS, k and the part excess g - 1 capped).""" top = r - 1 for k in range(0, top + 1 if kmax is None else min(top, kmax) + 1): base = flatten.pow_xm1(J + k) room = top - k stack = [(base, (), room, room if xmax is None else xmax)] while stack: p, parts, left, cap = stack.pop() yield flatten.height(p), k, parts if len(parts) >= J + k: continue for x in range(min(left, cap), 0, -1): stack.append((mul_block(p, x + 1), parts + (x + 1,), left - x, x)) def order(parts): return sum(1 for g in parts if g % 2 == 0) def champion(r): return (0, (2,) * (r - 1)) def read_cell(J, r, kmax=None, xmax=None): """ph with its shape, the champion's height, and per order s the least height with its shape.""" best = None by_s = {} hch = None for h, k, parts in pure_scan(J, r, kmax, xmax): if best is None or h < best[0]: best = (h, k, parts) s = order(parts) if s not in by_s or h < by_s[s][0]: by_s[s] = (h, k, parts) if (k, parts) == champion(r): hch = h return best, hch, by_s def theorem_e(r): l = r - 2 return math.isqrt(l // 2) if r % 2 else math.isqrt(l) // 2 def theorem_J(r): return (4 * theorem_e(r) - 1) * (r - 2) - 1 def name(k, parts): """The shape in words: k extra (x - 1)'s and the parts above 2.""" big = [g for g in parts if g != 2] twos = len(parts) - len(big) return "k=%d 2^%d%s" % (k, twos, "".join(" %d" % g for g in big)) def c1(kmax=None, xmax=None): """The restricted ph against flatten's full table, every M <= 24.""" table = flatten.pure_table(23) diff = 0 cells = 0 for M in range(2, 25): for J in range(1, M): r = M - J best, _, _ = read_cell(J, r, kmax, xmax) want, _ = flatten.ph_of(table, M, J) cells += 1 if best[0] != want: diff += 1 return cells, diff def walk_all(J, r): """Every admissible pure product's coefficients with its shape.""" for k in range(r): room = r - 1 - k stack = [(flatten.pow_xm1(J + k), (), room, room)] while stack: p, parts, left, cap = stack.pop() yield p, k, parts if len(parts) >= J + k: continue for x in range(min(left, cap), 0, -1): stack.append((mul_block(p, x + 1), parts + (x + 1,), left - x, x)) def test_family(c, pmax=16): """The best of ||/|w|_1 over w = (1 + x)^p (1 - x)^q, q = 0, 1, p <= pmax, k within one of centring: a lower bound on height(c).""" best = Fraction(0) n = len(c) for q in (0, 1): for p in range(pmax + 1): v = [math.comb(p, i) for i in range(p + 1)] if q: v = [a - b for a, b in zip(v + [0], [0] + v)] l1 = sum(abs(a) for a in v) k0 = (n - len(v)) // 2 for k in (k0 - 1, k0, k0 + 1): if 0 <= k and k + len(v) <= n: t = Fraction(abs(sum(c[k + i] * a for i, a in enumerate(v))), l1) best = max(best, t) return best def reads(): """READS, not verdicts: the middle coefficient and the test family as lower bounds over every non-champion pure product, against F_e's circle mean, at J*.""" for r in (11, 13, 15, 17, 18, 19, 21, 23, 24, 25): J, e = theorem_J(r), theorem_e(r) hch = circlemean.exact_Hch(r, J) mean = float(mp.exp(circlemean.log_mean(r, J, e))) / hch mid = tst = None o3 = [] for c, k, parts in walk_all(J, r): if order(parts) == r - 3 and r % 2: o3.append((Fraction(flatten.height(c), hch), name(k, parts))) if (k, parts) == champion(r): continue n = len(c) - 1 m = Fraction(max(abs(c[n // 2]), abs(c[(n + 1) // 2])), hch) t = test_family(c) / hch if mid is None or m < mid[0]: mid = (m, name(k, parts)) if tst is None or t < tst[0]: tst = (t, name(k, parts), flatten.height(c) / hch) print("READ r = %2d e = %d: mean(F_e)/H_ch %.4f; least middle" " coefficient/H_ch %.4f (%s); least test/H_ch %.4f at %s" " (its height/H_ch %.4f)" % (r, e, mean, mid[0], mid[1], tst[0], tst[1], tst[2]), flush=True) if o3: print(" order r - 3 heights/H_ch: %s" % ", ".join("%s %s" % (q, n) for q, n in sorted(o3))) def main(): t0 = time.time() cells, diff = c1() print("C1 restricted ph against pure_table(23): %d cells, %d differ" % (cells, diff), flush=True) pc, pdiff = c1(kmax=0) print("C1 PLANT (k pinned to 0): %d cells, %d differ [%.1f s]" % (pc, pdiff, time.time() - t0), flush=True) qc, qdiff = c1(xmax=2) print("C1 PLANT (parts at most 3): %d cells, %d differ [%.1f s]" % (qc, qdiff, time.time() - t0), flush=True) ok = diff == 0 and qdiff > 0 c2 = c3 = True worst = None rows = [] for r in range(5, 26): e = theorem_e(r) J = theorem_J(r) best, hch, by_s = read_cell(J, r) c2 &= hch == circlemean.exact_Hch(r, J) hf = circlemean.exact_height_F(r, J, e) fe = {t: circlemean.exact_height_F(r, J, t) for t in range(max(0, e - 2), e + 4)} eb = min(fe, key=fe.get) o2 = by_s.get(r - 2) nonch = min((v for v in by_s.values() if (v[1], v[2]) != champion(r)), key=lambda v: v[0]) print("r = %2d e = %d J* = %3d: ph/H_ch %.4f at %s (s = %d);" " F_e/H_ch %.4f, ph/F_e %.4f; best e %d (F/H_ch %.4f);" " least non-champion %s /H_ch %.4f; order r-2 %s /H_ch %.4f" " (Gauss %.4f) [%.1f s]" % (r, e, J, best[0] / hch, name(best[1], best[2]), order(best[2]), hf / hch, best[0] / hf, eb, fe[eb] / hch, name(nonch[1], nonch[2]), nonch[0] / hch, name(o2[1], o2[2]), o2[0] / hch, math.sqrt(J / (4 * (r - 1))), time.time() - t0), flush=True) print(" least per order s/H_ch: %s" % ", ".join("s=%d %.3f" % (s, by_s[s][0] / hch) for s in sorted(by_s, reverse=True)[:5]), flush=True) h, v, nodes = flatten.least_height(J + r, J) b, d, lam = flatten.lll(flatten.basis(J + r, J)) E = prunerr.err_bound(d, lam, r)[0] whole = 3 * E + 2 * E ** 2 + mp.mpf(2) ** -51 < mp.mpf(10) ** -9 c3 &= (flatten.clears(v, J) and flatten.height(v) == h and h <= hf and whole) print(" C3 lattice h/H_ch %.4f, h/ph %.4f, E %.2e (complete %s)," " %d nodes" % (h / hch, h / best[0], float(E), whole, nodes), flush=True) rows.append((r, e, J, best[0], hch, hf, h)) if e >= 2: q = best[0] / hf if worst is None or q < worst[0]: worst = (q, r) print("C2 the champion's height equals circlemean.exact_Hch: %s" % c2) print("C3 the lattice witness clears J, its height is h <= height(F_e)," " the enumeration complete by its bound, every rank: %s" % c3) fails = [r for r, e, J, ph, hch, hf, h in rows if h < ph] print("cells (J* + r, J*) with h < ph: r = %s" % fails) k1 = [r for r, e, J, ph, hch, hf, h in rows if e >= 2 and hf >= ph] print("K1 least ph/height(F_e) over the e >= 2 arm: %.4f at r = %d;" " K1 %s" % (worst[0], worst[1], "fires at %s" % k1 if k1 else "misses")) reads() print("%.1f s" % (time.time() - t0)) return not (ok and c2 and c3) or bool(k1) if __name__ == "__main__": sys.exit(1 if main() else 0)