lossdensity.py
Checks: THERMAL: The F₂ place loses at about 1 in 200 cubic fields theorem; rule.
Cited on: THERMAL.
To run it: save lossdensity.py, with cubic.py, cubrace.py, cubtop.py, fieldtop.py, growth.py, memory.py, quadtop.py, realrace.py, spectrum.py beside it; needs mpmath, numpy, scipy.
"""lossdensity.py -- the density of the cubic fields where the F_2 place
loses the top: one number for both signatures, read as a probability
under the product of the local splitting laws.
QUESTION. realrace.py's (5) proves the F_2 place loses at a positive
proportion of the cubic fields of each signature, ordered by |d_K|, and
tops at another. What is the proportion where it loses?
THE ARGUMENT.
(1) THE MEASURE. Bhargava, Shankar and Tsimerman (On the
Davenport-Heilbronn theorems and second order terms, Invent. Math.
193 (2013)), Lemma 18: the forms maximal at p of type 111, 12, 3,
1^2 1, 1^3 have p-adic density (p - 1)^2 (p + 1) p^-4 times p/6,
p/2, p/3, 1, 1/p, the proof uniform in p, 2 and 3 included; by
their Theorem 8 the number of cubic fields of a signature with
|d_K| < X and finitely many prescribed types is c X + o(X), c the
signature's constant times a product over p of local sums of
1/(Disc_p |Aut|), each proportional at p, by one factor for every
type, to these densities. So the types at distinct primes are
independent in the limit, with weights 1/6, 1/2, 1/3, 1/p, 1/p^2
over 1 + 1/p + 1/p^2 at p, the same in both signatures. At 2 the
types with an F_2 place weigh (1/6 + 1/2 + 1/2 + 1/4) / (7/4) =
17/21.
(2) EXISTENCE. The race of (3) by exact signs (no float guard, b* and
b0 any points with the stated signs) at a cut Y >= 10^6: zeta_K's
product, the column of norm 3's places of norm 2 3^a + 1 and the
F_2 column's places of Fermat norm all read below Y, and each part
past Y bounded: the tail at m = 0 and 3, missed_mass's sum taken
past Y, and exp(6 Y^-b) (the Fermat norms past Y are F_n = 2^(2^n)
+ 1, F_(n+1) - 1 = (F_n - 1)^2, so at b > 1 each N^-b is at most
half the one before, and they have at most three places each: the
places' sum of N^-b is below 6 Y^-b). Let LOSS_Y, TOP_Y be the
type sequences it decides. Each is a finite union of sets of
finitely many prescribed types, so has a density by (1), and the
lower and upper densities of the losing fields lie in [d(LOSS_Y),
d(LOSS_Y) + d(OPEN_Y)]. Every bound past Y tends to 0, so at a
sequence with beta_F2 != beta_3 the four ends close on the two
points uniformly on compacts, each Psi strictly decreasing, and
the race decides at every large Y, while a tie is OPEN at every Y:
by dominated convergence under the product measure, d(OPEN_Y) =
P(OPEN_Y) tends to the tie's mass. Fix the types at 2, 3, 5, 17,
257, 65537, at every Fermat prime past them and at every prime 2
3^a + 1, the set S the two sigmas read: then sigma_2 - sigma_3 is
a fixed Dirichlet series, analytic for b > 0 and not 0 (at every
failing configuration sigma_2 holds places of norm 5, 5 being of
type 111 or 1^2 1 there (cubtop.py's FINDINGS), and sigma_3 none,
since 4 | 5 - 1, and no product of norms 2, 3, 4 and 2 3^a + 1 is
5), so its zeros beta_0 in (1, oo) are countably many. A tie needs
zeta_K(beta_0) = 1 + sigma_2(beta_0) at one. log zeta_K(beta_0) is
the fixed part plus a sum over p not in S of independent laws on
five distinct atoms (a3 < a1 < a1 + a2 < 2 a1 < 3 a1, a_f = -log(1
- p^-f beta_0)), the likeliest of weight below 1/2. By Levy's
theorem (P. Levy, Sur les series dont les termes sont des
variables eventuelles independantes, Studia Math. 3 (1931),
119-153; stated as Theorem 12(a) of N. Bouleau, The lent
particle method, application to multiple Poisson integrals,
arXiv:1004.2819) an a.s. convergent sum of independent real
laws X_n is continuous when sum P(X_n != a_n) diverges for every
choice of constants a_n; here each term is above 1/2. So each
zero is hit with probability 0, the tie has mass 0, and the
density exists and is P(LOSS) under (1)'s product measure. The
engine below is the cut Y = 10^6, with the tighter Fermat bound
exp(6 2^-32b) (no Fermat number lies between 65537 and 2^32 + 1)
and the float tolerances of (3) and DESIGN, reached through the
cut of (5).
(3) THE RACE BY SIGNS. With Psi = zeta_K - sigma at a column's true
zeta_K and sigma, decreasing, its point where Psi = 1
(realrace.py's (6)). The four ends: F2 upper h = zeta_K(b, 3) -
sigma_2 - 1 >= Psi_2 - 1; F2 lower, zeta_K(b, 0) and sigma_2
times exp(6 2^-32b), the dropped Fermat places, <= Psi_2 - 1;
N=3 lower, zeta_K(b, 0) and sigma_3 times exp(missed_mass(b)),
<= Psi_3 - 1; N=3 upper, zeta_K(b, 3) and sigma_3, >= Psi_3 - 1.
zeta_K(b, m) is the product below the cut times the tail of
zeta's product to the m-th power, as cubrace.py's. LOSS where at
b* with the F2 upper h(b*) <= 0 the N=3 lower h(b*) > 0: there
Psi_2 <= 1 < Psi_3, so beta_F2 <= b* < beta_3. TOP where at b0
with the F2 lower h(b0) >= 0 the N=3 upper h(b0) < 0. Else OPEN,
read also at m = 1. A norm-3 end's |h| below 1e-9 counts as
OPEN.
(4) ONLY THE FAILING CONFIGURATIONS RACE. cubtop.py's bound decides
every column of every configuration but the column of norm 3 at
675 (its R3), all with 2 of type 1^2 1 or 1^3. So P(LOSS) = the
exact mass of the 675 times P(LOSS | failing), and only the
latter is drawn.
(5) THE CUT AT 10^3 IS EXACT. (3)'s ends at a cut Y < 10^6, zeta_K's
product below Y and the tail past it at m = 0 and 3, the sigmas
read whole as before: each local factor past Y lies in [1, (1 -
p^-b)^-3], so the F_2 upper end's h is still at least Psi_2 - 1
and the N=3 lower end's at most Psi_3 - 1, and a verdict at Y is
a verdict on the true points, as one at 10^6 is. Two decisions
therefore agree (but by the float tolerance), and the verdict at
Y where it decides, else at 10^6, is (3)'s at 10^6 wherever that
one decides. The primes past Y that the sigmas do not read are
independent of the rest, so they are drawn only where the cut is
OPEN.
DESIGN. The 675 from cubtop.py's sweep over the 6250 configurations
with 2 of type 1^2 1 or 1^3 (its R3 finds none elsewhere), their mass
as an exact fraction. A draw: a failing configuration in proportion to
its weight, then from (1)'s weights the types FIRST holds (every p <
10^3, and the primes past it the sigmas read: 1459, 39367, 65537),
raced at 10^3; where OPEN, every other p < 10^6, raced at 10^6 (5).
The verdict by (3), b* and b0 the ends of a root found by Brent's
method stepped 1e-10 outward and checked. The share decided at 10^3
is printed. Density = mass L/n, bracket
mass (L + O)/n, spread mass sqrt(q (1 - q)/n), q = L/n; among the
fields with an F_2 place, over 17/21.
CONTROLS, read before any verdict.
C1 (the race, answers known): -1420, 1940, -2764, 2804 and -2284
through the engine from their places (realrace.py's places): its
four ends equal cubrace.py's Race to 1e-8 and the verdicts LOSS,
LOSS, LOSS, LOSS, TOP; one random failing draw turned into places
and read by Race, the same. PLANT: sigma_3 without its held
places, C1 refuses at -1420.
C2 (the sampler): 20000 unconditional draws through draw_rows, the
run's row sampler, and its thresholds CUM, at 2, 3, 5, 7, 17, 257, 65537;
the type frequencies at 2, 3, 5, 7 within 4 sd of (1)'s weights,
the F_2 share printed beside 17/21, and the share in a failing
configuration within 4 sd of the exact mass; 20000 picks of
pick_config, the run's choice of configuration, with 2 of type 1^3
within 4 sd of that type's weight. PLANT: 1/p and 1/p^2 swapped at
2, C2 refuses.
C3 (the cut, (5)): 500 draws raced at 10^3 and, their types filled,
whole at 10^6: no draw both decide differently. PLANT: the cut's
F_2 upper end at m = 0 and N=3 lower end at m = 3, each bounding
the wrong side, C3 sees a disagreement. The controls draw from
their own stream (seed + 1).
PRINT (not a control): the types at 2, 3, 5 among the 515 cubic
fields with |d_K| <= 3000 against (1)'s weights; the second order
term is X^(5/6), 3000^(-1/6) = 0.26 relative, so no spread is
promised; and the count of them with an F_2 place, against which
the four losses the cubic races find (cubrace.py, realrace.py)
are read.
PREDICTIONS, fixed before code. (d1) failing mass
in [0.030, 0.045]. (d2) P(LOSS | failing) in [0.15, 0.50]. (d3) OPEN
below 3% of failing draws. (d4) P(LOSS | failing, 1^3) > P(LOSS |
failing, 1^2 1). (d5) at Y = 10^3, 30% to 80% of failing draws decided.
(d6) the density among fields with an F_2 place in [0.006, 0.020], the
fields' 4 of 388 within Poisson's two-sided 95% at 388 times it.
KILL, as printed: the "LOSS | failing" line; (d2) dies outside
[0.15, 0.50], and the hunch that the density rests on few primes dies
if the Y = 10^3 line prints decided below 30%. For the cut, fixed
before its code on the 10^4 draws' numbers: (s1) a draw costs 2 to 6
ms; (s2) decided at 10^3 in [0.912, 0.930]; (s3) LOSS | failing within
0.008 of the 10^4 draws' 0.1477. n = 34000 puts the spread at most
0.002 to q = 0.17.
FINDINGS. C0: 675 failing configurations, 250 at 1^2 1 and 425 at
1^3, their mass 86570/2598141 = 0.033320. C1: at -1420, 1940, -2764,
2804, -2284 the four ends equal Race's within 3.4e-10, verdicts LOSS,
LOSS, LOSS, LOSS, TOP; a random draw (2 = 1^2 1) within 3.3e-10, LOSS
both. C1p: the ends move by 4.5e-06. C2: largest z 2.47 (at 2), the
F_2 share 0.8085 beside 17/21 = 0.8095, the failing share 0.03385
beside 0.03332 (z 0.42), pick_config's 1^3 share 0.5135 beside 0.5124
(z 0.31). C2p: z 57.45 at 2, 6.33 on the failing share. C3: 463 of 500
decided at 10^3, no disagreement; C3p: 37 verdicts move, 24 against
the full race. PRINT, the
515 fields to 3000 against the weights: at 2, 0.041 0.280 0.247 0.256
0.177 against 0.095 0.286 0.190 0.286 0.143; at 3, 0.056 0.346 0.307
0.196 0.095 against 0.115 0.346 0.231 0.231 0.077: 111 at less than
half its weight at 2 and at 3, 1^2 1 under, 3 and 1^3 over, 12 on it.
RACE over 34000 failing draws: LOSS 5011, TOP 28971, OPEN 18 (at m =
1, 11 LOSS and 7 TOP). LOSS | failing 0.1474 +- 0.0019, bracket
[0.1474, 0.1479]; at 2 = 1^2 1, 2457 of 16533, 0.1486 +- 0.0028; at
1^3, 2554 of 17467, 0.1462 +- 0.0027.
DENSITY LOSS 0.00491 +- 0.00006 among all cubic fields, either
signature, bracket [0.00491, 0.00493]; 0.00607 +- 0.00008 among those
with an F_2 place. The fields to 3000, 4 of 388 with an F_2 place,
against 388 x 0.00607 = 2.35: P(X <= 4) 0.910, P(X >= 4) 0.212.
Y = 10^3: 31356 of 34000 decided, 0.922.
(d1) holds. (d2) dies: 0.1474 against the bound 0.15, 1.4 spreads
below it; the band was set by the four losses among the twelve fields
to 3000 in a failing configuration. (d3) holds: 18 OPEN in 34000. (d4)
fails: 1^2 1 at 0.1486, 1^3 at 0.1462, 0.6 combined spreads apart, the
other way from the 10^4 draws. (d5) fails upward: 92% decided, not 30%
to 80%, so the hunch that the density rests on few primes survives its
kill. (d6) holds at its edge: 0.00607 against the floor 0.006; it is
(d2)'s statistic times a constant, so not independent support. (s1)
holds, 2 ms a draw; (s2) holds, 0.922; (s3) holds, 0.0003 from 0.1477.
READING. The density exists by (2) and is about 1 in 200 cubic
fields, 1 in 165 of those with an F_2 place, in each signature: the
F_2 place loses at about 15% of the fields in a failing configuration,
the two ramified types of 2 alike within the spread. The fields to
3000 lose at 4 of 388, within chance of the limit, though their types
at 2 and 3 are not at the limit's weights. The Euler product below
10^3, the sigmas read whole, decides 92% of failing draws.
RUN RECORD. `python lossdensity.py` (n = 34000, seed 1, the controls
seed 2): 7 of 7 checks, 166.0 s (the race 81.0 s, 2 ms a draw), peak
commit 286.1 MB; its first record, C3p's plant then the F_2 upper end
at m = 0 alone, printed 13 verdicts moved, none against the full race,
the same race and 301.3 MB. Before the cut, every draw raced at
10^6 (26 ms a draw): 10^4 draws printed LOSS 1477 (0.1477 +- 0.0035),
the 1^3 rate above the 1^2 1 rate by 1.45 spreads, and 9209 decided at
10^3 with none against the full race. With the fields' count read off
the enumeration (388 of the 515 with an F_2 place, as before): 7 of 7,
168.7 s, peak commit 287.1 MB, every verdict line as above.
"""
import argparse
import os
import sys
import time
from fractions import Fraction
from itertools import product
from math import exp, log, sqrt
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import numpy as np # noqa: E402
from scipy.optimize import brentq # noqa: E402
from scipy.special import zeta # noqa: E402
from scipy.stats import poisson # noqa: E402
import cubrace # noqa: E402
import cubtop # noqa: E402
CODES = "sqirt" # 111, 12, 3, 1^2 1, 1^3
N1 = np.array([3, 1, 0, 2, 1], dtype=float) # places of norm p
N2 = np.array([0, 1, 0, 0, 0], dtype=float) # of norm p^2
N3 = np.array([0, 0, 1, 0, 0], dtype=float) # of norm p^3
PR = np.array(cubrace.PRIMES)
LOGP = np.log(PR.astype(float))
IDX = {int(p): i for i, p in enumerate(PR)}
REL = (2,) + tuple(cubtop.REL) # 2, 3, 5, 17, 257, 65537
HELD = [2 * 3 ** a + 1 for a in range(1, 13)
if 2 * 3 ** a + 1 < cubrace.X and int(2 * 3 ** a + 1) in IDX]
YCUT = 10 ** 3
NY = int(np.searchsorted(PR, YCUT))
# the first stage: every p below the cut, and past it the primes the
# sigmas read (the held 1459, 39367 and the Fermat 65537)
FIRST = np.array(sorted(set(range(NY)) | {IDX[p] for p in HELD + list(REL)}))
REST = np.setdiff1d(np.arange(len(PR)), FIRST)
CHECKS = []
def check(name, ok, detail=""):
CHECKS.append(bool(ok))
print(f" [{'PASS' if ok else 'FAIL'}] {name}" + (f" ({detail})"
if detail else ""))
def weights(p, plant=False):
w = [Fraction(1, 6), Fraction(1, 2), Fraction(1, 3), Fraction(1, p),
Fraction(1, p * p)]
if plant:
w[3], w[4] = w[4], w[3]
s = sum(w)
return [x / s for x in w]
W = np.array([[1 / 6, 1 / 2, 1 / 3, 1 / p, 1 / p ** 2] for p in
PR.astype(float)])
CUM = np.cumsum(W / W.sum(1, keepdims=True), axis=1)[:, :4]
def draw_rows(rng, cum, n):
"""n draws of a type index at each row's prime, cum its thresholds."""
return (rng.random((n, len(cum)))[:, :, None] > cum[None]).sum(2)
def draw_all(rng):
"""A type index at every p < 10^6 from (1)'s weights."""
return draw_rows(rng, CUM, 1)[0]
def pick_config(rng, cw):
return rng.choice(len(cw), p=cw)
def failing_configs():
"""cubtop.py's 675: {(t2, config)}."""
G = cubtop.Grid(1.0001 + cubtop.STEP * np.arange(29999))
top = cubtop.top_index(G, 3)
profs = [((t2, "".join(ty)), cubtop.cubic_profile(t2, ty))
for t2 in "rt" for ty in product(CODES, repeat=5)]
bad, _, _ = cubtop.sweep(G, profs, top)
return sorted({b[0] for b in bad})
def config_weight(key, plant=False):
t2, ty = key
w = weights(2, plant)[CODES.index(t2)]
for p, t in zip(cubtop.REL, ty):
w *= weights(p)[CODES.index(t)]
return w
def logtail(b, ls):
"""log of zeta(b) prod_(p below the cut) (1 - p^-b), ls the sum of
log(1 - p^-b) there."""
return log(zeta(b)) + ls
class Draw:
"""A cubic field's types at every p < 10^6, read by (3)."""
def __init__(self, types, plant=False):
self.set_types(types)
self.plant_cut = False # C3's plant
t = {p: CODES[self.types[IDX[p]]] for p in REL}
self.t2 = t[2]
n1 = lambda p: int(N1[self.types[IDX[p]]]) # noqa: E731
self.fermat = [] # norms of the Fermat places
for p in REL[1:]:
self.fermat += [p] * n1(p)
self.fermat += [9] * int(N2[self.types[IDX[3]]])
self.n3 = n1(3)
self.held = [] if plant else [M for M in HELD for _ in range(n1(M))]
self.g = {"s": 3, "q": 1, "r": 2, "t": 1, "i": 0}[self.t2]
self.twod = {"s": [3, 3, 3], "q": [3], "r": [2, 3], "t": [2],
"i": []}[self.t2]
self.four = self.t2 == "q"
def set_types(self, types):
self.types = np.asarray(types)
self.c1, self.c2, self.c3 = (N1[self.types], N2[self.types],
N3[self.types])
def fill(self, rng):
"""The second stage: a type at every prime outside FIRST."""
t = self.types.copy()
t[REST] = draw_rows(rng, CUM[REST], 1)[0]
self.set_types(t)
def logz(self, b, cut=None):
"""(log of the product below the cut, sum of log(1 - p^-b))."""
n = cut or len(PR)
x = np.exp(-b * LOGP[:n])
l1 = -np.log1p(-x)
lz = (self.c1[:n] @ l1 + self.c2[:n] @ -np.log1p(-x * x)
+ self.c3[:n] @ -np.log1p(-x * x * x))
return lz, -l1.sum()
def zk(self, b, m, cut=None):
lz, ls = self.logz(b, cut)
return exp(lz + m * logtail(b, ls))
def sigma2(self, b, lower=False):
x = 2.0 ** -b
s = (1 - x) ** -(self.g - 1)
for N in self.fermat:
s *= 1 + N ** -b
return s * exp(6 * 2.0 ** (-32 * b)) if lower else s
def sigma3(self, b, lower=False):
x = 2.0 ** -b
s = 1.0
for d in self.twod:
s *= sum(x ** j for j in range(d + 1))
if self.four:
s *= 1 + 4.0 ** -b + 16.0 ** -b
s /= (1 - 3.0 ** -b) ** (self.n3 - 1)
for M in self.held:
s *= 1 + M ** -b
return s * exp(cubrace.missed_mass(b)) if lower else s
def h(self, end, b, cut=None):
m = {"F2hi": 3, "F2lo": 0, "N3lo": 0, "N3hi": 3, "F2c": 1,
"N3c": 1}[end]
if cut and self.plant_cut and end in ("F2hi", "N3lo"):
m = 3 - m
if end[:2] == "F2":
sig = self.sigma2(b, lower=end == "F2lo")
else:
sig = self.sigma3(b, lower=end == "N3lo")
return self.zk(b, m, cut) - sig - 1
def root(self, end, cut=None, lo=1.05, hi=2.76):
f = lambda b: self.h(end, b, cut) # noqa: E731
assert f(lo) > 0 > f(hi), (end, lo, hi)
return brentq(f, lo, hi, xtol=1e-12, rtol=1e-14)
def verdict(self, cut=None):
"""LOSS, TOP or OPEN by (3); OPEN also read at m = 1."""
r = self.root("F2hi", cut)
bs = r + 1e-10
assert self.h("F2hi", bs, cut) <= 0
if self.h("N3lo", bs, cut) > 1e-9:
return "LOSS", None
r0 = self.root("F2lo", cut, lo=1.05, hi=r + 1e-9)
b0 = r0 - 1e-10
assert self.h("F2lo", b0, cut) >= 0
if self.h("N3hi", b0, cut) < -1e-9:
return "TOP", None
c = self.root("F2c", cut)
return "OPEN", ("LOSS" if self.h("N3c", c, cut) > 0 else "TOP")
def ends(self):
return (self.root("F2lo"), self.root("F2hi"), self.root("N3lo"),
self.root("N3hi"))
def places(self):
return [(int(p), fd, e) for p, t in zip(PR, self.types)
for fd, e in PLACES[CODES[t]]]
PLACES = {"s": [(1, 1)] * 3, "q": [(1, 1), (2, 1)], "i": [(3, 1)],
"r": [(1, 2), (1, 1)], "t": [(1, 3)]}
def from_places(pl):
by = {}
for p, fd, e in pl:
by.setdefault(p, []).append((fd, e))
key = {tuple(sorted(v)): c for c, v in PLACES.items()}
return np.array([CODES.index(key[tuple(sorted(by[int(p)]))])
for p in PR])
def race_ends(pl):
r = cubrace.Race(pl, 3).read()
return (r["F2"][0], r["F2"][2], r["N3"][0], r["N3"][2])
def control_race(rng, fails, cw):
import realrace
print("CONTROL C1: the engine against cubrace.py's Race")
known = {-1420: (cubrace.FIELDS[-1420], "LOSS"),
1940: ((-1, 0, 8, -2), "LOSS"),
-2764: (cubrace.FIELDS[-2764], "LOSS"),
2804: ((-1, -1, 9, -1), "LOSS"),
-2284: (cubrace.FIELDS[-2284], "TOP")}
ok, worst = True, 0.0
for dk, (form, want) in known.items():
pl = realrace.places(form, dk)
D = Draw(from_places(pl))
mine, theirs = D.ends(), race_ends(pl)
dev = max(abs(a - b) for a, b in zip(mine, theirs))
v, _ = D.verdict()
worst = max(worst, dev)
ok &= dev < 1e-8 and v == want
print(f" d_K {dk}: ends {' '.join(f'{x:.6f}' for x in mine)}, "
f"Race's within {dev:.1e}, {v} (want {want})")
if dk == -1420:
P = Draw(from_places(pl), plant=True)
pdev = max(abs(a - b) for a, b in zip(P.ends(), theirs))
pv, _ = P.verdict()
D = random_draw(rng, fails, cw)
mine, theirs = D.ends(), race_ends(D.places())
dev = max(abs(a - b) for a, b in zip(mine, theirs))
v, _ = D.verdict()
rv = ("LOSS" if theirs[2] > theirs[1] else "TOP" if theirs[3] <
theirs[0] else "OPEN")
print(f" a random failing draw (2 {D.t2}): ends within {dev:.1e}, "
f"{v}, Race's brackets {rv}")
check("C1 four ends equal Race's to 1e-8, verdicts LOSS x4, TOP, and "
"a random draw agrees", ok and dev < 1e-8 and v == rv,
f"worst {max(worst, dev):.1e}")
print(f" PLANT sigma_3 without held places at -1420: ends move by "
f"{pdev:.1e}, verdict {pv}")
check("C1p the plant is refused (an end moves past 1e-8)", pdev > 1e-8)
def control_sampler(rng, fails, cw, mass, plant=False):
tag = "PLANT 1/p and 1/p^2 swapped at 2" if plant else "C2"
print(f"CONTROL {tag}: 20000 unconditional draws through draw_rows")
n = 20000
ps = [2, 3, 5, 7, 17, 257, 65537]
cum = CUM[[IDX[p] for p in ps]].copy()
if plant:
cum[0] = np.cumsum([float(x) for x in weights(2, True)])[:4]
T = draw_rows(rng, cum, n)
worst = 0.0
for c, p in enumerate((2, 3, 5, 7)):
w = np.array([float(x) for x in weights(p)])
f = np.bincount(T[:, c], minlength=5) / n
z = np.abs(f - w) / np.sqrt(w * (1 - w) / n)
worst = max(worst, float(z.max()))
print(f" {p}: freq {' '.join(f'{x:.4f}' for x in f)} weights "
f"{' '.join(f'{x:.4f}' for x in w)} max z {z.max():.2f}")
f2 = float(np.isin(T[:, 0], [0, 1, 3, 4]).mean())
print(f" F2 share {f2:.4f} beside 17/21 = {17 / 21:.4f}")
fs = set(fails)
cols = [1, 2, 4, 5, 6] # 3, 5, 17, 257, 65537
hit = sum((CODES[T[i, 0]], "".join(CODES[T[i, c]] for c in cols))
in fs for i in range(n)) / n
zm = abs(hit - mass) / sqrt(mass * (1 - mass) / n)
print(f" in a failing configuration {hit:.5f} beside the mass "
f"{mass:.5f}, z {zm:.2f}")
ok = worst < 4 and zm < 4
if plant:
check("C2p the plant is refused (a z of 4 or more)", not ok)
return
tw = float(sum(cw[k] for k, key in enumerate(fails) if key[0] == "t"))
ks = [pick_config(rng, cw) for _ in range(n)]
ft = sum(fails[k][0] == "t" for k in ks) / n
zt = abs(ft - tw) / sqrt(tw * (1 - tw) / n)
print(f" pick_config: 2 of type 1^3 in {ft:.4f} of {n} picks beside "
f"its weight {tw:.4f}, z {zt:.2f}")
check("C2 frequencies within 4 sd, the failing share and the picks' "
"1^3 share within 4 sd", ok and zt < 4)
def print_fields():
import realrace
from cubic import enumerate_fields
print("PRINT: the types at 2, 3, 5 among the cubic fields to 3000")
fields = list(enumerate_fields(3000).items())
fields += list(realrace.enumerate_real(3000).items())
n = len(fields)
for p in (2, 3, 5):
c = np.zeros(5)
for form, dk in fields:
c[CODES.index(realrace.ptype(form, dk, p))] += 1
w = np.array([float(x) for x in weights(p)])
print(f" {p}: {n} fields, share {' '.join(f'{x:.3f}' for x in c / n)}"
f" weights {' '.join(f'{x:.3f}' for x in w)}")
nf2 = sum(N1[CODES.index(realrace.ptype(form, dk, 2))] > 0
for form, dk in fields)
print(f" {nf2} of the {n} with an F_2 place")
return nf2
def random_draw(rng, fails, cw, staged=False):
"""A failing configuration and the other types; staged, only those
in FIRST, the rest 0 until fill."""
k = pick_config(rng, cw)
t2, ty = fails[k]
if staged:
types = np.zeros(len(PR), dtype=int)
types[FIRST] = draw_rows(rng, CUM[FIRST], 1)[0]
else:
types = draw_all(rng)
for p, t in zip(REL, t2 + ty):
types[IDX[p]] = CODES.index(t)
return Draw(types)
def control_shortcut(rng, fails, cw, n=500):
"""C3: the cut's verdict against the full race on the same draws."""
print(f"CONTROL C3: {n} draws raced at {YCUT} and whole at 10^6")
bad = moved = pbad = dec = 0
for _ in range(n):
D = random_draw(rng, fails, cw, staged=True)
D.fill(rng)
v, _ = D.verdict()
vy, _ = D.verdict(cut=NY)
dec += vy != "OPEN"
bad += vy != "OPEN" and v != "OPEN" and vy != v
D.plant_cut = True
vp, _ = D.verdict(cut=NY)
moved += vp != vy
pbad += vp != "OPEN" and v != "OPEN" and vp != v
print(f" decided at the cut {dec} of {n}; disagreements {bad}")
check("C3 the cut never disagrees with the full race", bad == 0)
print(f" PLANT the cut's F_2 upper and N=3 lower tails swapped: {moved} "
f"verdicts move, {pbad} against the full race")
check("C3p the plant disagrees with the full race", pbad > 0)
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--n", type=int, default=34000)
ap.add_argument("--seed", type=int, default=1)
args = ap.parse_args()
rng = np.random.default_rng(args.seed)
crng = np.random.default_rng(args.seed + 1) # the controls
t0 = time.time()
fails = failing_configs()
mass = sum(config_weight(k) for k in fails)
print(f"the failing configurations: {len(fails)}, "
f"{sum(k[0] == 'r' for k in fails)} at 1^2 1, "
f"{sum(k[0] == 't' for k in fails)} at 1^3")
print(f" their mass {float(mass):.6f} = {mass.numerator}/"
f"{mass.denominator}"[:120])
cw = np.array([float(config_weight(k) / mass) for k in fails])
check("C0 cubtop.py's 675", len(fails) == 675)
control_race(crng, fails, cw)
control_sampler(crng, fails, cw, float(mass))
control_sampler(crng, fails, cw, float(mass), plant=True)
control_shortcut(crng, fails, cw)
nf2 = print_fields()
if not all(CHECKS):
print("a control failed: no verdict")
sys.exit(1)
print(f"\ncontrols read, {time.time() - t0:.1f} s; {args.n} draws")
tally = {"LOSS": 0, "TOP": 0, "OPEN": 0}
by2 = {"r": [0, 0], "t": [0, 0]}
central = {"LOSS": 0, "TOP": 0}
decY = 0
t1 = time.time()
for _ in range(args.n):
D = random_draw(rng, fails, cw, staged=True)
v, c = D.verdict(cut=NY)
if v == "OPEN":
D.fill(rng)
v, c = D.verdict()
else:
decY += 1
tally[v] += 1
by2[D.t2][0] += 1
by2[D.t2][1] += v == "LOSS"
if c:
central[c] += 1
dt = time.time() - t1
n = args.n
q = tally["LOSS"] / n
s = sqrt(q * (1 - q) / n)
m = float(mass)
print(f"RACE over {n} failing draws: {tally}, OPEN read at m = 1: "
f"{central} ({dt:.1f} s, {dt / n * 1e3:.0f} ms a draw)")
print(f"LOSS | failing {q:.4f} +- {s:.4f}, bracket [{q:.4f}, "
f"{(tally['LOSS'] + tally['OPEN']) / n:.4f}]")
for t2, (k, L) in by2.items():
if k:
print(f" 2 of type {'1^2 1' if t2 == 'r' else '1^3'}: "
f"{L} of {k}, {L / k:.4f} +- "
f"{sqrt(L / k * (1 - L / k) / k):.4f}")
d = m * q
hi = m * (tally["LOSS"] + tally["OPEN"]) / n
print(f"DENSITY LOSS {d:.5f} [{d:.5f}, {hi:.5f}] +- {m * s:.5f} among "
f"all cubic fields, either signature")
dF = d * 21 / 17
print(f" among those with an F_2 place {dF:.5f} +- {m * s * 21 / 17:.5f}")
lam = nf2 * dF
print(f" the fields to 3000: the races' 4 losses of {nf2} with an "
f"F_2 place; at {nf2} x "
f"{dF:.5f} = {lam:.2f}, P(X <= 4) {poisson.cdf(4, lam):.3f}, "
f"P(X >= 4) {poisson.sf(3, lam):.3f}")
print(f"Y = {YCUT}: decided {decY} of {n}, {decY / n:.3f}")
print(f"{sum(CHECKS)} of {len(CHECKS)} checks pass; "
f"{time.time() - t0:.1f} s")
if __name__ == "__main__":
main()