staircase.py
18.3 kB · python · 444 lines
1import math2import time3from fractions import Fraction45import mpmath as mp6import numpy as np7import sympy as sp89mp.mp.dps = 5010HALF = mp.mpf(3) / 21112def nine(x):13 x = mp.re(x)14 return mp.nstr(x, 9 + max(0, int(mp.floor(mp.log10(abs(x)))) + 1), strip_zeros=False)1516# CORNERS1718def corners(code, dim):19 return [i for i in range(2 ** dim) if code >> i & 1]2021def profile(code, dim):22 a = [0] * (dim + 1)23 for i in corners(code, dim):24 a[bin(i).count("1")] += 125 return tuple(a)2627def profiles(dim):28 seen = {}29 for code in range(1, 2 ** (2 ** dim)):30 seen.setdefault(profile(code, dim), []).append(code)31 return seen3233def coefficients(a):34 dim = len(a) - 135 c = [0] * (dim + 1)36 for j, aj in enumerate(a):37 for k in range(j + 1):38 c[k] += aj * math.comb(j, k) * (-1) ** k39 return c4041def fill(a, n):42 dim = len(a) - 143 return sum(aj * n ** (dim - j) * (n - 1) ** j for j, aj in enumerate(a))4445def shape(a):46 dim = len(a) - 147 w = sum(a)48 mean = Fraction(sum(j * aj for j, aj in enumerate(a)), w)49 var = Fraction(sum(j * j * aj for j, aj in enumerate(a)), w) - mean ** 250 drift = Fraction(dim, 2) - mean51 return w, drift, Fraction(dim, 8) + drift - var / 25253# RENDER5455def tile(code, dim, side):56 axes = np.indices((side,) * dim)57 index = sum((axes[j] % 2) << (dim - 1 - j) for j in range(dim))58 return ((code >> index) & 1).astype(np.int64)5960def section_render():61 bad = 062 checked = 063 for dim in (1, 2, 3):64 for code in range(1, 2 ** (2 ** dim)):65 a = profile(code, dim)66 w = sum(a)67 for side in (3, 5, 7, 9):68 bad += int(tile(code, dim, side).sum()) != fill(a, (side + 1) // 2)69 checked += 170 for side in (2, 4, 6, 8):71 bad += int(tile(code, dim, side).sum()) * 2 ** dim != w * side ** dim72 checked += 173 levels = {1: 5, 2: 3, 3: 2}[dim]74 word = tile(code, dim, 3)75 for level in range(2, levels + 1):76 word = np.kron(word, tile(code, dim, 2 * level + 1))77 want = math.prod(fill(a, k + 1) for k in range(1, level + 1))78 bad += int(word.sum()) != want79 checked += 180 print(f"render: {checked} counts over every code in dim 1, 2, 3 (odd sides 3..9, even sides 2..8, row words), mismatches {bad}")81 named = {(2, 7): (1, 2, 0), (2, 14): (0, 2, 1), (2, 6): (0, 2, 0), (2, 9): (1, 0, 1), (3, 23): (1, 3, 0, 0), (3, 232): (0, 0, 3, 1), (3, 150): (0, 3, 0, 1), (3, 105): (1, 0, 3, 0)}82 wrong = [k for k, v in named.items() if profile(k[1], k[0]) != v]83 print(f"render: crate corner order (corner i = binary digits of i, most significant first), named profiles wrong {len(wrong)}")8485# ROOTS8687def roots(a):88 n = sp.Symbol("n")89 p = sp.Poly(coefficients(a), n)90 out = []91 for factor, mult in sp.factor_list(p.as_expr())[1]:92 f = sp.Poly(factor, n)93 if f.degree() == 1:94 r = -f.all_coeffs()[1] / f.all_coeffs()[0]95 out += [mp.mpf(sp.Rational(r).p) / sp.Rational(r).q] * mult96 else:97 cs = [int(x) for x in f.all_coeffs()]98 out += list(mp.polyroots(cs, maxsteps=400, extraprec=400)) * mult99 return out100101def constant(rs, dim):102 value = HALF_GAMMA ** dim / mp.fprod(mp.gamma(2 - r) for r in rs)103 assert abs(mp.im(value)) < mp.mpf(10) ** -40 and mp.re(value) > 0104 return mp.re(value)105106HALF_GAMMA = mp.gamma(HALF)107108def log_exact_formula(a, rs, level):109 dim = len(a) - 1110 w = sum(a)111 s = level * mp.log(mp.mpf(w) / 2 ** dim)112 for r in rs:113 s += mp.loggamma(level + 2 - r) - mp.loggamma(2 - r)114 s -= dim * (mp.loggamma(level + HALF) - mp.loggamma(HALF))115 return mp.re(s)116117# EXACT118119def section_exact(table):120 worst = mp.mpf(0)121 for a, rs in table:122 dim = len(a) - 1123 top, bottom = 1, 1124 marks = {1, 2, 3, 10, 100, 1000}125 for k in range(1, 1001):126 top *= fill(a, k + 1)127 bottom *= (2 * k + 1) ** dim128 if k in marks:129 gap = abs(mp.log(top) - mp.log(bottom) - log_exact_formula(a, rs, k))130 worst = max(worst, gap)131 print(f"exact: R_L = (w/2^dim)^L prod Gamma(L+2-r)/Gamma(2-r) / (Gamma(L+3/2)/Gamma(3/2))^dim on {len(table)} profiles at L = 1, 2, 3, 10, 100, 1000, worst log gap {mp.nstr(worst, 3)}")132133# ASYMPTOTIC134135def direct_logs(a, top):136 dim = len(a) - 1137 w = sum(a)138 big = sp.Symbol("N")139 poly = sp.Poly(sum(aj * (big + 1) ** (dim - j) * (big - 1) ** j for j, aj in enumerate(a)) - w * big ** dim, big)140 d = [float(x) for x in poly.all_coeffs()] if not poly.is_zero else [0.0]141 side = 2.0 * np.arange(1, top + 1) + 1.0142 return np.log1p(np.polyval(d, side) / (w * side ** dim))143144def section_asymptotic(table, top=10 ** 6):145 marks = (10 ** 4, 10 ** 5, top)146 worst = {m: 0.0 for m in marks}147 for a, rs in table:148 dim = len(a) - 1149 w, drift, c1 = shape(a)150 c1_roots = mp.re(sum((mp.mpf(1) / 2 - r) * (mp.mpf(5) / 2 - r) for r in rs) / 2)151 assert abs(c1_roots - mp.mpf(c1.numerator) / c1.denominator) < mp.mpf(10) ** -30152 logs = direct_logs(a, top)153 log_c = float(mp.log(constant(rs, dim)))154 for m in marks:155 got = float(np.sum(logs[:m]))156 want = float(drift) * math.log(m) + log_c + float(c1) / m157 worst[m] = max(worst[m], abs(got - want))158 print(f"asymptotic: direct log sum against drift log L + log C + c1/L, c1 = dim/8 + drift - var/2, on {len(table)} profiles")159 for m in marks:160 print(f" L = {m:>8}: worst residual {worst[m]:.2e}, times L^2 {worst[m] * m * m:.2f}")161162# CLOSED FORMS163164def section_named(byprofile):165 pi = mp.pi166 sq = mp.sqrt167 forms = [168 (1, 1, "sqrt(pi)/2", sq(pi) / 2),169 (1, 2, "sqrt(pi)/2", sq(pi) / 2),170 (2, 7, "3 pi/(4 Gamma(1/3))", 3 * pi / (4 * mp.gamma(mp.mpf(1) / 3))),171 (2, 14, "3 pi/(8 Gamma(2/3))", 3 * pi / (8 * mp.gamma(mp.mpf(2) / 3))),172 (2, 6, "pi/4", pi / 4),173 (2, 9, "cosh(pi/2)/2", mp.cosh(pi / 2) / 2),174 (2, 11, "3 cosh(pi/(2 sqrt 3))/4", 3 * mp.cosh(pi / (2 * sq(3))) / 4),175 (2, 1, "pi/4", pi / 4),176 (2, 3, "sqrt(pi)/2", sq(pi) / 2),177 (3, 23, "pi^(3/2)/(2 Gamma(1/4))", pi ** HALF / (2 * mp.gamma(mp.mpf(1) / 4))),178 (3, 232, "pi^(3/2)/(6 Gamma(3/4))", pi ** HALF / (6 * mp.gamma(mp.mpf(3) / 4))),179 (3, 150, "pi^(3/2)/(8 |Gamma(7/4 - i sqrt3/4)|^2)", pi ** HALF / (8 * abs(mp.gamma(mp.mpc(7, -sq(3)) / 4)) ** 2)),180 (3, 105, "pi^(3/2)/(8 |Gamma(5/4 - i sqrt3/4)|^2)", pi ** HALF / (8 * abs(mp.gamma(mp.mpc(5, -sq(3)) / 4)) ** 2)),181 (3, 129, "cosh(pi sqrt3/2)/4", mp.cosh(pi * sq(3) / 2) / 4),182 (3, 126, "pi/4", pi / 4),183 (3, 24, "pi/4", pi / 4),184 ]185 print("named: dim, code, profile, w/2^dim, drift, C (Gamma at roots), closed form, gap")186 worst = mp.mpf(0)187 cs = {}188 for dim, code, text, value in forms:189 a = profile(code, dim)190 rs = byprofile[a]191 w, drift, _ = shape(a)192 c = constant(rs, dim)193 cs[(dim, code)] = c194 worst = max(worst, abs(c - value))195 print(f" dim {dim} code {code:>3} {str(a):<14} {str(Fraction(w, 2 ** dim)):>4} {str(drift):>5} {mp.nstr(c, 12):<15} {nine(c):<12} {text} {mp.nstr(abs(c - value), 3)}")196 print(f"named: worst closed-form gap {mp.nstr(worst, 3)}")197 mirrors = [((2, 7), (2, 14), "9 sqrt3 pi/64", 9 * sq(3) * pi / 64), ((3, 23), (3, 232), "sqrt2 pi^2/24", sq(2) * pi ** 2 / 24), ((1, 1), (1, 2), "pi/4", pi / 4)]198 for f, g, text, value in mirrors:199 print(f" mirror dim {f[0]} codes {f[1]} x {g[1]}: {mp.nstr(cs[f] * cs[g], 15)} against {text} gap {mp.nstr(abs(cs[f] * cs[g] - value), 3)}")200201# PAIRING202203def reflect(rs):204 out = mp.mpf(1)205 for r in rs:206 if abs(r) < mp.mpf(10) ** -40 or abs(r - 1) < mp.mpf(10) ** -40:207 out *= mp.pi / 4208 else:209 out *= mp.sin(mp.pi * r) / (4 * r * (1 - r))210 return mp.re(out)211212def section_pairing(byprofile):213 n = sp.Symbol("n")214 agree = 0215 total = 0216 for dim in (1, 2, 3, 4):217 for a in profiles(dim) if dim < 4 else four_profiles():218 p = sp.Poly(coefficients(a), n)219 q = sp.Poly(p.as_expr().subs(n, 1 - n) * (-1) ** dim, n)220 flipped = p == q221 palindrome = tuple(a) == tuple(reversed(a))222 agree += flipped == palindrome223 total += 1224 print(f"pairing: P(1-n) = (-1)^dim P(n) exactly when the profile is a palindrome, {agree} of {total} profiles at dim 1..4")225 worst = mp.mpf(0)226 for a, rs in byprofile.items():227 dim = len(a) - 1228 mirror = tuple(reversed(a))229 worst = max(worst, abs(constant(rs, dim) * constant(byprofile[mirror], dim) - reflect(rs)))230 print(f"pairing: C_F C_mirror = prod sin(pi r)/(4 r (1-r)) on all {len(byprofile)} profiles at dim 1..3, worst gap {mp.nstr(worst, 3)}")231 flips = 0232 palins = 0233 for dim in (2, 3):234 for code in range(1, 2 ** (2 ** dim)):235 cs = set(corners(code, dim))236 closed = all((2 ** dim - 1 - i) in cs for i in cs)237 a = profile(code, dim)238 palin = a == tuple(reversed(a))239 flips += closed240 palins += palin241 assert palin or not closed242 print(f"pairing: codes closed under the flip of every coordinate {flips}, codes with a palindromic profile {palins}, at dim 2 and 3 together")243 a = profile(11, 2)244 print(f"pairing: dim 2 code 11 corners {corners(11, 2)} profile {a}, not flip-closed, roots {[mp.nstr(r, 12) for r in byprofile[a]]}")245 tiny = mp.mpf(10) ** -30246 def core(rs):247 return [r for r in rs if min(abs(r), abs(r - mp.mpf(1) / 2), abs(r - 1)) > tiny]248 def paired(rs, shifts):249 left = list(rs)250 while left:251 r = left.pop()252 hit = next((i for i, s in enumerate(left) if any(abs(s - (1 - r + k)) < tiny for k in shifts)), None)253 if hit is None:254 return False255 left.pop(hit)256 return True257 palin = [a for a in byprofile if a == tuple(reversed(a))]258 exact = [a for a, rs in byprofile.items() if paired(core(rs), (0,))]259 shifted = [a for a, rs in byprofile.items() if paired(core(rs), range(-4, 5))]260 extra = [a for a in exact if a not in palin]261 print(f"pairing: of {len(byprofile)} profiles at dim 1..3, palindromic {len(palin)}, core roots (outside 0, 1/2, 1) closed under r -> 1-r {len(exact)}, closed up to an integer shift {len(shifted)}")262 print(f"pairing: reducible without a palindrome {len(extra)}: " + ", ".join(f"{a} code {byprofile_code(a)}" for a in extra))263 rest = [a for a in byprofile if a not in shifted]264 print(f"pairing: unreduced {len(rest)}, rational core roots among them {sum(all(abs(mp.im(r)) < tiny and abs(mp.re(r) * 840 - round(mp.re(r) * 840)) < tiny for r in core(byprofile[a])) for a in rest)}")265266def byprofile_code(a):267 dim = len(a) - 1268 return next(c for c in range(1, 2 ** (2 ** dim)) if profile(c, dim) == a)269270def four_profiles():271 out = []272 for a0 in range(2):273 for a1 in range(5):274 for a2 in range(7):275 for a3 in range(5):276 for a4 in range(2):277 if a0 + a1 + a2 + a3 + a4:278 out.append((a0, a1, a2, a3, a4))279 return out280281# PARITY282283def parity_constant(dim, sign):284 us = [mp.exp(2j * mp.pi * (k + (0 if sign < 0 else mp.mpf(1) / 2)) / dim) for k in range(dim)]285 return mp.re(HALF_GAMMA ** dim / mp.fprod(mp.gamma((3 - u) / 2) for u in us))286287def zeta_series(dim, sign):288 total = mp.mpf(0)289 for m in range(1, 400):290 tail = (1 - mp.mpf(2) ** (-dim * m)) * mp.zeta(dim * m) - 1291 term = tail / m * (-1 if sign < 0 else (-1) ** (m + 1))292 total += term293 if abs(tail) < mp.mpf(10) ** -60:294 break295 return mp.exp(total)296297def section_parity():298 pi = mp.pi299 forms = {300 (2, -1): ("pi/4", pi / 4),301 (2, 1): ("cosh(pi/2)/2", mp.cosh(pi / 2) / 2),302 (4, -1): ("pi cosh(pi/2)/8", pi * mp.cosh(pi / 2) / 8),303 (4, 1): ("(cosh(pi/sqrt2) + cos(pi/sqrt2))/4", (mp.cosh(pi / mp.sqrt(2)) + mp.cos(pi / mp.sqrt(2))) / 4),304 (6, -1): ("pi cosh(pi sqrt3/2)/24", pi * mp.cosh(pi * mp.sqrt(3) / 2) / 24),305 (6, 1): ("cosh(pi/2)(cosh(pi/2) + cos(pi sqrt3/2))/4", mp.cosh(pi / 2) * (mp.cosh(pi / 2) + mp.cos(pi * mp.sqrt(3) / 2)) / 4),306 }307 print("parity: dim, parity, code, Gamma at roots of unity, gap to the zeta log series exp(-+sum_m (+-1)^m (lambda(dim m) - 1)/m), direct log sum to L = 10^6 plus c1/L, closed form")308 for dim in (2, 3, 4, 5, 6):309 for sign, name in ((-1, "odd"), (1, "even")):310 want = 1 if name == "odd" else 0311 code = sum(1 << i for i in range(2 ** dim) if bin(i).count("1") % 2 == want)312 a = profile(code, dim)313 g = parity_constant(dim, sign)314 _, drift, c1 = shape(a)315 assert drift == 0316 logs = np.log1p(sign * (2.0 * np.arange(1, 10 ** 6 + 1) + 1.0) ** (-dim))317 direct = math.exp(float(np.sum(logs)) - float(c1) / 10 ** 6)318 series = zeta_series(dim, sign)319 text, value = forms.get((dim, sign), ("", None))320 gap = "" if value is None else f"closed-form gap {mp.nstr(abs(g - value), 3)}"321 print(f" dim {dim} {name:<4} code {code:<6} {nine(g):<12} series gap {mp.nstr(abs(g - series), 3):<9} direct {direct:.12f} {text} {gap}")322323# FAMILY324325def section_family():326 print("family: at most one odd coordinate, dim, code, w/2^dim, drift, C at the roots of the fill polynomial, gap to (sqrt(pi)/2)^dim (dim+1)/Gamma(1/(dim+1)), direct at L = 10^6")327 for dim in range(1, 9):328 a = tuple([1, dim] + [0] * (dim - 1))329 code = sum(1 << i for i in range(2 ** dim) if bin(i).count("1") <= 1)330 rs = roots(a)331 w, drift, c1 = shape(a)332 c = constant(rs, dim)333 form = (mp.sqrt(mp.pi) / 2) ** dim * (dim + 1) / mp.gamma(mp.mpf(1) / (dim + 1))334 logs = direct_logs(a, 10 ** 6)335 direct = math.exp(float(np.sum(logs)) - float(drift) * math.log(10 ** 6) - float(c1) / 10 ** 6)336 print(f" dim {dim} code {code:<20} {str(Fraction(w, 2 ** dim)):>6} {str(drift):>7} {mp.nstr(c, 12):<15} {nine(c):<12} gap {mp.nstr(abs(c - form), 3):<8} direct {direct:.10f}")337338# FENCE339340def section_fence(byprofile):341 print("fence: renormalised limit prod_k P((N_k+1)/2)/(w (N_k/2)^dim) along N_k = 3^k, k = 1..60")342 for dim, code in ((1, 1), (2, 7), (3, 23), (2, 6), (2, 9)):343 a = profile(code, dim)344 w = sum(a)345 out = mp.mpf(1)346 for k in range(1, 61):347 side = 3 ** k348 out *= mp.mpf(fill(a, (side + 1) // 2) * 2 ** dim) / (w * mp.mpf(side) ** dim)349 print(f" dim {dim} code {code}: {mp.nstr(out, 12)}, row-word constant {mp.nstr(constant(byprofile[a], dim), 12)}")350 a = profile(7, 2)351 for top in (10 ** 3, 10 ** 5):352 primes = [p for p in sp.primerange(3, top)]353 s = mp.mpf(1)354 for p in primes:355 s *= mp.mpf(fill(a, (p + 1) // 2) * 4) / (3 * mp.mpf(p) ** 2)356 print(f" dim 2 code 7 along the odd primes below {top}: renormalised product {mp.nstr(s, 8)}, sum 1/p {mp.nstr(sum(mp.mpf(1) / p for p in primes), 6)}")357358# STAIRCASE359360def stair_exact(a, rs, n):361 dim = len(a) - 1362 w = sum(a)363 g = mp.barnesg364 top = n * (n + 1) / 2 * mp.log(mp.mpf(w) / 2 ** dim)365 for r in rs:366 top += mp.log(g(n + 3 - r)) - mp.log(g(3 - r)) - n * mp.loggamma(2 - r)367 top -= dim * (mp.log(g(n + mp.mpf(5) / 2)) - mp.log(g(mp.mpf(5) / 2)) - n * mp.loggamma(HALF))368 side = n * (n + 1) / 2 * mp.log(2) + mp.log(g(n + mp.mpf(5) / 2)) - mp.log(g(mp.mpf(5) / 2)) - n * mp.loggamma(HALF)369 return dim + mp.re(top) / side370371def stair_direct(a, n):372 dim = len(a) - 1373 top = sum((n - j + 1) * mp.log(fill(a, j + 1)) for j in range(1, n + 1))374 side = sum((n - j + 1) * mp.log(2 * j + 1) for j in range(1, n + 1))375 return top / side376377def section_staircase(byprofile):378 a = profile(7, 2)379 print("staircase: dim 2 code 7, dimension(n) by Barnes G at the roots and by the direct fill sum")380 print(" " + ", ".join(f"{mp.nstr(stair_exact(a, byprofile[a], n), 10)}" for n in range(1, 6)))381 worst = mp.mpf(0)382 for b, rs in byprofile.items():383 for n in (1, 2, 3, 6, 20):384 worst = max(worst, abs(stair_exact(b, rs, n) - stair_direct(b, n)))385 print(f"staircase: Barnes G form against the direct sum on all {len(byprofile)} profiles at dim 1..3, n = 1, 2, 3, 6, 20, worst gap {mp.nstr(worst, 3)}")386 top = 2 * 10 ** 5387 marks = (10 ** 2, 10 ** 3, 10 ** 4, 10 ** 5, top)388 print("staircase: gap = dim - dimension, against log(2^dim/w)/(log 2n - 3/2) (staircase) and /(log 2L - 1) (row word); n x error")389 for dim, code in ((1, 1), (2, 7), (2, 9), (3, 23)):390 b = profile(code, dim)391 w = sum(b)392 lead = math.log(2 ** dim / w)393 logs = direct_logs(b, top)394 row_fill = np.cumsum(logs) + np.arange(1, top + 1) * math.log(w / 2 ** dim)395 row_side = np.cumsum(np.log(2.0 * np.arange(1, top + 1) + 1.0))396 stair_fill = np.cumsum(row_fill)397 stair_side = np.cumsum(row_side)398 cells = []399 for m in marks:400 stair_gap = -stair_fill[m - 1] / stair_side[m - 1]401 row_gap = -row_fill[m - 1] / row_side[m - 1]402 stair_err = m * (stair_gap - lead / (math.log(2 * m) - 1.5))403 row_err = m * (row_gap - lead / (math.log(2 * m) - 1))404 cells.append(f"n={m}: {stair_gap:.6f} ({stair_err:+.2f}), row {row_gap:.6f} ({row_err:+.2f})")405 print(f" dim {dim} code {code}:")406 for cell in cells:407 print(f" {cell}")408409# MAIN410411def main():412 clock = time.time()413 def lap(name):414 nonlocal clock415 print(f"[{name} {time.time() - clock:.1f}s]")416 clock = time.time()417 section_render()418 lap("render")419 byprofile = {}420 for dim in (1, 2, 3):421 for a in profiles(dim):422 byprofile[a] = roots(a)423 table = list(byprofile.items())424 print(f"roots: {len(table)} profiles (3 at dim 1, 11 at dim 2, 63 at dim 3) cover all {3 + 15 + 255} nonempty codes")425 lap("roots")426 section_exact(table)427 lap("exact")428 section_asymptotic(table)429 lap("asymptotic")430 section_named(byprofile)431 lap("named")432 section_pairing(byprofile)433 lap("pairing")434 section_parity()435 lap("parity")436 section_family()437 lap("family")438 section_fence(byprofile)439 lap("fence")440 section_staircase(byprofile)441 lap("staircase")442443if __name__ == "__main__":444 main()