origin.py
Cited on: PRINCIPAL.
To run it: save origin.py, with module_law.py, principal.py beside it; needs numpy, scipy.
"""origin.py -- the term at s = 0 in a count of prime powers, and the
real fields' remainder on the trivial class and k0.
QUESTION. K = Q(sqrt D), Cl+ its narrow class group of order h+ >= 2
(Cl over an imaginary field). principal.py counts PI_C(x), the prime-ideal
powers P^j of norm < x in the class C weighted 1/j, at six cuts x; read
over a population its level PI_C/(PI/h+) departs from 1 by a remainder
that over real fields sits on the trivial class (high) and on k0, the
narrow class of (alpha) with N(alpha) < 0 (low). Does the explicit
formula's term at s = 0 -- the zero's order and the leading coefficient
L*(0, chi) -- carry that remainder?
THE DERIVATION (before code). level_C - 1 = sum_{chi != 1} conj chi(C)
PI(x, chi) / PI(x, 1). Perron: PI(x, chi) = (1/2 pi i) int log L(s, chi)
x^s ds / s. With log L = m log s + log L*(0) + O(s) at s = 0, and s
int_1^oo loglog u u^(-s-1) du = -gamma - log s,
PI(x, chi) = -m (loglog x + gamma) + log |L*(0, chi)|
- sum_rho Li(x^rho) + O(1/(x log x)).
L(s, chi) is real on the real axis (conjugation sends C to C^-1, keeping
norms), and L*(0) > 0: it is Lambda(1) = Lambda(0) over the gamma
factor's positive leading coefficient, and L(1, chi) > 0; the modulus
taken below changes nothing. Over a real field chi with
chi(k0) = +1 carries Gamma_R(s)^2 (m = 2), chi(k0) = -1 carries
Gamma_R(s + 1)^2 (m = 0); no k0, every chi even. Over an imaginary field
Gamma_C(s), m = 1. With root number 1 (chi is trivial on the ideles of
Q, so Frohlich-Queyrut, Invent. Math. 20 (1973), known) and
Lambda(s) = Lambda(1 - s):
even L*(0) = Lambda(1)/4, odd and imaginary L*(0) = Lambda(1),
sum_n b_n W(n/A) = Lambda(1) (even, imaginary), pi Lambda(1) (odd),
b_n = sum_C chi(C) a_C(n),
a_C(n) the ideals of norm n in C, A = sqrt(D)/pi (real), sqrt|D|/(2 pi)
(imaginary), and from Theta(1/t) = t Theta(t),
W_e(y) = 4[(1/2y) int_2y^oo K0 + int_2y^oo K0(u)/u du],
W_o(y) = 2 K1(2y) + 2 int_2y^oo K0,
W_C(y) = e^-y / y + E1(y).
The prediction, every term but the nontrivial zeros:
E_C(x) = sum_{chi != 1} conj chi(C) [-m (loglog x + gamma)
+ log |L*(0, chi)|] / PI(x, 1).
SIZE by hand at D ~ 8000, x = 250: h+ = 2, E_T ~ +0.035, E_k0 ~ -0.035;
h+ = 4, E_T ~ +0.047, E_k0 ~ -0.095. Read: T 1.07, k0 0.82.
PREDICTIONS, frozen before the engine.
Q1 real fields, every cut: |mean(level - 1 - E)| < |mean(level - 1)|
on T and on K0.
Q2 printed, not judged: the residual by cell and cut with its se, real
and imaginary.
KILL CARRIES. "The term at s = 0 carries the real fields' remainder on
T and K0" only if the residual's mean on T and on K0 is within 0.01
of 0 at cut 250 and within 0.005 at 10^4.
C1 genus characters, the known answer: D = d1 d2, d1 the fundamental
factor of least |d1|. Real, d1, d2 < 0: L(0, chi) = L(0, chi_d1)
L(0, chi_d2); d1, d2 > 0: L*(0) = L'(0, chi_d1) L'(0, chi_d2);
imaginary, d1 > 0 > d2: L*(0) = L'(0, chi_d1) L(0, chi_d2). L(0,
chi_d) = -(1/|d|) sum a chi_d(a) (d < 0); L'(0, chi_d) = -(1/2)
sum chi_d(a) log sin(pi a/d) (d > 0). To 1e-8 relative, every field
with |D| <= C1_D.
C2 the characters are homomorphisms, h+ of them, orthogonal.
C3 Lambda(1) moves < 1e-10 when the sum runs to 1.5 N.
C4 THE BREAK: k0's parity flipped moves the K0 residual at 250 by more
than 0.05.
CHANGED BEFORE THE FULL RUN, after a rehearsal at --d 600 failed C1 and
C3 (nothing read): the odd sum carries pi Lambda (Gamma_R(s + 1)^2 =
pi^-(s+1) Gamma((s+1)/2)^2), so L(0) = sum/pi, C1's odd fields off by
exactly pi; scipy's it2i0k0 is 1.2e-7 off at z = 16 where the value is
1e-8, so the tails past z = 2 are Gauss-Laguerre sums of k0e (1e-14
against mpmath quadrature); W_C decays as e^-y, so the imaginary sum
runs to 2 NA A. Then C1 9.2e-14 over 255 fields, C3 7.8e-14.
ADDED AFTER THE FIRST FULL RUN, a print (Q3): the residual on T and K0
by |D| band and by h, at 250 and 10^4.
FINDINGS.
F0 CONTROLS (C1 to C4 held). Genus characters as their closed forms to
9.18e-14 over 671 fields with |D| <= 1500; 800 groups' characters
homomorphisms and orthogonal to 1.01e-11; Lambda(1) settled to
7.82e-14; k0's parity flipped moves K0's residual at 250 by
-0.2625 (C4 prints residual less flipped residual), a break any
change of m passes, unsaid at the freeze.
F1 THE TERM CARRIES PART (KILL CARRIES FIRED, Q1 failed at one cell).
Real, 4107 fields: T reads +0.0679 at 250, E +0.0321, residual
+0.0358 +- 0.0029; K0 (3352 fields) -0.1763, E -0.1405, residual
-0.0358 +- 0.0034. At 10^4 T +0.0050, E -0.0000, residual +0.0051
+- 0.0006 (Q1's one miss); K0 -0.0119, E -0.0082, residual -0.0037
+- 0.0007. The residual on T and K0 falls from 400 on (+0.0483,
-0.0593 at 400). E has the departures' signs at every cut but T
at 10^4, where it prints -0.0000.
F2 printed (Q3), real: at 250 the residual is small at h = 2 and
large at h >= 4 with or without complex characters, larger on K0
without; h is confounded with |D|, uncounted.
h = 2: T +0.0089 +-
0.0025, K0 -0.0082 +- 0.0029 at 250, +0.0011 and -0.0011 at 10^4;
h >= 4: T +0.0539, K0 -0.0506 at 250, +0.0077, -0.0051 at 10^4. On
T at 250 it grows with |D|: +0.0140, +0.0312, +0.0483 over |D| <=
4000, to 8000, beyond. h = 3 (88 fields, complex characters only,
the band added at the audit): T -0.0233 +- 0.0144 at 250. h >= 4
split at the audit by shape: elementary 2 (every chi genus, 1379
fields) T +0.0483 +- 0.0047, K0 -0.0304 at 250, T +0.0088 at
10^4; not (1144) T +0.0607 +- 0.0080, K0 -0.0799 at 250, T
+0.0064 at 10^4.
F3 printed (Q2), imaginary, 1208 fields: T at 10^4 -0.0127, E -0.0114,
residual -0.0014 +- 0.0022; at 250 E -0.0759 against the reading
+0.0128, residual +0.0887; that x = 250 below most conductors and
the floor's empty T carry it is a reading not tested.
TIERS. F1 to F3 are observations over this population.
RUN RECORD. 4/4 checks, 42.3 s, peak commit 117.0 MB under a memory
guard (512 MB); the first full run 42.7 s, the same prints without Q3.
At the audit C1 counts the fields whose genus character it cannot
find (0) and Q3 gains h = 3 and h >= 4 by shape: 4/4, 42.3 s, 116.8 MB.
"""
import os
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import sys
import time
from math import log, sqrt, pi, sin, ceil
HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, HERE)
import numpy as np # noqa: E402
from scipy.special import k0e, k1, iti0k0, it2i0k0, exp1 # noqa: E402
import principal as P # noqa: E402
from module_law import check, section, CHECKS # noqa: E402
GAMMA = 0.5772156649015329
C1_D = 1500 # the genus control's |D| range
NA = 15 # the sum runs to n <= NA * A
REAL_D, IMAG_D = P.REAL_D, P.IMAG_D
if "--d" in sys.argv: # a rehearsal: both populations cut to |D| <= this
REAL_D = IMAG_D = C1_D = int(sys.argv[sys.argv.index("--d") + 1])
def kronecker(d, n):
"""The Kronecker symbol (d/n), n >= 1."""
out = 1
while n % 2 == 0:
n //= 2
if d % 2 == 0:
return 0
out *= 1 if d % 8 in (1, 7) else -1
a, m = d % n if n > 1 else 0, n
if m == 1:
return out
while a:
while a % 2 == 0:
a //= 2
if m % 8 in (3, 5):
out = -out
a, m = m, a
if a % 4 == 3 and m % 4 == 3:
out = -out
a %= m
return out if m == 1 else 0
def characters(G):
"""The h characters of G as an (h, h) array, row = character, column =
class, from the joint eigenvectors of the regular representation."""
h = G.h
rng = np.random.default_rng(G.D % 1000003)
M = np.zeros((h, h), complex)
for g, c in zip(range(h), rng.normal(size=h) + 1j * rng.normal(size=h)):
for x in range(h):
M[G.mul(g, x), x] += c
_, V = np.linalg.eig(M)
X = (V / V[G.e]).T
return X
LAG_X, LAG_W = np.polynomial.laguerre.laggauss(40)
def tails(z):
"""int_z^oo K0 and int_z^oo K0(u)/u du. scipy's it2i0k0 loses all
accuracy past z ~ 10 (1.2e-7 off at z = 16, where the value is 1e-8),
so past z = 2 both are Gauss-Laguerre sums of the scaled k0e."""
z = np.asarray(z, float)
a = np.pi / 2 - iti0k0(z)[1]
b = it2i0k0(z)[1]
far = z >= 2
if far.any():
u = z[far][:, None] + LAG_X[None, :]
g = k0e(u)
a[far] = np.exp(-z[far]) * (g @ LAG_W)
b[far] = np.exp(-z[far]) * ((g / u) @ LAG_W)
return a, b
def weights(kind, y):
if kind == "e":
z = 2 * y
a, b = tails(z)
return 4 * (a / z + b)
if kind == "o":
z = 2 * y
return 2 * k1(z) + 2 * tails(z)[0]
return np.exp(-y) / y + exp1(y)
def local(G, chis, p, kr, e):
"""b(p^e) for each character."""
if kr == 1:
c = G.place(p)
v = chis[:, c]
return sum(v ** (2 * i - e) for i in range(e + 1))
if kr == -1:
return np.ones(len(chis)) if e % 2 == 0 else np.zeros(len(chis))
return chis[:, G.place(p)] ** e
def lam1(G, chis, kinds, A, N):
"""Lambda(1, chi) for each character, the sum to n <= N."""
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] *= local(G, chis, p, kr, e)[:, None]
e, pe = e + 1, pe * p
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]
w = weights(kind, n / A)
out[rows] = b[rows, 1:] @ w
return out
def lstar(G, chis, flip=False):
"""(m, L*(0)) per character, the trivial one included, and C3's move."""
D = G.D
if D > 0:
A = sqrt(D) / pi
if G.k0 != G.e:
odd = np.real(chis[:, G.k0]) < 0
else:
odd = np.zeros(len(chis), bool)
if flip:
odd = ~odd
kinds = ["o" if o else "e" for o in odd]
m = np.where(odd, 0, 2)
else:
A = sqrt(-D) / (2 * pi)
kinds = ["c"] * len(chis)
m = np.ones(len(chis), int)
N = int(ceil((NA if D > 0 else 2 * NA) * A)) + 5 # W_C decays as e^-y
lam = lam1(G, chis, kinds, A, N)
move = np.max(np.abs(lam1(G, chis, kinds, A, int(1.5 * N)) - lam)
/ np.maximum(np.abs(lam), 1e-300))
# even: Gamma(s/2)^2 ~ 4/s^2; odd: the sum carries pi Lambda, since
# Gamma_R(s + 1)^2 = pi^-(s+1) Gamma((s+1)/2)^2; imaginary: Gamma ~ 1/s
L = np.array([lam[i] / {"e": 4, "o": pi, "c": 1}[k]
for i, k in enumerate(kinds)])
return m, L, move
def l0(d):
return -sum(a * kronecker(d, a) for a in range(1, -d)) / -d
def l0p(d):
return -0.5 * sum(kronecker(d, a) * log(sin(pi * a / d))
for a in range(1, d))
def genus_split(D):
"""D = d1 d2, d1 the fundamental factor of least |d1| with d2
fundamental, or None."""
best = None
for d1 in range(-abs(D), abs(D) + 1):
if d1 in (0, 1) or D % d1 or not P.fundamental(d1):
continue
d2 = D // d1
if d2 == 1 or not P.fundamental(d2):
continue
if D < 0 and d1 < 0:
continue
if best is None or abs(d1) < abs(best[0]):
best = (d1, d2)
return best
def read_field(G, flip=False):
"""Per cut: level - 1 and E per class, and the field's L* data."""
raw, cor = P.counts(G)
chis = characters(G)
m, L, move = lstar(G, chis, flip)
nontriv = [i for i in range(G.h)
if np.max(np.abs(chis[i] - 1)) > 1e-6]
out = []
for j, x in enumerate(P.CUTS):
tot = sum(cor[j])
lx = log(log(x)) + GAMMA
term = np.array([-m[i] * lx + log(abs(L[i])) for i in nontriv])
E = np.real(np.conj(chis[nontriv]).T @ term) / tot
dev = [cor[j][c] / (tot / G.h) - 1 for c in range(G.h)]
out.append((dev, E))
return out, chis, m, L, move, nontriv
def cellmean(v, S):
return sum(v[c] for c in S) / len(S)
def main():
t0 = time.time()
section("A THE CONTROLS")
worst1, n1, worst2, nchar, worst3, skip = 0.0, 0, 0.0, 0, 0.0, 0
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 = characters(G)
nchar += 1
tab = np.array([[G.mul(a, b) for b in range(G.h)]
for a in range(G.h)])
hom = np.max(np.abs(chis[:, tab] - chis[:, :, None]
* chis[:, None, :]))
orth = np.max(np.abs(chis @ np.conj(chis).T - G.h * np.eye(G.h)))
worst2 = max(worst2, hom, orth)
sp = genus_split(D)
if sp is None:
continue
d1, d2 = sp
test = [p for p in P.SMALL[1:200] if P.kron(D, p) == 1]
want = np.array([P.kron(d1, p) for p in test])
cols = [G.place(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:
skip += 1
continue
m, L, move = lstar(G, chis)
worst3 = max(worst3, move)
known = ((l0(d1) if d1 < 0 else l0p(d1))
* (l0(d2) if d2 < 0 else l0p(d2)))
got = L[hit[0]]
worst1 = max(worst1, abs(abs(got) - abs(known)) / abs(known))
n1 += 1
check("C1 genus characters: |L*(0)| as the closed forms",
worst1 < 1e-8 and n1 > 100 and skip == 0,
f"{n1} fields with |D| <= {C1_D}, worst relative {worst1:.2e},"
f" genus character not found at {skip}")
check("C2 the characters are homomorphisms and orthogonal",
worst2 < 1e-8, f"{nchar} groups, worst {worst2:.2e}")
check("C3 Lambda(1) settled at N", worst3 < 1e-10,
f"worst relative move {worst3:.2e}")
for sign, name in ((+1, "real, narrow"), (-1, "imaginary")):
section(f"B {name.upper()}")
lo, hi = ((5, REAL_D + 1) if sign > 0 else (-IMAG_D, 0))
rows, flipped, meta = [], [], []
for D in range(lo, hi):
if not P.fundamental(D):
continue
G = P.Group(D)
if G.h < 2:
continue
cs = P.cells(G)
out = read_field(G)[0]
rows.append([{k: (cellmean(dv, S), cellmean(E, S))
for k, S in cs.items()} for dv, E in out])
meta.append((abs(D), G.h, len(G.two) == G.h))
if sign > 0 and "K0" in cs:
fo = read_field(G, flip=True)[0]
flipped.append((out[0][0][G.k0] - out[0][1][G.k0])
- (fo[0][0][G.k0] - fo[0][1][G.k0]))
print(f" {len(rows)} fields with h >= 2 [{time.time() - t0:.0f} s]")
names = ["T", "K0", "A", "S", "G"] if sign > 0 else ["T", "A", "S",
"G"]
verdict = {}
for k in names:
print(f" cell {k}: cut level-1 E"
" residual")
for j, x in enumerate(P.CUTS):
v = [r[j][k] for r in rows if k in r[j]]
if not v:
continue
d = np.array([a for a, _ in v])
e = np.array([b for _, b in v])
res = d - e
se = lambda a: a.std(ddof=1) / sqrt(len(a)) # noqa: E731
print(f" {x:6d} {d.mean():+.4f}+-{se(d):.4f}"
f" {e.mean():+.4f}+-{se(e):.4f}"
f" {res.mean():+.4f}+-{se(res):.4f} ({len(v)})")
verdict[(k, x)] = (d.mean(), res.mean())
print(" Q3 (added after the run, printed) the residual on T and K0"
" by |D| band and by h, at 250 and 10^4:")
EL = [False] # the field's group is elementary 2, every chi genus
bands = (("|D| <= 4000", lambda a, h: a <= 4000),
("4000 < |D| <= 8000", lambda a, h: 4000 < a <= 8000),
("|D| > 8000", lambda a, h: a > 8000),
("h = 2", lambda a, h: h == 2),
("h = 3", lambda a, h: h == 3),
("h >= 4", lambda a, h: h >= 4),
("h >= 4, all genus", lambda a, h: h >= 4 and EL[0]),
("h >= 4, not", lambda a, h: h >= 4 and not EL[0]))
for bn, test in bands:
line = f" {bn:20s}"
for k in (("T", "K0") if sign > 0 else ("T",)):
for j in (0, len(P.CUTS) - 1):
v = np.array([r[j][k][0] - r[j][k][1]
for r, (a, h, el) in zip(rows, meta)
if k in r[j] and not EL.__setitem__(0, el)
and test(a, h)])
if len(v) > 1:
line += (f" {k}@{P.CUTS[j]} {v.mean():+.4f}"
f"+-{v.std(ddof=1) / sqrt(len(v)):.4f}"
f" ({len(v)})")
print(line)
if sign > 0:
q1 = all(abs(verdict[(k, x)][1]) < abs(verdict[(k, x)][0])
for k in ("T", "K0") for x in P.CUTS)
print(f" Q1 the residual smaller than the departure on T and"
f" K0 at every cut: {q1}")
fl = np.array(flipped)
check("C4 THE BREAK: k0's parity flipped moves K0's residual"
" at 250 by > 0.05", abs(fl.mean()) > 0.05,
f"mean move {fl.mean():+.4f} over {len(fl)} fields")
ok = (all(abs(verdict[(k, 250)][1]) <= 0.01 for k in ("T", "K0"))
and all(abs(verdict[(k, 10000)][1]) <= 0.005
for k in ("T", "K0")))
print(" KILL CARRIES %s: residual T %+.4f, K0 %+.4f at 250;"
" T %+.4f, K0 %+.4f at 10^4"
% ("MISSED" if ok else "FIRED",
verdict[("T", 250)][1], verdict[("K0", 250)][1],
verdict[("T", 10000)][1], verdict[("K0", 10000)][1]))
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()