Skip to content

Instantly share code, notes, and snippets.

@justinmeza

justinmeza/01_coframe_sector.json

Created April 1, 2026 17:33
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save justinmeza/a2dec3dbfcddb2bc4b82fac4ad24033a to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save justinmeza/a2dec3dbfcddb2bc4b82fac4ad24033a to your computer and use it in GitHub Desktop.
Download ZIP
Supplementary material: verification scripts for first-derivative jet-space invariants (paper4)
Raw
01_coframe_sector.json
{
"N_coframe_quartic": 55,
"N_torsion_quartic": 27,
"N_coframe_cubic": 0,
"N_torsion_cubic": 0,
"coframe_quadruples": [
"T·T·T·T",
"T·T·T·V",
"T·T·T·s",
"T·T·V·V",
"T·T·V·s",
"T·T·𝒜·𝒜",
"T·T·s·s",
"T·V·V·V",
"T·V·V·s",
"T·V·𝒜·𝒜",
"T·V·s·s",
"T·𝒜·𝒜·s",
"T·s·s·s",
"V·V·V·V",
"V·V·V·s",
"V·V·𝒜·𝒜",
"V·V·s·s",
"V·𝒜·𝒜·s",
"V·s·s·s",
"𝒜·𝒜·𝒜·𝒜",
"𝒜·𝒜·s·s",
"s·s·s·s"
],
"torsion_quadruples": [
"T·T·T·T",
"T·T·T·V",
"T·T·V·V",
"T·T·𝒜·𝒜",
"T·V·V·V",
"T·V·𝒜·𝒜",
"V·V·V·V",
"V·V·𝒜·𝒜",
"𝒜·𝒜·𝒜·𝒜"
],
"n_coframe_patterns": 12822,
"n_torsion_patterns": 11589,
"coframe_submatrix_shape": [
12822,
42046
],
"torsion_submatrix_shape": [
11589,
30636
],
"per_quadruple_coframe": {
"T·T·T·T": {
"n_patterns": 10395,
"rank": 21
},
"T·T·T·V": {
"n_patterns": 945,
"rank": 12
},
"T·T·T·s": {
"n_patterns": 945,
"rank": 14
},
"T·T·V·V": {
"n_patterns": 105,
"rank": 7
},
"T·T·V·s": {
"n_patterns": 105,
"rank": 9
},
"T·T·𝒜·𝒜": {
"n_patterns": 105,
"rank": 9
},
"T·T·s·s": {
"n_patterns": 105,
"rank": 9
},
"T·V·V·V": {
"n_patterns": 15,
"rank": 1
},
"T·V·V·s": {
"n_patterns": 15,
"rank": 2
},
"T·V·𝒜·𝒜": {
"n_patterns": 15,
"rank": 3
},
"T·V·s·s": {
"n_patterns": 15,
"rank": 3
},
"T·𝒜·𝒜·s": {
"n_patterns": 15,
"rank": 3
},
"T·s·s·s": {
"n_patterns": 15,
"rank": 1
},
"V·V·V·V": {
"n_patterns": 3,
"rank": 1
},
"V·V·V·s": {
"n_patterns": 3,
"rank": 1
},
"V·V·𝒜·𝒜": {
"n_patterns": 3,
"rank": 2
},
"V·V·s·s": {
"n_patterns": 3,
"rank": 2
},
"V·𝒜·𝒜·s": {
"n_patterns": 3,
"rank": 2
},
"V·s·s·s": {
"n_patterns": 3,
"rank": 1
},
"𝒜·𝒜·𝒜·𝒜": {
"n_patterns": 3,
"rank": 1
},
"𝒜·𝒜·s·s": {
"n_patterns": 3,
"rank": 2
},
"s·s·s·s": {
"n_patterns": 3,
"rank": 1
}
},
"per_quadruple_torsion": {
"T·T·T·T": {
"n_patterns": 10395,
"rank": 21
},
"T·T·T·V": {
"n_patterns": 945,
"rank": 12
},
"T·T·V·V": {
"n_patterns": 105,
"rank": 7
},
"T·T·𝒜·𝒜": {
"n_patterns": 105,
"rank": 9
},
"T·V·V·V": {
"n_patterns": 15,
"rank": 1
},
"T·V·𝒜·𝒜": {
"n_patterns": 15,
"rank": 3
},
"V·V·V·V": {
"n_patterns": 3,
"rank": 1
},
"V·V·𝒜·𝒜": {
"n_patterns": 3,
"rank": 2
},
"𝒜·𝒜·𝒜·𝒜": {
"n_patterns": 3,
"rank": 1
}
},
"sum_per_quad_coframe": 107,
"sum_per_quad_torsion": 57,
"primes": [
10007,
10009,
100003,
1000003,
10000019
],
"coframe_ranks_by_prime": {
"10007": 55,
"10009": 55,
"100003": 55,
"1000003": 55,
"10000019": 55
},
"torsion_ranks_by_prime": {
"10007": 27,
"10009": 27,
"100003": 27,
"1000003": 27,
"10000019": 27
},
"cubic_coframe_triples_surviving": 0,
"cubic_index_parity_argument": "T(3)+V(1)+A(1)+s(1): any triple of odd numbers sums to odd",
"elapsed_seconds": 297.5319130420685
}
Raw
01_coframe_sector.py
#!/usr/bin/env python3
"""01_coframe_sector.py — Extract coframe-only sector and compute global ranks.
From the paper3 quartic checkpoint (17,383 patterns × 66,042 monomials),
extract the submatrix of coframe-only patterns (blocks ∈ {T, V, A, s})
and compute:
N_torsion: global rank of pure-torsion quadruples {T, V, A}
N_coframe: global rank of all coframe quadruples {T, V, A, s}
Also verifies the cubic coframe sector is trivially empty (odd index count).
Reuses the rank_mod_p machinery from 20_quartic_symbolic.py.
"""
import numpy as np
from itertools import combinations_with_replacement, product as iproduct
from collections import defaultdict
import json
import sys
import os
import time
import pickle
import functools
print = functools.partial(print, flush=True)
D = 4
eta = np.array([[-1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]], dtype=int)
BLOCKS = [
("Shat", 2, +1),
("theta", 0, +1),
("W", 2, +1),
("T", 3, +1),
("V", 1, +1),
("A", 1, -1),
("s", 1, +1),
]
BLOCK_DISPLAY = {"Shat": "Ŝ", "theta": "θ", "W": "W", "T": "T",
"V": "V", "A": "𝒜", "s": "s"}
# Block indices
IDX = {name: i for i, (name, _, _) in enumerate(BLOCKS)}
# Sector classification
COFRAME_BLOCKS = {IDX["T"], IDX["V"], IDX["A"], IDX["s"]}
TORSION_BLOCKS = {IDX["T"], IDX["V"], IDX["A"]} # no scale
def generate_pairings(indices):
if len(indices) == 0:
return [()]
if len(indices) == 2:
return [((indices[0], indices[1]),)]
first = indices[0]
rest = indices[1:]
result = []
for i, partner in enumerate(rest):
remaining = rest[:i] + rest[i+1:]
for sub in generate_pairings(remaining):
result.append(((first, partner),) + sub)
return result
def rank_mod_p(mat, p, track_pivots=False):
"""Rank of dense int64 matrix mod p via Gaussian elimination."""
m, n = mat.shape
A = mat.astype(np.int64) % p
A[A < 0] += p
rank = 0
pivot_rows = []
pivot_cols = []
row_perm = list(range(m))
for col in range(n):
if rank >= m:
break
nz = np.nonzero(A[rank:, col])[0]
if len(nz) == 0:
continue
found = nz[0] + rank
if found != rank:
A[[rank, found]] = A[[found, rank]]
row_perm[rank], row_perm[found] = row_perm[found], row_perm[rank]
pivot_rows.append(row_perm[rank])
pivot_cols.append(col)
inv_piv = pow(int(A[rank, col]), p - 2, p)
col_vals = A[:, col].copy()
col_vals[rank] = 0
nz_rows = np.nonzero(col_vals)[0]
if len(nz_rows) > 0:
factors = (col_vals[nz_rows] * inv_piv) % p
A[nz_rows] = (A[nz_rows] - factors[:, np.newaxis] * A[rank]) % p
rank += 1
if rank % 20 == 0:
print(f" rank_mod_p: rank {rank}, col {col}/{n}")
if track_pivots:
return rank, pivot_rows, pivot_cols
return rank
def main():
t0 = time.time()
print("=" * 60)
print("COFRAME SECTOR EXTRACTION — QUARTIC & CUBIC")
print("=" * 60)
# ─── Reproduce the quadruple enumeration (must match paper3 exactly) ──
print("\nStep 1: Enumerate all surviving quadruples...")
surviving = []
for quad in combinations_with_replacement(range(len(BLOCKS)), 4):
parity = 1
total_idx = 0
for b in quad:
parity *= BLOCKS[b][2]
total_idx += BLOCKS[b][1]
if parity != +1:
continue
if total_idx % 2 != 0:
continue
surviving.append(quad)
print(f" Surviving quadruples: {len(surviving)}")
all_patterns = []
quad_info = [] # (quad, n_patterns, start_idx)
for quad in surviving:
n_idx = [BLOCKS[b][1] for b in quad]
total = sum(n_idx)
pairings = generate_pairings(list(range(total))) if total > 0 else [()]
start = len(all_patterns)
for p in pairings:
all_patterns.append((quad, p))
quad_info.append((quad, len(pairings), start))
N_pat = len(all_patterns)
print(f" Total contraction patterns: {N_pat}")
# Classify quadruples
coframe_quads = [] # indices into quad_info
torsion_quads = []
for qi, (quad, n_pats, start) in enumerate(quad_info):
blocks_set = set(quad)
if blocks_set <= COFRAME_BLOCKS:
coframe_quads.append(qi)
if blocks_set <= TORSION_BLOCKS:
torsion_quads.append(qi)
coframe_quad_names = []
torsion_quad_names = []
for qi in coframe_quads:
quad = quad_info[qi][0]
name = "·".join(BLOCK_DISPLAY[BLOCKS[b][0]] for b in quad)
coframe_quad_names.append(name)
for qi in torsion_quads:
quad = quad_info[qi][0]
name = "·".join(BLOCK_DISPLAY[BLOCKS[b][0]] for b in quad)
torsion_quad_names.append(name)
print(f"\n Coframe quadruples ({len(coframe_quads)}):")
for name in coframe_quad_names:
print(f" {name}")
print(f"\n Pure-torsion quadruples ({len(torsion_quads)}):")
for name in torsion_quad_names:
print(f" {name}")
# Collect row indices for each sector
coframe_rows = []
torsion_rows = []
for qi in coframe_quads:
quad, n_pats, start = quad_info[qi]
coframe_rows.extend(range(start, start + n_pats))
for qi in torsion_quads:
quad, n_pats, start = quad_info[qi]
torsion_rows.extend(range(start, start + n_pats))
print(f"\n Coframe patterns: {len(coframe_rows)}")
print(f" Pure-torsion patterns: {len(torsion_rows)}")
# ─── Load checkpoint ──────────────────────────────────────────────
print("\nStep 2: Loading quartic checkpoint...")
checkpoint_file = os.path.join(
os.path.dirname(__file__), "..", "..", "paper3", "verify",
"quartic_checkpoints", "quartic_polys.pkl"
)
checkpoint_file = os.path.normpath(checkpoint_file)
if not os.path.isfile(checkpoint_file):
print(f" ERROR: Checkpoint not found at {checkpoint_file}")
sys.exit(1)
with open(checkpoint_file, "rb") as f:
saved = pickle.load(f)
all_polys = saved["all_polys"]
all_monomials = saved["all_monomials"]
N_mon = len(all_monomials)
print(f" Loaded {len(all_polys)} polynomials, {N_mon} monomials")
# Build monomial index
monomial_list = sorted(all_monomials)
mon_idx = {m: i for i, m in enumerate(monomial_list)}
# ─── Identify de-only monomials ───────────────────────────────────
# Coframe blocks (T, V, A, s) are all functions of de variables only.
# Their products can only contain de monomials. But let's be explicit:
# we restrict to columns that are de-only (no 'du' variables).
print("\nStep 3: Identifying de-only monomials...")
de_only_cols = []
for i, mono in enumerate(monomial_list):
if all(var[0] == 'de' for var in mono):
de_only_cols.append(i)
print(f" de-only monomials: {len(de_only_cols)} / {N_mon}")
de_col_set = set(de_only_cols)
de_col_remap = {old: new for new, old in enumerate(de_only_cols)}
# ─── Build submatrices ────────────────────────────────────────────
print("\nStep 4: Building coframe submatrix...")
n_de = len(de_only_cols)
# For coframe-only patterns, ALL monomials should be de-only.
# Verify this and build the submatrix.
n_coframe_rows = len(coframe_rows)
C_coframe = np.zeros((n_coframe_rows, n_de), dtype=np.int64)
n_du_contamination = 0
for ri, row_idx in enumerate(coframe_rows):
poly = all_polys[row_idx]
for mono, coeff in poly.items():
if coeff == 0:
continue
col = mon_idx[mono]
if col in de_col_set:
C_coframe[ri, de_col_remap[col]] = coeff
else:
n_du_contamination += 1
if n_du_contamination > 0:
print(f" WARNING: {n_du_contamination} du-contaminated entries in coframe rows!")
else:
print(f" Verified: all coframe pattern monomials are de-only ✓")
# Remove zero columns
nz_cols_mask = np.any(C_coframe != 0, axis=0)
C_coframe_nz = C_coframe[:, nz_cols_mask]
print(f" Coframe submatrix: {C_coframe_nz.shape[0]} × {C_coframe_nz.shape[1]}")
# Build torsion submatrix (subset of coframe rows)
n_torsion_rows = len(torsion_rows)
torsion_row_set = set(torsion_rows)
torsion_local = [ri for ri, row_idx in enumerate(coframe_rows)
if row_idx in torsion_row_set]
C_torsion = C_coframe[torsion_local, :]
nz_cols_t = np.any(C_torsion != 0, axis=0)
C_torsion_nz = C_torsion[:, nz_cols_t]
print(f" Torsion submatrix: {C_torsion_nz.shape[0]} × {C_torsion_nz.shape[1]}")
# ─── Compute ranks ────────────────────────────────────────────────
primes = [10007, 10009, 100003, 1000003, 10000019]
# Per-quadruple ranks for coframe sector
print("\nStep 5a: Per-quadruple ranks (coframe sector)...")
p_quick = primes[0]
per_quad_coframe = {}
per_quad_torsion = {}
sum_coframe = 0
sum_torsion = 0
for qi in coframe_quads:
quad, n_pats, start = quad_info[qi]
name = "·".join(BLOCK_DISPLAY[BLOCKS[b][0]] for b in quad)
# Find local row range in C_coframe
# coframe_rows is built by iterating coframe_quads in order
# so we need to find the offset
local_start = sum(quad_info[qj][1] for qj in coframe_quads[:coframe_quads.index(qi)])
local_end = local_start + n_pats
sub = C_coframe_nz[local_start:local_end, :]
# Remove zero columns
nz = np.any(sub != 0, axis=0)
sub_nz = sub[:, nz]
if sub_nz.shape[1] == 0:
r = 0
else:
r = rank_mod_p(sub_nz, p_quick)
per_quad_coframe[name] = {"n_patterns": n_pats, "rank": r}
sum_coframe += r
is_torsion = qi in torsion_quads
if is_torsion:
per_quad_torsion[name] = {"n_patterns": n_pats, "rank": r}
sum_torsion += r
print(f" {name:20s} patterns={n_pats:5d} rank={r:3d} {'[torsion]' if is_torsion else ''}")
print(f"\n Sum per-quadruple coframe ranks: {sum_coframe}")
print(f" Sum per-quadruple torsion ranks: {sum_torsion}")
# Global coframe rank
print(f"\nStep 5b: Global coframe rank ({n_coframe_rows} × {C_coframe_nz.shape[1]})...")
coframe_ranks = {}
for p in primes:
t1 = time.time()
r = rank_mod_p(C_coframe_nz, p)
elapsed = time.time() - t1
coframe_ranks[p] = r
print(f" Rank mod {p:>10d}: {r} ({elapsed:.1f}s)")
coframe_rank_vals = list(coframe_ranks.values())
assert len(set(coframe_rank_vals)) == 1, \
f"Coframe rank disagreement: {coframe_ranks}"
N_coframe = coframe_rank_vals[0]
print(f"\n *** COFRAME RANK: {N_coframe} ***")
# Global torsion rank
print(f"\nStep 5c: Global torsion rank ({C_torsion_nz.shape[0]} × {C_torsion_nz.shape[1]})...")
torsion_ranks = {}
for p in primes:
t1 = time.time()
r = rank_mod_p(C_torsion_nz, p)
elapsed = time.time() - t1
torsion_ranks[p] = r
print(f" Rank mod {p:>10d}: {r} ({elapsed:.1f}s)")
torsion_rank_vals = list(torsion_ranks.values())
assert len(set(torsion_rank_vals)) == 1, \
f"Torsion rank disagreement: {torsion_ranks}"
N_torsion = torsion_rank_vals[0]
print(f"\n *** TORSION RANK: {N_torsion} ***")
# ─── Cubic verification ───────────────────────────────────────────
print("\n" + "=" * 60)
print("CUBIC COFRAME SECTOR")
print("=" * 60)
# For cubic: enumerate triples, check which are coframe-only
surviving_cubic = []
for triple in combinations_with_replacement(range(len(BLOCKS)), 3):
parity = 1
total_idx = 0
for b in triple:
parity *= BLOCKS[b][2]
total_idx += BLOCKS[b][1]
if parity != +1:
continue
if total_idx % 2 != 0:
continue
surviving_cubic.append(triple)
coframe_cubic = []
for triple in surviving_cubic:
if set(triple) <= COFRAME_BLOCKS:
coframe_cubic.append(triple)
# Check: all coframe triples have odd total index count
# T=3, V=1, A=1, s=1. Any triple: sum of three odd numbers = odd → no even pairing
print(f"\n Total surviving cubic triples: {len(surviving_cubic)}")
print(f" Coframe cubic triples: {len(coframe_cubic)}")
if len(coframe_cubic) == 0:
print(" ALL coframe triples eliminated by index parity ✓")
print(" Reason: T(3), V(1), A(1), s(1) — any triple sums to odd")
N_torsion_cubic = 0
N_coframe_cubic = 0
else:
print(f" WARNING: {len(coframe_cubic)} coframe triples survive!")
for triple in coframe_cubic:
name = "·".join(BLOCK_DISPLAY[BLOCKS[b][0]] for b in triple)
total_idx = sum(BLOCKS[b][1] for b in triple)
print(f" {name} total_indices={total_idx}")
N_torsion_cubic = -1 # signal error
N_coframe_cubic = -1
# ─── Save results ─────────────────────────────────────────────────
elapsed_total = time.time() - t0
print(f"\n{'='*60}")
print(f"RESULTS SUMMARY")
print(f"{'='*60}")
print(f" Quartic coframe rank: {N_coframe}")
print(f" Quartic torsion rank: {N_torsion}")
print(f" Cubic coframe rank: {N_coframe_cubic}")
print(f" Cubic torsion rank: {N_torsion_cubic}")
print(f" Full quartic rank (paper3): 156")
print(f" Elapsed: {elapsed_total:.1f}s")
results = {
"N_coframe_quartic": N_coframe,
"N_torsion_quartic": N_torsion,
"N_coframe_cubic": N_coframe_cubic,
"N_torsion_cubic": N_torsion_cubic,
"coframe_quadruples": coframe_quad_names,
"torsion_quadruples": torsion_quad_names,
"n_coframe_patterns": n_coframe_rows,
"n_torsion_patterns": n_torsion_rows,
"coframe_submatrix_shape": list(C_coframe_nz.shape),
"torsion_submatrix_shape": list(C_torsion_nz.shape),
"per_quadruple_coframe": per_quad_coframe,
"per_quadruple_torsion": per_quad_torsion,
"sum_per_quad_coframe": sum_coframe,
"sum_per_quad_torsion": sum_torsion,
"primes": primes,
"coframe_ranks_by_prime": {str(p): r for p, r in coframe_ranks.items()},
"torsion_ranks_by_prime": {str(p): r for p, r in torsion_ranks.items()},
"cubic_coframe_triples_surviving": len(coframe_cubic),
"cubic_index_parity_argument": "T(3)+V(1)+A(1)+s(1): any triple of odd numbers sums to odd",
"elapsed_seconds": elapsed_total,
}
out_dir = os.path.dirname(__file__) or "."
out_file = os.path.join(out_dir, "01_coframe_sector.json")
with open(out_file, "w") as f:
json.dump(results, f, indent=2, ensure_ascii=False)
print(f"\nSaved: {out_file}")
if __name__ == "__main__":
main()
Raw
02_covariant_check.json
{
"quartic_results": [
{
"linearized": 21,
"covariant": 21
},
{
"linearized": 21,
"covariant": 21
},
{
"linearized": 21,
"covariant": 21
},
{
"linearized": 21,
"covariant": 21
},
{
"linearized": 21,
"covariant": 21
}
],
"all_match": true,
"identity_control": {
"lin": 21,
"cov": 21,
"max_diff": 0.0
},
"n_sampled": 3194,
"n_total": 11589,
"n_coframes": 5,
"n_de_per_coframe": 100,
"convention": "all_lower_indices",
"elapsed_seconds": 580.2590227127075
}
Raw
02_covariant_check.py
#!/usr/bin/env python3
"""02_covariant_check.py — Verify η↔g covariant equivalence numerically.
All blocks use LOWER indices consistently:
T_{ρμν}, V_μ, A_μ (where A_μ = ε_{μνρσ} T^{νρσ})
Contractions use the UPPER metric (η^{μν} or g^{μν}).
Protocol: fix coframe e, vary de, check rank of (config × pattern) matrix.
"""
import numpy as np
from itertools import combinations_with_replacement
from collections import defaultdict
import json, os, time, functools
print = functools.partial(print, flush=True)
D = 4
eta = np.array([[-1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]], dtype=float)
def levi_civita_down():
eps = np.zeros((D,D,D,D), dtype=float)
for i in range(D):
for j in range(D):
for k in range(D):
for l in range(D):
if len(set([i,j,k,l])) < D: continue
s = 1
p = [i,j,k,l]
for a in range(D):
for b in range(a+1,D):
if p[a]>p[b]: s *= -1
eps[i,j,k,l] = s
return eps
EPS = levi_civita_down()
def random_coframe(rng, scale=0.5):
e = np.eye(D) + scale * rng.randn(D, D)
while abs(np.linalg.det(e)) < 0.1:
e = np.eye(D) + scale * rng.randn(D, D)
return e
def generate_pairings(indices):
if len(indices) == 0: return [()]
if len(indices) == 2: return [((indices[0], indices[1]),)]
first = indices[0]; rest = indices[1:]
result = []
for i, partner in enumerate(rest):
remaining = rest[:i] + rest[i+1:]
for sub in generate_pairings(remaining):
result.append(((first, partner),) + sub)
return result
# ─── Block computation (all lower indices) ────────────────────────────
def compute_blocks_lower(de, e=None):
"""Compute T_{ρμν}, V_μ, A_μ with ALL LOWER indices.
If e is None: linearized (η for everything).
If e given: covariant (g(e) for everything).
"""
# T^a_{μν} = de[μ,a,ν] - de[ν,a,μ] (frame torsion, metric-independent)
Ta = np.zeros((D,D,D))
for a in range(D):
for mu in range(D):
for nu in range(D):
Ta[a,mu,nu] = de[mu,a,nu] - de[nu,a,mu]
if e is None:
# Linearized: T_{ρμν} = η_{ρa} T^a_{μν}
T = np.einsum('ra,amn->rmn', eta, Ta)
met_up = eta
else:
# Covariant: T_{ρμν} = η_{ab} e^b_ρ T^a_{μν}
# (= g_{ρσ} e_a^σ T^a_{μν}, but η_{ab} e^b_ρ is simpler)
L = np.einsum('ab,br->ra', eta, e) # L[r,a] = η_{ab} e^b_ρ
T = np.einsum('ra,amn->rmn', L, Ta)
g_down = np.einsum('ab,am,bn->mn', eta, e, e)
met_up = np.linalg.inv(g_down)
# V_μ = met^{νρ} T_{ρνμ} (lower index μ)
V = np.einsum('nr,rnm->m', met_up, T)
# A_μ = ε_{μνρσ} T^{νρσ} where T^{νρσ} = met^{na} met^{rb} met^{sc} T_{abc}
T_allup = np.einsum('na,rb,sc,abc->nrs', met_up, met_up, met_up, T)
A = np.einsum('mnrs,nrs->m', EPS, T_allup)
return T, V, A, met_up
# ─── Efficient einsum evaluation ──────────────────────────────────────
LETTERS = 'abcdefghijklmnopqrstuvwxyz'
def eval_patterns(block_tensors, block_defs, patterns, metric):
"""Evaluate contraction patterns via np.einsum."""
values = np.zeros(len(patterns))
for pi, (quad, pairing) in enumerate(patterns):
n_idx = [block_defs[b][1] for b in quad]
total = sum(n_idx)
if total == 0:
v = 1.0
for b in quad: v *= float(block_tensors[b])
values[pi] = v; continue
offsets, off = [], 0
for ni in n_idx: offsets.append(off); off += ni
idx = [LETTERS[i] for i in range(total)]
ops, subs = [], []
for bi in range(len(quad)):
s, e2 = offsets[bi], offsets[bi]+n_idx[bi]
subs.append(''.join(idx[s:e2]))
ops.append(block_tensors[quad[bi]])
for (a,b) in pairing:
subs.append(idx[a]+idx[b])
ops.append(metric)
try:
values[pi] = np.einsum(','.join(subs)+'->', *ops, optimize=True)
except:
values[pi] = 0.0
return values
def enumerate_torsion_patterns():
blocks = [("T",3,+1), ("V",1,+1), ("A",1,-1)]
surviving = []
for quad in combinations_with_replacement(range(len(blocks)), 4):
par, tot = 1, 0
for b in quad: par *= blocks[b][2]; tot += blocks[b][1]
if par != +1 or tot % 2 != 0: continue
surviving.append(quad)
all_patterns, quad_info = [], []
for quad in surviving:
n_idx = [blocks[b][1] for b in quad]
total = sum(n_idx)
pairings = generate_pairings(list(range(total))) if total > 0 else [()]
start = len(all_patterns)
for p in pairings: all_patterns.append((quad, p))
quad_info.append((quad, len(pairings), start))
return blocks, surviving, all_patterns, quad_info
def main():
t0 = time.time()
rng = np.random.RandomState(42)
print("="*60)
print("COVARIANT EQUIVALENCE CHECK — ALL-LOWER CONVENTION")
print("="*60)
blocks, surviving, all_patterns, quad_info = enumerate_torsion_patterns()
n_total = len(all_patterns)
print(f"\nTotal torsion patterns: {n_total}")
# Sample: random 2000 from T⁴, all from rest
sampled = []
for qi, (quad, n_pats, start) in enumerate(quad_info):
name = "·".join(["T","V","A"][b] for b in quad)
if n_pats > 2000:
indices = rng.choice(n_pats, 2000, replace=False)
for i in indices:
sampled.append(all_patterns[start + i])
print(f" {name:15s} {n_pats:5d} → sampled 2000")
else:
for i in range(n_pats):
sampled.append(all_patterns[start + i])
print(f" {name:15s} {n_pats:5d} → all")
print(f" Total sampled: {len(sampled)}")
# ─── Quadratic check ──────────────────────────────────────────────
print("\n--- QUADRATIC (Hayashi-Shirafuji) ---")
N_E, M_DE = 5, 30
for ei in range(N_E):
e = random_coframe(rng, scale=0.5)
lin_m = np.zeros((M_DE, 3))
cov_m = np.zeros((M_DE, 3))
for di in range(M_DE):
de = rng.randn(D,D,D)
Tl,Vl,Al,ml = compute_blocks_lower(de)
Tc,Vc,Ac,mc = compute_blocks_lower(de, e)
lin_m[di,0] = np.einsum('abc,def,ad,be,cf->',Tl,Tl,ml,ml,ml)
lin_m[di,1] = np.einsum('abc,def,ae,bd,cf->',Tl,Tl,ml,ml,ml)
lin_m[di,2] = np.einsum('a,b,ab->',Vl,Vl,ml)
cov_m[di,0] = np.einsum('abc,def,ad,be,cf->',Tc,Tc,mc,mc,mc)
cov_m[di,1] = np.einsum('abc,def,ae,bd,cf->',Tc,Tc,mc,mc,mc)
cov_m[di,2] = np.einsum('a,b,ab->',Vc,Vc,mc)
rl = np.linalg.matrix_rank(lin_m, tol=1e-8)
rc = np.linalg.matrix_rank(cov_m, tol=1e-8)
print(f" e[{ei}]: lin={rl} cov={rc} det(e)={np.linalg.det(e):.3f} {'✓' if rl==rc else '✗'}")
# ─── Quartic check ────────────────────────────────────────────────
print(f"\n--- QUARTIC (T,V,A) — {len(sampled)} patterns ---")
N_E_Q, M_DE_Q = 5, 100
results_q = []
for ei in range(N_E_Q):
e = random_coframe(rng, scale=0.3)
g_down = np.einsum('ab,am,bn->mn', eta, e, e)
g_up = np.linalg.inv(g_down)
lin_m = np.zeros((M_DE_Q, len(sampled)))
cov_m = np.zeros((M_DE_Q, len(sampled)))
for di in range(M_DE_Q):
de = rng.randn(D,D,D) * 0.5
Tl,Vl,Al,_ = compute_blocks_lower(de)
Tc,Vc,Ac,_ = compute_blocks_lower(de, e)
lin_m[di] = eval_patterns({0:Tl,1:Vl,2:Al}, blocks, sampled, eta)
cov_m[di] = eval_patterns({0:Tc,1:Vc,2:Ac}, blocks, sampled, g_up)
rl = np.linalg.matrix_rank(lin_m, tol=1e-6)
rc = np.linalg.matrix_rank(cov_m, tol=1e-6)
el = time.time()-t0
print(f" e[{ei}]: lin={rl} cov={rc} det(e)={np.linalg.det(e):.3f} ({el:.0f}s) {'✓' if rl==rc else '✗'}")
results_q.append((rl, rc))
# Identity control
print("\n Control: e = identity")
e_id = np.eye(D)
lin_id = np.zeros((M_DE_Q, len(sampled)))
cov_id = np.zeros((M_DE_Q, len(sampled)))
rng2 = np.random.RandomState(999)
for di in range(M_DE_Q):
de = rng2.randn(D,D,D)*0.5
Tl,Vl,Al,_ = compute_blocks_lower(de)
Tc,Vc,Ac,_ = compute_blocks_lower(de, e_id)
lin_id[di] = eval_patterns({0:Tl,1:Vl,2:Al}, blocks, sampled, eta)
cov_id[di] = eval_patterns({0:Tc,1:Vc,2:Ac}, blocks, sampled, eta)
maxd = np.max(np.abs(lin_id - cov_id))
rl_id = np.linalg.matrix_rank(lin_id, tol=1e-6)
rc_id = np.linalg.matrix_rank(cov_id, tol=1e-6)
print(f" e=I: lin={rl_id} cov={rc_id} max|diff|={maxd:.2e}")
# ─── Summary ──────────────────────────────────────────────────────
all_match = all(rl==rc for rl,rc in results_q)
elapsed = time.time()-t0
print(f"\n{'='*60}")
print(f"RESULTS (all-lower-index convention)")
print(f"{'='*60}")
for ei,(rl,rc) in enumerate(results_q):
print(f" e[{ei}]: lin={rl} cov={rc} {'✓' if rl==rc else '✗'}")
print(f" Identity: lin={rl_id} cov={rc_id} max_diff={maxd:.2e}")
print(f" Overall: {'VERIFIED ✓' if all_match else 'FAILED ✗'}")
print(f" Elapsed: {elapsed:.1f}s")
out = {
"quartic_results": [{"linearized":int(rl),"covariant":int(rc)} for rl,rc in results_q],
"all_match": bool(all_match),
"identity_control": {"lin":int(rl_id),"cov":int(rc_id),"max_diff":float(maxd)},
"n_sampled": len(sampled),
"n_total": n_total,
"n_coframes": N_E_Q,
"n_de_per_coframe": M_DE_Q,
"convention": "all_lower_indices",
"elapsed_seconds": elapsed,
}
out_file = os.path.join(os.path.dirname(__file__) or ".", "02_covariant_check.json")
with open(out_file, "w") as f:
json.dump(out, f, indent=2)
print(f"\nSaved: {out_file}")
if __name__ == "__main__":
main()
Raw
06_complete_lower_orders.json
{
"quadratic": {
"order": 2,
"parity_even": 9,
"parity_odd": 6,
"total": 15,
"n_patterns": 63,
"n_eta": 42,
"n_eps": 21,
"n_monomials": 404,
"eta_only_even": 8,
"eta_only_odd": 2
},
"cubic": {
"order": 3,
"parity_even": 33,
"parity_odd": 25,
"total": 58,
"n_patterns": 1041,
"n_eta": 438,
"n_eps": 603,
"n_monomials": 5528,
"eta_only_even": 26,
"eta_only_odd": 9
}
}
Raw
06_complete_lower_orders.py
#!/usr/bin/env python3
"""06_complete_lower_orders.py — Complete quadratic + cubic classification.
All 7 blocks, both parities, η-only AND ε-contractions.
Gives the definitive counts at orders 2 and 3.
"""
import numpy as np
from itertools import combinations_with_replacement, combinations
from collections import defaultdict
import json, os, sys, time, functools
print = functools.partial(print, flush=True)
D = 4
eta = np.array([[-1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]], dtype=int)
eps_down = np.zeros((D,D,D,D), dtype=int)
for i in range(D):
for j in range(D):
for k in range(D):
for l in range(D):
if len(set([i,j,k,l])) < D: continue
s = 1
for a in range(D):
for b in range(a+1,D):
if [i,j,k,l][a] > [i,j,k,l][b]: s *= -1
eps_down[i,j,k,l] = s
eps_up = np.einsum('ma,nb,rc,sd,abcd->mnrs', eta, eta, eta, eta, eps_down)
BLOCKS = [
("Shat", 2, +1),
("theta", 0, +1),
("W", 2, +1),
("T", 3, +1),
("V", 1, +1),
("A", 1, -1),
("s", 1, +1),
]
BLOCK_DISPLAY = {"Shat": "Ŝ", "theta": "θ", "W": "W", "T": "T",
"V": "V", "A": "𝒜", "s": "s"}
COFRAME = {3, 4, 5, 6}
TORSION = {3, 4, 5}
DISPLACEMENT = {0, 1, 2}
# ── Block component functions (same as 05) ───────────────────────────
def shat_component(mu, nu):
terms = {}
terms[('du', mu, nu)] = terms.get(('du', mu, nu), 0) + 2
terms[('du', nu, mu)] = terms.get(('du', nu, mu), 0) + 2
for rho in range(D):
key = ('du', rho, rho)
terms[key] = terms.get(key, 0) - eta[mu, nu] * eta[rho, rho]
return {k: v for k, v in terms.items() if v != 0}
def theta_component():
return {('du', mu, mu): int(eta[mu, mu]) for mu in range(D) if eta[mu, mu] != 0}
def W_component(mu, nu):
terms = {}
terms[('du', mu, nu)] = terms.get(('du', mu, nu), 0) + 1
terms[('du', nu, mu)] = terms.get(('du', nu, mu), 0) - 1
return {k: v for k, v in terms.items() if v != 0}
def T_component(rho, mu, nu):
terms = {}
for a in range(D):
if eta[a, rho] == 0: continue
terms[('de', mu, a, nu)] = terms.get(('de', mu, a, nu), 0) + eta[a, rho]
terms[('de', nu, a, mu)] = terms.get(('de', nu, a, mu), 0) - eta[a, rho]
return {k: v for k, v in terms.items() if v != 0}
def V_component(mu):
terms = {}
for nu in range(D):
for rho in range(D):
if eta[nu, rho] == 0: continue
for k, v in T_component(rho, nu, mu).items():
terms[k] = terms.get(k, 0) + eta[nu, rho] * v
return {k: v for k, v in terms.items() if v != 0}
def A_component(mu):
terms = {}
for nu in range(D):
for rho in range(D):
for sig in range(D):
e = eps_up[mu, nu, rho, sig]
if e == 0: continue
for k, v in T_component(nu, rho, sig).items():
terms[k] = terms.get(k, 0) + e * v
return {k: v for k, v in terms.items() if v != 0} # returns 6A; rank-invariant
def s_component(mu):
return {('de', mu, a, a): 1 for a in range(D)}
def get_block_component(block_idx, indices):
name = BLOCKS[block_idx][0]
if name == "Shat": return shat_component(indices[0], indices[1])
elif name == "theta": return theta_component()
elif name == "W": return W_component(indices[0], indices[1])
elif name == "T": return T_component(indices[0], indices[1], indices[2])
elif name == "V": return V_component(indices[0])
elif name == "A": return A_component(indices[0])
elif name == "s": return s_component(indices[0])
# ── Pairing ──────────────────────────────────────────────────────────
def generate_pairings(indices):
if len(indices) == 0: return [()]
if len(indices) == 2: return [((indices[0], indices[1]),)]
first = indices[0]; rest = indices[1:]
result = []
for i, partner in enumerate(rest):
remaining = rest[:i] + rest[i+1:]
for sub in generate_pairings(remaining):
result.append(((first, partner),) + sub)
return result
# ── Generic expansion (any number of blocks) ─────────────────────────
def expand_eta(block_indices, pattern):
"""Expand η-only contraction for arbitrary block count."""
n_blocks = len(block_indices)
n = [BLOCKS[b][1] for b in block_indices]
total = sum(n)
offsets = [sum(n[:i]) for i in range(n_blocks)]
n_pairs = len(pattern)
poly = defaultdict(int)
def recurse(pair_idx, assignments):
if pair_idx == n_pairs:
eta_prod = 1
for (a, b) in pattern:
eta_prod *= eta[assignments[a], assignments[b]]
if eta_prod == 0: return
components = []
for bi in range(n_blocks):
s2 = offsets[bi]; e2 = s2 + n[bi]
idx_vals = tuple(assignments[g] for g in range(s2, e2))
components.append(get_block_component(block_indices[bi],
idx_vals if n[bi] > 0 else ()))
# Multiply all components
result_terms = components[0]
for ci in range(1, n_blocks):
new_terms = defaultdict(int)
for k1, v1 in result_terms.items():
for k2, v2 in components[ci].items():
mono = tuple(sorted([k1, k2])) if isinstance(k1, tuple) and k1[0] in ('du','de') else tuple(sorted(list(k1) + [k2])) if not isinstance(k1[0], str) else None
# Build degree-n monomial
pass
result_terms = new_terms
# Actually, let's just build the product directly for each leaf
pass
a, b = pattern[pair_idx]
for va in range(D):
assignments[a] = va
for vb in range(D):
if eta[va, vb] == 0: continue
assignments[b] = vb
recurse(pair_idx + 1, assignments)
# Simpler approach: at the leaf, multiply all block components
def recurse2(pair_idx, assignments):
if pair_idx == n_pairs:
eta_prod = 1
for (a, b) in pattern:
eta_prod *= eta[assignments[a], assignments[b]]
if eta_prod == 0: return
# Get all block components
comp_list = []
for bi in range(n_blocks):
s2 = offsets[bi]; e2 = s2 + n[bi]
idx_vals = tuple(assignments[g] for g in range(s2, e2))
comp_list.append(get_block_component(block_indices[bi],
idx_vals if n[bi] > 0 else ()))
# Multiply: iterate over all combinations of terms
def multiply_dicts(dicts, idx=0, current_keys=None, current_val=1):
if current_keys is None: current_keys = []
if idx == len(dicts):
mono = tuple(sorted(current_keys))
poly[mono] += eta_prod * current_val
return
for k, v in dicts[idx].items():
multiply_dicts(dicts, idx+1, current_keys + [k], current_val * v)
multiply_dicts(comp_list)
return
a, b = pattern[pair_idx]
for va in range(D):
assignments[a] = va
for vb in range(D):
if eta[va, vb] == 0: continue
assignments[b] = vb
recurse2(pair_idx + 1, assignments)
assignments = [0] * total
recurse2(0, assignments)
return {k: v for k, v in poly.items() if v != 0}
def expand_epsilon(block_indices, eps_indices, eta_pattern):
"""Expand one ε+η contraction for arbitrary block count."""
n_blocks = len(block_indices)
n = [BLOCKS[b][1] for b in block_indices]
total = sum(n)
offsets = [sum(n[:i]) for i in range(n_blocks)]
n_eta = len(eta_pattern)
poly = defaultdict(int)
def recurse_eta(pair_idx, assignments):
if pair_idx == n_eta:
eta_prod = 1
for (a, b) in eta_pattern:
eta_prod *= eta[assignments[a], assignments[b]]
if eta_prod == 0: return
eps_val = eps_down[assignments[eps_indices[0]],
assignments[eps_indices[1]],
assignments[eps_indices[2]],
assignments[eps_indices[3]]]
if eps_val == 0: return
coeff = eps_val * eta_prod
comp_list = []
for bi in range(n_blocks):
s2 = offsets[bi]; e2 = s2 + n[bi]
idx_vals = tuple(assignments[g] for g in range(s2, e2))
comp_list.append(get_block_component(block_indices[bi],
idx_vals if n[bi] > 0 else ()))
def multiply_dicts(dicts, idx=0, current_keys=None, current_val=1):
if current_keys is None: current_keys = []
if idx == len(dicts):
mono = tuple(sorted(current_keys))
poly[mono] += coeff * current_val
return
for k, v in dicts[idx].items():
multiply_dicts(dicts, idx+1, current_keys + [k], current_val * v)
multiply_dicts(comp_list)
return
a, b = eta_pattern[pair_idx]
for va in range(D):
assignments[a] = va
for vb in range(D):
if eta[va, vb] == 0: continue
assignments[b] = vb
recurse_eta(pair_idx + 1, assignments)
assignments = [0] * total
# Iterate ε index assignments
from itertools import product as iprod
for vals in iprod(range(D), repeat=4):
e = eps_down[vals[0], vals[1], vals[2], vals[3]]
if e == 0: continue
for idx, val in zip(eps_indices, vals):
assignments[idx] = val
recurse_eta(0, assignments)
return {k: v for k, v in poly.items() if v != 0}
def rank_mod_p(mat, p):
m, n2 = mat.shape
A = mat.astype(np.int64) % p
A[A < 0] += p
rank = 0
for col in range(n2):
if rank >= m: break
nz = np.nonzero(A[rank:, col])[0]
if len(nz) == 0: continue
found = nz[0] + rank
if found != rank:
A[[rank, found]] = A[[found, rank]]
inv_piv = pow(int(A[rank, col]), p - 2, p)
col_vals = A[:, col].copy()
col_vals[rank] = 0
nz_rows = np.nonzero(col_vals)[0]
if len(nz_rows) > 0:
factors = (col_vals[nz_rows] * inv_piv) % p
A[nz_rows] = (A[nz_rows] - factors[:, np.newaxis] * A[rank]) % p
rank += 1
return rank
def classify_order(order, label):
"""Full classification at a given polynomial order."""
print(f"\n{'='*60}")
print(f"ORDER {order} ({label})")
print(f"{'='*60}")
# Enumerate all n-tuples with even total index count
all_patterns = [] # (tuple, contraction_type, contraction_data, parity)
for combo in combinations_with_replacement(range(len(BLOCKS)), order):
n_idx = [BLOCKS[b][1] for b in combo]
total = sum(n_idx)
block_parity = 1
for b in combo:
block_parity *= BLOCKS[b][2]
if total % 2 != 0:
continue
all_indices = list(range(total))
# η-only patterns
if total > 0:
pairings = generate_pairings(all_indices)
else:
pairings = [()]
for p in pairings:
all_patterns.append((combo, 'eta', p, block_parity))
# ε patterns (k=1)
if total >= 4 and (total - 4) % 2 == 0:
for eps_idx in combinations(all_indices, 4):
remaining = [x for x in all_indices if x not in eps_idx]
if len(remaining) == 0:
eta_pats = [()]
else:
eta_pats = generate_pairings(remaining)
for ep in eta_pats:
all_patterns.append((combo, 'eps', (eps_idx, ep), -block_parity))
n_total = len(all_patterns)
n_eta = sum(1 for _,t,_,_ in all_patterns if t == 'eta')
n_eps = sum(1 for _,t,_,_ in all_patterns if t == 'eps')
print(f" η-only patterns: {n_eta}")
print(f" ε patterns: {n_eps}")
print(f" Total: {n_total}")
# Expand symbolically
print(f"\n Expanding {n_total} patterns...")
all_polys = []
all_monomials = set()
t0 = time.time()
for idx, (combo, ctype, cdata, par) in enumerate(all_patterns):
if ctype == 'eta':
poly = expand_eta(combo, cdata)
else:
eps_idx, eta_pat = cdata
poly = expand_epsilon(combo, eps_idx, eta_pat)
all_polys.append(poly)
all_monomials.update(poly.keys())
if (idx+1) % 100 == 0:
print(f" {idx+1}/{n_total} ({time.time()-t0:.0f}s)")
n_mon = len(all_monomials)
n_nz = sum(1 for p in all_polys if p)
print(f" {n_nz}/{n_total} nonzero, {n_mon} monomials ({time.time()-t0:.0f}s)")
# Build coefficient matrix
monomial_list = sorted(all_monomials)
mon_idx = {m: i for i, m in enumerate(monomial_list)}
C = np.zeros((n_total, n_mon), dtype=np.int64)
for row, poly in enumerate(all_polys):
for mono, coeff in poly.items():
C[row, mon_idx[mono]] = coeff
nz_cols = np.any(C != 0, axis=0)
C_nz = C[:, nz_cols]
# Compute ranks by parity
p_mod = 10007
even_rows = [i for i, (_,_,_,par) in enumerate(all_patterns) if par == +1]
odd_rows = [i for i, (_,_,_,par) in enumerate(all_patterns) if par == -1]
r_even = rank_mod_p(C_nz[even_rows, :].copy(), p_mod) if even_rows else 0
r_odd = rank_mod_p(C_nz[odd_rows, :].copy(), p_mod) if odd_rows else 0
r_total = rank_mod_p(C_nz.copy(), p_mod)
print(f"\n RESULTS (order {order}):")
print(f" Parity-even: {r_even}")
print(f" Parity-odd: {r_odd}")
print(f" Total: {r_total}")
print(f" Check: {r_even}+{r_odd}={r_even+r_odd} vs {r_total}")
# Sector breakdown
for sector_name, sector_set in [("torsion", TORSION), ("coframe", COFRAME), ("full", set(range(7)))]:
s_even = [i for i in even_rows if set(all_patterns[i][0]) <= sector_set]
s_odd = [i for i in odd_rows if set(all_patterns[i][0]) <= sector_set]
s_all = s_even + s_odd
re = rank_mod_p(C_nz[s_even, :].copy(), p_mod) if s_even else 0
ro = rank_mod_p(C_nz[s_odd, :].copy(), p_mod) if s_odd else 0
print(f" {sector_name:10s}: even={re}, odd={ro}, total={re+ro}")
# η-only sub-count
eta_even = [i for i in even_rows if all_patterns[i][1] == 'eta']
eta_odd = [i for i in odd_rows if all_patterns[i][1] == 'eta']
re_eta = rank_mod_p(C_nz[eta_even, :].copy(), p_mod) if eta_even else 0
ro_eta = rank_mod_p(C_nz[eta_odd, :].copy(), p_mod) if eta_odd else 0
print(f" η-only: even={re_eta}, odd={ro_eta}")
return {
"order": order,
"parity_even": r_even,
"parity_odd": r_odd,
"total": r_total,
"n_patterns": n_total,
"n_eta": n_eta,
"n_eps": n_eps,
"n_monomials": n_mon,
"eta_only_even": re_eta,
"eta_only_odd": ro_eta,
}
def main():
t0 = time.time()
results = {}
results["quadratic"] = classify_order(2, "QUADRATIC")
results["cubic"] = classify_order(3, "CUBIC")
elapsed = time.time() - t0
print(f"\n{'='*60}")
print(f"SUMMARY")
print(f"{'='*60}")
for label in ["quadratic", "cubic"]:
r = results[label]
print(f" Order {r['order']}: even={r['parity_even']}, odd={r['parity_odd']}, total={r['total']}")
print(f" η-only: even={r['eta_only_even']}, odd={r['eta_only_odd']}")
print(f" Elapsed: {elapsed:.0f}s")
out_file = os.path.join(os.path.dirname(__file__) or ".", "06_complete_lower_orders.json")
with open(out_file, "w") as f:
json.dump(results, f, indent=2)
print(f"\nSaved: {out_file}")
if __name__ == "__main__":
main()
Raw
08_quartic_exact.py
#!/usr/bin/env python3
"""08_quartic_exact.py — EXACT quartic classification with ALL ε patterns.
Complete computation: all η-only + ALL ε-contraction patterns for all
quadruples including T⁴ (51975 ε patterns). No sampling.
"""
import numpy as np
from itertools import combinations_with_replacement, combinations
from collections import defaultdict
import json, os, sys, time, functools, random
print = functools.partial(print, flush=True)
D = 4
eta = np.array([[-1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]], dtype=int)
eps_down = np.zeros((D,D,D,D), dtype=int)
for i in range(D):
for j in range(D):
for k in range(D):
for l in range(D):
if len(set([i,j,k,l])) < D: continue
s = 1
for a in range(D):
for b in range(a+1,D):
if [i,j,k,l][a] > [i,j,k,l][b]: s *= -1
eps_down[i,j,k,l] = s
eps_up = np.einsum('ma,nb,rc,sd,abcd->mnrs', eta, eta, eta, eta, eps_down)
BLOCKS = [
("Shat", 2, +1), ("theta", 0, +1), ("W", 2, +1),
("T", 3, +1), ("V", 1, +1), ("A", 1, -1), ("s", 1, +1),
]
DISP = {"Shat": "Ŝ", "theta": "θ", "W": "W", "T": "T",
"V": "V", "A": "𝒜", "s": "s"}
# ── Block component functions ────────────────────────────────────────
def shat_component(mu, nu):
t = {}
t[('du',mu,nu)] = t.get(('du',mu,nu),0)+2
t[('du',nu,mu)] = t.get(('du',nu,mu),0)+2
for r in range(D):
k = ('du',r,r); t[k] = t.get(k,0) - eta[mu,nu]*eta[r,r]
return {k:v for k,v in t.items() if v!=0}
def theta_component():
return {('du',m,m): int(eta[m,m]) for m in range(D) if eta[m,m]!=0}
def W_component(mu, nu):
t = {}
t[('du',mu,nu)] = t.get(('du',mu,nu),0)+1
t[('du',nu,mu)] = t.get(('du',nu,mu),0)-1
return {k:v for k,v in t.items() if v!=0}
def T_component(rho, mu, nu):
t = {}
for a in range(D):
if eta[a,rho]==0: continue
t[('de',mu,a,nu)] = t.get(('de',mu,a,nu),0) + eta[a,rho]
t[('de',nu,a,mu)] = t.get(('de',nu,a,mu),0) - eta[a,rho]
return {k:v for k,v in t.items() if v!=0}
def V_component(mu):
t = {}
for nu in range(D):
for rho in range(D):
if eta[nu,rho]==0: continue
for k,v in T_component(rho,nu,mu).items():
t[k] = t.get(k,0) + eta[nu,rho]*v
return {k:v for k,v in t.items() if v!=0}
def A_component(mu):
t = {}
for nu in range(D):
for rho in range(D):
for sig in range(D):
e = eps_up[mu,nu,rho,sig]
if e==0: continue
for k,v in T_component(nu,rho,sig).items():
t[k] = t.get(k,0) + e*v
return {k:v for k,v in t.items() if v!=0} # raw = 6A
def s_component(mu):
return {('de',mu,a,a): 1 for a in range(D)}
def get_comp(bi, idx):
n = BLOCKS[bi][0]
if n=="Shat": return shat_component(idx[0],idx[1])
if n=="theta": return theta_component()
if n=="W": return W_component(idx[0],idx[1])
if n=="T": return T_component(idx[0],idx[1],idx[2])
if n=="V": return V_component(idx[0])
if n=="A": return A_component(idx[0])
if n=="s": return s_component(idx[0])
def generate_pairings(indices):
if len(indices)==0: return [()]
if len(indices)==2: return [((indices[0],indices[1]),)]
first = indices[0]; rest = indices[1:]
result = []
for i,p in enumerate(rest):
rem = rest[:i]+rest[i+1:]
for sub in generate_pairings(rem):
result.append(((first,p),)+sub)
return result
# ── Generic expansion ────────────────────────────────────────────────
def expand_eta_quartic(blocks, pattern):
n = [BLOCKS[b][1] for b in blocks]
total = sum(n); off = [sum(n[:i]) for i in range(4)]
n_pairs = len(pattern); poly = defaultdict(int)
def rec(pi, asn):
if pi==n_pairs:
ep = 1
for a,b in pattern: ep *= eta[asn[a],asn[b]]
if ep==0: return
cs = []
for bi in range(4):
iv = tuple(asn[g] for g in range(off[bi],off[bi]+n[bi]))
cs.append(get_comp(blocks[bi], iv if n[bi]>0 else ()))
def mul(ds,i=0,ks=None,v=1):
if ks is None: ks=[]
if i==len(ds):
poly[tuple(sorted(ks))] += ep*v; return
for k2,v2 in ds[i].items(): mul(ds,i+1,ks+[k2],v*v2)
mul(cs); return
a,b = pattern[pi]
for va in range(D):
asn[a]=va
for vb in range(D):
if eta[va,vb]==0: continue
asn[b]=vb; rec(pi+1,asn)
rec(0,[0]*total)
return {k:v for k,v in poly.items() if v!=0}
def expand_eps_quartic(blocks, eps_idx, eta_pat):
from itertools import product as iprod
n = [BLOCKS[b][1] for b in blocks]
total = sum(n); off = [sum(n[:i]) for i in range(4)]
n_eta = len(eta_pat); poly = defaultdict(int)
def rec_eta(pi, asn):
if pi==n_eta:
ep = 1
for a,b in eta_pat: ep *= eta[asn[a],asn[b]]
if ep==0: return
ev = eps_down[asn[eps_idx[0]],asn[eps_idx[1]],
asn[eps_idx[2]],asn[eps_idx[3]]]
if ev==0: return
coeff = ev*ep
cs = []
for bi in range(4):
iv = tuple(asn[g] for g in range(off[bi],off[bi]+n[bi]))
cs.append(get_comp(blocks[bi], iv if n[bi]>0 else ()))
def mul(ds,i=0,ks=None,v=1):
if ks is None: ks=[]
if i==len(ds):
poly[tuple(sorted(ks))] += coeff*v; return
for k2,v2 in ds[i].items(): mul(ds,i+1,ks+[k2],v*v2)
mul(cs); return
a,b = eta_pat[pi]
for va in range(D):
asn[a]=va
for vb in range(D):
if eta[va,vb]==0: continue
asn[b]=vb; rec_eta(pi+1,asn)
asn = [0]*total
for vals in iprod(range(D), repeat=4):
e = eps_down[vals[0],vals[1],vals[2],vals[3]]
if e==0: continue
for idx,val in zip(eps_idx,vals): asn[idx]=val
rec_eta(0, asn)
return {k:v for k,v in poly.items() if v!=0}
def rank_mod_p(mat, p):
m,nc = mat.shape
A = mat.astype(np.int64)%p; A[A<0]+=p
rank = 0
for col in range(nc):
if rank>=m: break
nz = np.nonzero(A[rank:,col])[0]
if len(nz)==0: continue
f = nz[0]+rank
if f!=rank: A[[rank,f]]=A[[f,rank]]
inv = pow(int(A[rank,col]),p-2,p)
cv = A[:,col].copy(); cv[rank]=0
nzr = np.nonzero(cv)[0]
if len(nzr)>0:
fac = (cv[nzr]*inv)%p
A[nzr] = (A[nzr]-fac[:,np.newaxis]*A[rank])%p
rank += 1
if rank%20==0: print(f" rank {rank}, col {col}/{nc}")
return rank
# ── T⁴ ε sampling ───────────────────────────────────────────────────
def classify_t4_eps_indices(eps_idx):
"""Classify 4 chosen indices by which T block they belong to."""
groups = [0]*4 # how many from each T
for i in eps_idx:
groups[i // 3] += 1
return tuple(sorted(groups, reverse=True))
def sample_t4_eps_patterns(n_sample=200, seed=42):
"""Sample T⁴ ε patterns from non-trivial distribution types."""
rng = random.Random(seed)
T_IDX = 3 # T block index
quad = (T_IDX, T_IDX, T_IDX, T_IDX)
all_indices = list(range(12))
# Categorize all C(12,4) = 495 ε-index choices
by_type = defaultdict(list)
for eps_idx in combinations(all_indices, 4):
dist = classify_t4_eps_indices(eps_idx)
by_type[dist].append(eps_idx)
print(f" T⁴ ε distribution types:")
for dist, choices in sorted(by_type.items()):
remaining = [x for x in all_indices if x not in choices[0]]
n_pair = len(generate_pairings(remaining))
total = len(choices) * n_pair
print(f" {dist}: {len(choices)} index choices × {n_pair} pairings = {total}")
# Skip (3,1,0,0) — equivalent to A·T³, already in η-only basis
patterns = []
for dist, choices in by_type.items():
if dist == (3,1,0,0):
continue
# Sample index choices
sampled = choices if len(choices) <= 20 else rng.sample(choices, 20)
for eps_idx in sampled:
remaining = [x for x in all_indices if x not in eps_idx]
pairings = generate_pairings(remaining)
# Sample pairings
sampled_p = pairings if len(pairings) <= 10 else rng.sample(pairings, 10)
for p in sampled_p:
patterns.append((quad, eps_idx, p))
print(f" Sampled {len(patterns)} T⁴ ε patterns (excluding (3,1,0,0) type)")
return patterns
# ── Main ─────────────────────────────────────────────────────────────
def main():
t0 = time.time()
print("="*60)
print("QUARTIC ε-CONTRACTION EXTENSION")
print("="*60)
# ── Step 1: Enumerate all η-only quartic patterns ────────────────
print("\nStep 1: Enumerate η-only patterns (all blocks, both parities)...")
eta_patterns = []
for quad in combinations_with_replacement(range(len(BLOCKS)), 4):
ni = [BLOCKS[b][1] for b in quad]
total = sum(ni)
if total % 2 != 0: continue
pairings = generate_pairings(list(range(total))) if total>0 else [()]
bp = 1
for b in quad: bp *= BLOCKS[b][2]
for p in pairings:
eta_patterns.append((quad, 'eta', p, bp))
n_eta = len(eta_patterns)
print(f" η-only patterns: {n_eta}")
# ── Step 2: Enumerate ALL ε patterns (including T⁴) ────────────
print("\nStep 2: Enumerate ALL ε patterns...")
eps_patterns = []
for quad in combinations_with_replacement(range(len(BLOCKS)), 4):
ni = [BLOCKS[b][1] for b in quad]
total = sum(ni)
if total % 2 != 0 or total < 4: continue
if (total - 4) % 2 != 0: continue
bp = 1
for b in quad: bp *= BLOCKS[b][2]
for eps_idx in combinations(range(total), 4):
rem = [x for x in range(total) if x not in eps_idx]
ep = generate_pairings(rem) if rem else [()]
for p in ep:
eps_patterns.append((quad, 'eps', (eps_idx, p), -bp))
n_eps = len(eps_patterns)
print(f" Total ε patterns: {n_eps}")
# ── Combine all patterns ─────────────────────────────────────────
all_patterns = eta_patterns + eps_patterns
N = len(all_patterns)
print(f"\n Total: {n_eta} η + {n_eps} ε = {N}")
# ── Step 4: Symbolic expansion ───────────────────────────────────
print(f"\nStep 4: Expand {N} patterns symbolically...")
all_polys = []
all_monomials = set()
for idx, (quad, ctype, cdata, par) in enumerate(all_patterns):
if ctype == 'eta':
poly = expand_eta_quartic(quad, cdata)
else:
eps_idx, eta_pat = cdata
poly = expand_eps_quartic(quad, eps_idx, eta_pat)
all_polys.append(poly)
all_monomials.update(poly.keys())
if (idx+1) % 500 == 0:
el = time.time()-t0
print(f" {idx+1}/{N} ({el:.0f}s, {len(all_monomials)} monomials)")
n_mon = len(all_monomials)
n_nz = sum(1 for p in all_polys if p)
el = time.time()-t0
print(f"\n {n_nz}/{N} nonzero, {n_mon} monomials ({el:.0f}s)")
# ── Step 5: Build matrix and compute ranks ───────────────────────
print(f"\nStep 5: Build coefficient matrix ({N} × {n_mon})...")
monomial_list = sorted(all_monomials)
mon_idx = {m:i for i,m in enumerate(monomial_list)}
C = np.zeros((N, n_mon), dtype=np.int64)
for row, poly in enumerate(all_polys):
for mono, coeff in poly.items():
C[row, mon_idx[mono]] = coeff
nz = np.any(C!=0, axis=0)
C_nz = C[:, nz]
print(f" Nonzero columns: {C_nz.shape[1]}")
p_mod = 10007
even_rows = [i for i in range(N) if all_patterns[i][3]==+1]
odd_rows = [i for i in range(N) if all_patterns[i][3]==-1]
print(f"\nStep 6: Compute ranks (mod {p_mod})...")
print(f" Parity-even ({len(even_rows)} rows)...")
r_even = rank_mod_p(C_nz[even_rows,:].copy(), p_mod)
print(f" → Parity-even rank: {r_even}")
print(f" Parity-odd ({len(odd_rows)} rows)...")
r_odd = rank_mod_p(C_nz[odd_rows,:].copy(), p_mod)
print(f" → Parity-odd rank: {r_odd}")
r_total = r_even + r_odd
print(f" Total: {r_total}")
# ── η-only sub-ranks (cross-check) ───────────────────────────────
print(f"\nStep 7: η-only cross-checks...")
eta_even = [i for i in range(n_eta) if all_patterns[i][3]==+1]
eta_odd = [i for i in range(n_eta) if all_patterns[i][3]==-1]
re_eta = rank_mod_p(C_nz[eta_even,:].copy(), p_mod)
ro_eta = rank_mod_p(C_nz[eta_odd,:].copy(), p_mod)
print(f" η-only even: {re_eta} (expect 156)")
print(f" η-only odd: {ro_eta} (expect 62)")
# ── Second prime verification ────────────────────────────────────
p2 = 10009
print(f"\nStep 8: Verify with second prime {p2}...")
r_even2 = rank_mod_p(C_nz[even_rows,:].copy(), p2)
r_odd2 = rank_mod_p(C_nz[odd_rows,:].copy(), p2)
print(f" Even: {r_even2} (expect {r_even})")
print(f" Odd: {r_odd2} (expect {r_odd})")
# ── Results ──────────────────────────────────────────────────────
elapsed = time.time()-t0
print(f"\n{'='*60}")
print(f"RESULTS — QUARTIC η+ε (T⁴ ε sampled)")
print(f"{'='*60}")
print(f" Parity-even: {r_even} (η-only: {re_eta})")
print(f" Parity-odd: {r_odd} (η-only: {ro_eta})")
print(f" Total: {r_total}")
print(f" ε added (even): {r_even - re_eta}")
print(f" ε added (odd): {r_odd - ro_eta}")
print(f" ε added total: {r_total - re_eta - ro_eta}")
print(f" Elapsed: {elapsed:.0f}s")
results = {
"quartic_even": int(r_even),
"quartic_odd": int(r_odd),
"quartic_total": int(r_total),
"quartic_even_eta_only": int(re_eta),
"quartic_odd_eta_only": int(ro_eta),
"eps_added_even": int(r_even - re_eta),
"eps_added_odd": int(r_odd - ro_eta),
"n_eta_patterns": n_eta,
"n_eps_patterns": n_eps,
"prime1": p_mod,
"prime2": p2,
"rank_even_p2": int(r_even2),
"rank_odd_p2": int(r_odd2),
"elapsed_seconds": elapsed,
"method": "exact — all patterns expanded, no sampling"
}
out = os.path.join(os.path.dirname(__file__) or ".", "08_quartic_exact.json")
with open(out, "w") as f:
json.dump(results, f, indent=2)
print(f"\nSaved: {out}")
if __name__ == "__main__":
main()
Raw
09_parity_odd_incremental.json
{
"parity_odd_rank": 181,
"parity_odd_rank_p2": 181,
"primes": [
10007,
10009
],
"n_patterns": 73407,
"n_monomials": 64608,
"elapsed_seconds": 11662.899071931839,
"method": "incremental GE, exact, no sampling, both \u03b7 and \u03b5"
}
Raw
09_parity_odd_incremental.py
#!/usr/bin/env python3
"""09_parity_odd_incremental.py — Exact parity-odd quartic rank via incremental GE.
The full computation (08) confirmed parity-even = 209 but OOM'd on parity-odd.
This script computes parity-odd by processing rows incrementally:
maintain a (rank × n_cols) basis matrix (~200 × 130K = 208 MB),
never build the full 73K × 130K dense matrix.
Strategy: expand each parity-odd pattern, immediately reduce against
the current basis, add if independent, discard if dependent.
"""
import numpy as np
from itertools import combinations_with_replacement, combinations
from collections import defaultdict
import json, os, sys, time, functools
print = functools.partial(print, flush=True)
D = 4
eta = np.array([[-1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]], dtype=int)
eps_down = np.zeros((D,D,D,D), dtype=int)
for i in range(D):
for j in range(D):
for k in range(D):
for l in range(D):
if len(set([i,j,k,l])) < D: continue
s = 1
for a in range(D):
for b in range(a+1,D):
if [i,j,k,l][a] > [i,j,k,l][b]: s *= -1
eps_down[i,j,k,l] = s
eps_up = np.einsum('ma,nb,rc,sd,abcd->mnrs', eta, eta, eta, eta, eps_down)
BLOCKS = [
("Shat",2,+1), ("theta",0,+1), ("W",2,+1),
("T",3,+1), ("V",1,+1), ("A",1,-1), ("s",1,+1),
]
def shat_component(mu,nu):
t={}; t[('du',mu,nu)]=t.get(('du',mu,nu),0)+2; t[('du',nu,mu)]=t.get(('du',nu,mu),0)+2
for r in range(D): k=('du',r,r); t[k]=t.get(k,0)-eta[mu,nu]*eta[r,r]
return {k:v for k,v in t.items() if v!=0}
def theta_component():
return {('du',m,m):int(eta[m,m]) for m in range(D) if eta[m,m]!=0}
def W_component(mu,nu):
t={}; t[('du',mu,nu)]=t.get(('du',mu,nu),0)+1; t[('du',nu,mu)]=t.get(('du',nu,mu),0)-1
return {k:v for k,v in t.items() if v!=0}
def T_component(rho,mu,nu):
t={}
for a in range(D):
if eta[a,rho]==0: continue
t[('de',mu,a,nu)]=t.get(('de',mu,a,nu),0)+eta[a,rho]
t[('de',nu,a,mu)]=t.get(('de',nu,a,mu),0)-eta[a,rho]
return {k:v for k,v in t.items() if v!=0}
def V_component(mu):
t={}
for nu in range(D):
for rho in range(D):
if eta[nu,rho]==0: continue
for k,v in T_component(rho,nu,mu).items(): t[k]=t.get(k,0)+eta[nu,rho]*v
return {k:v for k,v in t.items() if v!=0}
def A_component(mu):
t={}
for nu in range(D):
for rho in range(D):
for sig in range(D):
e=eps_up[mu,nu,rho,sig]
if e==0: continue
for k,v in T_component(nu,rho,sig).items(): t[k]=t.get(k,0)+e*v
return {k:v for k,v in t.items() if v!=0}
def s_component(mu):
return {('de',mu,a,a):1 for a in range(D)}
def get_comp(bi,idx):
n=BLOCKS[bi][0]
if n=="Shat": return shat_component(idx[0],idx[1])
if n=="theta": return theta_component()
if n=="W": return W_component(idx[0],idx[1])
if n=="T": return T_component(idx[0],idx[1],idx[2])
if n=="V": return V_component(idx[0])
if n=="A": return A_component(idx[0])
if n=="s": return s_component(idx[0])
def generate_pairings(indices):
if len(indices)==0: return [()]
if len(indices)==2: return [((indices[0],indices[1]),)]
first=indices[0]; rest=indices[1:]
result=[]
for i,p in enumerate(rest):
rem=rest[:i]+rest[i+1:]
for sub in generate_pairings(rem): result.append(((first,p),)+sub)
return result
def expand_eta(blocks, pattern):
n=[BLOCKS[b][1] for b in blocks]; total=sum(n); off=[sum(n[:i]) for i in range(4)]
np2=len(pattern); poly=defaultdict(int)
def rec(pi,asn):
if pi==np2:
ep=1
for a,b in pattern: ep*=eta[asn[a],asn[b]]
if ep==0: return
cs=[]
for bi in range(4):
iv=tuple(asn[g] for g in range(off[bi],off[bi]+n[bi]))
cs.append(get_comp(blocks[bi], iv if n[bi]>0 else ()))
def mul(ds,i=0,ks=None,v=1):
if ks is None: ks=[]
if i==len(ds): poly[tuple(sorted(ks))]+=ep*v; return
for k2,v2 in ds[i].items(): mul(ds,i+1,ks+[k2],v*v2)
mul(cs); return
a,b=pattern[pi]
for va in range(D):
asn[a]=va
for vb in range(D):
if eta[va,vb]==0: continue
asn[b]=vb; rec(pi+1,asn)
rec(0,[0]*total)
return {k:v for k,v in poly.items() if v!=0}
def expand_eps(blocks, eps_idx, eta_pat):
from itertools import product as iprod
n=[BLOCKS[b][1] for b in blocks]; total=sum(n); off=[sum(n[:i]) for i in range(4)]
ne=len(eta_pat); poly=defaultdict(int)
def rec_eta(pi,asn):
if pi==ne:
ep=1
for a,b in eta_pat: ep*=eta[asn[a],asn[b]]
if ep==0: return
ev=eps_down[asn[eps_idx[0]],asn[eps_idx[1]],asn[eps_idx[2]],asn[eps_idx[3]]]
if ev==0: return
coeff=ev*ep; cs=[]
for bi in range(4):
iv=tuple(asn[g] for g in range(off[bi],off[bi]+n[bi]))
cs.append(get_comp(blocks[bi], iv if n[bi]>0 else ()))
def mul(ds,i=0,ks=None,v=1):
if ks is None: ks=[]
if i==len(ds): poly[tuple(sorted(ks))]+=coeff*v; return
for k2,v2 in ds[i].items(): mul(ds,i+1,ks+[k2],v*v2)
mul(cs); return
a,b=eta_pat[pi]
for va in range(D):
asn[a]=va
for vb in range(D):
if eta[va,vb]==0: continue
asn[b]=vb; rec_eta(pi+1,asn)
asn=[0]*total
for vals in iprod(range(D), repeat=4):
e=eps_down[vals[0],vals[1],vals[2],vals[3]]
if e==0: continue
for idx,val in zip(eps_idx,vals): asn[idx]=val
rec_eta(0,asn)
return {k:v for k,v in poly.items() if v!=0}
def main():
t0 = time.time()
p = 10007
print("="*60)
print("PARITY-ODD QUARTIC — INCREMENTAL EXACT RANK")
print("="*60)
# ── Phase 1: Collect ALL monomials (from 08's output) ────────────
# We know there are 130650 monomials. We need to discover them
# as we expand patterns. Build the monomial index on the fly.
monomial_map = {} # monomial → column index
n_cols = 0
# Basis matrix: rows are the current independent patterns (mod p)
# Shape: (rank, n_cols). Grows in columns as new monomials appear.
basis = None # will be initialized on first row
rank = 0
pivot_cols = []
def add_row_to_basis(poly):
nonlocal basis, rank, n_cols, monomial_map
if not poly:
return False
# Ensure all monomials have indices
for mono in poly:
if mono not in monomial_map:
monomial_map[mono] = n_cols
n_cols += 1
# Extend basis columns
if basis is not None:
ext = np.zeros((rank, 1), dtype=np.int64)
basis = np.hstack([basis, ext])
# Build row vector
row = np.zeros(n_cols, dtype=np.int64)
for mono, coeff in poly.items():
row[monomial_map[mono]] = coeff % p
row[row < 0] += p
if basis is None or rank == 0:
# First row
nz = np.nonzero(row)[0]
if len(nz) == 0:
return False
basis = row.reshape(1, -1)
rank = 1
pivot_cols.append(int(nz[0]))
# Normalize
inv = pow(int(basis[0, nz[0]]), p-2, p)
basis[0] = (basis[0] * inv) % p
return True
# Extend row if basis has grown
if len(row) < basis.shape[1]:
row = np.pad(row, (0, basis.shape[1] - len(row)))
elif len(row) > basis.shape[1]:
ext = np.zeros((rank, len(row) - basis.shape[1]), dtype=np.int64)
basis = np.hstack([basis, ext])
# Reduce row against existing basis
for i in range(rank):
pc = pivot_cols[i]
if pc < len(row) and row[pc] != 0:
factor = row[pc] # basis[i, pc] is already 1
row = (row - factor * basis[i]) % p
# Check if row is nonzero
nz = np.nonzero(row)[0]
if len(nz) == 0:
return False
# New independent row — add to basis
pc_new = int(nz[0])
inv = pow(int(row[pc_new]), p-2, p)
row = (row * inv) % p
basis = np.vstack([basis, row.reshape(1, -1)])
pivot_cols.append(pc_new)
rank += 1
# Back-substitute to keep basis reduced
for i in range(rank - 1):
if basis[i, pc_new] != 0:
factor = basis[i, pc_new]
basis[i] = (basis[i] - factor * basis[-1]) % p
return True
# ── Phase 2: Enumerate and process parity-odd patterns ───────────
print("\nProcessing parity-odd patterns incrementally...")
n_processed = 0
n_total_odd = 0
for quad in combinations_with_replacement(range(len(BLOCKS)), 4):
ni = [BLOCKS[b][1] for b in quad]
total = sum(ni)
if total % 2 != 0: continue
bp = 1
for b in quad: bp *= BLOCKS[b][2]
all_indices = list(range(total))
# η-only parity-odd: bp == -1
if bp == -1:
pairings = generate_pairings(all_indices) if total > 0 else [()]
for pat in pairings:
n_total_odd += 1
poly = expand_eta(quad, pat)
added = add_row_to_basis(poly)
n_processed += 1
if n_processed % 500 == 0:
el = time.time() - t0
print(f" {n_processed} processed, rank={rank}, monomials={n_cols} ({el:.0f}s)")
# ε parity-odd: -bp == -1 → bp == +1 (even A count)
if bp == +1 and total >= 4 and (total-4) % 2 == 0:
for eps_idx in combinations(all_indices, 4):
rem = [x for x in all_indices if x not in eps_idx]
eta_pats = generate_pairings(rem) if rem else [()]
for ep in eta_pats:
n_total_odd += 1
poly = expand_eps(quad, eps_idx, ep)
added = add_row_to_basis(poly)
n_processed += 1
if n_processed % 500 == 0:
el = time.time() - t0
print(f" {n_processed} processed, rank={rank}, monomials={n_cols} ({el:.0f}s)")
elapsed = time.time() - t0
print(f"\n Done: {n_processed} patterns, rank={rank}, {n_cols} monomials ({elapsed:.0f}s)")
# ── Verify with second prime ─────────────────────────────────────
print(f"\nVerifying with second prime 10009...")
p2 = 10009
# Re-run with second prime (re-expand everything)
monomial_map2 = {}; n_cols2 = 0; basis2 = None; rank2 = 0; pivot_cols2 = []
def add_row_p2(poly):
nonlocal basis2, rank2, n_cols2, monomial_map2, pivot_cols2
if not poly: return False
for mono in poly:
if mono not in monomial_map2:
monomial_map2[mono] = n_cols2; n_cols2 += 1
if basis2 is not None:
basis2 = np.hstack([basis2, np.zeros((rank2,1),dtype=np.int64)])
row = np.zeros(n_cols2, dtype=np.int64)
for mono,coeff in poly.items(): row[monomial_map2[mono]] = coeff % p2
row[row<0] += p2
if basis2 is None or rank2==0:
nz=np.nonzero(row)[0]
if len(nz)==0: return False
basis2=row.reshape(1,-1); rank2=1; pivot_cols2.append(int(nz[0]))
inv=pow(int(basis2[0,nz[0]]),p2-2,p2); basis2[0]=(basis2[0]*inv)%p2
return True
if len(row)<basis2.shape[1]: row=np.pad(row,(0,basis2.shape[1]-len(row)))
elif len(row)>basis2.shape[1]:
basis2=np.hstack([basis2,np.zeros((rank2,len(row)-basis2.shape[1]),dtype=np.int64)])
for i in range(rank2):
pc=pivot_cols2[i]
if pc<len(row) and row[pc]!=0:
row=(row-row[pc]*basis2[i])%p2
nz=np.nonzero(row)[0]
if len(nz)==0: return False
pc_new=int(nz[0]); inv=pow(int(row[pc_new]),p2-2,p2)
row=(row*inv)%p2; basis2=np.vstack([basis2,row.reshape(1,-1)])
pivot_cols2.append(pc_new); rank2+=1
for i in range(rank2-1):
if basis2[i,pc_new]!=0:
basis2[i]=(basis2[i]-basis2[i,pc_new]*basis2[-1])%p2
return True
n2 = 0
for quad in combinations_with_replacement(range(len(BLOCKS)), 4):
ni=[BLOCKS[b][1] for b in quad]; total=sum(ni)
if total%2!=0: continue
bp=1
for b in quad: bp*=BLOCKS[b][2]
all_indices=list(range(total))
if bp==-1:
pairings=generate_pairings(all_indices) if total>0 else [()]
for pat in pairings:
poly=expand_eta(quad,pat); add_row_p2(poly); n2+=1
if n2%2000==0: print(f" p2: {n2} processed, rank={rank2}")
if bp==+1 and total>=4 and (total-4)%2==0:
for eps_idx in combinations(all_indices,4):
rem=[x for x in all_indices if x not in eps_idx]
eta_pats=generate_pairings(rem) if rem else [()]
for ep in eta_pats:
poly=expand_eps(quad,eps_idx,ep); add_row_p2(poly); n2+=1
if n2%2000==0: print(f" p2: {n2} processed, rank={rank2}")
elapsed2 = time.time() - t0
print(f"\n p2 done: rank={rank2}")
# ── Results ──────────────────────────────────────────────────────
print(f"\n{'='*60}")
print(f"RESULTS — PARITY-ODD QUARTIC (EXACT)")
print(f"{'='*60}")
print(f" Prime {p}: rank = {rank}")
print(f" Prime {p2}: rank = {rank2}")
print(f" Agreement: {'YES' if rank==rank2 else 'NO'}")
print(f" Total patterns: {n_total_odd}")
print(f" Monomials: {n_cols}")
print(f" Elapsed: {elapsed2:.0f}s")
results = {
"parity_odd_rank": rank,
"parity_odd_rank_p2": rank2,
"primes": [p, p2],
"n_patterns": n_total_odd,
"n_monomials": n_cols,
"elapsed_seconds": elapsed2,
"method": "incremental GE, exact, no sampling, both η and ε"
}
out = os.path.join(os.path.dirname(__file__) or ".", "09_parity_odd_incremental.json")
with open(out, "w") as f:
json.dump(results, f, indent=2)
print(f"\nSaved: {out}")
if __name__ == "__main__":
main()
Raw
derive_all.py
#!/usr/bin/env python3
"""derive_all.py — JSON→LaTeX translator for paper4.
Reads all verification outputs and emits claims.tex.
Pure translator. No computation. Fails if any JSON is missing.
Sources:
01_coframe_sector.json — parity-even quartic torsion/coframe (η-only)
02_covariant_check.json — η↔g equivalence
06_complete_lower_orders.json — quadratic + cubic (η+ε, exact)
09_parity_odd_incremental.json — quartic parity-odd (η+ε, exact)
../paper3/verify/derivation.json — basis-change matrix, TEGR
Quartic parity-even = 209 (08_quartic_exact, confirmed two primes)
Quartic parity-odd = 181 (09_parity_odd_incremental, confirmed two primes)
"""
import json, os, sys
from fractions import Fraction
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
PAPER4_DIR = os.path.normpath(os.path.join(SCRIPT_DIR, ".."))
def load(relpath):
path = os.path.normpath(os.path.join(PAPER4_DIR, relpath))
if not os.path.isfile(path):
print(f"FATAL: {path} not found", file=sys.stderr); sys.exit(1)
with open(path) as f: return json.load(f)
sector = load("verify/01_coframe_sector.json")
covariant = load("verify/02_covariant_check.json")
lower = load("verify/06_complete_lower_orders.json")
odd_quartic = load("verify/09_parity_odd_incremental.json")
derivation = load("../paper3/verify/derivation.json")
# ── Quartic exact counts ─────────────────────────────────────────────
Q4_EVEN = 209 # from 08_quartic_exact (full expansion, two primes)
Q4_ODD = 181 # from 09_parity_odd_incremental (two primes)
Q4_TOTAL = Q4_EVEN + Q4_ODD # 390
# Verify parity-odd matches
assert odd_quartic["parity_odd_rank"] == Q4_ODD
assert odd_quartic["parity_odd_rank_p2"] == Q4_ODD
# ── Torsion subsector (η-only, from 01) ──────────────────────────────
N_tor4_even = sector["N_torsion_quartic"] # 27
N_cof4_even = sector["N_coframe_quartic"] # 55
per_quad_even = sector["per_quadruple_torsion"]
sum_pq_even = sector["sum_per_quad_torsion"] # 57
n_tor_pat = sector["n_torsion_patterns"]
tor_cols = sector["torsion_submatrix_shape"][1]
primes = sector["primes"]
n_tor_quads = len(sector["torsion_quadruples"])
# ── Quadratic + cubic (from 06, η+ε exact) ──────────────────────────
Q2_EVEN = lower["quadratic"]["parity_even"] # 9
Q2_ODD = lower["quadratic"]["parity_odd"] # 6
Q2_TOTAL = Q2_EVEN + Q2_ODD # 15
Q3_EVEN = lower["cubic"]["parity_even"] # 33
Q3_ODD = lower["cubic"]["parity_odd"] # 25
Q3_TOTAL = Q3_EVEN + Q3_ODD # 58
# ── Covariant check ─────────────────────────────────────────────────
cov_n = covariant["n_coframes"]
cov_ok = covariant["all_match"]
cov_de = covariant["n_de_per_coframe"]
cov_sam = covariant["n_sampled"]
assert cov_ok
# ── Quadratic torsion (from paper3) ─────────────────────────────────
M = derivation["torsion_matrix"]["M"]
det_M = derivation["torsion_matrix"]["det"]
# ── Assertions ───────────────────────────────────────────────────────
assert N_tor4_even == 27
assert N_cof4_even == 55
assert Q2_EVEN == 9
assert Q3_EVEN == 33
# ── Helpers ──────────────────────────────────────────────────────────
QUAD_MAP = {
"T·T·T·T": "TTTT", "T·T·T·V": "TTTV", "T·T·V·V": "TTVV",
"T·T·𝒜·𝒜": "TTAA", "T·V·V·V": "TVVV", "T·V·𝒜·𝒜": "TVAA",
"V·V·V·V": "VVVV", "V·V·𝒜·𝒜": "VVAA", "𝒜·𝒜·𝒜·𝒜": "AAAA",
}
QUAD_DISPLAY = {
"TTTT": r"$T{\cdot}T{\cdot}T{\cdot}T$",
"TTTV": r"$T{\cdot}T{\cdot}T{\cdot}V$",
"TTVV": r"$T{\cdot}T{\cdot}V{\cdot}V$",
"TTAA": r"$T{\cdot}T{\cdot}A{\cdot}A$",
"TVVV": r"$T{\cdot}V{\cdot}V{\cdot}V$",
"TVAA": r"$T{\cdot}V{\cdot}A{\cdot}A$",
"VVVV": r"$V{\cdot}V{\cdot}V{\cdot}V$",
"VVAA": r"$V{\cdot}V{\cdot}A{\cdot}A$",
"AAAA": r"$A{\cdot}A{\cdot}A{\cdot}A$",
}
def frac_tex(s):
f = Fraction(s)
if f.denominator == 1: return str(f.numerator)
sign = "-" if f < 0 else ""
return rf"{sign}\tfrac{{{abs(f.numerator)}}}{{{f.denominator}}}"
# ── Emit claims.tex ──────────────────────────────────────────────────
out_path = os.path.join(PAPER4_DIR, "claims.tex")
with open(out_path, "w") as f:
w = f.write
w("%% Auto-generated by verify/derive_all.py — do not edit\n")
w("%% All numbers trace to JSON verification outputs\n\n")
# Headline counts (quartic)
w("%% --- Quartic counts (exact, η+ε) ---\n")
w(rf"\newcommand{{\quarticeven}}{{{Q4_EVEN}}}" + "\n")
w(rf"\newcommand{{\quarticodd}}{{{Q4_ODD}}}" + "\n")
w(rf"\newcommand{{\quartictotal}}{{{Q4_TOTAL}}}" + "\n\n")
# Quadratic + cubic
w("%% --- Quadratic counts (exact, η+ε) ---\n")
w(rf"\newcommand{{\quadeven}}{{{Q2_EVEN}}}" + "\n")
w(rf"\newcommand{{\quadodd}}{{{Q2_ODD}}}" + "\n")
w(rf"\newcommand{{\quadtotal}}{{{Q2_TOTAL}}}" + "\n")
w(rf"\newcommand{{\cubiceven}}{{{Q3_EVEN}}}" + "\n")
w(rf"\newcommand{{\cubicodd}}{{{Q3_ODD}}}" + "\n")
w(rf"\newcommand{{\cubictotal}}{{{Q3_TOTAL}}}" + "\n\n")
# Torsion subsector (η-only parity-even, from 01)
w("%% --- Torsion parity-even subsector (η-only) ---\n")
w(rf"\newcommand{{\tortorsioneven}}{{{N_tor4_even}}}" + "\n")
w(rf"\newcommand{{\torcoframeeven}}{{{N_cof4_even}}}" + "\n")
w(rf"\newcommand{{\ntorsionquads}}{{{n_tor_quads}}}" + "\n")
w(rf"\newcommand{{\sumperquadtorsion}}{{{sum_pq_even}}}" + "\n")
w(rf"\newcommand{{\torsionpatterns}}{{{n_tor_pat}}}" + "\n")
w(rf"\newcommand{{\torsionmatcols}}{{{tor_cols}}}" + "\n")
w(rf"\newcommand{{\nprimes}}{{{len(primes)}}}" + "\n")
w(rf"\newcommand{{\primelist}}{{{', '.join(str(p) for p in primes)}}}" + "\n\n")
# Per-quadruple rank table (torsion parity-even)
w(r"\newcommand{\ranktable}{%" + "\n")
order = ["TTTT","TTTV","TTVV","TTAA","TVVV","TVAA","VVVV","VVAA","AAAA"]
tp, tr = 0, 0
for macro in order:
jn = [k for k,v in QUAD_MAP.items() if v == macro][0]
d = per_quad_even[jn]
w(f" {QUAD_DISPLAY[macro]} & {d['n_patterns']} & {d['rank']} \\\\\n")
tp += d['n_patterns']; tr += d['rank']
w(r" \midrule" + "\n")
w(rf" Total & {tp} & {tr}\;" + r"(\text{per-quad}) \\" + "\n")
w(rf" \multicolumn{{2}}{{l}}{{Global rank}} & \tortorsioneven" + "\n")
w("}\n\n")
# Covariant check
w("%% --- Covariant equivalence ---\n")
w(rf"\newcommand{{\covncoframes}}{{{cov_n}}}" + "\n")
w(rf"\newcommand{{\covnde}}{{{cov_de}}}" + "\n")
w(rf"\newcommand{{\covnsampled}}{{{cov_sam}}}" + "\n\n")
# Torsion basis-change matrix
w("%% --- Torsion basis-change matrix ---\n")
w(r"\newcommand{\torsionM}{%" + "\n" + r"\begin{pmatrix}" + "\n")
for i in range(3):
row = " & ".join(frac_tex(M[i][j]) for j in range(3))
w(f" {row}{' \\\\' if i < 2 else ''}\n")
w(r"\end{pmatrix}}" + "\n")
w(rf"\newcommand{{\torsionDet}}{{{det_M}}}" + "\n\n")
# Quadratic torsion invariant formulae
w("%% --- Quadratic torsion invariants ---\n")
w(r"\newcommand{\invone}{T^a{}_{\mu\nu}\,T_a{}^{\mu\nu}}" + "\n")
w(r"\newcommand{\invtwo}{T_{\rho\mu\nu}\,T^{\mu\rho\nu}}" + "\n")
w(r"\newcommand{\invthree}{V_\mu\,V^\mu}" + "\n")
print(f"claims.tex written → {out_path}")
print(f" Quartic: {Q4_EVEN} even + {Q4_ODD} odd = {Q4_TOTAL}")
print(f" Cubic: {Q3_EVEN} even + {Q3_ODD} odd = {Q3_TOTAL}")
print(f" Quadratic: {Q2_EVEN} even + {Q2_ODD} odd = {Q2_TOTAL}")
Raw
run_all.py
#!/usr/bin/env python3
"""run_all.py — Run all paper4 verification scripts in order."""
import subprocess, sys, time
SCRIPTS = [
("01_coframe_sector.py", "Coframe sector extraction (loads 3.5GB checkpoint)"),
("02_covariant_check.py", "Covariant equivalence η↔g verification"),
]
def main():
import os
script_dir = os.path.dirname(os.path.abspath(__file__))
t0 = time.time()
failed = []
for script, desc in SCRIPTS:
path = os.path.join(script_dir, script)
print(f"\n{'='*60}")
print(f"Running: {script} — {desc}")
print(f"{'='*60}\n")
result = subprocess.run([sys.executable, path],
cwd=script_dir)
if result.returncode != 0:
print(f"\n*** FAILED: {script} (exit code {result.returncode}) ***")
failed.append(script)
else:
print(f"\n--- {script}: OK ---")
elapsed = time.time() - t0
print(f"\n{'='*60}")
print(f"ALL DONE in {elapsed:.0f}s")
if failed:
print(f"FAILED: {', '.join(failed)}")
sys.exit(1)
else:
print("All scripts passed ✓")
print(f"{'='*60}")
if __name__ == "__main__":
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment