The script src/21-gravity/validate_counting.py of the validation package of the paper 21-grav, rendered with an anchor per line. ▶ run the predicates on Google Colab.
1#!/usr/bin/env python32"""Unified counting checks (round-01, W-J / W-N): exact integer arithmetic only.341. Christoffel registration pattern: at rational rate a/b the drive registers5 exactly a advances per recurrence b, with gaps taking exactly two values6 floor(b/a), ceil(b/a) of mean b/a (Lemma isopair's tally identity).72. Censored geometric identity (the [approx] chart of the gap spectrum),8 checked symbolically over the rationals.93. The readout worked instance of the quantum paper (Omega=641, g=3):10 s_4 = 114, R_4 = 23/54 and 2/27, largest-remainder allocation 35+34+6+6 = 81.114. Negative control: naive prefix counts violate unit discrepancy already at12 B = 108 (Kesten), so the theorem is stated over complete recurrences.135. Per-stratum split at B = 2, 3 on Omega = 641: over the complete dial14 recurrence the outcome-string tallies equal |dial| times the15 Galois-symmetric multinomial coefficients of (W+ + W-)^B = 4^B,16 with W(j) = 2 +/- r_j the dial-shifted Bell weights, exactly in Z[sqrt(2)]17 and as residues mod 641.18"""19from fractions import Fraction20from itertools import product2122ok = True23def chk(name, cond):24 global ok25 ok = ok and bool(cond)26 print(('PASS ' if cond else 'FAIL ') + name)2728# ---- 1. Christoffel words: five (a,b) pairs incl. 23/54 ----29for a, b in ((23, 54), (2, 27), (1, 5), (3, 8), (7, 12)):30 marks = [((k + 1) * a) // b - (k * a) // b for k in range(b)] # 0/1 registration word31 regs = sum(marks)32 gaps, run = [], 033 idx = [k for k in range(b) if marks[k] == 1]34 gaps = [((idx[(i + 1) % len(idx)] - idx[i]) % b) or b for i in range(len(idx))]35 gapset = sorted(set(gaps))36 chk(f"christoffel a/b={a}/{b}: a regs per recurrence, two-valued gaps",37 regs == a and set(gapset) <= {b // a, -(-b // a)} and38 Fraction(sum(gaps), len(gaps)) == Fraction(b, a))3940# ---- 2. censored identity E[min(G,T)] = (1-(1-q)^T)/q over Q ----41q = Fraction(1, 7); T = 1142lhs = sum(k * q * (1 - q) ** (k - 1) for k in range(1, T)) + T * (1 - q) ** (T - 1)43rhs = (1 - (1 - q) ** T) / q44chk("censored geometric identity (exact rational, q=1/7, T=11)", lhs == rhs)4546# ---- 3. the readout worked instance ----47from math import isqrt48s4 = isqrt(2 * 3 ** 8)49R4 = Fraction(2 * 81 + s4, 8 * 81)50chk("s_4 = 114 and R_4 = 23/54", s4 == 114 and R4 == Fraction(23, 54))51weights = [Fraction(23, 54), Fraction(23, 54), Fraction(2, 27), Fraction(2, 27)]52B = 8153prov = [int(B * w) for w in weights]54resid = B - sum(prov)55rema = sorted(range(4), key=lambda j: (B * weights[j]) - prov[j], reverse=True)56for j in rema[:resid]:57 prov[j] += 158chk("largest-remainder allocation 35+34+6+6 = 81, discrepancy <= 1",59 sorted(prov, reverse=True) == [35, 34, 6, 6] and60 all(abs(prov[j] - B * weights[j]) <= 1 for j in range(4)))6162# ---- 4. negative control: prefix counts against out-of-lattice thresholds ----63# Threshold counts along the raw drive orbit (3^k mod 641 below the grid64# rational's scaled threshold) are NOT discrepancy-bounded (Kesten: bounded65# prefix discrepancy only for thresholds in the orbit's own lattice); the66# theorem is therefore stated over complete recurrences, where counts are forced.67Om_, g_ = 641, 368th = Fraction(23, 54)69x, c, viol_B = 1, 0, None70for B in range(1, 200):71 if x < th * Om_:72 c += 173 x = x * g_ % Om_74 if abs(c - B * th) >= 2 and viol_B is None:75 viol_B = B76chk(f"negative control: raw-orbit prefix count breaks unit discrepancy (first at B={viol_B})",77 viol_B is not None)7879# ---- 5. per-stratum split, B = 2, 3, Omega = 641 ----80Om = 64181r = 67 # sqrt(2) mod 64182assert (r * r) % Om == 283g = 384z8 = pow(g, (Om - 1) // 8, Om) # primitive 8th root85dial = [(pow(z8, 2 * j + 1, Om) + pow(z8, -(2 * j + 1), Om)) % Om for j in range(8)]86chk("dial values are +/- sqrt2 mod 641 and sum to zero",87 set(dial) == {r % Om, (-r) % Om} and sum(dial) % Om == 0)8889# exact Z[sqrt2] arithmetic: elements as (rational part, sqrt2 part)90def mul(x, y):91 return (x[0] * y[0] + 2 * x[1] * y[1], x[0] * y[1] + x[1] * y[0])92for B in (2, 3):93 total = [0] * (2 ** B)94 for j_idx, rj_sign in enumerate([+1, -1] * 4): # dial recurrence: r_j alternates +r,-r95 Wp = (2, rj_sign) # 2 + r_j96 Wm = (2, -rj_sign) # 2 - r_j97 for si, string in enumerate(product((0, 1), repeat=B)):98 t = (1, 0)99 for bit in string:100 t = mul(t, Wp if bit == 0 else Wm)101 # accumulate rational part; sqrt2 parts must cancel over the dial102 total[si] = total[si] + t[0] if isinstance(total[si], int) else total[si]103 # recompute properly with both parts104 tallies = []105 for string in product((0, 1), repeat=B):106 acc = (0, 0)107 for rj_sign in [+1, -1] * 4:108 t = (1, 0)109 for bit in string:110 t = mul(t, (2, rj_sign) if bit == 0 else (2, -rj_sign))111 acc = (acc[0] + t[0], acc[1] + t[1])112 tallies.append((string, acc))113 dialN = 8114 # Galois-symmetric coefficient = rational part of prod W_{s_i} at r_j=+sqrt2115 def sym_coeff(string):116 t = (1, 0)117 for bit in string:118 t = mul(t, (2, 1) if bit == 0 else (2, -1))119 return t[0]120 all_ok = all(acc[1] == 0 and acc[0] == dialN * sym_coeff(st) for st, acc in tallies)121 tot = sum(acc[0] for _, acc in tallies)122 chk(f"B={B}: string tallies = |dial| x Galois-symmetric coefficients; total = 4^B x |dial| ({tot} = {4**B*8})",123 all_ok and tot == (4 ** B) * dialN)124 # residue check mod 641125 res_ok = True126 for st, acc in tallies:127 acc_res = 0128 for rj in dial:129 t = 1130 for bit in st:131 t = t * ((2 + rj) if bit == 0 else (2 - rj)) % Om132 acc_res = (acc_res + t) % Om133 if acc_res != (dialN * sym_coeff(st)) % Om:134 res_ok = False135 chk(f"B={B}: residue tallies mod 641 agree", res_ok)136137print("\nALL PASS" if ok else "\nFAILURES PRESENT")138raise SystemExit(0 if ok else 1)139