restricted_franel.py
23.6 kB · python · 715 lines
1import sys2import time3from fractions import Fraction4from math import gcd, pi56import numpy as np78BASE3 = (3, frozenset({0, 1}), "base 3, digits {0,1}")9BASE10 = (10, frozenset(range(9)), "base 10, digit 9 missing")1011TWO_PI_I = 2j * pi1213# THE DESIGN141516def holds(base, digits, n):17 if n == 0:18 return True19 while n > 0:20 if n % base not in digits:21 return False22 n //= base23 return True242526def flags(base, digits, q):27 keep = np.zeros(q + 1, dtype=bool)28 for n in range(q + 1):29 keep[n] = holds(base, digits, n)30 return keep313233def mobius_sieve(n):34 mu = np.ones(n + 1, dtype=np.int64)35 primes = np.ones(n + 1, dtype=bool)36 primes[:2] = False37 for p in range(2, n + 1):38 if not primes[p]:39 continue40 primes[p * p :: p] = False41 mu[p::p] *= -142 mu[p * p :: p * p] = 043 mu[0] = 044 return mu454647def totients(n):48 phi = np.arange(n + 1, dtype=np.int64)49 for p in range(2, n + 1):50 if phi[p] == p:51 phi[p::p] -= phi[p::p] // p52 return phi535455# THE DILATED MERTENS SUMS565758def mertens_dilated(keep, mu, q, d):59 c = np.arange(1, q // d + 1, dtype=np.int64)60 return int(mu[c][keep[d * c]].sum())616263def mertens_design(keep, mu, q):64 return mertens_dilated(keep, mu, q, 1)656667# THE LITERAL FOURIER SIDE, DENOMINATOR SET686970def literal_denominator(keep, q, freqs):71 out = {m: 0j for m in freqs}72 for b in range(1, q + 1):73 if not keep[b]:74 continue75 a = np.arange(1, b + 1, dtype=np.int64)76 a = a[np.gcd(a, b) == 1]77 for m in freqs:78 r = (m * a) % b79 out[m] += complex(np.exp(TWO_PI_I * (r / b)).sum())80 return out818283def card_denominator(keep, phi, q):84 return int(phi[1 : q + 1][keep[1 : q + 1]].sum())858687def ladder_check(keep, mu, q):88 total = 0j89 worst = 0.090 at = 091 bad = 092 for b in range(1, q + 1):93 if not keep[b]:94 continue95 a = np.arange(1, b + 1, dtype=np.int64)96 a = a[np.gcd(a, b) == 1]97 v = complex(np.exp(TWO_PI_I * (a % b) / b).sum())98 total += v99 e = abs(v - int(mu[b]))100 if e > worst:101 worst, at = e, b102 if round(v.real) != int(mu[b]) or abs(v.imag) > 0.5:103 bad += 1104 return total, worst, at, bad105106107# THE STRICT SET108109110def design_list(keep, q):111 return np.nonzero(keep[1 : q + 1])[0] + 1112113114def strict_literal(keep, q):115 sf = design_list(keep, q)116 total = 0j117 card = 0118 for i, b in enumerate(sf):119 b = int(b)120 a = sf[: i + 1]121 a = a[np.gcd(a, b) == 1]122 card += a.size123 total += complex(np.exp(TWO_PI_I * (a % b) / b).sum())124 return total, card125126127def strict_ramanujan_literal(keep, b):128 a = np.arange(1, b + 1, dtype=np.int64)129 a = a[keep[a] & (np.gcd(a, b) == 1)]130 return complex(np.exp(TWO_PI_I * (a % b) / b).sum())131132133def strict_ramanujan_divisor(keep, mu, b):134 total = 0j135 for d in range(1, b + 1):136 if b % d or mu[d] == 0:137 continue138 n = b // d139 ap = np.arange(1, n + 1, dtype=np.int64)140 ap = ap[keep[d * ap]]141 total += int(mu[d]) * complex(np.exp(TWO_PI_I * (ap % n) / n).sum())142 return total143144145def strict_weights(keep, mu, phi, q):146 sf = design_list(keep, q)147 plain = 0148 counted = 0149 ratio = 0.0150 for i, b in enumerate(sf):151 b = int(b)152 a = sf[: i + 1]153 pf = int((np.gcd(a, b) == 1).sum())154 plain += int(mu[b])155 counted += int(mu[b]) * pf156 ratio += int(mu[b]) * pf / int(phi[b])157 return plain, counted, ratio158159160# THE GCD KERNEL161162163def kernel_sum(keep, mu, q):164 ms = {}165 for d in range(1, q + 1):166 v = mertens_dilated(keep, mu, q, d)167 if v:168 ms[d] = v169 total = Fraction(0)170 ks = sorted(ms)171 for d in ks:172 for e in ks:173 total += Fraction(gcd(d, e) ** 2 * ms[d] * ms[e], d * e)174 return total, ms175176177def fourier_side(keep, q, cut):178 nodes = []179 for b in range(1, q + 1):180 if not keep[b]:181 continue182 for a in range(1, b + 1):183 if gcd(a, b) == 1:184 nodes.append((a, b))185 m = len(nodes)186 num = np.array([a for a, _ in nodes], dtype=np.int64)187 den = np.array([b for _, b in nodes], dtype=np.int64)188 total = 0.0189 step = max(1, 2_000_000 // max(m, 1))190 k = 1191 while k <= cut:192 hi = min(cut, k + step - 1)193 ks = np.arange(k, hi + 1, dtype=np.int64)194 r = (ks[:, None] * num[None, :]) % den[None, :]195 s = np.exp(TWO_PI_I * (r / den[None, :])).sum(axis=1)196 total += float((np.abs(s) ** 2 / ks.astype(float) ** 2).sum())197 k = hi + 1198 return 2.0 * total, m199200201# THE BACKWARD SCAN202203204def kernel_float(keep, mu, w, q):205 ms = np.zeros(q + 1)206 for d in range(1, q + 1):207 ms[d] = mertens_dilated(keep, mu, q, d)208 v = ms[1 : q + 1]209 return float(v @ w[:q, :q] @ v)210211212def gcd_weights(qmax):213 d = np.arange(1, qmax + 1, dtype=np.int64)214 g = np.gcd.outer(d, d).astype(np.float64)215 return g * g / np.outer(d.astype(np.float64), d.astype(np.float64))216217218def scan_backward(design, qmax):219 keep = keep_of(design, qmax)220 mu = mobius_sieve(qmax)221 w = gcd_weights(qmax)222 best = 0.0223 at = 0224 tail = 0.0225 tail_at = 0226 seen = 0227 bad = 0228 for q in range(1, qmax + 1):229 if not keep[q]:230 continue231 seen += 1232 g = kernel_float(keep, mu, w, q)233 mf = mertens_design(keep, mu, q)234 ident = g * pi * pi / 3.0235 r = 2.0 * mf * mf / ident236 if r > 1.0 + 1e-9:237 bad += 1238 if r > best:239 best, at = r, q240 if q >= SCAN_FLOOR and r > tail:241 tail, tail_at = r, q242 return seen, bad, best, at, tail, tail_at243244245# THE CONTROL246247248def farey_delta_square(keep, q):249 nodes = []250 for b in range(1, q + 1):251 if not keep[b]:252 continue253 for a in range(1, b + 1):254 if gcd(a, b) == 1:255 nodes.append(Fraction(a, b))256 nodes.sort()257 m = len(nodes)258 s2 = sum((r - Fraction(j + 1, m)) ** 2 for j, r in enumerate(nodes))259 return m, s2260261262def reflection_closed(keep, q):263 for b in range(1, q + 1):264 if not keep[b]:265 continue266 for a in range(1, b):267 if gcd(a, b) == 1 and not (keep[a] == keep[b - a]):268 return False269 return True270271272# THE VERBS273274FREQS = [1, 2, 3, 4, 5, 6, 12]275276DEN_SMALL = [(BASE3, 2187), (BASE10, 1000), ((0, None, "full set, control"), 300)]277278DEN_LADDER = [279 (BASE3, [243, 2187, 19683, 100000]),280 (BASE10, [100, 1000, 10000]),281 ((0, None, "full set, control"), [100, 1000, 10000]),282]283284285def keep_of(design, q):286 base, digits, _ = design287 if digits is None:288 k = np.ones(q + 1, dtype=bool)289 return k290 return flags(base, digits, q)291292293def verb_denominator():294 print("DENOMINATOR SET, THE FREQUENCY-m IDENTITY")295 for design, q in DEN_SMALL:296 keep = keep_of(design, q)297 mu = mobius_sieve(q)298 lit = literal_denominator(keep, q, FREQS)299 print(f"\n{design[2]}, Q = {q}")300 print(" m | sum_(d|m) d M_F(Q/d; d) | Re literal | Im literal | err")301 for m in FREQS:302 exact = sum(d * mertens_dilated(keep, mu, q, d) for d in range(1, m + 1) if m % d == 0)303 v = lit[m]304 err = abs(v - exact)305 print(f"{m:5d} | {exact:23d} | {v.real:13.9f} | {v.imag:13.9f} | {err:.3e}")306307 print("\nDENOMINATOR SET, FREQUENCY 1 AGAINST THE DESIGN MERTENS SUM")308 print(309 " set | Q | M_F(Q) | Re literal | Im literal | agg err |"310 " per-b err | at b | bad | s"311 )312 for design, ladder in DEN_LADDER:313 for q in ladder:314 t0 = time.time()315 keep = keep_of(design, q)316 mu = mobius_sieve(q)317 exact = mertens_design(keep, mu, q)318 v, worst, at, bad = ladder_check(keep, mu, q)319 dt = time.time() - t0320 print(321 f" {design[2]:26s} | {q:6d} | {exact:7d} | {v.real:13.9f} | "322 f"{v.imag:13.9f} | {abs(v - exact):.2e} | {worst:9.2e} | {at:5d} | {bad:3d} | {dt:4.1f}"323 )324325326SCAN_FLOOR = 100327328SCAN = [(BASE3, 2187), (BASE10, 400), ((0, None, "full set, control"), 400)]329330ID_CASES = [(BASE3, 81, 200000), (BASE3, 243, 200000), (BASE10, 40, 200000), ((0, None, "full set, control"), 40, 200000)]331332333def verb_identity():334 print("\nTHE DIGIT-RESTRICTED FRANEL IDENTITY, GCD-KERNEL SIDE AGAINST THE FOURIER SIDE")335 print(" set | Q | m | G_F(Q) | (pi^2/3) G_F | truncated | gap | tail bound")336 for design, q, cut in ID_CASES:337 keep = keep_of(design, q)338 mu = mobius_sieve(q)339 g, _ = kernel_sum(keep, mu, q)340 ident = float(g) * pi * pi / 3.0341 lhs, m = fourier_side(keep, q, cut)342 tail = 2.0 * m * m / cut343 print(344 f" {design[2]:26s} | {q:4d} | {m:4d} | {float(g):11.6f} | {ident:12.6f} | "345 f"{lhs:10.6f} | {abs(ident - lhs):.2e} | {tail:.4f}"346 )347348 print("\n THE RANK FORM, G_F(Q) - 1 = 12 m sum delta^2, exact rationals")349 print(" set | Q | m | sum delta^2 | G_F - 1 == 12 m sum delta^2")350 for design, q in [(BASE3, 81), (BASE3, 243), (BASE10, 40), ((0, None, "full set, control"), 40)]:351 keep = keep_of(design, q)352 mu = mobius_sieve(q)353 g, _ = kernel_sum(keep, mu, q)354 m, s2 = farey_delta_square(keep, q)355 print(f" {design[2]:26s} | {q:4d} | {m:4d} | {float(s2):13.10f} | {g - 1 == 12 * m * s2}")356357 print("\n the backward inequality 2 M_F(Q)^2 <= (pi^2/3) G_F(Q), sampled")358 print(" set | Q | M_F(Q) | (pi^2/3) G_F | ratio")359 for design, ladder in [(BASE3, [243, 2187]), (BASE10, [100, 1000]), ((0, None, "full set, control"), [100, 1000])]:360 for q in ladder:361 keep = keep_of(design, q)362 mu = mobius_sieve(q)363 g, _ = kernel_sum(keep, mu, q)364 mf = mertens_design(keep, mu, q)365 ident = float(g) * pi * pi / 3.0366 print(f" {design[2]:26s} | {q:6d} | {mf:7d} | {ident:14.6f} | {2.0 * mf * mf / ident:.6f}")367368 print("\n the same inequality asserted at EVERY integer Q up to the bound")369 print(" both sides step only at Q in S_F, so scanning S_F covers every integer Q")370 print(" set | Qmax | jumps | violations | max ratio | at Q | max over Q >= 100 | at Q")371 for design, qmax in SCAN:372 seen, bad, best, at, tail, tail_at = scan_backward(design, qmax)373 assert bad == 0374 print(375 f" {design[2]:26s} | {qmax:6d} | {seen:5d} | {bad:10d} | {best:10.6f} | {at:5d} | "376 f"{tail:18.6f} | {tail_at:5d}"377 )378379380STRICT_LADDER = [(BASE3, [243, 2187, 19683, 177147]), (BASE10, [100, 1000, 10000])]381382383def verb_strict():384 print("\nTHE STRICT SET, THE DIVISOR-SUM IDENTITY")385 print(" set | Q | denominators | max abs diff")386 for design, q in [(BASE3, 729), (BASE10, 200)]:387 keep = keep_of(design, q)388 mu = mobius_sieve(q)389 worst = 0.0390 n = 0391 for b in range(1, q + 1):392 if not keep[b]:393 continue394 n += 1395 worst = max(worst, abs(strict_ramanujan_literal(keep, b) - strict_ramanujan_divisor(keep, mu, b)))396 print(f" {design[2]:26s} | {q:4d} | {n:12d} | {worst:.3e}")397398 print("\nTHE STRICT SET, THE FIRST Q WITH A NONZERO IMAGINARY PART")399 for design, q in [(BASE3, 100), (BASE10, 100)]:400 keep = keep_of(design, q)401 run = 0j402 hit = None403 for b in range(1, q + 1):404 if not keep[b]:405 continue406 run += strict_ramanujan_literal(keep, b)407 if hit is None and abs(run.imag) > 1e-9:408 hit = (b, run)409 b, run = hit410 print(f" {design[2]:26s}: Q = {b}, sum = {run.real:.9f} + {run.imag:.9f} i")411412 print("\nTHE STRICT SET, FREQUENCY 1 AGAINST EVERY MERTENS-TYPE CANDIDATE")413 print(414 " set | Q | card | Re T | Im T | abs T | abs T/card |"415 " M_F(Q) | sum mu phi_F | sum mu phi_F/phi"416 )417 for design, ladder in STRICT_LADDER:418 for q in ladder:419 keep = keep_of(design, q)420 mu = mobius_sieve(q)421 phi = totients(q)422 t, card = strict_literal(keep, q)423 plain, counted, ratio = strict_weights(keep, mu, phi, q)424 print(425 f" {design[2]:26s} | {q:6d} | {card:8d} | {t.real:8.3f} | {t.imag:8.3f} | "426 f"{abs(t):9.3f} | {abs(t) / card:10.6f} | {plain:7d} | {counted:13d} | {ratio:17.6f}"427 )428429430# THE DILATE AUTOMATON431432433def dilate_matrix(base, digits, d):434 t = [[0] * d for _ in range(d)]435 for r in range(d):436 for e in range(base):437 v = d * e + r438 if v % base in digits:439 t[r][v // base] += 1440 return t441442443def dilate_accept(base, digits, d):444 return [1 if (r == 0 or holds(base, digits, r)) else 0 for r in range(d)]445446447def qfree(d, base):448 while True:449 g = gcd(d, base)450 if g == 1:451 return d452 d //= g453454455def automaton_count(base, digits, d, power):456 t = dilate_matrix(base, digits, d)457 v = [0] * d458 v[0] = 1459 for _ in range(power):460 w = [0] * d461 for i, x in enumerate(v):462 if x:463 for j, y in enumerate(t[i]):464 if y:465 w[j] += x * y466 v = w467 return sum(x * y for x, y in zip(v, dilate_accept(base, digits, d)))468469470def dilate_mask(base, digits, d, x):471 lut = np.zeros(base, dtype=bool)472 for f in digits:473 lut[f] = True474 t = d * np.arange(1, x + 1, dtype=np.int64)475 ok = np.ones(x, dtype=bool)476 while t.any():477 ok &= lut[t % base] | (t == 0)478 t //= base479 return ok480481482def divisor_counts(base, digits, cut):483 m = np.nonzero(flags(base, digits, cut))[0]484 m = m[m > 0]485 n = np.zeros(cut + 1, dtype=np.int64)486 for d in range(1, cut + 1):487 n[d] = int((m % d == 0).sum())488 return n489490491def smith_bilinear(n, cut, block=400):492 u = np.sqrt(n[1 : cut + 1].astype(np.float64))493 idx = np.arange(1, cut + 1, dtype=np.int64)494 total = 0.0495 for lo in range(0, cut, block):496 hi = min(lo + block, cut)497 g = np.gcd.outer(idx[lo:hi], idx).astype(np.float64)498 w = g * g / (idx[lo:hi, None] * idx[None, :])499 total += float((w * u[lo:hi, None] * u[None, :]).sum())500 return total501502503DILATE_D3 = [1, 2, 3, 4, 5, 7, 8, 11, 13, 16, 22, 31]504DILATE_D10 = [1, 2, 3, 4, 5, 7, 11, 13, 121, 243, 729]505506507def verb_dilate():508 base, digits, name = BASE3509 print(f"\nTHE DILATE AUTOMATON, {name}, d = 2, states are the carries")510 t = dilate_matrix(base, digits, 2)511 acc = dilate_accept(base, digits, 2)512 for r, row in enumerate(t):513 print(f" carry {r} -> {row} accept {acc[r]}")514 print(f" row sums {[sum(r) for r in t]}, column sums {[sum(c) for c in zip(*t)]}, |F| = {len(digits)}")515 print("\n L | 3^L | transfer matrix | brute force | 2^L - 1")516 for lvl in range(1, 13):517 x = base**lvl518 auto = automaton_count(base, digits, 2, lvl) - 1519 brute = int(dilate_mask(base, digits, 2, x - 1).sum())520 print(f" {lvl:2d} | {x:6d} | {auto:15d} | {brute:11d} | {2**lvl - 1:10d} {'ok' if auto == brute else 'MISMATCH'}")521522 print("\nTHE COLUMN-SUM LAW, every dilate d <= 64, both designs")523 for bb, dd, nm in (BASE3, BASE10):524 bad, off = 0, []525 for d in range(1, 65):526 m = dilate_matrix(bb, dd, d)527 if any(sum(c) != len(dd) for c in zip(*m)):528 bad += 1529 if any(sum(r) != len(dd) for r in m):530 off.append(d)531 print(f" {nm:26s} columns off |F| at {bad} of 64; rows off |F| at {len(off)} of 64, every one with gcd(d,q) > 1: {all(gcd(d, bb) > 1 for d in off)}")532533 for bb, dd, nm, ds in ((BASE3[0], BASE3[1], BASE3[2], DILATE_D3), (BASE10[0], BASE10[1], BASE10[2], DILATE_D10)):534 al = np.log(len(dd)) / np.log(bb)535 print(f"\nTHE ACCEPT MASS, {nm}, K_d = A_d(q^24)/|F|^24 against |Acc_d|/d")536 print(" d | gcd(d,q) | K_d | |Acc_d|/d | d^(alpha-1)")537 for d in ds + ([8, 16, 32] if bb == 10 else []):538 k = automaton_count(bb, dd, d, 24) / len(dd) ** 24539 a = sum(dilate_accept(bb, dd, d))540 print(f" {d:6d} | {gcd(d, bb):8d} | {k:7.4f} | {a / d:9.6f} | {d ** (al - 1):11.6f}")541542 for bb, dd, nm, top, ds in ((3, BASE3[1], BASE3[2], 12, DILATE_D3), (10, BASE10[1], BASE10[2], 7, DILATE_D10[:7])):543 x = bb**top544 mu = mobius_sieve(x)545 al = np.log(len(dd)) / np.log(bb)546 root = x ** (al / 2)547 print(f"\nTHE DILATED METER, {nm}, x = {bb}^{top} = {x}")548 print(" d | A_d(x) | M_F(x;d) | max |M| | log max/log x | same at x/q | (U) ratio | (U\u0027) ratio | local exponents")549 for d in ds:550 ok = dilate_mask(bb, dd, d, x)551 run = np.cumsum(np.where(ok, mu[1 : x + 1], 0))552 peaks = [int(np.abs(run[: bb**lvl]).max()) for lvl in range(top - 4, top + 1)]553 loc = [round(float(np.log(peaks[i + 1] / peaks[i]) / np.log(bb)), 3) for i in range(len(peaks) - 1)]554 mass, peak = int(ok.sum()), peaks[-1]555 ru = peak / (d ** ((al - 1) / 2) * root)556 rv = peak / (qfree(d, bb) ** ((al - 1) / 2) * root)557 under = np.log(peaks[-2]) / np.log(x / bb)558 print(559 f" {d:6d} | {mass:7d} | {int(run[-1]):9d} | {peak:7d} | {np.log(peak) / np.log(x):13.6f} | "560 f"{under:11.6f} | {ru:9.4f} | {rv:10.4f} | {loc}"561 )562563 print(f"\nTHE LADDER, {BASE3[2]}, the meter against its dilate at d = 2")564 x = 3**12565 mu = mobius_sieve(x)566 r1 = np.cumsum(np.where(dilate_mask(3, BASE3[1], 1, x), mu[1 : x + 1], 0))567 r2 = np.cumsum(np.where(dilate_mask(3, BASE3[1], 2, x), mu[1 : x + 1], 0))568 print(" L | M_F(3^L) | M_F(3^L;2) | max |M_F| | max |M_F(;2)|")569 for lvl in range(1, 13):570 n = 3**lvl571 print(572 f" {lvl:2d} | {int(r1[n - 1]):9d} | {int(r2[n - 1]):11d} | {int(np.abs(r1[:n]).max()):9d} | "573 f"{int(np.abs(r2[:n]).max()):13d}"574 )575576 al = np.log(2) / np.log(3)577 print(f"\nTHE SMITH REDUCTION, B(Q) = sum gcd(d,e)^2/(d e) sqrt(N_F(Q;d) N_F(Q;e)), {BASE3[2]}")578 print(" Q | A_F(Q) | B(Q) | B/Q^alpha | B/(Q^alpha ln Q) | B/(Q^alpha ln^2 Q) | local exponent")579 prev = None580 for lvl in range(4, 9):581 cut = 3**lvl582 b = smith_bilinear(divisor_counts(3, BASE3[1], cut), cut)583 qa = cut**al584 loc = "" if prev is None else f"{np.log(b / prev) / np.log(3):14.3f}"585 print(586 f" {cut:8d} | {2**lvl:7d} | {b:11.3f} | {b / qa:9.4f} | {b / (qa * np.log(cut)):16.4f} | "587 f"{b / (qa * np.log(cut) ** 2):18.4f} | {loc}"588 )589 prev = b590591592# THE CONVERSE593594595def digit_gap(digits):596 f = sorted(digits)597 g = 0598 for x in f[1:]:599 g = gcd(g, x - f[0])600 return g601602603def brute_dilate(base, digits, d, power):604 return sum(1 for c in range(1, base**power) if holds(base, digits, d * c))605606607def positive_count(base, digits, d, power):608 return automaton_count(base, digits, d, power) - (1 if 0 in digits else 0)609610611def subsets(base):612 for m in range(1, 1 << base):613 yield {e for e in range(base) if m >> e & 1}614615616SCOPE_BASES = (3, 4, 5)617GAP_BASES = range(2, 8)618GAP_POWER = 400619620621def verb_converse():622 print("\nTHE SCOPE OF D1, D2, D5, automaton against brute force, q = 3, 4, 5, every F, d <= 6, L <= 5")623 tally = {True: [0, 0], False: [0, 0]}624 for base in SCOPE_BASES:625 for digits in subsets(base):626 for d in range(1, 7):627 for lvl in range(1, 6):628 row = tally[0 in digits]629 row[0] += 1630 row[1] += positive_count(base, digits, d, lvl) != brute_dilate(base, digits, d, lvl)631 for key, label in ((True, "0 in F"), (False, "0 not in F")):632 n, bad = tally[key]633 print(f" {label:12s} {bad:4d} mismatches of {n:4d}")634 base, digits, d, lvl = 3, {1, 2}, 1, 3635 print(f" witness q = 3, F = {{1,2}}, d = 1, L = 3: automaton {positive_count(base, digits, d, lvl)}, "636 f"true {brute_dilate(base, digits, d, lvl)}, |F|^L {len(digits) ** lvl}")637638 print(f"\nTHE MASS CONSTANT, K_d against |Acc_d|/d at L = {GAP_POWER}, every F with 0 in F, gcd(d,q) = 1, 2 <= d <= 24")639 split = {True: [0, 0], False: [0, 0]}640 odd = []641 for base in GAP_BASES:642 for digits in subsets(base):643 if 0 not in digits or len(digits) < 2:644 continue645 gap = digit_gap(digits)646 for d in range(2, 25):647 if gcd(d, base) != 1:648 continue649 k = automaton_count(base, digits, d, GAP_POWER) / len(digits) ** GAP_POWER650 row = split[gcd(d, gap) == 1]651 row[0] += 1652 good = abs(k - sum(dilate_accept(base, digits, d)) / d) < 1e-9653 row[1] += good654 if good == (gcd(d, gap) > 1):655 odd.append((base, sorted(digits), d, round(k, 6), sum(dilate_accept(base, digits, d)) / d))656 for key, label in ((True, "gcd(d, gap) = 1"), (False, "gcd(d, gap) > 1")):657 n, ok = split[key]658 print(f" {label:16s} {ok:5d} agree, {n - ok:5d} fail, of {n:5d}")659 print(f" off the split: {odd}")660 base, digits, d = 3, {0, 2}, 2661 print(f" witness q = 3, F = {{0,2}}, gap 2, d = 2: counts "662 f"{[brute_dilate(base, digits, d, lvl) + 1 for lvl in range(1, 9)]}, |F|^L, so K_2 = 1 against "663 f"|Acc_2|/2 = {sum(dilate_accept(base, digits, d)) / d}")664665 base, digits, name = BASE3666 al = np.log(len(digits)) / np.log(base)667 top = 12668 x = base**top669 mu = mobius_sieve(x)670 run = np.cumsum(np.where(dilate_mask(base, digits, 1, x), mu[1 : x + 1], 0))671 peak, tail = int(np.abs(run).max()), int(run[-1])672 print(f"\n(U) IS UNSATISFIABLE, {name}, x = {base}^{top}, M_F(x; {base}^j) = M_F(x) = {tail} at every j by D4")673 print(" j | d = q^j | d^((alpha-1)/2) x^(alpha/2) | |M_F| / bound | max |M| / bound")674 for j in range(top + 1):675 d = base**j676 bound = d ** ((al - 1) / 2) * x ** (al / 2)677 print(f" {j:4d} | {d:14d} | {bound:27.4f} | {abs(tail) / bound:13.3f} | {peak / bound:15.3f}")678679 base, digits, name = BASE10680 print(f"\nTHE BASE-SMOOTH MASS, {name}, K_d at L = 24 over the base-smooth d <= 1000")681 smooth = sorted(d for d in range(1, 1001) if qfree(d, base) == 1)682 masses = {d: automaton_count(base, digits, d, 24) / len(digits) ** 24 for d in smooth}683 lo = min(masses, key=masses.get)684 hi = max(masses, key=masses.get)685 print(f" {len(smooth)} such d; K_d in [{masses[lo]:.4f} at d = {lo}, {masses[hi]:.4f} at d = {hi}]")686 print(f" sample {[(d, round(masses[d], 4)) for d in (2, 4, 5, 8, 16, 32, 64, 128, 256, 512, 625)]}")687 ratio = {}688 for d in range(1, 201):689 f = qfree(d, base)690 k = automaton_count(base, digits, d, 24) / len(digits) ** 24691 kf = masses.get(f) or automaton_count(base, digits, f, 24) / len(digits) ** 24692 masses[f] = kf693 ratio[d] = k / kf694 lo = min(ratio, key=ratio.get)695 hi = max(ratio, key=ratio.get)696 print(f" K_d / K_(q-free part) over d <= 200 in [{ratio[lo]:.4f} at d = {lo}, {ratio[hi]:.4f} at d = {hi}]")697698699def main():700 verbs = {701 "denominator": verb_denominator,702 "identity": verb_identity,703 "strict": verb_strict,704 "dilate": verb_dilate,705 "converse": verb_converse,706 }707 want = sys.argv[1:] or list(verbs)708 t0 = time.time()709 for v in want:710 verbs[v]()711 print(f"\ntotal runtime {time.time() - t0:.1f} s")712713714if __name__ == "__main__":715 main()