#!/usr/bin/env python3 """verify_hv4.py -- checker for hv4.json (Hulek-Verrill fourfold, S_6-invariant sector; main text Section 4, SM Section S7). Python 3 standard library only, exact integer/rational arithmetic. Run in this directory: python3 verify_hv4.py # full run: includes the re-enumeration of all fluxes in the cut and of their # classes under the residual monodromy (about half a minute) python3 verify_hv4.py --quick # everything except the enumeration (seconds) Exit status 0 iff every check passes. What is checked: Gram5 (det 6480, signature (3,2)), W5, the polynomials N_flux = Q(g)/2, u = Gram5 g, n_2 and the on-cone factorization; M_d preserves Gram5 and is unipotent of index 5; the Hodge-number identities and chi/24 = 30; the cut 751 = floor(60/mu_lo) from the stored dyadic interval and the component caps (27,11,7,11,27); the box R; N_flux, n_2, u and r = g5/g1 of every named flux (the charge-4, 3, 2, 1 fluxes, the r = 5/2 edge flux, the 28 class representatives, the three vanishing classes); that the certified positions lie inside/above R as stated and their distances to 1/64 and to the box top; the frame identities (C integral, det C = 6, C^T Sigma_th C = Gram5; the printed monodromies preserve Sigma_th and multiply to 1; their pull-backs are integral, preserve Gram5, M_g[0] = M_d, and at each conifold point M_g is the reflection in the vanishing class with (2/Q) Gram5 g0 integral); the L5 data (symbol factorization, indicial polynomial theta^5 at 0, the holomorphic series c_n = sum (n!/prod n_i!)^2 and its annihilation by L5 through n = 40); the count identities of Table S12, and (full mode) an independent re-enumeration of the 3,509,207 / 59,360 / 55,516 / 3,872 / 230 / 12,082 signed flux vectors and of the 58,388 / 626 / 1,994 / 1,428 / 4,048 = 3,818 + 230 classes under ; the special-point catalogue (eleven points, only 1/64 in R). The Krawczyk certificates (positions, radii, |W| bounds), mu_lo, the field-sign table and the block margins are period-level results quoted from the computation and are not re-derived here. """ import json, os, sys, math, time from fractions import Fraction as Fr HERE = os.path.dirname(os.path.abspath(__file__)) D = json.load(open(os.path.join(HERE, "hv4.json"))) QUICK = "--quick" in sys.argv[1:] FAILS = []; NP = 0 def check(name, ok, detail=""): global NP NP += bool(ok) if not ok: FAILS.append(name) print(("PASS " if ok else "FAIL ") + name + (f" [{detail}]" if detail != "" else ""), flush=True) def T(M): return [list(r) for r in zip(*M)] def mul(A, B): Bt = T(B); return [[sum(a * b for a, b in zip(r, c)) for c in Bt] for r in A] def mv(A, v): return [sum(a * b for a, b in zip(r, v)) for r in A] def eye(n): return [[int(i == j) for j in range(n)] for i in range(n)] def madd(A, B, s=1): return [[a + s * b for a, b in zip(r, q)] for r, q in zip(A, B)] def mpow(M, e): R = eye(len(M)) while e: if e & 1: R = mul(R, M) M = mul(M, M); e >>= 1 return R def det(M): A = [[Fr(x) for x in r] for r in M]; n = len(A); d = Fr(1) for c in range(n): p = next((r for r in range(c, n) if A[r][c] != 0), None) if p is None: return Fr(0) if p != c: A[c], A[p] = A[p], A[c]; d = -d d *= A[c][c] for r in range(c + 1, n): f = A[r][c] / A[c][c] if f: A[r] = [a - f * b for a, b in zip(A[r], A[c])] return d def inertia(Q): n = len(Q); A = [[Fr(x) for x in r] for r in Q]; pos = neg = zero = 0; active = list(range(n)) while active: i = next((k for k in active if A[k][k] != 0), None) if i is None: pair = next(((a, b) for a in active for b in active if a < b and A[a][b] != 0), None) if pair is None: zero += len(active); break a, b = pair for k in range(n): A[a][k] += A[b][k] for k in range(n): A[k][a] += A[k][b] continue p = A[i][i]; pos += p > 0; neg += p < 0; row = A[i][:] for r in active: if r != i and A[r][i] != 0: f = A[r][i] / p for c in active: A[r][c] -= f * row[c] for c in active: A[i][c] = Fr(0) for r in active: A[r][i] = Fr(0) active.remove(i) return (pos, neg, zero) t0 = time.time() IL = D["invariant_lattice"]; G5 = IL["Gram5"]; W5 = IL["W5_diag"]; Md = IL["M_d"] def Q(g): return sum(G5[i][j] * g[i] * g[j] for i in range(5) for j in range(5)) def n2(g): return sum(W5[i] * g[i] * g[i] for i in range(5)) def u_of(g): return mv(G5, g) # ---- lattice check("Gram5 symmetric, det 6480 = 2^4 3^4 5, signature (3,2)", G5 == T(G5) and det(G5) == 6480 == IL["det_Gram5"] and inertia(G5) == (3, 2, 0) and IL["signature_Gram5"] == [3, 2]) check("W5 = diag(1,6,15,6,1) = orbit sizes (B^T B for the orbit-indicator basis)", W5 == [1, 6, 15, 6, 1] == IL["orbit_sizes"]) # polynomial identities: compare coefficient dictionaries of Q/2, u, n2 with the printed strings by evaluation on enough integer points import itertools, re def evalpoly(s, g): env = {f"g{i+1}": g[i] for i in range(5)}; env["N"] = g[0] return eval(s.replace("^", "**"), {"__builtins__": {}}, env) pts = [list(p) for p in itertools.product([-2, 1, 3], repeat=5)] check("N_flux = Q(g)/2 equals the printed polynomial " + IL["N_flux_polynomial"], all(Fr(Q(g), 2) == evalpoly(IL["N_flux_polynomial"], g) for g in pts)) ulist = [x.strip() for x in IL["u_of_g"].strip("[]").split(",")] check("u = Gram5 g equals the printed vector " + IL["u_of_g"], all(u_of(g) == [evalpoly(c, g) for c in ulist] for g in pts)) check("n_2(g) = g^T W5 g equals " + IL["n2_polynomial"], all(n2(g) == evalpoly(IL["n2_polynomial"], g) for g in pts)) check("on the cone g2 = g3 = 0: N_flux = g1 (g1 + g5)", all(Fr(Q([a, 0, 0, b, c]), 2) == a * (a + c) for a in range(-3, 4) for b in range(-2, 3) for c in range(-3, 4))) fam = IL["family"] check("family (N,1,0,0,3N): N_flux = 4N^2 - 30, u = N(5,0,30,0,1) + (0,-60,0,6,0)", all(Fr(Q([N, 1, 0, 0, 3 * N]), 2) == 4 * N * N - 30 and u_of([N, 1, 0, 0, 3 * N]) == [5 * N, -60, 30 * N, 6, N] for N in range(-5, 6))) I5 = eye(5); Nm = madd(Md, I5, -1) check("M_d integral, det 1, M_d^T Gram5 M_d = Gram5, (M_d - 1)^5 = 0 != (M_d - 1)^4", det(Md) == 1 and mul(mul(T(Md), G5), Md) == G5 and mpow(Nm, 5) == [[0] * 5] * 5 and mpow(Nm, 4) != [[0] * 5] * 5) # ---- Hodge, cut, box H = D["fourfold"]["hodge_numbers"] check("chi = 4 + 2 h11 - 4 h21 + 2 h31 + h22 = 720; h22 = 44 + 4 h11 - 2 h21 + 4 h31 = 492; chi/24 = 30", 4 + 2 * H["h11"] - 4 * H["h21"] + 2 * H["h31"] + H["h22"] == H["chi"] == 720 and 44 + 4 * H["h11"] - 2 * H["h21"] + 4 * H["h31"] == H["h22"] and H["chi"] // 24 == 30 == D["fourfold"]["chi_over_24"] and H["chi"] % 24 == 0) check("horizontal lattice rank 1+6+15+6+1 = 29, det 2^16 3^6 = 47775744", sum(IL["orbit_sizes"]) == 29 == D["fourfold"]["horizontal_lattice"]["rank"] and 2 ** 16 * 3 ** 6 == D["fourfold"]["horizontal_lattice"]["det_value"] == 47775744) CB = D["cut_and_box"]; lo = Fr(CB["mu_lo_dyadic_interval"]["lower"]); hi = Fr(CB["mu_lo_dyadic_interval"]["upper"]) check("cut: floor(60/mu) = 751 at both ends of the stored dyadic interval for mu_lo (60 = 2 chi/24)", (60 / lo).__floor__() == 751 == (60 / hi).__floor__() == CB["n2_cut"] and lo.denominator == 2 ** 56 and 60 == 2 * H["chi"] // 24) check("component caps floor(sqrt(751/w_i)) = (27,11,7,11,27)", [math.isqrt(751 // w) for w in W5] == CB["component_caps"] == [27, 11, 7, 11, 27]) Re_lo, Re_hi = map(Fr, CB["box_R"]["Re_phi"]); Im_lo, Im_hi = map(Fr, CB["box_R"]["Im_phi"]) check("box R: Re phi in [13/1024, 19/1024], |Im phi| <= 1/256; contains 1/64", Re_lo == Fr(13, 1024) and Re_hi == Fr(19, 1024) and -Im_lo == Im_hi == Fr(1, 256) and Re_lo < Fr(1, 64) < Re_hi and CB["N_flux_range"] == [1, 30]) # ---- L5 PF = D["picard_fuchs_L5"]; Qs = PF["Q_theta"] def polymul(a, b): r = [0] * (len(a) + len(b) - 1) for i, x in enumerate(a): for j, y in enumerate(b): r[i + j] += x * y return r sym = [q[5] for q in Qs] check("symbol (theta^5 coefficients) = [1,-56,784,-2304] = (1-4phi)(1-16phi)(1-36phi); leading term phi^0 Q_0 = theta^5 (MUM at 0); singular points {0,1/36,1/16,1/4,inf}", sym == PF["symbol"] == [1, -56, 784, -2304] and polymul(polymul([1, -4], [1, -16]), [1, -36]) == sym and Qs[0] == [0, 0, 0, 0, 0, 1] and PF["singular_points"] == ["0", "1/36", "1/16", "1/4", "inf"]) # holomorphic series c_n = sum_{n1+..+n6=n} (n!/prod n_i!)^2, computed by adding one letter at a time: c^(k)_n = sum_j C(n,j)^2 c^(k-1)_(n-j) NMAX = 40 ck = [1] * (NMAX + 1) # one letter: (n!/n!)^2 = 1 for k in range(5): # add the remaining five letters: c'_n = sum_j C(n,j)^2 c_(n-j) ck = [sum(math.comb(n, j) ** 2 * ck[n - j] for j in range(n + 1)) for n in range(NMAX + 1)] c = ck check("holomorphic series c_0..c_40 as printed (1, 6, 66, 996, 18306, ...)", [str(x) for x in c] == PF["holomorphic_series_c0_c40"]) def Lc(n): # coefficient of phi^n in L5 applied to sum_m c_m phi^m: sum_j Q_j(n-j) c_{n-j} (phi^j Q_j(theta) phi^m = Q_j(m) phi^(m+j)) return sum(sum(Qs[j][i] * (n - j) ** i for i in range(6)) * c[n - j] for j in range(4) if n - j >= 0) check("L5 annihilates the holomorphic series through order 40 (recurrence sum_j Q_j(n-j) c_(n-j) = 0, n = 0..40)", all(Lc(n) == 0 for n in range(NMAX + 1))) # local exponents at a conifold point s: indicial polynomial of L5 at phi = s from the leading D-form; we check instead the stated sets are as printed and consistent (sum of exponents rule not applied) check("local exponents as printed: {0,0,0,0,0} at 0; {0,1,3/2,2,3} at 1/36, 1/16, 1/4; {1,1,3/2,2,2} at infinity", PF["local_exponents"]["0"] == ["0"] * 5 and all(sorted(PF["local_exponents"][p], key=Fr) == ["0", "1", "3/2", "2", "3"] for p in ("1/36", "1/16", "1/4")) and sorted(PF["local_exponents"]["inf"], key=Fr) == ["1", "1", "3/2", "2", "2"]) # ---- frame F = D["frame"]; C = F["C"]; Sg = F["Sigma_th"] check("C integral, det C = 6, C^T Sigma_th C = Gram5 (index-6 sublattice: det Sigma_th * 36 = det Gram5)", det(C) == 6 == F["det_C"] and mul(mul(T(C), Sg), C) == G5 and det(Sg) * 36 == det(G5) and F["det_Sigma_th"] == det(Sg)) Mth = F["M_th"] prod = eye(5) for k in ("0", "1/36", "1/16", "1/4", "inf"): prod = mul(prod, Mth[k]) check("printed monodromies preserve Sigma_th (M^T Sigma M = Sigma), det +-1, and M_0 M_1/36 M_1/16 M_1/4 M_inf = 1", all(mul(mul(T(Mth[k]), Sg), Mth[k]) == Sg and abs(det(Mth[k])) == 1 for k in Mth) and prod == I5) Cinv = [[Fr(x) for x in r] for r in C] # invert C exactly def inv(M): n = len(M); A = [[Fr(x) for x in r] + [Fr(int(i == j)) for j in range(n)] for i, r in enumerate(M)] for cc in range(n): p = next(r for r in range(cc, n) if A[r][cc] != 0); A[cc], A[p] = A[p], A[cc] piv = A[cc][cc]; A[cc] = [x / piv for x in A[cc]] for r in range(n): if r != cc and A[r][cc] != 0: f = A[r][cc]; A[r] = [a - f * b for a, b in zip(A[r], A[cc])] return [r[n:] for r in A] Ci = inv(C) okp = True for k in Mth: Mg = mul(mul(Ci, [[Fr(x) for x in r] for r in Mth[k]]), [[Fr(x) for x in r] for r in C]) okp = okp and Mg == [[Fr(x) for x in r] for r in F["M_g"][k]] and mul(mul(T(F["M_g"][k]), G5), F["M_g"][k]) == G5 check("pulled-back monodromies M_g = C^-1 M_th C are integral, equal the printed ones, preserve Gram5; M_g[0] = M_d", okp and F["M_g"]["0"] == Md) okr = True for vc in F["vanishing_classes"]: g0 = vc["g0"]; Qg = Q(g0); w = [Fr(2 * x, Qg) for x in u_of(g0)] refl = [[Fr(int(i == j)) - w[j] * g0[i] for j in range(5)] for i in range(5)] # x -> x - (2/Q)(g0^T Gram5 x) g0 ... as matrix acting on columns: R = I - g0 (2/Q) (Gram5 g0)^T okr = okr and Qg == vc["Q"] and Fr(Qg, 2) == vc["N_flux"] and all(x.denominator == 1 for x in w) and [int(x) for x in w] == vc["two_over_Q_Gram5_g0"] and refl == [[Fr(x) for x in r] for r in F["M_g"][vc["point"]]] and mv(F["M_g"][vc["point"]], g0) == [-x for x in g0] check("at 1/36, 1/16, 1/4: M_g is the reflection in the vanishing class g0 = (1,0,0,0,0), (6,1,0,0,0), (15,5,1,0,0) with Q = 2, 12, 30 (N_flux 1, 6, 15), (2/Q) Gram5 g0 integral, M_g g0 = -g0", okr) # ---- named fluxes and vacua def r_of(g): return Fr(g[4], g[0]) cv = {c["name"]: c for c in D["certified_vacua"]} okv = True for c in D["certified_vacua"]: g = c["g"]; okv = okv and Fr(Q(g), 2) == c["N_flux"] and n2(g) == c["n2"] <= 751 and u_of(g) == c["u"] and g[1] == g[2] == 0 check("certified vacua: g = (1,0,0,0,3), (1,0,0,0,2), (1,0,0,0,1) have N_flux 4, 3, 2, n2 10, 5, 2 (inside the cut), u = (5,0,30,0,1), (4,0,30,0,1), (3,0,30,0,1), on the cone with r = 3, 2, 1", okv and [c["N_flux"] for c in D["certified_vacua"]] == [4, 3, 2] and [r_of(c["g"]) for c in D["certified_vacua"]] == [3, 2, 1]) okz = True for c in D["certified_vacua"]: z = Fr(c["Re_z_mid"]); okz = okz and abs(float(z - Fr(19, 1024)) - float(c["z_minus_box_top"])) < 1e-15 and abs(float(abs(z - Fr(1, 64))) - float(c["dist_to_1_64"])) < 1e-15 inside = Re_lo <= z <= Re_hi okz = okz and (inside == (not c["above_R_certified"])) and c["Im_enclosure_contains_0"] and c["W_nonzero"] and float(c["absW_lower_bound"]) > 0 and float(c["e_minus_K"]) > 0 okz = okz and float(c["Re_z_rad"]) < 1e-22 and float(c["Im_z_rad"]) < 1e-22 and (abs(float(z - Fr(19, 1024))) > 100 * float(c["Re_z_rad"])) check("certified positions: charge 4 inside R (z* - 19/1024 = -0.00127...), charges 3 and 2 above R by 0.00202... and 0.00608... (far outside the radii < 1e-22); Im-enclosures contain 0; |W| > 0 and e^-K > 0 as stored; distances to 1/64 and to the box top consistent with the stored midpoints", okz) V = D["vacua_in_R"]; rows = V["rows"] okr = all(Fr(Q(s["g"]), 2) == s["N_flux"] and Q(s["g"]) == s["Q"] and n2(s["g"]) == s["n2"] <= 751 and 1 <= s["N_flux"] <= 30 and s["g"][1] == s["g"][2] == 0 and s["W_nonzero_certified"] for s in rows) inR = all(Re_lo <= Fr(s["Re_z"]) <= Re_hi and Im_lo <= Fr(s["Im_z"]) <= Im_hi for s in rows) from collections import Counter per = Counter(s["N_flux"] for s in rows) check("the 28 vacua in R: 28 rows, all on the cone, N_flux and n2 of each representative as stated, all positions inside R, per-charge counts {4:2, 5:2, 16:6, 18:6, 20:6, 22:6}, 12 distinct positions, all W != 0", len(rows) == 28 == V["n_vacua"] and okr and inR and dict(per) == {4: 2, 5: 2, 16: 6, 18: 6, 20: 6, 22: 6} and {str(k): v for k, v in per.items()} == V["per_Nflux"] and len(V["positions"]) == 12 == V["n_positions"] and len({(round(float(s["Re_z"]), 12), round(abs(float(s["Im_z"])), 12) if abs(float(s["Im_z"])) > 1e-20 else 0.0, float(s["Im_z"]) > 1e-20) for s in rows}) == 12) check("max ball radius %.4g <= 2.2e-22; min |z* - 1/64| = %.4g > 0 (no vacuum at the special point 1/64); |W| lower bounds in [%s, %s]" % (float(V["max_ball_radius"]), float(V["min_dist_to_1_64"]), V["absW_lower_bound_range"][0][:8], V["absW_lower_bound_range"][1][:9]), float(V["max_ball_radius"]) == max(float(s["ball_radius"]) for s in rows) <= 2.2e-22 and float(V["min_dist_to_1_64"]) == min(float(s["dist_to_1_64"]) for s in rows) > 1e-4 and all(abs(abs(float(Fr(s["Re_z"]) - Fr(1, 64)) if abs(float(s["Im_z"])) < 1e-20 else math.hypot(float(s["Re_z"]) - 1 / 64, float(s["Im_z"]))) - float(s["dist_to_1_64"])) < 1e-12 for s in rows)) c1 = D["charge_1"] check("charge 1: g = (1,0,0,0,0) has N_flux 1, n2 1, u = (2,0,30,0,1), r = 0", Fr(Q(c1["g"]), 2) == 1 == c1["N_flux"] and n2(c1["g"]) == 1 and u_of(c1["g"]) == c1["u"] == [2, 0, 30, 0, 1]) e = D["r52_edge_zero"] check("r = 5/2 edge flux (2,0,0,0,5): N_flux 14, n2 29, u = (9,0,60,0,2); stored enclosure lies above 19/1024", Fr(Q(e["g"]), 2) == 14 and n2(e["g"]) == 29 and u_of(e["g"]) == e["u"] == [9, 0, 60, 0, 2] and Fr(e["Re_z_interval"][0]) - Fr(19, 1024) > 0 and abs(float(Fr(e["Re_z_interval"][0]) - Fr(19, 1024)) - float(e["z_minus_box_top_interval"][0])) < 1e-20) check("r = 2 ray: u = (4,0,30,0,1) = Gram5 (1,0,0,0,2) = (1/2) Gram5 (2,0,0,0,4)", u_of([1, 0, 0, 0, 2]) == D["r2_ray"]["u"] == [4, 0, 30, 0, 1] and u_of([2, 0, 0, 0, 4]) == [8, 0, 60, 0, 2]) # the 230 on-cone list onc = [] for g4 in range(-11, 12): for s in (1, -1): for g5 in (0, 1, 2): onc.append((s, 0, 0, g4, s * g5)) onc.append((2 * s, 0, 0, g4, -s)); onc.append((3 * s, 0, 0, g4, -2 * s)) onc = set(onc) check("the predicted on-cone list at N_flux <= 3 has 230 members = 2*23*3 + 2*23 + 2*23, all with 1 <= N_flux <= 3, n2 <= 751, r <= 2", len(onc) == 230 == D["on_cone_Nflux_le_3"]["n"] and all(1 <= Fr(Q(g), 2) <= 3 and n2(g) <= 751 and r_of(g) <= 2 for g in onc)) # ---- counts CA = D["counts"]["unit_A_signed_flux_vectors"]; CBc = D["counts"]["unit_B_classes_under_Md"] check("count identities: 4102 = 3872 + 230; per-charge 626 + 1996 + 1480 = 4102; classes 4048 = 3818 + 230 = 626 + 1994 + 1428; 58388 = 28 + 58360; sum of per-charge classes = 58388", CA["Nflux_le_3"] == CA["offcone_Nflux_le_3"] + CA["oncone_Nflux_le_3"] == sum(CA["per_Nflux_1_2_3"]) == 4102 and CBc["Nflux_le_3"] == sum(CBc["Nflux_le_3_oncone_offcone"]) == 4048 == sum(CBc["per_Nflux_1_to_30"][k] for k in "123") and CBc["in_cut"] == CBc["with_vacuum_in_R"] + CBc["without_vacuum_in_R"] == 58388 == sum(CBc["per_Nflux_1_to_30"].values()) and CBc["Nflux_le_3_oncone_offcone"][0] == CA["oncone_Nflux_le_3"] == 230) if not QUICK: tt = time.time(); caps = CB["component_caps"] n_ell = 0; adm = [] r1 = range(-caps[0], caps[0] + 1) for g2 in range(-caps[1], caps[1] + 1): a2 = 6 * g2 * g2 for g3 in range(-caps[2], caps[2] + 1): a3 = a2 + 15 * g3 * g3 if a3 > 751: continue for g4 in range(-caps[3], caps[3] + 1): a4 = a3 + 6 * g4 * g4 if a4 > 751: continue rem = 751 - a4 for g1 in r1: a1 = g1 * g1 if a1 > rem: continue m5 = math.isqrt(rem - a1) n_ell += 2 * m5 + 1 base = 2 * g1 * g1 - 60 * g2 * g2 + 180 * g3 * g3 + 60 * g1 * g3 + 12 * g2 * g4 # Q without the g5 term 2 g1 g5 for g5 in range(-m5, m5 + 1): q = base + 2 * g1 * g5 if 2 <= q <= 60: adm.append((g1, g2, g3, g4, g5)) A = set(adm) Bv = lambda g: (2 * g[0] * g[0] - 60 * g[1] * g[1] + 180 * g[2] * g[2] + 60 * g[0] * g[2] + 12 * g[1] * g[3] + 2 * g[0] * g[4]) // 2 off = [g for g in adm if (g[1], g[2]) != (0, 0)] per3 = Counter(Bv(g) for g in adm if Bv(g) <= 3) check("enumeration of signed flux vectors: %d with n2 <= 751; %d in the cut; off-cone %d (N_flux <= 30), %d (N_flux <= 3), %d (|g1| <= 2); on-cone N_flux <= 3: %d (= the predicted list); per-charge %s; stated in hv4.json: %s (%.0fs)" % (n_ell, len(adm), len(off), sum(1 for g in off if Bv(g) <= 3), sum(1 for g in off if abs(g[0]) <= 2), sum(1 for g in adm if Bv(g) <= 3 and (g[1], g[2]) == (0, 0)), [per3[1], per3[2], per3[3]], [CA["ellipsoid_points_n2_le_751"], CA["in_cut"], CA["offcone_Nflux_le_30"], CA["offcone_Nflux_le_3"], CA["offcone_Nflux_le_30_g1_abs_le_2"], CA["oncone_Nflux_le_3"], CA["per_Nflux_1_2_3"]], time.time() - tt), n_ell == CA["ellipsoid_points_n2_le_751"] == 3509207 and len(adm) == CA["in_cut"] == 59360 and len(off) == CA["offcone_Nflux_le_30"] == 55516 and sum(1 for g in off if Bv(g) <= 3) == CA["offcone_Nflux_le_3"] == 3872 and {g for g in adm if Bv(g) <= 3 and (g[1], g[2]) == (0, 0)} == onc and sum(1 for g in off if abs(g[0]) <= 2) == CA["offcone_Nflux_le_30_g1_abs_le_2"] == 12082 and [per3[1], per3[2], per3[3]] == CA["per_Nflux_1_2_3"]) # classes under : union g ~ M_d^k g whenever both are in the cut. All k with n2(M_d^k g) <= 751 are found from the exact polynomial p(k) = n2(M_d^k g) (degree <= 8) and a root bound. tt = time.time() Nm2 = mul(Nm, Nm); Nm3 = mul(Nm2, Nm); Nm4 = mul(Nm3, Nm) Mdinv = [[int(x) for x in r] for r in inv(Md)] parent = {g: g for g in adm} def find(x): while parent[x] != x: parent[x] = parent[parent[x]]; x = parent[x] return x maxk = 0 for g in adm: v1 = mv(Nm, g) if not any(v1): continue # fixed by M_d: singleton v = [list(g), v1, mv(Nm2, g), mv(Nm3, g), mv(Nm4, g)] # M_d^k g = sum_j binom(k,j) v_j # p(k) = n2(sum_j binom(k,j) v_j): even-degree polynomial with positive leading coefficient; bound its real roots (of p(k) - 751) by Cauchy's bound on the expansion in powers of k # compute coefficients of p(k) - 751 in the monomial basis via exact interpolation on k = 0..8 vals = [] for k in range(9): w = [sum(math.comb(k, j) * v[j][i] for j in range(5)) for i in range(5)]; vals.append(n2(w) - 751) # Newton forward differences -> monomial coefficients coef = [Fr(0)] * 9; diffs = [Fr(x) for x in vals]; fall = [Fr(1)] # falling factorial basis newton = [] d_ = list(diffs) for i in range(9): newton.append(d_[0] / math.factorial(i)); d_ = [d_[j + 1] - d_[j] for j in range(len(d_) - 1)] # convert sum_i newton[i] * k(k-1)...(k-i+1) to monomials poly = [Fr(0)] * 9; basis = [Fr(1)] for i in range(9): for j, b in enumerate(basis): poly[j] += newton[i] * b nb = [Fr(0)] * (len(basis) + 1) # basis *= (k - i) for j, b in enumerate(basis): nb[j + 1] += b; nb[j] -= i * b basis = nb while poly and poly[-1] == 0: poly.pop() lead = poly[-1]; bound = 1 + max(abs(c_ / lead) for c_ in poly[:-1]) if len(poly) > 1 else 1 K = int(bound) + 1; maxk = max(maxk, K) w = list(g) for k in range(1, K + 1): w = mv(Md, w) if n2(w) <= 751: tw = tuple(w) if tw in parent: parent[find(tw)] = find(g) w = list(g) for k in range(1, K + 1): w = mv(Mdinv, w) if n2(w) <= 751: tw = tuple(w) if tw in parent: parent[find(tw)] = find(g) roots = Counter(find(g) for g in adm) ncls = len(roots) cls_bin = Counter(Bv(r) for r in roots) cls3 = [r for r in roots if Bv(r) <= 3] on3 = 0; off3 = 0 members = {} for g in adm: members.setdefault(find(g), []).append(g) for r in cls3: if any((m[1], m[2]) == (0, 0) for m in members[r]): on3 += 1 else: off3 += 1 single_on = all(len(members[find(g)]) == 1 for g in onc) check("classes under : %d in the cut; per charge N_flux = 1..30 as printed (1: %d, 2: %d, 3: %d); N_flux <= 3: %d = %d on-cone (all singletons: %s) + %d off-cone; stated in hv4.json: %d in the cut, [on, off] = %s (largest |k| scanned %d, %.0fs)" % (ncls, cls_bin[1], cls_bin[2], cls_bin[3], on3 + off3, on3, single_on, off3, CBc["in_cut"], CBc["Nflux_le_3_oncone_offcone"], maxk, time.time() - tt), ncls == CBc["in_cut"] == 58388 and all(cls_bin[int(k)] == v for k, v in CBc["per_Nflux_1_to_30"].items()) and [on3, off3] == CBc["Nflux_le_3_oncone_offcone"] == [230, 3818] and single_on) reps = {tuple(s["g"]) for s in rows} check("the 28 class representatives of the vacua in R are in the cut and lie in 28 distinct classes", reps <= A and len({find(g) for g in reps}) == 28) # ---- special points SPc = D["special_points"] def algval(s): s = s.replace("sqrt3", "*1.7320508075688772").replace("sqrt6", "*2.449489742783178").replace("(*", "(").replace(" ", "") s = re.sub(r'(?