determinant.py
20.1 kB · python · 503 lines
1import argparse2import importlib.util3import math4import os5import time6from collections import Counter7from fractions import Fraction89import numpy as np10from mpmath import iv, mp1112HERE = os.path.dirname(os.path.abspath(__file__))13JP_LO = "0.53128050627720514162446864736847178549305910901839877988839780392752953564383134591810957018118523987"14JP_HI = "0.53128050627720514162446864736847178549305910901839877988839780392752953564383134591810957018118523989"15PV = {(1, 3): "0.454489077661828743845", (2, 3): "0.337436780806063636304", (1, 2, 3): "0.705660908028738230607"}16S1_NOTE = ("0.457015235231", "6.958882679527")17WALKS = 4 * 10 ** 6181920def sibling():21 path = os.path.join(HERE, "..", "question-mark", "question.py")22 spec = importlib.util.spec_from_file_location("question", path)23 mod = importlib.util.module_from_spec(spec)24 spec.loader.exec_module(mod)25 return mod262728def ivq(q):29 return iv.mpf(q.numerator) / iv.mpf(q.denominator)303132def up(x):33 return mp.mpf(x._mpi_[1])343536def low(x):37 return mp.mpf(x._mpi_[0])383940def exact(x):41 man, ex = x.man, x.exp42 return Fraction(int(man)) * (Fraction(2) ** int(ex))434445def fixed(q, d):46 sign = "-" if q < 0 else ""47 q = abs(q)48 whole = q.numerator // q.denominator * 10 ** d + (q - q.numerator // q.denominator) * 10 ** d49 n = int(whole)50 text = str(n).rjust(d + 1, "0")51 return f"{sign}{text[:-d]}.{text[-d:]}"525354def bracket(lo, hi, d=None):55 lo, hi = exact(lo), exact(hi)56 if d is None:57 d = 158 while Fraction(1, 10 ** (d + 1)) > (hi - lo) / 10 or d < 3:59 d += 160 down = Fraction(math.floor(lo * 10 ** d), 10 ** d)61 upq = Fraction(math.ceil(hi * 10 ** d), 10 ** d)62 return f"[{fixed(down, d)}, {fixed(upq, d)}]"636465def ceil_sig(x, k):66 q = x if isinstance(x, Fraction) else exact(mp.mpf(x))67 assert q > 068 e = int(mp.floor(mp.log10(mp.mpf(q.numerator) / q.denominator)))69 while q >= Fraction(10) ** (e + 1):70 e += 171 while q < Fraction(10) ** e:72 e -= 173 n = math.ceil(q * Fraction(10) ** (k - 1 - e))74 if n >= 10 ** k:75 n, e = math.ceil(Fraction(n, 10)), e + 176 t = str(n)77 if -4 <= e <= 5:78 d = max(0, k - 1 - e)79 return fixed(Fraction(n) * Fraction(10) ** (e - k + 1), d)80 return f"{t[0]}.{t[1:]}e{e}" if k > 1 else f"{t}e{e}"818283def trunc(x, digits):84 return math.floor(x * 10 ** digits) / 10 ** digits858687# GRAPH8889def closed_walks(edges, period):90 states = sorted({u for u, _, _ in edges} | {v for _, v, _ in edges})91 count = Counter()92 for u0 in states:93 st = np.array([u0])94 mat = np.array([[1, 0, 0, 1]], dtype=np.int64)95 for n in range(1, period + 1):96 ns, nm = [], []97 for u, v, a in edges:98 pick = st == u99 if pick.any():100 x = mat[pick]101 nm.append(np.stack([x[:, 1], x[:, 0] + a * x[:, 1], x[:, 3], x[:, 2] + a * x[:, 3]], axis=1))102 ns.append(np.full(len(x), v))103 st, mat = np.concatenate(ns), np.concatenate(nm)104 assert len(st) < WALKS and mat.max() < 2 ** 60105 shut = st == u0106 t, c = np.unique(mat[shut, 0] + mat[shut, 3], return_counts=True)107 for x, y in zip(t, c):108 count[(n, int(x))] += int(y)109 return sorted(count.items())110111112def minimize(states, edges):113 out = {u: {} for u in states}114 for u, v, a in edges:115 out[u][a] = v116 part = {u: tuple(sorted(out[u])) for u in states}117 while True:118 new = {u: (part[u], tuple((a, part[out[u][a]]) for a in sorted(out[u]))) for u in states}119 if len(set(new.values())) == len(set(part.values())):120 break121 part = new122 idx = {b: i for i, b in enumerate(sorted(set(part.values()), key=str))}123 return len(idx), sorted({(idx[part[u]], idx[part[out[u][a]]], a) for u in states for a in out[u]})124125126def period_of(n, edges):127 level = {0: 0}128 todo = [0]129 while todo:130 u = todo.pop()131 for x, y, _ in edges:132 if x == u and y not in level:133 level[y] = level[u] + 1134 todo.append(y)135 assert len(level) == n136 g = 0137 for x, y, _ in edges:138 g = math.gcd(g, level[x] + 1 - level[y])139 return g140141142def pieces(q, code, k):143 n = 1 << (k - 1)144 edges, tails = q.graph_of(code, k)145 comp = q.components(n, edges, tails)146 out = []147 for c in sorted(set(comp.values()), key=sorted):148 e = [(u, v, a) for (u, v), ls in edges.items() for a in ls if u in c and v in c]149 cof = any(u in c and v in c for u, v in tails.items())150 if cof:151 out.append(("tail", len(c), None))152 elif len(e) > len(c):153 s, me = minimize(sorted(c), e)154 out.append(("finite", len(c), (s, me)))155 return out156157158# DISC159160def disc_numbers(edges, c, r, beta):161 labels = sorted({a for _, _, a in edges})162 ends = [abs(1 / (c - r + a) - c) for a in labels] + [abs(c - 1 / (c + r + a)) for a in labels]163 return labels, max(ends) / r164165166def choose_disc(edges, s, states, period):167 th = np.linspace(0, 2 * np.pi, 721)168 best = None169 for c in np.arange(0.30, 0.92, 0.04):170 for r in np.arange(0.30, 1.30, 0.05):171 if c - r <= -0.9:172 continue173 labels, h = disc_numbers(edges, c, r, None)174 if h >= 0.95:175 continue176 z = c + r * np.exp(1j * th)177 for beta in np.arange(0.6, 4.01, 0.2):178 if c - r <= -beta + 0.05:179 continue180 w = {a: np.exp((2 * s * (np.log(z + beta) - np.log(beta * (z + a) + 1))).real).max() for a in labels}181 sig = sum(sum(w[a] for u, v, a in edges if u == x and v == y) ** 2 for x in range(states) for y in range(states))182 cc = math.sqrt(sig / (1 - h * h))183 m = np.arange(1, 60)184 e = np.concatenate([[1.0], np.cumprod(cc * h ** (m - 1) / (1 - h ** m))])185 f = e186 for _ in range(states - 1):187 f = np.convolve(f, e)[:60]188 tail = f[period + 1:].sum()189 if best is None or tail < best[0]:190 best = (tail, round(c, 2), round(r, 2), round(beta, 2))191 return tuple(Fraction(str(x)) for x in best[1:])192193194def weight_sup(labels, c, r, beta, sig, tau, pad, arcs):195 out = {}196 for a in labels:197 top = None198 for j in range(arcs):199 th = iv.mpf([2 * j, 2 * j + 2]) * iv.pi / arcs200 x = ivq(c) + ivq(r) * iv.cos(th)201 y = ivq(r) * iv.sin(th)202 if beta is None:203 u = x + a204 assert low(u) > 0205 re = -iv.log(u * u + y * y) / 2206 im = -iv.atan2(y, u)207 else:208 u1 = x + ivq(beta)209 u2 = ivq(beta) * x + ivq(beta * a + 1)210 v2 = ivq(beta) * y211 assert beta > 0 and low(u1) > 0 and low(u2) > 0212 re = (iv.log(u1 * u1 + y * y) - iv.log(u2 * u2 + v2 * v2)) / 2213 im = iv.atan2(y, u1) - iv.atan2(v2, u2)214 val = (2 * (sig * re - tau * im) + 2 * pad * iv.sqrt(re * re + im * im)).b215 top = val if top is None or val > top else top216 out[a] = iv.exp(iv.mpf(top))217 return out218219220def euler_tail(cc, h, states, period, radius):221 x = cc * radius222 e = [iv.mpf(1)]223 m = 0224 while True:225 m += 1226 e.append(e[-1] * x * h ** (m - 1) / (1 - h ** m))227 if m > period + 2 and (x * h ** m / (1 - h ** (m + 1))).b < 0.5:228 break229 rest = e[-1]230 f = e231 for _ in range(states - 1):232 f = [sum((f[i] * e[j - i] for i in range(max(0, j - len(e) + 1), min(j, len(f) - 1) + 1)), iv.mpf(0)) for j in range(len(f) + len(e) - 1)]233 total = sum(e, iv.mpf(0)) + rest234 return up(sum(f[period + 1:], iv.mpf(0)) + states * rest * total ** (states - 1))235236237def walk_cap(edges, states, period):238 adj = np.zeros((states, states))239 for u, v, _ in edges:240 adj[u, v] += 1241 growth = max(abs(np.linalg.eigvals(adj)))242 return min(period, int(math.log(WALKS / 2) / math.log(growth)))243244245class Certificate:246 def __init__(self, edges, states, period, point):247 period = walk_cap(edges, states, period)248 self.edges, self.states, self.period = edges, states, period249 self.labels = sorted({a for _, _, a in edges})250 self.c, self.r, self.beta = choose_disc(edges, point, states, period)251 assert self.c - self.r > -min(self.labels) and self.c - self.r > -self.beta252 self.hq = max(max(abs(1 / (self.c - self.r + a) - self.c), abs(self.c - 1 / (self.c + self.r + a))) for a in self.labels) / self.r253 assert self.hq < 1254 self.h = ivq(self.hq)255 self.groups = closed_walks(edges, period)256 self.iv_rows = []257 self.mp_rows = []258 for (n, t), m in self.groups:259 e = (-1) ** n260 mu = (t + iv.sqrt(iv.mpf(t * t - 4 * e))) / 2261 self.iv_rows.append((n, iv.log(mu), m / (1 - e / (mu * mu))))262 mum = (t + mp.sqrt(t * t - 4 * e)) / 2263 self.mp_rows.append((n, mp.log(mum), m / (1 - e / (mum * mum))))264265 def constant(self, sig, tau, pad, arcs):266 w = weight_sup(self.labels, self.c, self.r, self.beta, sig, tau, pad, arcs)267 sig2 = iv.mpf(0)268 for x in range(self.states):269 for y in range(self.states):270 sig2 += sum((w[a] for u, v, a in self.edges if u == x and v == y), iv.mpf(0)) ** 2271 return iv.sqrt(sig2 / (1 - self.h * self.h))272273 def tail(self, sig, tau, pad, radius=1, arcs=720):274 return euler_tail(self.constant(sig, tau, pad, arcs), self.h, self.states, self.period, radius)275276 def coefficients(self, s, ctx, deriv=False):277 rows = self.iv_rows if ctx is iv else self.mp_rows278 zero = ctx.mpc(0) if isinstance(s, (ctx.mpc,)) else ctx.mpf(0)279 tr = [zero] * (self.period + 1)280 dtr = [zero] * (self.period + 1)281 for n, lg, f in rows:282 x = f * ctx.exp(-2 * s * lg)283 tr[n] += x284 if deriv:285 dtr[n] += -2 * lg * x286 d, dd = [zero + 1], [zero]287 for n in range(1, self.period + 1):288 d.append(-sum((tr[k] * d[n - k] for k in range(1, n + 1)), zero) / n)289 if deriv:290 dd.append(-sum((dtr[k] * d[n - k] + tr[k] * dd[n - k] for k in range(1, n + 1)), zero) / n)291 return d, dd292293 def newton(self, s, steps=8):294 for _ in range(steps):295 d, dd = self.coefficients(s, mp, True)296 s = s - sum(d) / sum(dd)297 return s298299 def winding(self, d, radius, bound, arcs=360):300 def horner(z):301 v = d[-1]302 for coef in reversed(d[:-1]):303 v = v * z + coef304 return v305 planes = []306 pts = []307 for j in range(arcs):308 th = iv.mpf([2 * j, 2 * j + 2]) * iv.pi / arcs309 v = horner(radius * iv.mpc(iv.cos(th), iv.sin(th)))310 tests = [(0, 1, low(v.real) - bound > 0), (0, -1, up(v.real) + bound < 0), (1, 1, low(v.imag) - bound > 0), (1, -1, up(v.imag) + bound < 0)]311 plane = next((a, b) for a, b, ok in tests if ok)312 planes.append(plane)313 t = iv.mpf(2 * j) * iv.pi / arcs314 u = horner(radius * iv.mpc(iv.cos(t), iv.sin(t)))315 pts.append(complex(float(u.real.mid), float(u.imag.mid)))316 turn = 0.0317 for j in range(arcs):318 a, b = planes[j]319 for p in (pts[j], pts[(j + 1) % arcs]):320 assert b * (p.real if a == 0 else p.imag) > 0321 q = pts[(j + 1) % arcs] / pts[j]322 turn += math.atan2(q.imag, q.real)323 count = round(turn / (2 * math.pi))324 assert abs(turn / (2 * math.pi) - count) < 0.01325 return count326327328def top_zero(cert):329 f = lambda x: sum(cert.coefficients(mp.mpf(x), mp)[0])330 hi = 1.0331 fh = f(hi)332 assert fh > 0333 while True:334 lo = hi - 0.01335 assert lo > 0336 fl = f(lo)337 if fl < 0:338 break339 hi, fh = lo, fl340 for _ in range(30):341 mid = (lo + hi) / 2342 if f(mid) < 0:343 lo = mid344 else:345 hi = mid346 return cert.newton(mp.mpf((lo + hi) / 2), 4)347348349def certify_real(cert, guess, width_floor):350 s = cert.newton(mp.mpf(guess), 4) if guess is not None else top_zero(cert)351 cc = cert.constant(iv.mpf(s), iv.mpf(0), iv.mpf("1e-6"), 720)352 t1 = euler_tail(cc, cert.h, cert.states, cert.period, 1)353 t2 = euler_tail(cc, cert.h, cert.states, cert.period, 2)354 _, dd = cert.coefficients(s, mp, True)355 r = mp.mpf(mp.nstr(max(10 * t1 / abs(sum(dd)), width_floor), 2))356 assert r < mp.mpf("1e-7")357 vals = []358 for x in (s - r, s + r):359 d, _ = cert.coefficients(iv.mpf(x), iv)360 assert cert.winding(d, 2, t2) == cert.period_graph361 vals.append(sum(d, iv.mpf(0)))362 assert (vals[0] + t1).b < 0 < (vals[1] - t1).a, (vals, t1)363 return s - r, s + r, t1364365366# ZERO367368def zero(period, arcs):369 t0 = time.time()370 mp.dps = iv.dps = 90371 cert = Certificate([(0, 0, 1), (0, 0, 2)], 1, period, complex(0.457015235231, 6.958882679527))372 cert.period_graph = 1373 h = float(cert.hq)374 print(f"zero: A = {{1,2}}, periods 1..{period}, {len(cert.groups)} (period, trace) classes, disc centre {cert.c} radius {cert.r}, conjugation beta {cert.beta}, contraction h = {cert.hq} = {h:.6f}")375 s0 = cert.newton(mp.mpc(*S1_NOTE))376 S0 = iv.mpc(s0.real, s0.imag)377 kk = cert.constant(S0.real, S0.imag, iv.mpf(0), arcs) / cert.h378 tail0 = cert.tail(S0.real, S0.imag, iv.mpf(0), arcs=arcs)379 f0, fp0 = [sum(x, iv.mpc(0)) for x in cert.coefficients(S0, iv, True)]380 r = mp.mpf(mp.nstr(10 * tail0 / abs(mp.mpc(fp0.real.mid, fp0.imag.mid)), 2))381 rho2 = mp.mpf("0.001")382 tail1 = up(iv.mpf(cert.tail(S0.real, S0.imag, iv.mpf(rho2) + 4 * r, arcs=arcs)) / rho2)383 box = iv.mpf([-r, r])384 loose = S0 + iv.mpc(box, box)385 re_text = bracket(low(loose.real), up(loose.real))386 im_text = bracket(low(loose.imag), up(loose.imag))387 X = iv.mpc(iv.mpf(re_text[1:-1].split(", ")), iv.mpf(im_text[1:-1].split(", ")))388 assert up(X.real) - low(X.real) < 3 * r and up(X.imag) - low(X.imag) < 3 * r389 half_re, half_im = [(Fraction(t[1:-1].split(", ")[1]) - Fraction(t[1:-1].split(", ")[0])) / 2 for t in (re_text, im_text)]390 _, fpX = [sum(x, iv.mpc(0)) for x in cert.coefficients(X, iv, True)]391 y = 1 / mp.mpc(fp0.real.mid, fp0.imag.mid)392 Y = iv.mpc(y.real, y.imag)393 b0 = iv.mpf([-tail0, tail0])394 b1 = iv.mpf([-tail1, tail1])395 K = S0 - Y * (f0 + iv.mpc(b0, b0)) + (1 - Y * (fpX + iv.mpc(b1, b1))) * (X - S0)396 inside = K.real.a > X.real.a and K.real.b < X.real.b and K.imag.a > X.imag.a and K.imag.b < X.imag.b397 assert inside398 print(f"zero: weight constant K <= {ceil_sig(up(kk), 6)} on {arcs} arcs, Euler tail past period {period} at s_1 <= {ceil_sig(tail0, 3)}, derivative tail on the box <= {ceil_sig(tail1, 3)}")399 print(f"zero: |D_{period}(s_0)| = {mp.nstr(abs(mp.mpc(f0.real.mid, f0.imag.mid)), 3)}, D_{period}'(s_0) = {mp.nstr(mp.mpc(fp0.real.mid, fp0.imag.mid), 12)}, radius before rounding {mp.nstr(r, 2)}, printed box half-widths {ceil_sig(half_re, 2)} (Re) and {ceil_sig(half_im, 2)} (Im)")400 print(f"zero: Krawczyk image, on the printed box, inside it: {inside}; exactly one zero of D in")401 print(f"zero: Re s in {re_text}")402 print(f"zero: Im s in {im_text}")403 re_t, im_t = trunc(float(X.real.a), 12), trunc(float(X.imag.a), 12)404 note = (float(S1_NOTE[0]), float(S1_NOTE[1]))405 rounds = abs(note[0] - float(s0.real)) < 5e-13 and abs(note[1] - float(s0.imag)) < 5e-13406 assert rounds407 print(f"zero: the note reads {S1_NOTE[0]} + {S1_NOTE[1]} i, the box rounded to twelve digits: {rounds}")408 jc, jr = Fraction("0.758687144013554292899790137015621955739"), Fraction("0.957589818521375342814351002388265920293")409 jh = ivq(max(max(abs(1 / (jc - jr + a) - jc), abs(jc - 1 / (jc + jr + a))) for a in (1, 2)) / jr)410 for sig, tau, name in ((iv.mpf(JP_LO), iv.mpf(0), "dim E_2"), (S0.real, S0.imag, "s_1")):411 ks = []412 for beta in (None, Fraction(33, 20), cert.beta):413 ws = weight_sup([1, 2], jc, jr, beta, sig, tau, iv.mpf(0), arcs)414 ks.append(ceil_sig(up((ws[1] + ws[2]) / (jh * iv.sqrt(1 - jh * jh))), 7))415 print(f"zero: on the Jenkinson-Pollicott disc, h = {mp.nstr(up(jh), 6)}, upper bounds on the constant K at {name}, rounded up: {ks[0]} unconjugated, {ks[1]} at beta = 33/20, {ks[2]} at beta = {cert.beta}")416 print(f"zero: {time.time() - t0:.1f} s")417 return X, re_t, im_t418419420def control(period, arcs):421 t0 = time.time()422 mp.dps = iv.dps = 90423 rows = [((1, 2), JP_LO, JP_HI, "Jenkinson-Pollicott 2018 Theorem 1")]424 rows += [(a, PV[a], PV[a], "Pollicott-Vytnova 2022 Table 3, +- 1e-20") for a in ((1, 3), (2, 3), (1, 2, 3))]425 for alpha, lo, hi, src in rows:426 cert = Certificate([(0, 0, a) for a in alpha], 1, period, 0.5)427 cert.period_graph = 1428 guess = mp.mpf(lo)429 slo, shi, tail = certify_real(cert, guess, mp.mpf("1e-60"))430 if src.startswith("Pollicott"):431 ok = slo < mp.mpf(lo) + mp.mpf("1e-20") and mp.mpf(lo) - mp.mpf("1e-20") < shi432 else:433 ok = slo < mp.mpf(lo) and mp.mpf(hi) < shi434 print(f"control: A = {{{','.join(map(str, alpha))}}}, periods 1..{cert.period}, disc ({cert.c}, {cert.r}, beta {cert.beta}), h = {float(cert.hq):.6f}, tail {mp.nstr(tail, 3)}: dim in {bracket(slo, shi)}, {src} inside: {ok}")435 assert ok436 print(f"control: {time.time() - t0:.1f} s")437438439# ORPHANS440441def orphans(period, modes):442 t0 = time.time()443 mp.dps = iv.dps = 60444 q = sibling()445 x, w = q.nodes(modes)446 B = q.branch_table(x, w, 2000)447 rowsq = q.taylor_rows(x, w, 4)448 _, _, _, live = q.symmetric_orphans(4)449 order = [11892] + [c for c in live if c != 11892]450 print("| code | states | minimal | period | P | h at most | tail at most | dim_CF in | collocation |")451 print("|---|---|---|---|---|---|---|---|---|")452 walls = []453 kinds = {}454 for k, code in [(4, c) for c in order] + [(3, 54)]:455 parts = pieces(q, code, k)456 if any(p[0] == "tail" for p in parts):457 walls.append(code)458 continue459 colloc, _ = q.graph_zero(code, k, x, B, rowsq, 2000)460 found = []461 seen = set()462 for _, n, (s, me) in parts:463 if (s, tuple(me)) in seen:464 continue465 seen.add((s, tuple(me)))466 cert = Certificate(me, s, period, colloc)467 cert.period_graph = period_of(s, me)468 try:469 slo, shi, tail = certify_real(cert, colloc, mp.mpf("1e-45"))470 except AssertionError:471 slo, shi, tail = certify_real(cert, None, mp.mpf("1e-45"))472 found.append((slo, shi, tail, n, s, cert))473 slo, shi, tail, n, s, cert = max(found, key=lambda f: f[0])474 assert all(f[1] < slo for f in found if f[5] is not cert)475 reading = Fraction(math.floor(Fraction(colloc) * 10 ** 12), 10 ** 12)476 cut = [Fraction(math.floor(exact(x) * 10 ** 12), 10 ** 12) for x in (slo, shi)]477 kind = "truncates" if cut[0] == cut[1] == reading else ("contains" if exact(slo) <= reading <= exact(shi) else None)478 ok = kind is not None479 kinds[kind] = kinds.get(kind, 0) + 1480 print(f"| {code} | {n} | {s} | {cert.period_graph} | {cert.period} | {fixed(Fraction(math.ceil(cert.hq * 10 ** 4), 10 ** 4), 4)} | {ceil_sig(tail, 2)} | {bracket(slo, shi)} | {trunc(colloc, 12):.12f} |")481 assert ok, (code, slo, shi, colloc)482 print(f"orphans: the printed twelve-digit reading is the bracket truncated at twelve digits on {kinds.get('truncates', 0)} rows and lies inside the bracket on {kinds.get('contains', 0)}")483 print(f"orphans: a cofinite edge on a cycle, outside the periodic-point determinant: {walls}")484 print(f"orphans: {time.time() - t0:.1f} s")485486487def main():488 ap = argparse.ArgumentParser()489 ap.add_argument("verb", choices=["zero", "control", "orphans", "all"])490 ap.add_argument("--period", type=int, default=18)491 ap.add_argument("--arcs", type=int, default=2000)492 ap.add_argument("--modes", type=int, default=40)493 a = ap.parse_args()494 if a.verb in ("zero", "all"):495 zero(a.period, a.arcs)496 if a.verb in ("control", "all"):497 control(a.period, a.arcs)498 if a.verb in ("orphans", "all"):499 orphans(a.period + 4, a.modes)500501502if __name__ == "__main__":503 main()