verify.py
16.1 kB · python · 474 lines
1from fractions import Fraction2from itertools import product3from math import comb, log, cos, pi45# SLICE67def choose(n, k):8 if k < 0 or k > n or n < 0:9 return 010 return comb(n, k)1112def polymul(a, b):13 out = [0] * (len(a) + len(b) - 1)14 for i, x in enumerate(a):15 if x:16 for j, y in enumerate(b):17 if y:18 out[i + j] += x * y19 return out2021def digit_poly(D):22 p = [1]23 for _ in range(D - 1):24 p = polymul(p, [1, 0, 1])25 return polymul(p, [1, D, 1])2627def brute_poly(D):28 out = [0] * (2 * D + 1)29 for v in product((0, 1, 2), repeat=D):30 if sum(1 for x in v if x == 1) <= 1:31 out[sum(v)] += 132 return out3334def bcoef(D, s):35 if s < 0 or s > 2 * D:36 return 037 if s % 2 == 0:38 return choose(D - 1, s // 2) + choose(D - 1, (s - 2) // 2)39 return D * choose(D - 1, (s - 1) // 2)4041def transfer(D):42 P = digit_poly(D)43 size = 2 * D + 144 M = [[0] * size for _ in range(size)]45 for c in range(-D, D + 1):46 for s in range(0, 2 * D + 1):47 if (c + D - s) % 3 == 0:48 nc = (c + D - s) // 349 if abs(nc) <= D:50 M[nc + D][c + D] += P[s]51 return M5253def m_even(D):54 n = (D + 1) // 255 M = [[0] * n for _ in range(n)]56 for i in range(n):57 M[i][0] = bcoef(D, D - 3 * i)58 for j in range(1, n):59 M[i][j] = bcoef(D, D + j - 3 * i) + bcoef(D, D - j - 3 * i)60 return M6162def matmul(A, B):63 n, m, p = len(A), len(B), len(B[0])64 C = [[0] * p for _ in range(n)]65 for i in range(n):66 Ai = A[i]67 Ci = C[i]68 for k in range(m):69 a = Ai[k]70 if a:71 Bk = B[k]72 for j in range(p):73 if Bk[j]:74 Ci[j] += a * Bk[j]75 return C7677def matpow(A, L):78 n = len(A)79 R = [[1 if i == j else 0 for j in range(n)] for i in range(n)]80 for _ in range(L):81 R = matmul(R, A)82 return R8384def ladder(D, top):85 M = m_even(D)86 R = [[1 if i == j else 0 for j in range(len(M))] for i in range(len(M))]87 out = [1]88 for _ in range(top):89 R = matmul(R, M)90 out.append(R[0][0])91 return out9293def ladder_auto(D, top):94 M = transfer(D)95 R = [[1 if i == j else 0 for j in range(len(M))] for i in range(len(M))]96 out = [1]97 for _ in range(top):98 R = matmul(R, M)99 out.append(R[D][D])100 return out101102def substitution(D, L):103 P = digit_poly(D)104 out = [1]105 for j in range(L):106 step = 3 ** j107 new = [0] * (len(out) + 2 * D * step)108 for i, x in enumerate(out):109 if x:110 for s, y in enumerate(P):111 if y:112 new[i + s * step] += x * y113 out = new114 return out115116def reachable(D):117 M = transfer(D)118 seen = {0}119 frontier = [0]120 while frontier:121 c = frontier.pop()122 for nc in range(-D, D + 1):123 if M[nc + D][c + D] and nc not in seen:124 seen.add(nc)125 frontier.append(nc)126 return sorted(seen)127128def det(rows):129 A = [[Fraction(x) for x in row] for row in rows]130 n = len(A)131 d = Fraction(1)132 for i in range(n):133 piv = None134 for r in range(i, n):135 if A[r][i]:136 piv = r137 break138 if piv is None:139 return Fraction(0)140 if piv != i:141 A[i], A[piv] = A[piv], A[i]142 d = -d143 d *= A[i][i]144 inv = Fraction(1) / A[i][i]145 for r in range(i + 1, n):146 f = A[r][i] * inv147 if f:148 for c in range(i, n):149 A[r][c] -= f * A[i][c]150 return d151152def hankel(D):153 n = (D + 1) // 2154 a = ladder(D, 2 * n)155 return det([[a[i + j + 1] for j in range(n)] for i in range(n)])156157def charpoly(A):158 n = len(A)159 I = [[Fraction(1 if i == j else 0) for j in range(n)] for i in range(n)]160 F = [[Fraction(x) for x in row] for row in A]161 M = I162 coeffs = [Fraction(1)]163 for k in range(1, n + 1):164 AM = matmul(F, M)165 c = -Fraction(sum(AM[i][i] for i in range(n)), k)166 coeffs.append(c)167 M = [[AM[i][j] + (c if i == j else 0) for j in range(n)] for i in range(n)]168 return [int(x) for x in coeffs]169170def roots(coeffs):171 n = len(coeffs) - 1172 z = [complex(0.4, 0.9) ** k for k in range(n)]173 for _ in range(500):174 moved = 0.0175 for i in range(n):176 num = 0j177 for c in coeffs:178 num = num * z[i] + c179 den = 1 + 0j180 for j in range(n):181 if j != i:182 den *= z[i] - z[j]183 step = num / den184 z[i] -= step185 moved = max(moved, abs(step))186 if moved < 1e-14:187 break188 return z189190def gap_ratio(D):191 n = (D + 1) // 2192 if n < 2:193 return None194 fill = 2 ** (D - 1) * (D + 2)195 s = Fraction(fill, 3)196 cp = charpoly(m_even(D))197 scaled = [float(Fraction(cp[k]) / s ** k) for k in range(len(cp))]198 z = sorted((abs(r) for r in roots(scaled)), reverse=True)199 return z[0] / z[1]200201202# SIGN LAW203204def core_transfer(D):205 r = (D - 1) // 2206 states = list(range(-r, r + 1))207 return [[bcoef(D, c + D - 3 * cp) for c in states] for cp in states], r208209def sheaf(D, top):210 M, r = core_transfer(D)211 n = 2 * r + 1212 v = [1 if i == r else 0 for i in range(n)]213 out = [1]214 for _ in range(top):215 v = [sum(M[i][j] * v[j] for j in range(n)) for i in range(n)]216 out.append(sum(v))217 return out218219def sheaf_from_product(D, L):220 coeffs = substitution(D, L)221 NL = D * (3 ** L - 1) // 2222 step = 3 ** L223 return sum(x for T, x in enumerate(coeffs) if (T - NL) % step == 0)224225def phi_circle(D, psi):226 c2 = 2.0 * cos(psi)227 return c2 ** (D - 1) * (D + c2)228229def sheaf_dft(D, L):230 Q = 3 ** L231 tot = 0.0232 for m in range(Q):233 p = 1.0234 for j in range(L):235 p *= phi_circle(D, 2 * pi * ((m * 3 ** j) % Q) / Q)236 tot += p237 return tot / Q238239def main():240 for D in range(2, 9):241 got = digit_poly(D)242 want = brute_poly(D)243 assert got == want, f"D={D}: digit polynomial got {got} want {want}"244 closed = [bcoef(D, s) for s in range(2 * D + 1)]245 assert closed == want, f"D={D}: entry form got {closed} want {want}"246 print(f"D={D}: P(t) factorisation and entry form match enumeration of 3^{D} tuples")247248 for D in range(2, 8):249 auto = ladder_auto(D, 4)250 even = ladder(D, 4)251 subs = [substitution(D, L)[D * (3 ** L - 1) // 2] for L in range(0, 5)]252 assert auto == subs, f"D={D}: automaton got {auto} want {subs}"253 assert even == subs, f"D={D}: even block got {even} want {subs}"254 print(f"D={D}: three generators agree, a(0..4) = {subs}")255256 for D in range(2, 25):257 r = (D - 1) // 2258 got = reachable(D)259 want = list(range(-r, r + 1))260 assert got == want, f"D={D}: reachable carries got {got} want {want}"261 print("D=2..24: reachable carries are exactly {|c| <= floor((D-1)/2)}")262263 for D in range(2, 25):264 n = (D + 1) // 2265 h = hankel(D)266 assert h != 0, f"D={D}: Hankel determinant got 0 want nonzero"267 if D <= 6:268 print(f"D={D}: order {n}, Hankel determinant {h}")269 print("D=2..24: Hankel determinant nonzero, so the order is exactly ceil(D/2)")270271 for D in range(2, 25):272 M = m_even(D)273 got = sum(M[i][i] for i in range(len(M)))274 want = 3 * 2 ** (D - 2) - 1 if D % 2 == 0 else 3 * D * 2 ** (D - 3)275 assert got == want, f"D={D}: trace got {got} want {want}"276 print("D=2..24: trace is 3*2^(D-2)-1 at even D and 3*D*2^(D-3) at odd D")277278 for D in range(2, 13):279 P = digit_poly(D)280 fill = sum(P)281 assert fill == 2 ** (D - 1) * (D + 2), f"D={D}: fill got {fill} want {2 ** (D - 1) * (D + 2)}"282 eps = Fraction((D - 1) * (-1) ** (D - 1), 3)283 got = [sum(P[s] for s in range(len(P)) if s % 3 == j) for j in range(3)]284 want = [Fraction(fill, 3) + (2 * eps if j == D % 3 else -eps) for j in range(3)]285 assert got == want, f"D={D}: class sums got {got} want {want}"286 print("D=2..12: fill = 2^(D-1)(D+2) and the three class sums are fill/3+2eps, fill/3-eps, fill/3-eps")287288 for D in range(3, 8):289 n = (D + 1) // 2290 r = (D - 1) // 2291 cp = charpoly(m_even(D))292 for m in range(1, r + 1):293 b = [substitution(D, L)[D * (3 ** L - 1) // 2 + m] for L in range(1, 2 * n + 4)]294 for start in range(len(b) - n):295 got = sum(cp[k] * b[start + n - k] for k in range(n + 1))296 assert got == 0, f"D={D}, m={m}, L={start}: off-centre residual got {got} want 0"297 print(f"D={D}: off-centre censuses at offsets 1..{r} obey the central recurrence")298299 cp4 = charpoly(m_even(4))300 b4 = [substitution(4, L)[4 * (3 ** L - 1) // 2 + 2] for L in range(1, 9)]301 stray = [sum(cp4[k] * b4[start + 2 - k] for k in range(3)) for start in range(6)]302 assert all(x != 0 for x in stray), f"D=4, m=2: residuals got {stray} want all nonzero"303 print(f"D=4: offset 2 sits on the stalling carry D/2 and breaks the recurrence, residuals {stray}")304305 anchor = m_even(3)306 assert anchor == [[6, 6], [1, 3]], f"D=3: even block got {anchor} want [[6, 6], [1, 3]]"307 cp3 = charpoly(anchor)308 assert cp3 == [1, -9, 12], f"D=3: characteristic polynomial got {cp3} want [1, -9, 12]"309 a3 = ladder(3, 6)310 want3 = [1, 6, 42, 306, 2250, 16578, 122202]311 assert a3 == want3, f"D=3: ladder got {a3} want {want3}"312 rho3 = (9 + 33 ** 0.5) / 2313 dim3 = log(rho3) / log(3)314 assert abs(dim3 - 1.818410) < 1e-6, f"D=3: slice dimension got {dim3} want 1.818410"315 print(f"D=3: [[6,6],[1,3]], x^2-9x+12, ladder {want3}, dim {dim3:.6f}")316317 a2 = ladder(2, 8)318 assert a2 == [2 ** L for L in range(9)], f"D=2: ladder got {a2} want powers of two"319 cp4 = charpoly(m_even(4))320 assert cp4 == [1, -11, -66], f"D=4: characteristic polynomial got {cp4} want [1, -11, -66]"321 a4 = ladder(4, 6)322 want4 = [1, 6, 132, 1848, 29040, 441408, 6772128]323 assert a4 == want4, f"D=4: ladder got {a4} want {want4}"324 a5 = ladder(5, 4)325 want5 = [1, 30, 1000, 35700, 1321600]326 assert a5 == want5, f"D=5: ladder got {a5} want {want5}"327 a6 = ladder(6, 4)328 want6 = [1, 20, 4030, 242300, 24642700]329 assert a6 == want6, f"D=6: ladder got {a6} want {want6}"330 print(f"D=4: ladder {want4}; D=5: {want5}; D=6: {want6}")331332 prev = None333 for D in range(4, 21):334 got = gap_ratio(D)335 assert got > 1.0, f"D={D}: spectral ratio got {got} want above 1"336 if prev is not None and D >= 6:337 assert got < prev, f"D={D}: spectral ratio got {got} want below {prev}"338 if D >= 6:339 free = (D + 2) / (D - 2)340 assert abs(got - free) < 0.05, f"D={D}: spectral ratio got {got} want near {free}"341 prev = got342 print(f"D={D}: rho/|lambda_2| = {got:.6f}")343344 for D in range(2, 21):345 fill = 2 ** (D - 1) * (D + 2)346 cp = charpoly(m_even(D))347 s = Fraction(fill, 3)348 scaled = [float(Fraction(cp[k]) / s ** k) for k in range(len(cp))]349 rho = max(abs(r) for r in roots(scaled))350 got = 1 if rho > 1.0 else -1351 want = (-1) ** (D + 1)352 assert got == want, f"D={D}: sign of rho - fill/3 got {got} want {want}"353 print("D=2..20: sign(rho - fill/3) alternates as (-1)^(D+1)")354355 for D in range(2, 31):356 fill = Fraction(2 ** (D - 1) * (D + 2), 3)357 cp = charpoly(m_even(D))358 val = Fraction(0)359 for c in cp:360 val = val * fill + c361 got = 1 if val > 0 else -1362 want = (-1) ** D363 assert got == want, f"D={D}: sign of chi(fill/3) got {got} want {want}"364 print("D=2..30: exact rational sign of chi(fill/3) alternates as (-1)^D")365366 for D in range(2, 8):367 b = sheaf(D, 5)368 for L in range(1, 5):369 got = sheaf_from_product(D, L)370 assert got == b[L], f"D={D} L={L}: class sum {got} want {b[L]}"371 for L in range(1, 6):372 approx = sheaf_dft(D, L)373 assert abs(approx - b[L]) < 1e-6 * max(1.0, b[L]), f"D={D} L={L}: dft {approx} want {b[L]}"374 print("D=2..7: sheaf census = coefficient class sums = trigonometric product formula")375376 for D in range(2, 8):377 fill = 2 ** (D - 1) * (D + 2)378 b = sheaf(D, 5)379 for k in range(1, 6):380 W = 3 ** (k - 1) * (3 * b[k] - fill * b[k - 1])381 U = [u for u in range(1, 3 ** k) if u % 3]382 direct = 0.0383 for u in U:384 p = 1.0385 for j in range(k):386 p *= phi_circle(D, 2 * pi * ((u * 3 ** j) % 3 ** k) / 3 ** k)387 direct += p388 assert abs(direct - W) < 1e-6 * max(1.0, abs(W)), f"D={D} k={k}: W {W} vs unit sum {direct}"389 par = phi_circle(D, 2 * pi / 3)390 assert abs(par - (-1) ** (D - 1) * (D - 1)) < 1e-9, f"D={D}: parity factor {par}"391 print("D=2..7: step identity W_k matches unit sums, parity factor (-1)^(D-1)(D-1)")392393 for D in range(3, 26, 2):394 fill = 2 ** (D - 1) * (D + 2)395 b = sheaf(D, 8)396 for k in range(1, 9):397 W = 3 ** (k - 1) * (3 * b[k] - fill * b[k - 1])398 assert W > 0, f"D={D} k={k}: W_k = {W} want positive"399 assert 3 ** 8 * b[8] >= fill ** 8, f"D={D}: sheaf census below (fill/3)^8"400 print("odd D=3..25: every W_k > 0 for k<=8 and b(8) >= (fill/3)^8, the theorem's mechanism")401402 for D in range(6, 25, 2):403 fill = 2 ** (D - 1) * (D + 2)404 b = sheaf(D, 2)405 W2 = 3 * (3 * b[2] - fill * b[1])406 assert W2 > 0, f"D={D}: W_2 = {W2} want positive, i.e. V_2 < 0"407 print("even D=6..24: W_2 > 0 so V_2 < 0, the even-side obstruction")408409 for D in range(2, 81):410 fill = 2 ** (D - 1) * (D + 2)411 Me = m_even(D)412 n = len(Me)413 J = [[fill * (i == j) - 3 * Me[i][j] for j in range(n)] for i in range(n)]414 dj = det(J)415 assert dj.denominator == 1 and dj != 0, f"D={D}: det J = {dj} want nonzero integer"416 dji = int(dj)417 assert (dji - fill ** n) % 3 == 0, f"D={D}: det J != fill^n mod 3"418 if D % 3 != 1:419 assert dji % 3 != 0, f"D={D}: det J divisible by 3 with D != 1 mod 3"420 print("D=2..80: det(fill I - 3 M_even) is a nonzero integer, == fill^n mod 3")421422 D = 61423 fill = 2 ** (D - 1) * (D + 2)424 cp = charpoly(m_even(D))425 def chi_at(x):426 v = Fraction(0)427 for c in cp:428 v = v * x + c429 return v430 lo, hi = Fraction(fill, 3), Fraction(fill, 3) + 1431 assert chi_at(lo) < 0 and chi_at(hi) > 0, "D=61: bracketing of rho failed"432 for _ in range(100):433 mid = (lo + hi) / 2434 if chi_at(mid) < 0:435 lo = mid436 else:437 hi = mid438 gap61 = float(lo - Fraction(fill, 3))439 pred = 2 * (D - 1) / 3440 k = 2441 while 3 ** k < 10 ** 13:442 c2 = 2 * cos(2 * pi / 3 ** k)443 pred *= (c2 / 2) ** (D - 1) * (D + c2) / (D + 2)444 k += 1445 assert abs(gap61 / pred - 1) < 1e-8, f"D=61: gap {gap61} vs tower product {pred}"446 print("D=61: exact-bisected rho - fill/3 matches the tower product within 1e-8")447448 for D in range(2, 21):449 fill = 2 ** (D - 1) * (D + 2)450 M, r = core_transfer(D)451 n = 2 * r + 1452 v = [1.0] * n453 for _ in range(300):454 w = [sum(M[i][j] * v[j] for j in range(n)) for i in range(n)]455 top = max(w)456 v = [x / top for x in w]457 w = [sum(M[i][j] * v[j] for j in range(n)) for i in range(n)]458 rho = sum(w) / sum(v)459 eps = (D - 1) * (-1) ** (D - 1) / 3.0460 p = sum(v[i] for i in range(n) if (i - r) % 3 == 0) / sum(v)461 lhs = rho - fill / 3.0462 rhs = eps * (3 * p - 1)463 assert abs(lhs - rhs) < 1e-6 * max(1.0, abs(eps) * 3), \464 "D=%d: mass identity lhs %r rhs %r" % (D, lhs, rhs)465 assert abs(lhs) <= 2 * (D - 1) / 3.0 + 1e-9, "D=%d: pinning violated: %r" % (D, lhs)466 lo = fill / 3.0 - (2 if D % 2 == 0 else 1) * (D - 1) / 3.0 - 1e-9467 hi = fill / 3.0 + (1 if D % 2 == 0 else 2) * (D - 1) / 3.0 + 1e-9468 assert lo <= rho <= hi, "D=%d: parity-refined pinning violated: %r" % (D, rho)469 print("D=2..20: mass identity 3 rho = fill + (-1)^(D-1)(D-1)(3p-1) and both pinning brackets")470471 print("all green")472473if __name__ == "__main__":474 main()