ladder.py

1from fractions import Fraction2from math import comb, log3import mpmath4import numpy as np5import sympy as sp67DIGITS = ((0, 0), (1, 0), (0, 1))8LOG3 = log(3)9KAPPA = 3 - log(5) / LOG310WALL = 2 / (3 + log(5) / LOG3)1112def log3(x):13    return log(x) / LOG31415def beta(kappa):16    return kappa / (2 * kappa + 2 - KAPPA)1718def delta_counts(K):19    types = [(b, c, comb(K, b) * comb(K - b, c)) for b in range(K + 1) for c in range(K + 1 - b)]20    out = {}21    for b, c, m in types:22        for B, C, M in types:23            out[(b - B, c - C)] = out.get((b - B, c - C), 0) + m * M24    return out2526def carry_matrix(K):27    radius = (K - 1) // 228    states = [(i, j) for i in range(-radius, radius + 1) for j in range(-radius, radius + 1)]29    index = {s: i for i, s in enumerate(states)}30    M = [[0] * len(states) for _ in states]31    for s in states:32        for (d1, d2), w in delta_counts(K).items():33            z1, z2 = s[0] + d1, s[1] + d234            if z1 % 3 == 0 and z2 % 3 == 0:35                M[index[s]][index[(z1 // 3, z2 // 3)]] += w36    return M, index[(0, 0)]3738def energies(K, levels):39    M, zero = carry_matrix(K)40    v = [0] * len(M)41    v[zero] = 142    out = [1]43    for _ in range(levels):44        v = [sum(v[i] * M[i][j] for i in range(len(v))) for j in range(len(v))]45        out.append(v[zero])46    return out4748def gasket(a):49    pts = [(0, 0)]50    for l in range(a):51        pts = [(x + dx * 3**l, y + dy * 3**l) for x, y in pts for dx, dy in DIGITS]52    return pts5354def energy_direct(K, a):55    pts = gasket(a)56    top = K * (3**a - 1) + 157    sums = np.zeros((top, top), dtype=np.int64)58    sums[0, 0] = 159    for _ in range(K):60        fresh = np.zeros_like(sums)61        for x, y in pts:62            fresh[x:, y:] += sums[: top - x, : top - y]63        sums = fresh64    return sum(int(v) * int(v) for v in sums.flat if v)6566print("kappa = 3 - log_3 5 = %.6f, beta_0^(4) = %.6f, wall 2/(3 + log_3 5) = %.12f" % (KAPPA, beta(KAPPA), WALL))6768print("energies E_2K(G_a), carry matrix exact integers, a = 0..6")69rows = {}70for K in (2, 3, 4, 5):71    rows[K] = energies(K, 6)72    print("2K = %d:" % (2 * K), rows[K])73print("E_4 = 15^a for a <= 6:", all(rows[2][a] == 15**a for a in range(7)))74for K in (2, 3, 4, 5):75    top = 4 if K < 5 else 376    agree = all(energy_direct(K, a) == rows[K][a] for a in range(1, top + 1))77    print("2K = %d direct convolution agrees at a = 1..%d: %s" % (2 * K, top, agree))7879print("characteristic polynomials and Perron roots")80x = sp.symbols("x")81mpmath.mp.dps = 3082perron = {}83for K in (3, 4, 5):84    M, _ = carry_matrix(K)85    poly = sp.Matrix(M).charpoly(x).as_expr()86    print("2K = %d states %d factors: %s" % (2 * K, len(M), sp.factor(poly)))87    perron[K] = max(mpmath.mpf(str(r)) for f, _ in sp.factor_list(poly)[1] for r in sp.Poly(f, x).nroots(n=25) if r.is_real)88lam6 = 57 + 6 * mpmath.sqrt(46)89lam8 = 456 + 3 * mpmath.sqrt(11017)90lam10 = perron[5]91print("lambda_6 = 57 + 6 sqrt 46 = %s, matches largest factor root %s" % (mpmath.nstr(lam6, 12), abs(lam6 - perron[3]) < 1e-15))92print("lambda_8 = 456 + 3 sqrt 11017 = %s, matches largest factor root %s" % (mpmath.nstr(lam8, 12), abs(lam8 - perron[4]) < 1e-15))93print("lambda_10 largest root of the quartic factor = %s" % mpmath.nstr(lam10, 13))9495def reachable(A, start):96    seen, stack = {start}, [start]97    while stack:98        i = stack.pop()99        for j in range(len(A)):100            if A[i][j] and j not in seen:101                seen.add(j)102                stack.append(j)103    return seen104105def scc_zero(M, zero):106    T = [[M[j][i] for j in range(len(M))] for i in range(len(M))]107    S = sorted(reachable(M, zero) & reachable(T, zero))108    sub = [[M[i][j] for j in S] for i in S]109    return S, sub, all(len(reachable(sub, i)) == len(S) for i in range(len(S)))110111print("energy cap: strongly connected component of the zero state, Perron root, E_2K(G_a) <= lambda_2K^a")112for K in (2, 3, 4, 5):113    M, zero = carry_matrix(K)114    radius = (K - 1) // 2115    spread = max(max(abs(d1), abs(d2)) for d1, d2 in delta_counts(K))116    closed = (radius + spread) // 3 <= radius117    S, sub, irreducible = scc_zero(M, zero)118    root = max(mpmath.mpf(str(r)) for f, _ in sp.factor_list(sp.Matrix(sub).charpoly(x).as_expr())[1] for r in sp.Poly(f, x).nroots(n=25) if r.is_real)119    lam = mpmath.mpf(15) if K == 2 else perron[K]120    print("2K = %2d box radius %d digit spread %d closed %s scc %2d of %2d irreducible %s self loop %s perron equals lambda_2K %s cap holds a = 0..6 %s ratio at a = 6 %s" % (2 * K, radius, spread, closed, len(S), len(M), irreducible, M[zero][zero] > 0, abs(root - lam) < mpmath.mpf(10) ** -18, all(mpmath.mpf(rows[K][a]) <= lam**a for a in range(7)), mpmath.nstr(mpmath.mpf(rows[K][6]) / lam**6, 6)))121122QUARTIC = sp.Poly(x**4 - 7833 * x**3 + 7916949 * x**2 - 850684437 * x + 13054946580, x)123LO, HI = sp.Rational(66641136625, 10**7), sp.Rational(66641136626, 10**7)124print("lambda_10 certified by Sturm on the exact quartic: roots above %s %d, roots in the bracket %d, width %s" % (HI, QUARTIC.count_roots(HI, sp.oo), QUARTIC.count_roots(LO, HI), sp.nsimplify(HI - LO)))125KAP10 = 10 - mpmath.log(mpmath.mpf(HI.p) / HI.q) / mpmath.log(3)126127def truncate(v, d):128    q = int(mpmath.floor(v * mpmath.mpf(10) ** d))129    return "%d.%0*d" % (q // 10**d, d, q % 10**d)130131print("certified kappa_10 >= %s, beta_0^(10) >= %s, both truncated down at 12 places" % (truncate(KAP10, 12), truncate(beta(KAP10), 12)))132print("safe short edge, truncated down at 7 places: %s" % truncate(beta(KAP10), 7))133134print("ladder rungs: moments, kappa_2K, beta_0^(2K)")135kappas = {}136for K, lam in ((2, 15), (3, lam6), (4, lam8), (5, lam10)):137    kappas[K] = 2 * K - mpmath.log(lam) / mpmath.log(3)138    print("2K = %2d kappa %s beta %s" % (2 * K, mpmath.nstr(kappas[K], 13), mpmath.nstr(beta(kappas[K]), 12)))139print("tenth rung above eighth by %s" % mpmath.nstr(beta(kappas[5]) - beta(kappas[4]), 9))140print("Holder block exponents 2K = 6, 8, 10:", ", ".join(str(Fraction(2, 1) / (1 - Fraction(1, 2 * K))) for K in (3, 4, 5)))141142print("rows above the tenth, floating Perron roots of exact matrices")143for K in range(6, 11):144    M, _ = carry_matrix(K)145    lam = max(np.linalg.eigvals(np.array(M, dtype=float)).real)146    kappa = 2 * K - log3(lam)147    print("2K = %2d states %3d lambda %.6f beta %.9f wall - beta %.2e" % (2 * K, len(M), lam, beta(kappa), WALL - beta(kappa)))148149def primes(lo, hi):150    return list(sp.primerange(lo, hi + 1))151152def fourier(p, n):153    t = np.arange(p)154    e = np.exp(2j * np.pi * t / p)155    W = np.abs(1 + e[:, None] + e[None, :]) / 3156    F = np.ones((p, p))157    out = []158    for l in range(n):159        I = (t * pow(3, l, p)) % p160        F = F * W[np.ix_(I, I)]161        out.append(F.copy())162    return out163164print("moment identities on (Z/p)^2, primes 5..199")165checks, worst = 0, 0.0166for p in primes(5, 199):167    F = fourier(p, 6)168    for a in range(1, 6):169        targets = []170        if 3**a <= p:171            targets.append((2, p * p * 3.0 ** (-a)))172        if 2 * 3**a <= p:173            targets.append((4, p * p * (5 / 27) ** a))174        for K in (3, 4, 5):175            if K * (3**a - 1) < p and a < len(rows[K]):176                targets.append((2 * K, p * p * rows[K][a] / 3 ** (2 * K * a)))177        for power, target in targets:178            worst = max(worst, abs((F[a - 1] ** power).sum() - target) / target)179            checks += 1180print("identities checked %d, worst relative error %.1e" % (checks, worst))181182def half(p):183    a = 0184    while 2 * 3 ** (a + 1) <= p:185        a += 1186    return a187188KAPS = {K: float(kappas[K]) for K in (3, 4, 5)}189EDGE10 = KAPS[5] / (2 - KAPPA + 2 * KAPS[5])190191def master10(p, n):192    a = half(p)193    b = min(a - 1, n - 2 * a)194    if n < 2 * a or b < 0:195        return None196    return p * p * 3.0 ** (-((8 + KAPPA) * a + KAPS[5] * b) / 10)197198def master(p, n):199    a = half(p)200    bounds = [float(p * p)]201    if 3 ** (n // 2) <= p:202        bounds.append(p * p * 3.0 ** (-(n // 2)))203    if n >= 4 * a:204        bounds.append(p * p * 3.0 ** (-KAPPA * a))205    if n >= 3 * a:206        bounds.append(p * p * 3.0 ** (-(1 + KAPPA) * a / 2))207    if n >= 2 * a:208        b4 = min(a, n - 2 * a)209        bounds.append(p * p * 3.0 ** (-a / 2 - KAPPA * (a + b4) / 4))210        b = min(a - 1, n - 2 * a)211        if b >= 0:212            for K in (3, 4, 5):213                bounds.append(p * p * 3.0 ** (-((2 * K - 2 + KAPPA) * a + KAPS[K] * b) / (2 * K)))214                if n >= 2 * K * b:215                    bounds.append(p * p * 3.0 ** (-KAPS[K] * b))216    return min(bounds)217218print("master bound against exact L_n(p), primes 5..199, n = 2..24, the min over all orders and the order-10 block alone")219solo, solo_worst, solo_where, solo_cases = 0, 0.0, None, 0220violations, worst, where = 0, 0.0, None221for p in primes(5, 199):222    for n, F in enumerate(fourier(p, 24), start=1):223        if n < 2:224            continue225        only = master10(p, n)226        if only is not None:227            solo_cases += 1228            r10 = F.sum() / only229            if r10 > 1 + 1e-9:230                solo += 1231            if r10 > solo_worst:232                solo_worst, solo_where = r10, (p, n)233        ratio = F.sum() / master(p, n)234        if ratio > 1 + 1e-9:235            violations += 1236        if ratio > worst:237            worst, where = ratio, (p, n)238print("min over all orders: violations %d, worst ratio L/bound %.4f at (p, n) = %s" % (violations, worst, where))239print("order-10 block alone: %d cases, violations %d, worst ratio L/bound %.4f at (p, n) = %s" % (solo_cases, solo, solo_worst, solo_where))240241def regime(p, n):242    a = half(p)243    if 2 * a > n:244        return None245    return "I" if n >= 4 * a else "II" if n >= 3 * a else "III"246247def tail_bound(p, n):248    a = half(p)249    r = regime(p, n)250    if r == "I":251        return p * p * (5 / 27) ** a252    if r == "II":253        return p * p * (5 / 81) ** (a / 2)254    return p * p * 3.0 ** (-a / 2) * (5 / 27) ** ((n - a) / 4)255256print("dyadic tail: exponents (kappa - 1)/8 = %.4f, (2 + kappa)/4 = %.4f" % ((KAPPA - 1) / 8, (2 + KAPPA) / 4))257print("tail constants 16/(kappa - 1) = %.2f, 10 (1 + kappa)/(kappa - 1) 2^((1 - kappa)/2) = %.2f, 8/(2 + kappa) = %.2f" % (16 / (KAPPA - 1), 10 * (1 + KAPPA) / (KAPPA - 1) * 2 ** ((1 - KAPPA) / 2), 8 / (2 + KAPPA)))258uncovered = [(p, n) for n in (8, 10, 12, 14, 16, 20, 24) for p in primes(5, 199) if p <= 3 ** (EDGE10 * n) and regime(p, n) is None]259print("regime cover below the order-10 edge 3^(%.12f n) for n in 8..24: uncovered %d" % (EDGE10, len(uncovered)))260fails, cases = 0, 0261for n in (8, 10, 12):262    for p in primes(5, 199):263        if regime(p, n) is None:264            continue265        cases += 1266        R = fourier(p, n)[-1].sum() - 1267        if R > tail_bound(p, n) * (1 + 1e-9):268            fails += 1269print("per-prime tail bounds at n = 8, 10, 12: %d cases, %d failures" % (cases, fails))270271LAMBDA = {K: 2 - KAPPA + 2 * KAPS[K] for K in (3, 4, 5)}272print("master-bound constants per order: Lambda_2K = 2 - kappa + 2 kappa_2K, edge kappa_2K/Lambda_2K, decay Lambda_2K/2K, geometric constant 2/(1 - 3^(-Lambda_2K/2K))")273for K in (3, 4, 5):274    print("2K = %2d Lambda %.9f edge %.12f decay %.9f constant %.4f" % (2 * K, LAMBDA[K], KAPS[K] / LAMBDA[K], LAMBDA[K] / (2 * K), 2 / (1 - 3.0 ** (-LAMBDA[K] / (2 * K)))))275276def gain10(a, n):277    return (KAPS[5] * n - LAMBDA[5] * a) / 10278279def four_term(n, z, eta):280    return 2 / z + 35 * z ** (1 - KAPPA) + 40 * 3.0 ** (-(KAPPA - 1) * n / 8) + 6 * 3.0 ** (-(LAMBDA[5] / 10) * eta * n)281282edge10 = EDGE10283bad, ratio = 0, 0.0284for eta in (0.001, 0.01, 0.05, 0.1):285    for n in range(6, 401):286        lo, hi = n // 3 + 1, int((edge10 - eta) * n)287        if hi < lo:288            continue289        block = [gain10(a, n) for a in range(lo, hi + 1)]290        total = sum(2 * 3.0 ** (-g) for g in block)291        cap = 6 * 3.0 ** (-(LAMBDA[5] / 10) * eta * n)292        if min(block) <= 0 or total > cap:293            bad += 1294        ratio = max(ratio, total / cap)295print("order-10 main range, primes with 3a > n: gains positive and geometric cap holds at eta in 0.001..0.1, n = 6..400, failures %d, worst sum/cap %.4f" % (bad, ratio))296297def energy_sum(p, n):298    t = np.arange(p)299    e = np.exp(2j * np.pi * t / p)300    W = np.abs(1 + e[:, None] + e[None, :]) / 3301    F = np.ones((p, p))302    for l in range(n):303        I = (t * pow(3, l, p)) % p304        F = F * W[np.ix_(I, I)]305    return F.sum()306307print("four-term dyadic bound against the exact prime sum, eta = 0.02")308for n in (6, 8, 10, 12):309    for z in (5, 11):310        top = 3.0 ** ((edge10 - 0.02) * n)311        total = sum(energy_sum(p, n) / (p * p) for p in primes(z + 1, int(top)))312        print("n = %2d z = %2d exact %.6f bound %.6f holds %s" % (n, z, total, four_term(n, z, 0.02), total <= four_term(n, z, 0.02)))