"""room.py -- the room the explicit formula of L(s, rho) leaves at a complex cubic field's least odd partial prime, against the side of that place in Cl/2Cl. QUESTION. K a complex cubic field of 2-rank 1, L/K its unramified quadratic extension, chi the character of Cl/2Cl. The least odd partially split place lies outside 2Cl at 0.837 (|d_K| <= 24000) and 0.800 (24000 < |d_K| <= 48000) where the limit share is 1/2 (seat.py). Does the explicit formula of L(s, rho) = L(s, chi; L/K) leave too little room at that prime p0 for a(p0) = +1, the place in 2Cl, and does the room predict the side? THE ARGUMENT (written before the engine). (1) THE REPRESENTATION. Gal(M/Q) = S4, K fixed by D4 = <(1234), (13)>, whose index-2 subgroups are C4, the normal V and V' = <(13), (24)>. Complex conjugation is a transposition and the transpositions of D4 lie in V' alone, so L, where the real place of K splits, is the fixed field of V'. Ind chi has character 3, 1, -1, -1, 0 at e, (12), (12)(34), (1234), (123): the standard rho. L(s, rho) is a Hecke L-function of a nontrivial character, entire. (2) CONDUCTOR AND GAMMA. f(Ind chi) = |d_K|^chi(1) N f(chi) = |d_K|, chi unramified at every place. rho(c) has trace 1 in dimension 3, so gamma(s) = Gamma_R(s)^2 Gamma_R(s + 1), the gamma of zeta_K as well (Gamma_R(s) Gamma_R(s + 1) = Gamma_C(s)). (3) THE FORMULA. F even, supported in [-T, T], F^(t) = int F(u) e^{itu} du. Over the strip, Lambda'/Lambda(s) = -Lambda'/Lambda(1 - s) (rho real): H(T) = sum_gamma F^(gamma) = F(0) log q + G(T) - 2 sum_n Lambda(n) c(n) F(log n) / sqrt(n), G(T) = (1/2 pi) int F^(t) sum_j [-log pi + Re psi((1/2 + kappa_j + it)/2)] dt, kappa = (0, 0, 1), c(p^k) = tr(Frob^k | V^I), the coefficient of -L'/L. A Dedekind zeta adds Phi(0) + Phi(1) = 2 int F(u) cosh(u/2) du, its poles. Under GRH, with F^ >= 0, the HEADROOM H(T) is >= 0. (4) THE TEST FUNCTION. Odlyzko's f(x) = (1 - |x|) cos(pi x) + sin(pi |x|)/pi on |x| <= 1, the autocorrelation of cos(pi x) on |x| <= 1/2: f >= 0 and f^(t) = 8 pi^2 cos^2(t/2) / (pi^2 - t^2)^2 >= 0, f^(pi) = 1/2. F_T(u) = f(u/T), F_T^(t) = T f^(tT), F(0) = 1. (5) THE COEFFICIENTS, each read off the local factor of zeta_L/zeta_K from how p and its places split in K and in L, so wild primes are included: split, three places in 2Cl 3 at every k (e) split, two outside -1 at odd k, 3 even (12)(34) partial, in 2Cl 1, 3, 1, 3, ... (12) partial, outside: 1/((1 + p^-s)(1 + p^-2s)) -1, -1, -1, 3, ... (1234) inert 0, 0, 3, ... (123) P^2 Q, [P] in 2Cl ([Q] = -2[P]) 2 at every k P^2 Q, [P] outside 0 at odd k, 2 even P^3 (3[P] = 0, so [P] in 2Cl) 1 at every k and for zeta_K, c_K(p^k) = sum over P | p with f_P | k of f_P. (6) THE ROOM. p_r the r-th odd unramified partial prime, term_+-(T) = 2 log p_r sum_k c_+-(p_r^k) F_T(k log p_r) / p_r^(k/2) with the partial rows of (5), and R(T) = log q + G(T) - 2 (the prime sum without p_r's powers), so H = R - term(a). FORCED: R(T) < term_+(T) at some width, a GRH theorem that a(p_r) = -1 (forced by the formula, not by generation as in seat.py's (5)). SLACK s(T) = (R - term_+)/Delta, Delta = term_+ - term_- > 0, read only where Delta >= 0.25 (s blows up as Delta -> 0); s* its least over the widths. (7) THE CONFOUND. R = H + term(a) with H >= 0: were the zero sum H indifferent to a, R would still sit lower by Delta on the fields reading -1, a rise in R's bins that the identity alone makes. So the kill leans on the CLASS-BLIND room R_K, R with every class-dependent coefficient replaced by its mean over Cl/2Cl uniform: split 0 at odd k, 3 at even; partial 0 at odd k, 1 at k = 2 mod 4, 3 at k = 0 mod 4; P^2 Q 1 at odd k, 2 at even. R_K is a function of K's splitting and conductor alone, not built from a. (The frozen record had the partial mean 1 at every even k; the mean at k = 0 mod 4 is (3 + 3)/2, corrected before the first run. It reaches only 2^4 and 3^4.) s_K* as s*. 2's splitting type grades the side and moves R_K through 2's own terms, so the bar sits within 2-types. (8) THE WAIT. A prime is OUTSIDE when a degree-1 place over it lies outside 2Cl; the split and partial primes, 2 included, in order of size. A split prime's three classes sum to [(p)] = 0, so mod 2Cl its pattern is 000, 110, 101 or 011, outside at 3/4 under Chebotarev, partial at 1/2; the third place over a split prime is fixed by the first two, so the wait is read prime by prime. At h = 2, 2Cl = {0} and the principal level of a place is 2 (places in 2Cl)/places, null 1. (9) THE WORST CASE. Every class-dependent coefficient of (5) is least, at every k, on its outside value, and the weights are >= 0. So U_A(T), the headroom with the primes of A in 2Cl and every other class-dependent prime outside, bounds H(T) over every configuration with A in: U_A < 0 at some width forbids A under GRH from K's splitting and |d_K| alone, not from the sides (the method of Belabas, Diaz y Diaz and Friedman, Math. Comp. 77 (2008), for small generators of Cl). THE ENGINE. cubic.py: the fields by Hunter's theorem, the certified class group (always_cert), the class of every degree-1 place; a class is outside 2Cl when odd at the even invariant. Every prime p <= e^6 = 403 is read, checkpointed. Widths T in {1.5, 2, 2.5, 3, 3.5, 4, 5, 6}. G(T) and the pole term by mpmath quadrature at 20 digits, cubic.py's 60 held everywhere else. PREDICTIONS, frozen before the engine. C1 f^ closed form against a numerical transform, within 1e-8. C2 the formula on Riemann's zeta (q = 1, kappa = 0, two poles) against its first 100 zeros, T = 2, 4, 6: within 1e-5. C3 the formula on L(s, chi_-3) (q = 3, kappa = 1) against its zeros found on the critical line to t = 60: within 1e-4. C4 every field: [Q] = -2[P] at P^2 Q, 3[P] = 0 at P^3, an even number of places outside 2Cl at every split prime, every place placed. C5 GRH positivity: H(T) >= 0 for L(s, rho) and for zeta_K at every field and width, and no forced field reads +1 (the strong arm, read as an error before anything else). C6 zeta_K's bound with the prime sum dropped, the largest over T of -(G + Phi(0) + Phi(1)), at most log 23. C7 the break: a(p) flipped at every partial prime makes some H(T) negative, or C5 cannot catch a wrong coefficient. C8 the kill's control: a planted side, outside with chance 0.95 on each 2-type's lower half of s_K* and 0.6 on its upper, reads D within 3 se of 0.35. KILL. D = the share outside 2Cl over the half of each 2-type with the smaller s_K* minus the half with the larger, pooled over the 2-types weighted by size (a type of fewer than 40 fields dropped), se binomial. The class-blind room explains the lean only if D > 3 se and D > 0.05. Q2 the forced share (full R, position 1) between 0.1 and 0.6. Q3 the KILL misses, D near 0.1 to 0.2. (Doubt: seat.py reads the share flat over p0 = 3, 5, 7, 11, and p0 is most of s_K*'s spread within a 2-type.) Q4 printed, not judged: the share by quintile of s* and of s_K*, the forced share at positions 1 to 4, corr(R - R_K, side). Frozen before section G (the wait): C9 a memory plant (outside at 0.83 partial, 0.93 split until the first outside prime, then at the null) reads the weld's gap > 0.3 and 3 se and holds the memoryless kill. C10 a fixed plant (each (p, type) at its observed share) reads the weld's gap within 3 se of 0. Q5 odd partial position 3 after two in 2Cl: 0.75 to 0.9; after exactly one earlier outside, 0.6 to 0.7 at positions 2 and 3. Q6 the prime-level all-in share above 0.7 (partial) and 0.85 (split), flat in r. Q7 at h = 2, p <= 7: L0 <= 0.5, L1 >= 0.8. KILL MEMORYLESS. The partial all-in share pooled over r = 2..4 minus r = 1; memoryless dies past 3 se. KILL WELD. At h = 2, |d_K| <= 12000, p <= 7, the level after an outside prime minus before (L1 - L0), within (p, type) since the level rises with p, pooled; the shortfall is the wait's only if it exceeds 0.3 and 3 se. Frozen after section G, before section H (the worst case): C11 U_empty >= H at every field and width; no field's own first-two pattern forbidden. C12 U at the actual sides reproduces H to 1e-9. Q8 the first prime in 2Cl forbidden at 0.3 to 0.7 of the fields. KILL CONDUCTOR. "The first two primes are never both in 2Cl" is the conductor's only if the worst case forbids it at >= 0.95 of the fields. CHANGED FROM THE FROZEN RECORD while the engine was written, before any run: C3's zeros read to t = 60, not 100 (the tail past 60 is inside the agreement C3 prints); C8's plant set on each 2-type's halves of s_K*, not the global median, its band 3 se; C4 also asks every place placed. The frozen record of section G worded KILL WELD inverted ("the wait's unless" the gap passes); its bar's direction, a gap a fixed weight cannot make, is the code's, written before the run: alive only past it. ADDED AFTER THE FIRST RUN, prints and not verdicts: D at each fixed width (s_K at one T, not its least), and position 2's side split by position 1's. After section H's first run: the count of fields below the least free discriminant, and the first prime forbidden principal at h = 2, |d_K| <= 12000. FINDINGS. F0 CONTROLS (C1 to C8 held). f^ to 3.3e-17; the formula against zeta's first 100 zeros to 1.5e-7 and against L(s, chi_-3)'s 22 zeros to t = 60 to 1.1e-5; at 1532 fields of 2-rank 1, no split parity, [Q] = -2[P] or 3[P] = 0 broken, no place unplaced; the least H(T) over fields and widths 0.0013, zeta_K's 0.0067; zeta_K's prime-free bound 3.0875 against log 23 = 3.1355; the flipped partial sign drives H negative at 1471 fields; the planted side reads D = 0.363 +- 0.020. The population reproduces seat.py's shares at positions 1 to 4: 0.817, 0.691, 0.654, 0.624. F1 THE ZERO SUM VANISHES AT WIDE SUPPORT. H's median falls 4.586, 3.573, 2.803, 2.133, 1.581, 1.079, 0.352, 0.064 over T = 1.5 to 6, while log 24000 + G(T) falls from 4.53 to 0.11. On a field reading -1 with H below Delta the formula, an identity, leaves no room for +1; on one reading +1 both values fit, so only a = -1 is ever forced: 1206 of 1532 fields (0.787) are FORCED at some width, a(p0) = -1 under GRH, every one reading outside; of the 326 unforced, 45 (0.138) read outside, the rest in 2Cl. Q2 FAILED (0.787, not 0.1 to 0.6): the full room reads the side back and explains nothing, the confound (7) at full strength. Through position 4 the forced share tracks the side (0.646, 0.587, 0.550 against 0.691, 0.654, 0.624). F2 THE KILL FIRED (Q3 FAILED). Within 2-types D = 0.030 +- 0.019 (1.6 se): inert 0.062 +- 0.041 (357), partial -0.020 +- 0.041 (508), P^2 P' 0.047 +- 0.032 (340), P^3 -0.009 +- 0.026 (222), split 0.192 +- 0.068 (105). Pooled over the types the quintiles of s_K* fall 0.863, 0.846, 0.860, 0.788, 0.726, a fall the within-type halves show only as D = 0.030; 2's grading, the confound (7) places the bar against, is the likely carrier (s_K* by type not printed). s_K* is below 0 at every field, so the class-blind room never leaves room for +1 at its least width. The class-dependent part of the room carries what R reads: corr(R - R_K, side) = -0.345 at T = 3 over the 1516 fields with p0 < e^3. F3 printed, not judged: D at one width, each over the fields with Delta >= 0.25 at that width (no 2-type reaches 40 such fields at T = 1.5), reads 0.032, 0.039, 0.059, 0.085, 0.070, 0.035, 0.030 at T = 2 to 6 (+- 0.019 to 0.031). At T = 3.5 (4.2 se) and T = 4 (3.7 se) it passes the frozen bar at that width, two of a scan of seven; the frozen kill, on s_K*, fired. F4 printed, not judged: THE NEXT PLACE PAYS. Position 2 lies outside 2Cl at 0.840 +- 0.022 (281) where position 1 lies in 2Cl and at 0.658 +- 0.013 (1251) where it lies outside. A mixture of field types biasing both positions alike would correlate the two sides positively; they anticorrelate (a reading, not tested: a budget paid by whichever early place reads -1 would). A position-2 place after a square reads as position 1 does (0.840 against 0.817). F5 THE WAIT IS NOT MEMORYLESS (C9, C10 held; KILL MEMORYLESS FIRED). The first prime is outside at 0.857 +- 0.010 (1235, partial) and 0.997 +- 0.003 (297, split, null 3/4); a second prime after a first in 2Cl is outside at 1.000 (142 partial, 35 split): the first two primes are never both in 2Cl over the 1532 fields. After exactly one outside prime the partial share rises with r, 0.711, 0.831, 0.864 at r = 2, 3, 4 (1051, 384, 44). Q6's flatness failed on the all-in share itself (0.857, then 1.000; the pool at 1.000 has no spread, its se 0, and 0.857^142 is 3e-10); both rises are the shape a budget makes. Odd partial places (W1): position 3 after two in 2Cl 0.778 +- 0.062 (45), Q5's band held; after exactly one outside 0.658 and 0.726 at positions 2 and 3, Q5's band held at 2 and missed at 3. F6 THE WELD (KILL WELD MISSED, Q7 half failed). At h = 2, |d_K| <= 12000 (220 fields), p <= 7: L0 0.176 (224 field-primes), L1 0.541 (238), within (p, type) L1 - L0 = 0.432 +- 0.064; to |d_K| <= 48000 (1016 fields) 0.329, 0.644, 0.452 +- 0.036. The history carries part of the shortfall at the smallest primes, but a place after one paid still reads 0.54, not 1. By 11 <= p <= 89 only 4 field-primes have no earlier outside prime. F7 THE CONDUCTOR FORBIDS IT (C11, C12 held; KILL CONDUCTOR MISSED). U_empty - H least -8.9e-16, U at the actual sides H to 2.7e-15. Under GRH, from the splitting and |d_K| alone: the first prime in 2Cl is forbidden at 925 of 1532 fields (0.604, Q8 held), the first two both in at 1511 (0.986), at all 601 fields with |d_K| < 21767. The 21 left free all have two partial primes first and |d_K| >= 21767; none reads both in (6 of 21 read the first in, 6 the second; independent at those shares, about 1.7 would). At h = 2, |d_K| <= 12000 the first prime principal is forbidden at 207 of 220. TIERS. The coefficient table and the explicit formula are property; C5's positivity is GRH read at 1532 fields; H's smallness at T = 6 and F1 to F6 are observations over this population; F7's forbidding is a rule under GRH at each field read (every field to |d_K| < 21767). RUN RECORD. 8/8 controls, 343.0 s with the population read (274.8 s, 7975 fields, checkpointed), 68 s from the checkpoint; peak commit 111.9 MB and 112.1 MB under a memory guard (512 MB). The first rehearsal at --cap 6000 failed C4 at 12 fields: a global mpmath.mp.dps = 20 had cut cubic.py's 60 digits under its regulator, now scoped to the formula section. The F3 and F4 prints were added after the first full run, which reproduced every earlier line. Sections G and H: 12/12 controls, 68.3 s from the checkpoint, peak commit 135.6 MB under a memory guard; the first G run stopped on a helper named as a local list, renamed; H's added prints and W3's relabelled counts reproduced every earlier line. """ import hashlib import inspect import json import os import random import sys import tempfile import time from math import cos, exp, log, pi, prod, sin, sqrt os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") HERE = os.path.dirname(os.path.abspath(__file__)) sys.path.insert(0, HERE) import mpmath # noqa: E402 import cubic as C # noqa: E402 RESULTS = [] def _flag(name, default): if name in sys.argv: return sys.argv[sys.argv.index(name) + 1] return default CAP = int(_flag("--cap", 48000)) # the population: |d_K| <= CAP WIDTHS = [1.5, 2, 2.5, 3, 3.5, 4, 5, 6] PMAX = int(exp(max(WIDTHS))) # 403: every p^k <= e^T is read DELTA_MIN = 0.25 # s is read only where Delta >= this POSITIONS = 4 MIN_TYPE = 40 # a 2-type below this leaves the kill SEED = 1655 def check(name, ok, detail=""): RESULTS.append(bool(ok)) print(f" [{'PASS' if ok else 'FAIL'}] {name}" + (f" -- {detail}" if detail else "")) def section(title): print(f"\n{title}") # --------------------------------------------------------------------------- # the test function and the archimedean terms def f(x): x = abs(float(x)) if x >= 1: return 0.0 return (1 - x) * cos(pi * x) + sin(pi * x) / pi def fhat(x): x = mpmath.mpf(x) if abs(abs(x) - mpmath.pi) < mpmath.mpf(10) ** -8: return mpmath.mpf(1) / 2 return 8 * mpmath.pi ** 2 * mpmath.cos(x / 2) ** 2 \ / (mpmath.pi ** 2 - x * x) ** 2 def gamma_piece(T, kappa): """(1/2 pi) int F_T^(t) [-log pi + Re psi((1/2 + kappa + it)/2)] dt for one Gamma_R(s + kappa).""" half = mpmath.mpf(1) / 2 def g(t): return T * fhat(t * T) * mpmath.re( mpmath.digamma((half + kappa + 1j * t) / 2)) pts = [k * mpmath.pi / T for k in range(201)] + [mpmath.inf] return float(-mpmath.log(mpmath.pi) + mpmath.quad(g, pts) / mpmath.pi) def pole(T): """Phi(0) + Phi(1) = 2 int F_T(u) cosh(u/2) du over [-T, T].""" return float(4 * mpmath.quad(lambda u: f(u / T) * mpmath.cosh(u / 2), [0, T])) PRIMES = [p for p in C.PRIMES if p <= PMAX] def weights(T): """{(p, k): 2 log p F_T(k log p) / p^(k/2)} over p^k < e^T.""" W = {} for p in PRIMES: k = 1 while k * log(p) < T: W[(p, k)] = 2 * log(p) * f(k * log(p) / T) / p ** (k / 2) k += 1 return W W = {T: weights(T) for T in WIDTHS} G = {} # G(T) for L(s, rho) and zeta_K: kappa = 0, 0, 1 POLE = {} # --------------------------------------------------------------------------- # coefficients def coeff(typ, v, k, flip=False): """c(p^k) for L(s, rho), (5); v is the side datum of the prime.""" if typ == "s": # v: places outside 2Cl, 0 or 2 return 3 if (v == 0 or k % 2 == 0) else -1 if typ == "p": # v: 1 when outside 2Cl if v ^ flip: return 3 if k % 4 == 0 else -1 return 1 if k % 2 else 3 if typ == "n": return 3 if k % 3 == 0 else 0 if typ == "r2": # v: 1 when [P] is outside 2Cl return 2 if (not v or k % 2 == 0) else 0 return 1 # r3 def blind(typ, k): """The mean of coeff over Cl/2Cl uniform, (7).""" if typ == "s": return 0 if k % 2 else 3 if typ == "p": return 0 if k % 2 else (3 if k % 4 == 0 else 1) if typ == "r2": return 1 if k % 2 else 2 return coeff(typ, 0, k) def coeff_K(typ, k): """c_K(p^k) for zeta_K.""" return {"s": 3, "p": 1 + 2 * (k % 2 == 0), "n": 3 * (k % 3 == 0), "r2": 2, "r3": 1}[typ] def partial_term(T, p, out): return sum(w * coeff("p", out, k) for (q, k), w in W[T].items() if q == p) # --------------------------------------------------------------------------- # the population def read_field(key, dk): K = C.Field(key, dk) if not K.class_group(always_cert=True): return {"d": dk, "ok": False} inv = list(K.cl) rec = {"d": dk, "ok": True, "inv": inv} if sum(1 for n in inv if n % 2 == 0) != 1: return rec def outside(c): return any(a % 2 for a, n in zip(c, inv) if n % 2 == 0) def zero(c): return all(a % n == 0 for a, n in zip(c, inv)) primes, bad, unplaced = [], 0, 0 for p in PRIMES: pls = K.places(p) typ = K.ptype(p) if typ == ((3, 1),): primes.append([p, "n", 0]) continue cs = {pl: K.place_class(pl) for pl in pls if pl.deg == 1} if any(c is None for c in cs.values()): unplaced += 1 continue if typ == ((1, 1), (1, 1), (1, 1)): n_out = sum(outside(c) for c in cs.values()) bad += n_out % 2 primes.append([p, "s", n_out]) elif typ == ((1, 1), (2, 1)): primes.append([p, "p", int(outside(list(cs.values())[0]))]) elif typ == ((1, 1), (1, 2)): cP = [c for pl, c in cs.items() if pl.e == 2][0] cQ = [c for pl, c in cs.items() if pl.e == 1][0] bad += not zero([a + 2 * b for a, b in zip(cQ, cP)]) primes.append([p, "r2", int(outside(cP))]) elif typ == ((1, 3),): cP = list(cs.values())[0] bad += not zero([3 * a for a in cP]) primes.append([p, "r3", 0]) else: raise ValueError("splitting %r" % (typ,)) rec.update(primes=primes, bad=bad, unplaced=unplaced) return rec _SRC = hashlib.sha1( open(os.path.join(HERE, "..", "code", "cubic.py"), "rb").read() .split(b'"""', 2)[2] + inspect.getsource(read_field).encode() + str(PMAX).encode()).hexdigest()[:12] CKPT = os.path.join(tempfile.gettempdir(), "room_%d_%s.json" % (CAP, _SRC)) def population(): if "--fresh" not in sys.argv and os.path.exists(CKPT): print(" read from the checkpoint %s" % CKPT) with open(CKPT) as fh: return json.load(fh) t0 = time.time() F = C.enumerate_fields(CAP) print(" %d complex cubic fields with |d_K| <= %d" % (len(F), CAP)) recs = [] for i, (key, dk) in enumerate(sorted(F.items(), key=lambda kv: (-kv[1], kv[0]))): recs.append(read_field(key, dk)) if (i + 1) % 1000 == 0: print(" %d fields read, %.0f s" % (i + 1, time.time() - t0), flush=True) with open(CKPT, "w") as fh: json.dump(recs, fh) print(" read in %.1f s, checkpointed to %s" % (time.time() - t0, CKPT)) return recs # --------------------------------------------------------------------------- # one field's room class Room: """The headroom, the room at the first POSITIONS odd partial primes (full and, at position 1, class-blind), and zeta_K's headroom.""" def __init__(self, rec, flip=False): self.d = rec["d"] logq = log(-self.d) pr = {p: (t, v) for p, t, v in rec["primes"]} self.two = pr[2][0] self.parts = [(p, v) for p, t, v in rec["primes"] if t == "p" and p > 2][:POSITIONS] self.H, self.HK, self.R, self.RK = {}, {}, {}, {} for T in WIDTHS: byp, byp_b = {}, {} hk = logq + G[T] + POLE[T] for (p, k), w in W[T].items(): t, v = pr[p] byp[p] = byp.get(p, 0.0) + w * coeff(t, v, k, flip) byp_b[p] = byp_b.get(p, 0.0) + w * blind(t, k) hk -= w * coeff_K(t, k) full, full_b = sum(byp.values()), sum(byp_b.values()) self.H[T] = logq + G[T] - full self.HK[T] = hk self.R[T] = [logq + G[T] - full + byp.get(p, 0.0) for p, _ in self.parts] self.RK[T] = (logq + G[T] - full_b + byp_b.get(self.parts[0][0], 0.0) if self.parts else None) def forced(self, r): p = self.parts[r][0] return any(self.R[T][r] < partial_term(T, p, 0) for T in WIDTHS) def slack(self, r, blind_room=False): p = self.parts[r][0] best = None for T in WIDTHS: tp, tm = partial_term(T, p, 0), partial_term(T, p, 1) if tp - tm < DELTA_MIN: continue room = self.RK[T] if blind_room else self.R[T][r] s = (room - tp) / (tp - tm) best = s if best is None else min(best, s) return best # --------------------------------------------------------------------------- # controls on the formula def control_transform(): worst = 0.0 for t in (0.0, 1.0, float(pi), 5.0, 10.0): num = 2 * mpmath.quad(lambda x: f(x) * mpmath.cos(t * x), [0, 1]) worst = max(worst, abs(float(num - fhat(t)))) check("C1 f^ closed form against the numerical transform", worst < 1e-8, "worst %.1e" % worst) def control_zeta(): zs = [float(mpmath.im(mpmath.zetazero(n))) for n in range(1, 101)] worst = 0.0 for T in (2, 4, 6): lhs = 2 * sum(T * float(fhat(g * T)) for g in zs) rhs = pole(T) + gamma_piece(T, 0) - sum(W[T].values()) print(" zeta, T = %d: zeros %.8f, formula %.8f" % (T, lhs, rhs)) worst = max(worst, abs(lhs - rhs)) check("C2 the formula on zeta against its first 100 zeros", worst < 1e-5, "worst %.1e" % worst) def lchi3(t): s = mpmath.mpf(1) / 2 + 1j * mpmath.mpf(t) return mpmath.re((mpmath.mpf(3) / mpmath.pi) ** (s / 2) * mpmath.gamma((s + 1) / 2) * mpmath.dirichlet(s, [0, 1, -1])) def control_chi3(tmax=60.0, step=0.05): zs = [] a, va = 0.01, lchi3(0.01) while a < tmax: b = a + step vb = lchi3(b) if va * vb < 0: lo, hi, vlo = a, b, va for _ in range(34): mid = (lo + hi) / 2 vm = lchi3(mid) if vm * vlo < 0: hi = mid else: lo, vlo = mid, vm zs.append((lo + hi) / 2) a, va = b, vb chi = {1: 1, 2: -1, 0: 0} worst = 0.0 for T in (2, 4, 6): lhs = 2 * sum(T * float(fhat(g * T)) for g in zs) rhs = log(3) + gamma_piece(T, 1) - sum( w * chi[p % 3] ** k for (p, k), w in W[T].items()) print(" L(s, chi_-3), T = %d: %d zeros %.8f, formula %.8f" % (T, len(zs), lhs, rhs)) worst = max(worst, abs(lhs - rhs)) check("C3 the formula on L(s, chi_-3) against its zeros to t = %g" % tmax, worst < 1e-4, "worst %.1e" % worst) # --------------------------------------------------------------------------- # statistics def share(v): n = len(v) if n == 0: return float("nan"), float("nan"), 0 s = sum(v) / n return s, sqrt(s * (1 - s) / n), n def halves(rooms, side, key): """D over the 2-types (KILL): lower half of key minus upper, pooled.""" by = {} for r, o in zip(rooms, side): k = key(r) if k is not None: by.setdefault(r.two, []).append((k, o)) num = den = var = 0.0 rows = [] for two in sorted(by): g = sorted(by[two]) n = len(g) if n < MIN_TYPE: rows.append((two, n, None)) continue m = n // 2 lo, hi = share([o for _, o in g[:m]]), share([o for _, o in g[n - m:]]) d, se = lo[0] - hi[0], sqrt(lo[1] ** 2 + hi[1] ** 2) rows.append((two, n, (lo[0], hi[0], d, se))) num += n * d var += (n * se) ** 2 den += n if den == 0: # no 2-type of MIN_TYPE fields return float("nan"), float("nan"), rows return num / den, sqrt(var) / den, rows def quintiles(name, rooms, side, key): pairs = sorted((key(r), o) for r, o in zip(rooms, side) if key(r) is not None) n = len(pairs) print(" %s, by quintile (%d fields):" % (name, n)) for i in range(5): g = pairs[i * n // 5:(i + 1) * n // 5] s, se, m = share([o for _, o in g]) print(" %.3f to %.3f %.3f +- %.3f (%d)" % (g[0][0], g[-1][0], s, se, m)) # --------------------------------------------------------------------------- # the wait: primes in order of size, each OUTSIDE when a degree-1 place # over it lies outside 2Cl; unramified primes with a degree-1 place only NULL = {"p": 0.5, "s": 0.75} def sequence(rec): """[p, type, outside, places in 2Cl, places] for the split and partial primes, in order.""" out = [] for p, t, v in rec["primes"]: if t == "s": out.append([p, "s", int(v > 0), 3 - v, 3]) elif t == "p": out.append([p, "p", v, 1 - v, 1]) return out def planted(seq, rng, before, after=None): """A planted side: outside with chance before(p, type) until the first outside prime and after(p, type) from then on (after None: before).""" new, paid = [], False for p, t, _, _, _ in seq: o = int(rng.random() < (after if paid and after else before)(p, t)) paid = paid or o n_in = (1 - o) if t == "p" else (1 if o else 3) new.append([p, t, o, n_in, 1 if t == "p" else 3]) return new def cluster(ys, ns): """Share sum(ys)/sum(ns) and its standard error over fields.""" N = sum(ns) if N == 0: return float("nan"), float("nan"), 0 s = sum(ys) / N return s, sqrt(sum((y - s * n) ** 2 for y, n in zip(ys, ns))) / N, N def wait(seqs, typ, rmax=6): """Prime r of type typ outside, by history: rows[r][j] = the share over the fields whose primes 1..r-1 hold j outside (j = 0, 1).""" rows = {} for seq in seqs: for r, (p, t, o, _, _) in enumerate(seq[:rmax]): j = sum(x[2] for x in seq[:r]) if t == typ and j <= 1: rows.setdefault((r + 1, j), []).append(o) return rows def memoryless(seqs): """The MEMORYLESS kill on partial primes: the all-in share pooled over r = 2..4 minus r = 1, the pool's error over fields.""" first = [s[0][2] for s in seqs if s and s[0][1] == "p"] ys, ns = [], [] for seq in seqs: y = n = 0 for r, (p, t, o, _, _) in enumerate(seq[:4]): if r == 0 or any(x[2] for x in seq[:r]): continue if t == "p": y, n = y + o, n + 1 if n: ys.append(y) ns.append(n) a, b = share(first), cluster(ys, ns) return a, b, b[0] - a[0], sqrt(a[1] ** 2 + b[1] ** 2) OUT = {"s": 2, "p": 1, "r2": 1} # each type's side datum outside 2Cl def upper(rec, held, actual=False): """U_A(T): the headroom with the primes of A in 2Cl and every other class-dependent prime outside (actual: at its own side).""" logq, U = log(-rec["d"]), {} pr = {p: (t, v) for p, t, v in rec["primes"]} for T in WIDTHS: s = 0.0 for (p, k), w in W[T].items(): t, v = pr[p] if p in held: v = 0 elif not actual and t in OUT: v = OUT[t] s += w * coeff(t, v, k) U[T] = logq + G[T] - s return U def forbids(rec, held): return min(upper(rec, held).values()) < 0 def weld(seqs, pmin, pmax): """The principal level at h = 2, 2 (places in 2Cl) / places, with no earlier outside prime (L0) and with one (L1), within (p, type).""" st = {} for seq in seqs: paid = 0 for p, t, o, n_in, n in seq: if pmin <= p <= pmax: st.setdefault((p, t), ([], []))[paid].append(2 * n_in / n) paid = paid or o num = den = var = 0.0 lv = ([], []) for key in sorted(st): for h in (0, 1): lv[h].extend(st[key][h]) g0, g1 = st[key] if len(g0) < 3 or len(g1) < 3: continue m0, m1 = sum(g0) / len(g0), sum(g1) / len(g1) v0 = sum((x - m0) ** 2 for x in g0) / (len(g0) - 1) / len(g0) v1 = sum((x - m1) ** 2 for x in g1) / (len(g1) - 1) / len(g1) n = len(g0) + len(g1) num += n * (m1 - m0) var += n * n * (v0 + v1) den += n L = [(sum(g) / len(g), len(g)) if g else (float("nan"), 0) for g in lv] if den == 0: return float("nan"), float("nan"), L return num / den, sqrt(var) / den, L # --------------------------------------------------------------------------- def main(): t0 = time.time() section("A THE FORMULA") with mpmath.workdps(20): # cubic.py's 60 digits hold elsewhere for T in WIDTHS: G[T] = 2 * gamma_piece(T, 0) + gamma_piece(T, 1) POLE[T] = pole(T) print(" T G(T) pole log 24000 + G") for T in WIDTHS: print(" %-4g %9.4f %9.4f %9.4f" % (T, G[T], POLE[T], log(24000) + G[T])) control_transform() control_zeta() control_chi3() section("B THE POPULATION") recs = population() seat = [r for r in recs if r["ok"] and "primes" in r] print(" fields %d, uncertified %d, 2-rank 1 %d" % (len(recs), sum(not r["ok"] for r in recs), len(seat))) bad = sum(r["bad"] for r in seat) unpl = sum(r["unplaced"] for r in seat) check("C4 [Q] = -2[P], 3[P] = 0, split parity, every place placed", bad == 0 and unpl == 0 and all(r["ok"] for r in recs), "bad %d, unplaced %d" % (bad, unpl)) rooms = [Room(r) for r in seat] rooms = [r for r in rooms if r.parts] side = [r.parts[0][1] for r in rooms] s, se, n = share(side) print(" fields with an odd partial prime <= %d: %d, position 1 " "outside 2Cl %.3f +- %.3f" % (PMAX, n, s, se)) section("C POSITIVITY") hmin = min(min(r.H.values()) for r in rooms) hkmin = min(min(r.HK.values()) for r in rooms) forced = [r.forced(0) for r in rooms] strong = sum(1 for fo, o in zip(forced, side) if fo and not o) print(" least H(T) over fields and widths %.4f; zeta_K's %.4f" % (hmin, hkmin)) for T in WIDTHS: hs = sorted(r.H[T] for r in rooms) print(" T = %-4g H median %.3f, lower decile %.3f" % (T, hs[len(hs) // 2], hs[len(hs) // 10])) check("C5 H >= 0 for L(s, rho) and zeta_K; no forced field reads +1", hmin >= -1e-9 and hkmin >= -1e-9 and strong == 0, "forced and in 2Cl: %d" % strong) odl = max(-(G[T] + POLE[T]) for T in WIDTHS) check("C6 zeta_K's bound with the primes dropped is at most log 23", odl <= log(23), "%.4f against %.4f" % (odl, log(23))) flipped = [Room(r, flip=True) for r in seat] neg = sum(1 for r in flipped if min(r.H.values()) < 0) check("C7 the partial sign flipped makes some H negative", neg > 0, "%d of %d fields" % (neg, len(flipped))) section("D THE ROOM AT POSITION 1") nf = sum(forced) print(" forced (GRH, full R): %d of %d = %.3f; outside among the " "unforced %.3f" % (nf, len(rooms), nf / len(rooms), share([o for fo, o in zip(forced, side) if not fo])[0])) quintiles("full slack s*", rooms, side, lambda r: r.slack(0)) quintiles("class-blind slack s_K*", rooms, side, lambda r: r.slack(0, True)) xs, ys = [], [] for r, o in zip(rooms, side): if r.RK[3] is not None and r.parts[0][0] < exp(3): xs.append(r.R[3][0] - r.RK[3]) ys.append(o) mx, my = sum(xs) / len(xs), sum(ys) / len(ys) cov = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) cor = cov / sqrt(sum((x - mx) ** 2 for x in xs) * sum((y - my) ** 2 for y in ys)) print(" corr(R - R_K, outside) at T = 3, p0 < e^3: %.3f over %d" % (cor, len(xs))) nk = sum(1 for r in rooms if r.slack(0, True) is None) print(" fields with no width where Delta >= %.2f: %d" % (DELTA_MIN, nk)) section("E THE KILL") rng = random.Random(SEED) plant = [] D, se, rows = halves(rooms, side, lambda r: r.slack(0, True)) for two, n, v in rows: if v is None: print(" 2-type %-3s %4d fields, dropped" % (two, n)) else: print(" 2-type %-3s %4d fields: lower half %.3f, upper %.3f," " D %.3f +- %.3f" % ((two, n) + v)) # C8: the planted side, by each type's halves of s_K* by = {} for i, r in enumerate(rooms): k = r.slack(0, True) if k is not None: by.setdefault(r.two, []).append((k, i)) psid = list(side) for two, g in by.items(): g.sort() m = len(g) // 2 for j, (_, i) in enumerate(g): psid[i] = int(rng.random() < (0.95 if j < m else 0.6)) Dp, sep, _ = halves(rooms, psid, lambda r: r.slack(0, True)) check("C8 the planted side reads D within 3 se of 0.35", abs(Dp - 0.35) <= 3 * sep, "%.3f +- %.3f" % (Dp, sep)) print(" D = %.3f +- %.3f (%.1f se)" % (D, se, D / se)) alive = D > 3 * se and D > 0.05 print(" KILL %s: the class-blind room %s the lean" % ( "MISSED" if alive else "FIRED", "tracks" if alive else "does not track")) print(" added after the first run, prints: D at each fixed width") for T in WIDTHS: def at(r, T=T): p0 = r.parts[0][0] tp, tm = partial_term(T, p0, 0), partial_term(T, p0, 1) if tp - tm < DELTA_MIN: return None return (r.RK[T] - tp) / (tp - tm) Dt, st, _ = halves(rooms, side, at) print(" T = %-4g D %.3f +- %.3f" % (T, Dt, st)) section("F POSITIONS") for r_ in range(POSITIONS): g = [r for r in rooms if len(r.parts) > r_] sd = share([r.parts[r_][1] for r in g]) fo = [r.forced(r_) for r in g] print(" position %d: outside %.3f +- %.3f (%d), forced %.3f, " "outside among unforced %.3f" % ( r_ + 1, sd[0], sd[1], sd[2], sum(fo) / len(g), share([r.parts[r_][1] for r, x in zip(g, fo) if not x])[0])) print(" added after the first run, a print: position 2 by position 1") for o1 in (0, 1): g = [r.parts[1][1] for r in rooms if len(r.parts) > 1 and r.parts[0][1] == o1] sd = share(g) print(" position 1 %s: position 2 outside %.3f +- %.3f (%d)" % ("outside" if o1 else "in 2Cl", sd[0], sd[1], sd[2])) section("G THE WAIT") print(" W1 odd partial places by position, outside given j earlier " "outside:") parts = [[v for p, t, v in r["primes"] if t == "p" and p > 2][:POSITIONS] for r in seat] for r_ in range(POSITIONS): cells = [] for j in range(r_ + 1): sd = share([x[r_] for x in parts if len(x) > r_ and sum(x[:r_]) == j]) cells.append("j=%d %.3f +- %.3f (%d)" % ((j,) + sd)) print(" r = %d: %s" % (r_ + 1, ", ".join(cells))) seqs = [sequence(r) for r in seat] for typ, name in (("p", "partial"), ("s", "split")): rows = wait(seqs, typ) print(" W2 %s primes (null %.2f), prime r outside given j earlier " "outside primes:" % (name, NULL[typ])) for r_ in range(1, 7): cells = ["j=%d %.3f +- %.3f (%d)" % ((j,) + share(rows[(r_, j)])) for j in (0, 1) if rows.get((r_, j))] print(" r = %d: %s" % (r_, ", ".join(cells))) h2 = [r for r in seat if prod(r["inv"]) == 2] small = [sequence(r) for r in h2 if -r["d"] <= 12000] print(" W3 THE WELD at h = 2: %d fields with |d_K| <= 12000, %d in all" % (len(small), len(h2))) for lo, hi in ((2, 7), (11, 89)): g, se_, L = weld(small, lo, hi) print(" %d <= p <= %d: L0 %.3f (%d primes), L1 %.3f (%d), " "L1 - L0 within (p, type) %.3f +- %.3f" % (lo, hi, L[0][0], L[0][1], L[1][0], L[1][1], g, se_)) g, se_, L = weld([sequence(r) for r in h2], 2, 7) print(" |d_K| <= %d, 2 <= p <= 7: L0 %.3f (%d), L1 %.3f (%d), " "L1 - L0 %.3f +- %.3f" % (CAP, L[0][0], L[0][1], L[1][0], L[1][1], g, se_)) G7, S7, _ = weld(small, 2, 7) rng = random.Random(SEED + 1) mem = [planted(s, rng, lambda p, t: 0.83 if t == "p" else 0.93, lambda p, t: NULL[t]) for s in seqs] msm = [planted(s, rng, lambda p, t: 0.83 if t == "p" else 0.93, lambda p, t: NULL[t]) for s in small] gm, sm, _ = weld(msm, 2, 7) _, _, dmm, smm = memoryless(mem) check("C9 the memory plant: weld gap > 0.3 and 3 se, memoryless holds", gm > 0.3 and gm > 3 * sm and abs(dmm) <= 3 * smm, "gap %.3f +- %.3f, diff %.3f +- %.3f" % (gm, sm, dmm, smm)) obs = {} for s in small: for p, t, o, _, _ in s: obs.setdefault((p, t), []).append(o) fx = [planted(s, rng, lambda p, t: sum(obs[(p, t)]) / len(obs[(p, t)])) for s in small] gf, sf, _ = weld(fx, 2, 7) check("C10 the fixed plant: weld gap within 3 se of 0", abs(gf) <= 3 * sf, "%.3f +- %.3f" % (gf, sf)) a, b, dm, sem = memoryless(seqs) print(" MEMORYLESS: partial all-in share r = 1 %.3f +- %.3f (%d), " "r = 2..4 %.3f +- %.3f (%d), diff %.3f +- %.3f (%.1f se)" % (a + b + (dm, sem, dm / sem))) print(" KILL MEMORYLESS %s" % ("FIRED" if abs(dm) > 3 * sem else "MISSED")) print(" KILL WELD %s: L1 - L0 at p <= 7 = %.3f +- %.3f" % ( "MISSED" if G7 > 0.3 and G7 > 3 * S7 else "FIRED", G7, S7)) section("H THE WORST-CASE ROOM") gap = err = 0.0 wc1 = wc2 = clash = 0 free = [] for rec in seat: H = Room(rec).H U0, Ua = upper(rec, set()), upper(rec, set(), actual=True) gap = min(gap, min(U0[T] - H[T] for T in WIDTHS)) err = max(err, max(abs(Ua[T] - H[T]) for T in WIDTHS)) seq = sequence(rec) p1, p2 = seq[0], seq[1] f1 = forbids(rec, {p1[0]}) f2 = forbids(rec, {p1[0], p2[0]}) wc1 += f1 wc2 += f2 clash += (f1 and not p1[2]) + (f2 and not p1[2] and not p2[2]) if not f2: free.append((rec["d"], p1[:3], p2[:3])) check("C11 U_empty >= H everywhere; no actual pattern forbidden", gap >= -1e-9 and clash == 0, "least U - H %.2e, forbidden yet seen %d" % (gap, clash)) check("C12 U at the actual sides reproduces H", err < 1e-9, "worst %.1e" % err) n = len(seat) print(" WC1 first prime in 2Cl forbidden: %d of %d = %.3f" % (wc1, n, wc1 / n)) print(" WC2 first two in 2Cl forbidden: %d of %d = %.3f" % (wc2, n, wc2 / n)) for d, a, b in free: print(" free: d %d, first %s, second %s" % (d, a, b)) least = min(-d for d, _, _ in free) if free else CAP print(" added after the first run, prints: %d fields with |d_K| < %d," " every one forbidden" % (sum(-r["d"] < least for r in seat), least)) sm2 = [r for r in h2 if -r["d"] <= 12000] s1 = [forbids(r, {sequence(r)[0][0]}) for r in sm2] print(" h = 2, |d_K| <= 12000: first prime principal forbidden at " "%d of %d" % (sum(s1), len(sm2))) print(" KILL CONDUCTOR %s" % ("MISSED" if wc2 >= 0.95 * n else "FIRED")) print(f"\n{sum(RESULTS)}/{len(RESULTS)} controls PASS, " f"{time.time() - t0:.1f} s") raise SystemExit(0 if all(RESULTS) else 1) if __name__ == "__main__": main()