"""residual.py -- the real fields' remainder after the term at s = 0, read one class character at a time. QUESTION. origin.py reads, over a population of quadratic fields with h+ >= 2, the level of a class C counted in prime powers, PI_C/(PI/h+), against E_C, the explicit formula's term at s = 0. What is left on a cell C is, per field, an identity of orthogonality: level_C - 1 - E_C = sum_{chi != 1} conj chi(C) r_chi / PI(x, 1), r_chi = PI(x, chi) - [-m (loglog x + gamma) + log L*(0, chi)], the nontrivial zeros' sum of chi up to O(1/(x log x)). Over real fields it is small at h+ = 2 and large at h+ >= 4. Is that growth the count of nontrivial characters, each kind of character carrying a fixed mean rho = r_chi/PI(x, 1), or does a character's own mean move with h+ or with |D|? KINDS. Real even (chi(k0) = +1, or no k0), real odd (chi(k0) = -1), imaginary; each split genus (real-valued) or complex, the split printed and not fitted. HAND, from origin.py's F2 at 250. h+ = 2 with k0: one odd character, rho_o = T = -K0, about +0.008; h+ = 3 (no k0, two complex even characters): rho_e about -0.012 +- 0.007; h+ = 4 elementary with k0, T = rho_e + 2 rho_o, K0 = rho_e - 2 rho_o: rho_e about +0.009, rho_o about +0.019. So rho_o at h+ >= 4, at larger |D|, is about twice h+ = 2's: a fixed mean pooled over |D| should fail, and |D| is where it moves. PREDICTIONS, frozen before the engine. R1 printed: rho by kind x (genus, complex) x |D| band x h band, at 250 and 10^4; and the fields by h band x |D| band. R2 printed: the field residual on T regressed on (h+ - 1)/PI(x, 1) at every cut, real and imaginary. R3 predicted: rho_o rises with log|D| at 250, slope > 0 by 3 se. KILL FIXED. "A fixed mean per character of each kind carries the field residual's dependence on h+": the kind means fitted over every real field, a field's prediction the sum over its characters of conj chi(C) times its kind's mean; dies if at 250 the T or K0 residual of any h band (2, 3, 4, >= 5, and h >= 4 elementary 2 or not) departs from its prediction by more than 3 of its own se. KILL FIXED|D|: the same with the kind means fitted within each |D| band. C1w every genus character of every field to |D| = C1_D, D = d1 d2, signed: L*(0) real and positive, as the closed form to 1e-8. C5 complex characters, the known answer: Lambda(2) from the kernels W(y; 2) = int_1^oo phi(yt) (t + t^-2) dt, phi = 4 K0(2t) (even), 4t K0(2t) (odd), e^-t (imaginary), against L(2, chi)'s Euler product over p < 10^4, every nontrivial character to |D| = C1_D, relative gap below 2.5/(P ln P). Its break: the Euler product with chi(P bar) = chi(P) fails at the complex characters (> 1e-3) and passes at the genus ones. C6 the identity: the characters' sum reproduces origin.py's cell residual per field to 1e-10. C7 THE KILL'S CONTROL: each rho replaced by its kind's mean plus noise of its kind's spread, KILL FIXED misses; by the kind's mean times h+/4, it fires. THE KERNELS AT s = 2. With the sum Lambda~(s) = sum_n b_n int_0^oo phi(nt/A) t^s dt/t = A^s G(s) L(s, chi), G = Gamma(s/2)^2, Gamma((s+1)/2)^2, Gamma(s), and Theta(1/t) = t Theta(t), Lambda~(s) = sum_n b_n W(n/A; s), W(y; s) = int_1^oo phi(yt)(t^(s-1) + t^-s) dt, which at s = 1 is origin.py's W. At s = 2, z = 2y: W_e = 4 K1(z)/z + 4z int_z^oo K0(v) v^-2 dv, W_o = 2 I2(z)/z^2 + 2z int_z^oo K0(v)/v dv, I2 = int_z^oo v^2 K0 = z^2 K1(z) + z K0(z) + int_z^oo K0, W_c = (1 + y) e^-y / y^2 + e^-y - y E1(y), and L(2) = sum/A^2 (even, imaginary), sum/(A^2 pi/4) (odd: Gamma(3/2)^2 = pi/4). CHANGED AFTER THE FIRST FULL RUN, labelled so. C5's break failed its frozen bar (complex least gap 6.53e-05 against > 1e-3, a hand guess at the break's size), yet sits above C5's own bar at every complex character, so C5 catches the break wherever it can act: the break's bar is now C5's. Its genus half cannot fail, the break being the identity on real characters. The pooled kind means, computed, became a print. ADDED AT THE AUDIT, a check (C8): L*(0, chi) real and positive at every nontrivial character of both populations, complex ones included (C1w reads genus characters only); the trivial row is a positive theta sum by construction, not zeta_K's L*(0), and is left out. ADDED AFTER THE FIRST READING, frozen before its code (section C, real fields, cut 250): the step in rho_o from h+ = 2 to h+ >= 4. S0 printed: the step's targets, rho_o over the odd characters at h = 2 and at h >= 4; each kind's sd per h band. C9u every complex character's conjugate is in its field, every genus character stands alone: a conjugate pair is one unit, rho_chi-bar = conj rho_chi. C9a the seeded kill reproduces KILL FIXED on the data. C9 THE STEP'S CONTROL, 200 seeds, one draw per unit: plant FIXED (each kind its pooled mean plus noise of its sd) and STEP (odd at the step's targets, even at its pooled mean). STEP's z median above 1.5 under STEP, |median| below 0.5 under FIXED, FIXED's STEP sd within 25% of STEP's se, KILL FIXED firing below 0.10 under FIXED. Predicted: KILL FIXED fires at 0.10 to 0.45 of the STEP seeds. Z2 printed: the sum of z^2 over the T and K0 cells of h = 2, 3, 4, >= 5 (seven), and over the four T cells; p against the FIXED seeds, power at their 95th. KILL STEP. STEP = mean rho_o over units at h >= 4 minus at h = 2, against the FIXED seeds' sd: above 3 STANDS, below 1 KILLED, between UNDECIDED; if UNDECIDED, the real population widened to |D| <= 48000 (--dreal 48000), the same kill. Not in the freeze, written with the code: C9u, C9a, and C9's clause that KILL FIXED fires below 0.10 under FIXED. CHANGED AFTER THE WIDENED RUN, labelled so. The STEP plant took the population's own targets, so at 48000 it planted the widened step and C9 failed (median 0.89); it now plants NAMED, the targets at 16000, at either population. The prediction 0.10 to 0.45 is the 16000 one. The line "the named step: observed ... sd from it" was added with it. NAMED is read from the 16000 fields, which the 48000 population contains, so its distance at 48000 is in sample, a bias toward the step. ADDED AT THE AUDIT: R3u (R3 over units, a conjugate pair once, R3 having counted each pair twice); each |D| band's STEP printed as its distance from the named step; a band with fewer than two odd characters a side skipped. THE PARITY IDENTITY, by hand, read by parity.py. Over a real field with k0, the even characters sum chi(C) to (h+/2)[C in {1, k0}], the odd ones to (h+/2)([C = 1] - [C = k0]), so with PI_C the count in C: odd mean of PI(x, chi) = PI_1 - PI_k0, even nontrivial mean of PI(x, chi) = ((h+/2)(PI_1 + PI_k0) - PI)/(h+/2 - 1). The even line holds at h+ >= 4, there being no even nontrivial character at h+ = 2. The odd mean of rho is then (PI_1 - PI_k0 - the odd mean of log L*(0, chi))/PI(x, 1): a fixed rho_o asks that difference of absolute counts, the small primes the trivial class holds over k0, not to scale with h+. And by the functional equation log L*(0, chi) is (1/2) log D + log L(1, chi) plus a constant set by the parity, so every rho carries -(1/2) log D/PI(x, 1) against whatever the counts do with |D|. FINDINGS, the real fields at the default |D| <= 48000 (12599 fields, 63511 nontrivial characters, 47189 units), the imaginary at |D| <= 4000 (1208 fields, 22385 characters). F0 CONTROLS. C1w: 1177 genus characters to |D| <= 1500, L*(0) real and positive at every one, the closed form to 9.18e-14. C5: L(2, chi) from the s = 2 kernels against the Euler product, genus 1177 worst 4.43e-07, complex 4632 worst 4.13e-07, bar 2.72e-05; the break moves every complex character by 6.53e-05 or more. C6 to 6.33e-15 (724 real fields), 3.50e-14 (902 imaginary). C7: the fixed plant misses (worst +2.0 se), the h+/4 plant fires (-28.6 se). C8: L*(0) real and positive at all 85896 nontrivial characters. F1 KILL FIXED MISSED, real: pooled rho_e +0.0133, rho_o +0.0091 at 250 put every h band within 3 se, worst -2.8 se (h >= 4 elementary 2, K0 -0.0118 against -0.0025; its T +2.6 se); h = 2 T +0.0077 against +0.0101, h >= 5 T +0.1247 against +0.1245. KILL FIXED|D| MISSED, worst +2.8 se (h >= 4 elementary 2, T). A fixed mean per parity is not killed at its bar, and carries most of the residual's growth with h+; the kill's control never planted the step the hand reading named (rho_o from h = 2 to h >= 4, flat past 4): F6 prices it. F2 R3 MISSED, the other way: rho_o on log|D| at 250, slope -0.0035 +- 0.0007 (-5.1 se, 32585 odd characters taken as independent); R3u, over units, -0.0030 +- 0.0008 (-3.9 se, 25484). The kind means by |D| band at 250, odd +0.0151, +0.0162, +0.0083, even -0.0084, -0.0016, +0.0153: |D| moves both, the odd falling by its slope, the even rising band on band, and pooled the even mean sits above the odd. At 10^4 every real R1 cell lies within 0.005 of 0. F3 R2, the mean of r_chi in counts: real slope +0.63 +- 0.02 at 250, +1.14, +1.30, +1.62, +1.67 at 400 to 2500, +1.86 +- 0.09 at 10^4, intercepts +0.0043 +- 0.0026, +0.0099 +- 0.0022, +0.0082 +- 0.0018 at 250, 400, 630, within 0.003 of 0 from 1000 on; imaginary +0.88 +- 0.04 at 250 with intercept -0.172 +- 0.013, slope -0.16 +- 0.19 and intercept +0.001 at 10^4. F4 printed: the imaginary arm of KILL FIXED fires at h = 4 (T -0.0669 +- 0.0237 against +0.0144, -3.4 se, 54 fields) and at h >= 4 elementary 2 (-3.3 se, 37 fields), the other bands within 0.5 se. F5 the confound, counted: 3793 of the 4271 real fields with h+ >= 5, and 3338 of the 4145 with h+ = 2, have |D| > 8000. F6 THE STEP: targets +0.0075 (3219 odd characters at h = 2), +0.0093 (29366 at h >= 4). KILL STEP UNDECIDED at its frozen bar: +0.0023 +- 0.0020 over units (h >= 4 +0.0098, 22265 units; h = 2 +0.0075, 3219), +1.15 against the FIXED seeds' sd 0.0020, and 2.54 sd below the named step +0.0073, which fires KILL FIXED at 0.735 of its seeds and passes Z2's 95th at 0.915 (the print's "MISSED (0.735)" is the prediction frozen at 16000). Z2 reads 7.12, p 0.215; Z2_T 5.59, p 0.135. NAMED is read from the |D| <= 16000 fields, which this population contains: in sample, a bias toward the step. By |D| band the step reads +3.2, -0.3, +0.8 se, from the named step +1.8, -1.5, -2.5 se: past 8000, +0.0017 +- 0.0023 (2584 odd units at h = 2), the band sits furthest below it, where at 16000 the band 8000 to 16000 read +0.0070 +- 0.0048 (537), on it; the shortfall is the fields past 16000's. TIERS. F1 to F6 are observations over these populations. RUN RECORD. A rehearsal at --d 600: 6/7 (C7's h+/4 plant at 134 fields lacks power, -1.9 se). The first full run 36.7 s, peak commit 131.9 MB, 6/7 (C5's break as frozen); after the change 7/7, 36.7 s, 132.3 MB, under a memory guard (512 MB). At the audit, with C8: 8/8, 36.9 s, 132.0 MB; C8 then narrowed to the nontrivial characters: 8/8, 37.1 s, 132.2 MB. With section C: a rehearsal at --d 600 9/11 (its cells too small for C7 and C9); 11/11, 38.5 s, 135.2 MB; --dreal 48000 10/11, 114.2 s, 189.0 MB (C9 as frozen). After the change: 11/11, 38.2 s, 134.6 MB; --dreal 48000 11/11, 113.4 s, 189.0 MB. At the audit (R3u, the band distances): 11/11, 38.3 s, 135.4 MB; --dreal 48000 11/11, 114.1 s, 190.5 MB. The real default then moved from 16000 to 48000: 11/11, 112.5 s, 190.4 MB. The reading at |D| <= 16000 (--dreal 16000; 4107 real fields, 16640 characters): KILL FIXED missed, pooled rho_e +0.0025, rho_o +0.0146, worst +2.3 se (h >= 4 elementary 2, T); R3 +0.0009 +- 0.0012 (8751 characters); R2 real +0.42 +- 0.04 at 250, +1.43 +- 0.18 at 10^4. The step's targets +0.0082 (1172), +0.0155 (7579), named from here; KILL FIXED fired at 0.175 of the STEP seeds (predicted 0.10 to 0.45: HELD), 0.005 of the FIXED; KILL STEP UNDECIDED, +0.0073 +- 0.0032, +2.38 against the FIXED sd 0.0031; Z2 15.21, p 0.015 (power 0.53); by |D| band +3.2, -0.3, +1.5 se; R3u +0.9 se. """ import os os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import sys import time from math import log, sqrt, pi, ceil HERE = os.path.dirname(os.path.abspath(__file__)) sys.path.insert(0, HERE) import numpy as np # noqa: E402 from scipy.special import k0, k1, k0e, exp1 # noqa: E402 from scipy.integrate import quad # noqa: E402 import principal as P # noqa: E402 import origin as O # noqa: E402 from module_law import check, section, CHECKS # noqa: E402 C1_D = 1500 CUTS = P.CUTS J250, J10K = 0, len(CUTS) - 1 REAL_D, IMAG_D = 48000, P.IMAG_D # the real population widened (F6) if "--d" in sys.argv: # a rehearsal: both populations cut to |D| <= this REAL_D = IMAG_D = int(sys.argv[sys.argv.index("--d") + 1]) C1_D = min(C1_D, REAL_D) if "--dreal" in sys.argv: # the real population alone set to |D| <= this REAL_D = int(sys.argv[sys.argv.index("--dreal") + 1]) SEEDS = 200 NAMED = (0.0082, 0.0155) # the step's targets at |D| <= 16000 (S0), planted DBANDS = (("|D|<=4000", 0, 4000), ("4000-8000", 4000, 8000), ("|D|>8000", 8000, 10 ** 9)) HBANDS = (("h=2", lambda h, el: h == 2), ("h=3", lambda h, el: h == 3), ("h=4", lambda h, el: h == 4), ("h>=5", lambda h, el: h >= 5), ("h>=4 el", lambda h, el: h >= 4 and el), ("h>=4 not", lambda h, el: h >= 4 and not el)) def im2(z): """int_z^oo K0(v) v^-2 dv: Gauss-Laguerre on k0e past 2, quad below.""" z = np.asarray(z, float) out = np.empty_like(z) far = z >= 2 if far.any(): u = z[far][:, None] + O.LAG_X[None, :] out[far] = np.exp(-z[far]) * ((k0e(u) / u ** 2) @ O.LAG_W) at2 = (np.exp(-2.0) * ((k0e(2 + O.LAG_X) / (2 + O.LAG_X) ** 2) @ O.LAG_W)) for i in np.flatnonzero(~far): out[i] = quad(lambda v: k0(v) / v ** 2, z[i], 2.0, epsabs=0, epsrel=1e-13, limit=200)[0] + at2 return out def weights2(kind, y): z = 2 * y if kind == "e": return 4 * k1(z) / z + 4 * z * im2(z) if kind == "o": a, b = O.tails(z) i2 = z ** 2 * k1(z) + z * k0(z) + a return 2 * i2 / z ** 2 + 2 * z * b return (1 + y) * np.exp(-y) / y ** 2 + np.exp(-y) - y * exp1(y) def coeffs(G, chis, N): """b_n for n <= N, each character, as origin.lam1 builds them.""" b = np.ones((len(chis), N + 1), complex) b[:, 0] = 0 for p in P.SMALL: if p > N: break kr = P.kron(G.D, p) e, pe = 1, p while pe <= N: idx = np.arange(pe, N + 1, pe) idx = idx[idx % (pe * p) != 0] b[:, idx] *= O.local(G, chis, p, kr, e)[:, None] e, pe = e + 1, pe * p return b def parity(G, chis): if G.D > 0 and G.k0 != G.e: return ["o" if np.real(chis[i, G.k0]) < 0 else "e" for i in range(len(chis))] return ["e" if G.D > 0 else "c"] * len(chis) def l2_afe(G, chis, kinds): D = G.D A = sqrt(D) / pi if D > 0 else sqrt(-D) / (2 * pi) N = int(ceil((O.NA if D > 0 else 2 * O.NA) * A)) + 5 b = coeffs(G, chis, N) n = np.arange(1, N + 1) out = np.zeros(len(chis), complex) for kind in set(kinds): rows = [i for i, k in enumerate(kinds) if k == kind] out[rows] = b[rows, 1:] @ weights2(kind, n / A) norm = np.array([A * A * (pi / 4 if k == "o" else 1) for k in kinds]) return out / norm def l2_euler(G, chis, broken=False): out = np.ones(len(chis), complex) for p in P.SMALL: kr = P.kron(G.D, p) q = 1.0 / (p * p) if kr == 1: v = chis[:, G.place(p)] w = v if broken else np.conj(v) out /= (1 - v * q) * (1 - w * q) elif kr == -1: out /= 1 - q * q else: out /= 1 - chis[:, G.place(p)] * q return out def genus_pairs(D): """Every unordered D = d1 d2, both fundamental, neither 1.""" out = set() for d1 in range(-abs(D), abs(D) + 1): if d1 in (0, 1, D) or D % d1 or not P.fundamental(d1): continue d2 = D // d1 if d2 != 1 and P.fundamental(d2): out.add(tuple(sorted((d1, d2)))) return sorted(out) def nontrivial(chis): return [i for i in range(len(chis)) if np.max(np.abs(chis[i] - 1)) > 1e-6] def section_controls(): section("A THE CONTROLS") w1, n1, bad1, miss = 0.0, 0, 0, 0 w5g, w5c, n5g, n5c, brk_g, brk_c = 0.0, 0.0, 0, 0, 0.0, np.inf P5 = P.SMALL[-1] bar5 = 2.5 / (P5 * log(P5)) for D in list(range(5, C1_D + 1)) + list(range(-C1_D, 0)): if not P.fundamental(D): continue G = P.Group(D) if G.h < 2: continue chis = O.characters(G) m, L, _ = O.lstar(G, chis) test = [p for p in P.SMALL[1:200] if P.kron(D, p) == 1] cols = [G.place(p) for p in test] pairs = genus_pairs(D) if len(pairs) != len(G.two) - 1: miss += 1 for d1, d2 in pairs: want = np.array([P.kron(d1, p) for p in test]) hit = [i for i in range(G.h) if np.max(np.abs(chis[i, cols] - want)) < 1e-6] if len(hit) != 1: miss += 1 continue known = ((O.l0(d1) if d1 < 0 else O.l0p(d1)) * (O.l0(d2) if d2 < 0 else O.l0p(d2))) got = L[hit[0]] if not (known > 0 and got.real > 0 and abs(got.imag) < 1e-10 * abs(got)): bad1 += 1 w1 = max(w1, abs(got - known) / known) n1 += 1 nt = nontrivial(chis) kinds = parity(G, chis) afe = l2_afe(G, chis[nt], [kinds[i] for i in nt]) eu = l2_euler(G, chis[nt]) br = l2_euler(G, chis[nt], broken=True) for j, i in enumerate(nt): rel = abs(afe[j] - eu[j]) / abs(eu[j]) relb = abs(afe[j] - br[j]) / abs(br[j]) if np.max(np.abs(chis[i].imag)) < 1e-9: w5g, n5g, brk_g = max(w5g, rel), n5g + 1, max(brk_g, relb) else: w5c, n5c, brk_c = max(w5c, rel), n5c + 1, min(brk_c, relb) check("C1w every genus character: L*(0) real, positive, the closed form", w1 < 1e-8 and bad1 == 0 and miss == 0 and n1 > 100, f"{n1} characters, |D| <= {C1_D}, worst relative {w1:.2e}," f" sign or reality off at {bad1}, unmatched {miss}") check("C5 L(2, chi) from the s = 2 kernels against the Euler product", max(w5g, w5c) < bar5 and n5c > 100, f"genus {n5g} worst {w5g:.2e}, complex {n5c} worst {w5c:.2e}," f" bar {bar5:.2e}") check("C5 break: chi(P bar) = chi(P) fails every complex character," " passes every genus one", brk_c > bar5 and brk_g < bar5, f"complex least gap {brk_c:.2e}, genus worst {brk_g:.2e}") def read_population(sign): """Per field (D, h, elementary, k0 present, residual T and K0 per cut) and per character (field index, kind, genus, chi(k0), rho per cut).""" lo, hi = (5, REAL_D + 1) if sign > 0 else (-IMAG_D, 0) fields, chars, worst6, n6 = [], [], 0.0, 0 UNITS.clear() UNIT_N[0] = 0 global POS_BAD, POS_N for D in range(lo, hi): if not P.fundamental(D): continue G = P.Group(D) if G.h < 2: continue raw, cor = P.counts(G) chis = O.characters(G) m, L, _ = O.lstar(G, chis) nt = nontrivial(chis) Ln = L[nt] POS_BAD += int(np.sum(~(Ln.real > 0) | ~(np.abs(Ln.imag) <= 1e-9 * np.abs(Ln)))) POS_N += len(Ln) kinds = parity(G, chis) cor = np.array(cor) # (cuts, h) tot = cor.sum(axis=1) PIchi = cor @ chis.T # (cuts, h): PI(x, chi) lx = np.log(np.log(np.array(CUTS, float))) + O.GAMMA term = -np.outer(lx, m) + np.log(np.abs(L))[None, :] rho = (PIchi - term) / tot[:, None] # (cuts, h) hask0 = D > 0 and G.k0 != G.e resT = np.real(rho[:, nt].sum(axis=1)) resK = (np.real((np.conj(chis[nt, G.k0])[None, :] * rho[:, nt]) .sum(axis=1)) if hask0 else None) fi = len(fields) fields.append((abs(D), G.h, len(G.two) == G.h, hask0, resT, resK, tot)) local = {} for i in nt: chars.append((fi, kinds[i], np.max(np.abs(chis[i].imag)) < 1e-9, np.real(chis[i, G.k0]), np.real(rho[:, i]))) gap = np.max(np.abs(chis[nt] - np.conj(chis[i])), axis=1) key = min(i, nt[int(np.argmin(gap))]) if gap.min() > 1e-6: key = -1 - i # no conjugate: C9u fails if key not in local: local[key] = UNIT_N[0] UNIT_N[0] += 1 UNITS.append(local[key]) if abs(D) <= 3000: # C6 on the small fields out = O.read_field(G)[0] for j in range(len(CUTS)): dv, E = out[j] worst6 = max(worst6, abs(dv[G.e] - E[G.e] - resT[j])) if hask0: worst6 = max(worst6, abs(dv[G.k0] - E[G.k0] - resK[j])) n6 += 1 return fields, chars, worst6, n6 POS_BAD, POS_N = 0, 0 UNITS, UNIT_N = [], [0] # per character its unit: a conjugate pair is one def mean_se(v): v = np.asarray(v, float) if len(v) < 2: return (v.mean() if len(v) else np.nan), np.nan return v.mean(), v.std(ddof=1) / sqrt(len(v)) def fit_means(chars, j, fsel=None): """Kind means of rho at cut j over the characters of the selected fields.""" acc = {} for fi, kind, _, _, r in chars: if fsel is None or fsel(fi): acc.setdefault(kind, []).append(r[j]) return {k: np.mean(v) for k, v in acc.items()} def predict(fields, chars, means_of): """Per field the predicted T and K0 residual from kind means; means_of(fi) gives the means that field uses.""" pT = np.zeros(len(fields)) pK = np.zeros(len(fields)) for fi, kind, _, ck0, _ in chars: mu = means_of(fi)[kind] pT[fi] += mu pK[fi] += ck0 * mu return pT, pK def kill(fields, chars, j, means_of, label, quiet=False): """The KILL FIXED test at cut j: True when it fires.""" pT, pK = predict(fields, chars, means_of) fired, worst = False, (0.0, "") for bn, test in HBANDS: for cell in ("T", "K0"): obs, pre = [], [] for fi, (a, h, el, hk, rT, rK, _) in enumerate(fields): if not test(h, el): continue if cell == "T": obs.append(rT[j]), pre.append(pT[fi]) elif hk: obs.append(rK[j]), pre.append(pK[fi]) if len(obs) < 2: continue mo, se = mean_se(obs) z = (mo - np.mean(pre)) / se if abs(z) > abs(worst[0]): worst = (z, f"{bn} {cell}") fired |= abs(z) > 3 if not quiet: print(f" {label:9s} {bn:9s} {cell:2s} observed " f"{mo:+.4f}+-{se:.4f} predicted {np.mean(pre):+.4f}" f" ({z:+.1f} se, {len(obs)})") return fired, worst def dband(a): return next(i for i, (_, lo, hi) in enumerate(DBANDS) if lo < a <= hi) def section_population(sign, name): t0 = time.time() section(f"B {name.upper()}") fields, chars, worst6, n6 = read_population(sign) print(f" {len(fields)} fields, {len(chars)} nontrivial characters" f" [{time.time() - t0:.0f} s]") check(f"C6 {name}: the characters' sum is origin.py's cell residual", worst6 < 1e-10, f"{n6} fields with |D| <= 3000, worst {worst6:.2e}") print(" the fields by h band x |D| band:") for bn, test in HBANDS: cnt = [sum(1 for a, h, el, *_ in fields if test(h, el) and dband(a) == b) for b in range(3)] print(f" {bn:9s} " + " ".join(f"{DBANDS[b][0]} {cnt[b]:5d}" for b in range(3))) print(" R1 rho per character, mean+-se (count), by kind, |D| band and" " h band:") for j in (J250, J10K): print(f" cut {CUTS[j]}") for kind in sorted({c[1] for c in chars}): for gen in (True, False): for b in range(3): line = (f" {kind} {'genus ' if gen else 'complex'}" f" {DBANDS[b][0]:10s}") for bn, test in HBANDS[:4]: v = [r[j] for fi, k, g, _, r in chars if k == kind and g == gen and dband(fields[fi][0]) == b and test(fields[fi][1], fields[fi][2])] if len(v) >= 2: mo, se = mean_se(v) line += f" {bn} {mo:+.4f}+-{se:.4f} ({len(v)})" if "(" in line: print(line) print(" R2 the field residual on T against (h+ - 1)/PI(x, 1):") for j, x in enumerate(CUTS): y = np.array([f[4][j] for f in fields]) xs = np.array([(f[1] - 1) / f[6][j] for f in fields]) X = np.column_stack([np.ones(len(y)), xs]) beta, res, *_ = np.linalg.lstsq(X, y, rcond=None) resid = y - X @ beta s2 = resid @ resid / (len(y) - 2) cov = s2 * np.linalg.inv(X.T @ X) print(f" {x:6d} slope {beta[1]:+.4f}+-{sqrt(cov[1, 1]):.4f}" f" intercept {beta[0]:+.5f}+-{sqrt(cov[0, 0]):.5f}") if sign > 0: lg = np.array([log(fields[fi][0]) for fi, k, *_ in chars if k == "o"]) ro = np.array([r[J250] for fi, k, g, ck, r in chars if k == "o"]) X = np.column_stack([np.ones(len(lg)), lg]) beta = np.linalg.lstsq(X, ro, rcond=None)[0] resid = ro - X @ beta cov = resid @ resid / (len(ro) - 2) * np.linalg.inv(X.T @ X) z = beta[1] / sqrt(cov[1, 1]) print(f" R3 rho_o on log|D| at 250: slope {beta[1]:+.4f}" f"+-{sqrt(cov[1, 1]):.4f} ({z:+.1f} se, {len(ro)} characters," f" taken as independent): {'HELD' if z > 3 else 'MISSED'}") un = np.array(UNITS) uo = np.array([k == "o" for _, k, *_ in chars]) _, first = np.unique(un[uo], return_index=True) lgu, rou = lg[first], ro[first] Xu = np.column_stack([np.ones(len(lgu)), lgu]) bu = np.linalg.lstsq(Xu, rou, rcond=None)[0] ru = rou - Xu @ bu cu = ru @ ru / (len(rou) - 2) * np.linalg.inv(Xu.T @ Xu) print(f" R3u the same over units, a conjugate pair once (added at" f" the audit): slope {bu[1]:+.4f}+-{sqrt(cu[1, 1]):.4f}" f" ({bu[1] / sqrt(cu[1, 1]):+.1f} se, {len(rou)} units)") # C7 first: the kill's control on planted rho rng = np.random.default_rng(1660) mu = fit_means(chars, J250) sd = {k: np.std([r[J250] for _, kk, _, _, r in chars if kk == k], ddof=1) for k in mu} for plant, fn in (("fixed", lambda h, k: mu[k] + rng.normal(0, sd[k])), ("h+/4", lambda h, k: mu[k] * h / 4 + rng.normal(0, sd[k]))): fake = [] for fi, kind, g, ck0, r in chars: v = r.copy() v[J250] = fn(fields[fi][1], kind) fake.append((fi, kind, g, ck0, v)) ff = [] for fi, f in enumerate(fields): ff.append(list(f[:4]) + [np.zeros(len(CUTS)), np.zeros(len(CUTS)) if f[3] else None, f[6]]) for fi, kind, g, ck0, v in fake: ff[fi][4][J250] += v[J250] if ff[fi][3]: ff[fi][5][J250] += ck0 * v[J250] mf = fit_means(fake, J250) fired, worst = kill(ff, fake, J250, lambda fi: mf, "", quiet=True) want = plant == "h+/4" check(f"C7 the kill on planted rho ({plant}) " f"{'fires' if want else 'misses'}", fired == want, f"worst {worst[0]:+.1f} se at {worst[1]}") print(" KILL FIXED at 250, kind means pooled: " + ", ".join( f"{k} {v:+.4f}" for k, v in sorted(mu.items()))) fired, worst = kill(fields, chars, J250, lambda fi: mu, "pooled") print(f" KILL FIXED {'FIRED' if fired else 'MISSED'}: worst" f" {worst[0]:+.1f} se at {worst[1]}") mb = [fit_means(chars, J250, lambda fi, b=b: dband(fields[fi][0]) == b) for b in range(3)] print(" kind means by |D| band at 250: " + "; ".join( f"{DBANDS[b][0]} " + ", ".join(f"{k} {v:+.4f}" for k, v in sorted(mb[b].items())) for b in range(3))) fired2, worst2 = kill(fields, chars, J250, lambda fi: mb[dband(fields[fi][0])], "by |D|") print(f" KILL FIXED|D| {'FIRED' if fired2 else 'MISSED'}: worst" f" {worst2[0]:+.1f} se at {worst2[1]}") section_step(fields, chars, fired, worst) else: mu = fit_means(chars, J250) print(" the imaginary arm of KILL FIXED at 250, printed, mean" f" {mu['c']:+.4f}:") fired, worst = kill(fields, chars, J250, lambda fi: mu, "pooled") print(f" imaginary FIXED {'fires' if fired else 'misses'}: worst" f" {worst[0]:+.1f} se at {worst[1]}") def section_step(fields, chars, fired0, worst0): """C: the step in rho_o from h+ = 2 to h+ >= 4, priced by planted seeds with conjugate pairs tied.""" section("C THE STEP IN RHO_O, REAL, CUT 250") F = len(fields) fi = np.array([c[0] for c in chars]) odd = np.array([c[1] == "o" for c in chars]) ck0 = np.array([c[3] for c in chars]) r = np.array([c[4][J250] for c in chars]) unit = np.array(UNITS) h = np.array([f[1] for f in fields]) el = np.array([f[2] for f in fields]) hk = np.array([f[3] for f in fields]) a = np.array([f[0] for f in fields]) hc = h[fi] sizes = np.bincount(unit) cplx = np.array([not c[2] for c in chars]) check("C9u every complex character has its conjugate, every genus one" " stands alone", np.all(sizes[unit] == np.where(cplx, 2, 1)), f"{len(sizes)} units over {len(chars)} characters") bands = [(bn, np.array([t(hh, ee) for hh, ee in zip(h, el)])) for bn, t in HBANDS] disj = bands[:4] def cells(v, which): """z per (band, cell), pooled kind means as KILL FIXED fits.""" mo, me = v[odd].mean(), v[~odd].mean() pre = np.where(odd, mo, me) T = np.bincount(fi, v, minlength=F) K = np.bincount(fi, ck0 * v, minlength=F) pT = np.bincount(fi, pre, minlength=F) pK = np.bincount(fi, ck0 * pre, minlength=F) out = [] for bn, m in which: for cell, o, p_, mm in (("T", T, pT, m), ("K0", K, pK, m & hk)): if mm.sum() < 2: continue x = o[mm] out.append((f"{bn} {cell}", (x.mean() - p_[mm].mean()) / (x.std(ddof=1) / sqrt(len(x))))) return out def step(v, sel=None): """mean rho_o over units at h >= 4 minus at h = 2, its se.""" out = [] for m in (odd & (hc >= 4), odd & (hc == 2)): if sel is not None: m = m & sel u, first = np.unique(unit[m], return_index=True) x = v[m][first] out.append((x.mean(), x.var(ddof=1) / len(x), len(x))) return out[0][0] - out[1][0], sqrt(out[0][1] + out[1][1]), out zs = cells(r, bands) wz = max(zs, key=lambda t: abs(t[1])) check("C9a the seeded kill reproduces KILL FIXED on the data", abs(wz[1] - worst0[0]) < 1e-9 and (abs(wz[1]) > 3) == fired0, f"worst {wz[1]:+.6f} at {wz[0]} against {worst0[0]:+.6f}") t2 = r[odd & (hc == 2)].mean() t4 = r[odd & (hc >= 4)].mean() print(f" S0 targets, count-weighted: rho_o at h=2 {t2:+.4f}" f" ({np.sum(odd & (hc == 2))}), at h>=4 {t4:+.4f}" f" ({np.sum(odd & (hc >= 4))})") mu = {True: r[odd].mean(), False: r[~odd].mean()} sd = {True: r[odd].std(ddof=1), False: r[~odd].std(ddof=1)} for k, nm in ((False, "even"), (True, "odd")): line = f" S0 sd of rho {nm:4s} pooled {sd[k]:.4f}; by band" for bn, m in disj: sel = (odd == k) & m[fi] if sel.sum() >= 2: line += f" {bn} {r[sel].std(ddof=1):.4f} ({sel.sum()})" print(line) S, se, parts = step(r) print(f" STEP observed {S:+.4f} +- {se:.4f} ({S / se:+.1f} se by its" f" own se): h>=4 {parts[0][0]:+.4f} ({parts[0][2]} units)," f" h=2 {parts[1][0]:+.4f} ({parts[1][2]})") ab = np.array([dband(x) for x in a])[fi] for b in range(3): sel = ab == b if min(np.sum(odd & sel & (hc >= 4)), np.sum(odd & sel & (hc == 2))) < 2: continue d, s_, pb = step(r, sel) print(f" {DBANDS[b][0]:10s} STEP {d:+.4f} +- {s_:.4f}" f" ({d / s_:+.1f} se; units {pb[0][2]}, {pb[1][2]});" f" from the named step {(d - (NAMED[1] - NAMED[0])) / s_:+.1f}" f" se (added at the audit)") zd = cells(r, disj) Z2 = sum(z * z for _, z in zd) Z2T = sum(z * z for c, z in zd if c.endswith(" T")) print(" the disjoint cells: " + ", ".join(f"{c} {z:+.2f}" for c, z in zd)) rng = np.random.default_rng(1661) sdk = np.where(odd, sd[True], sd[False]) base = {"FIXED": np.where(odd, mu[True], mu[False]), "STEP": np.where(odd, np.where(hc == 2, NAMED[0], NAMED[1]), mu[False])} res = {} for plant in ("FIXED", "STEP"): fires, steps, z2s, z2ts = 0, [], [], [] for _ in range(SEEDS): v = base[plant] + rng.normal(0, 1, len(sizes))[unit] * sdk fires += any(abs(z) > 3 for _, z in cells(v, bands)) steps.append(step(v)[0]) zz = cells(v, disj) z2s.append(sum(z * z for _, z in zz)) z2ts.append(sum(z * z for c, z in zz if c.endswith(" T"))) res[plant] = (fires / SEEDS, np.array(steps), np.array(z2s), np.array(z2ts)) sdn = res["FIXED"][1].std(ddof=1) for plant in ("FIXED", "STEP"): fr, st, z2, z2t = res[plant] print(f" C9 plant {plant:5s}: KILL FIXED fires {fr:.3f} of" f" {SEEDS}; STEP z median {np.median(st) / sdn:+.2f}" f" (sd {st.std(ddof=1):.4f}); Z2 median {np.median(z2):.1f}," f" Z2_T {np.median(z2t):.1f}") pw = res["STEP"][0] print(f" prediction: KILL FIXED fires 0.10 to 0.45 under STEP:" f" {'HELD' if 0.10 <= pw <= 0.45 else 'MISSED'} ({pw:.3f})") check("C9 the plants: STEP z median above 1.5 under STEP, |median|" " below 0.5 under FIXED; FIXED's STEP sd within 25% of the se;" " KILL FIXED fires below 0.10 under FIXED", np.median(res["STEP"][1]) / sdn > 1.5 and abs(np.median(res["FIXED"][1]) / sdn) < 0.5 and abs(sdn / se - 1) < 0.25 and res["FIXED"][0] < 0.10, f"FIXED seeds' STEP sd {sdn:.4f} against se {se:.4f}") for nm, obs, k in (("Z2 (T, K0; 7 cells)", Z2, 2), ("Z2_T (T; 4 cells)", Z2T, 3)): null = res["FIXED"][k] q95 = np.quantile(null, 0.95) print(f" {nm}: observed {obs:.2f}, FIXED null median" f" {np.median(null):.2f}, 95th {q95:.2f}, p" f" {np.mean(null >= obs):.3f}; power at the 95th" f" {np.mean(res['STEP'][k] > q95):.3f}") z = S / sdn verdict = "STANDS" if z > 3 else "KILLED" if z < 1 else "UNDECIDED" print(f" KILL STEP: STEP {S:+.4f} against the FIXED seeds' sd" f" {sdn:.4f}, {z:+.2f}: {verdict}") d = NAMED[1] - NAMED[0] print(f" the named step {d:+.4f}: observed {(S - d) / sdn:+.2f} sd" f" from it") def main(): t0 = time.time() section_controls() print(f" [{time.time() - t0:.0f} s]") for sign, name in ((+1, "real, narrow"), (-1, "imaginary")): section_population(sign, name) check("C8 (added at the audit) L*(0, chi) real and positive at every" " nontrivial character", POS_BAD == 0 and POS_N > 0, f"{POS_N} characters, off at {POS_BAD}") print(f"\n{sum(CHECKS)}/{len(CHECKS)} checks PASS, " f"{time.time() - t0:.1f} s") raise SystemExit(0 if all(CHECKS) else 1) if __name__ == "__main__": main()