room.py
Checks: TRIPLE: The least partial place leans non-square observation; rule.
Cited on: TRIPLE.
To run it: save room.py, with cubic.py beside it; needs mpmath.
"""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()