ubound.py
13.2 kB · python · 363 lines
1import itertools2import math3import sys4import time56import numpy as np7from mpmath import iv89FAILS = []1011# CHECKS1213def chk(value, want, tol):14 ok = abs(value - want) <= tol15 if not ok:16 FAILS.append((value, want))17 return "OK" if ok else f"MISMATCH want {want}"1819def chk_le(value, cap, name):20 ok = value <= cap21 if not ok:22 FAILS.append((name, value, cap))23 return "OK" if ok else f"OVER {cap}"2425def close(t0):26 if FAILS:27 raise SystemExit(f"{len(FAILS)} rows off: {FAILS[:4]}")28 print(f"all rows agree ({time.time() - t0:.1f}s)")2930# FLOAT TRANSFORM3132def hat_missing(q, a0, t):33 th = np.asarray(t) % 1.034 s = np.sin(np.pi * th)35 d = np.where(np.abs(s) < 1e-14, float(q), np.sin(q * np.pi * th) / np.where(np.abs(s) < 1e-14, 1.0, s))36 ph = 2 * np.pi * (a0 - (q - 1) / 2.0) * th37 return np.sqrt(np.maximum(d * d - 2 * d * np.cos(ph) + 1.0, 0.0)) / (q - 1)3839def dirichlet(q, t):40 th = np.asarray(t) % 1.041 s = np.sin(np.pi * th)42 return np.where(np.abs(s) < 1e-14, float(q), np.abs(np.sin(q * np.pi * th) / np.where(np.abs(s) < 1e-14, 1.0, s)))4344def hat_excl(q, E, t):45 th = np.asarray(t) % 1.046 acc = np.zeros(th.shape, dtype=complex)47 for d in range(q):48 if d not in E:49 acc = acc + np.exp(2j * np.pi * d * th)50 return np.abs(acc) / (q - len(E))5152def domination(bases, m=1, pts=4001):53 t0 = time.time()54 t = (np.arange(pts) + 0.5) / pts55 for q in bases:56 u = np.minimum(1.0, (dirichlet(q, t) + float(m)) / (q - m))57 worst = 0.058 gap = 1.059 n = 060 for E in itertools.combinations(range(q), m):61 h = hat_excl(q, E, t)62 worst = max(worst, float((h - u).max()))63 gap = min(gap, float((u - h).max()))64 n += 165 print(f"base {q} at {m} excluded digits: all {n} sets, max(|hat F| - u_q) = {worst:.3e} {chk_le(worst, 1e-12, f'dom {q}')}, every set slack at least {gap:.6f}")66 close(t0)6768# THE SUBSET IDENTITY6970def runs(E, N):71 out = []72 s = None73 for j in range(N + 1):74 if j < N and j in E:75 if s is None:76 s = j77 elif s is not None:78 out.append((s, j - s))79 s = None80 return out8182def sigma_u(q, N, x, m=1):83 i = np.arange(q ** N)84 t = x + i / q ** N85 p = np.ones(len(i))86 for j in range(N):87 p *= (dirichlet(q, (q ** j) * t) + float(m)) / (q - m)88 return float(p.sum())8990def subset_sum(q, N, x, m=1):91 i = np.arange(q ** N)92 t = x + i / q ** N93 tot = 0.094 for msk in range(1 << N):95 E = {j for j in range(N) if msk >> j & 1}96 p = np.ones(len(i))97 for s, l in runs(E, N):98 p *= dirichlet(q ** l, (q ** s) * t)99 tot += float(m) ** (N - len(E)) * float(p.sum())100 return tot / (q - m) ** N101102def identity(rows, m=1):103 t0 = time.time()104 for q, N in rows:105 for x in (0.0, 0.1234567, 0.5, 1.0 / (3 * q ** N), 0.4999):106 a, b = sigma_u(q, N, x, m), subset_sum(q, N, x, m)107 if abs(a - b) > 1e-9 * max(1.0, a):108 FAILS.append((q, N, x, a, b))109 print(f"base {q} levels {N} at {m} excluded digits: Sigma_N^u against the subset sum over {2 ** N} runs decompositions weighted {m}^(N - |E|) {chk(0.0, 0.0, 0.0)}, Sigma_N^u(0) = {sigma_u(q, N, 0.0, m):.6f}")110 close(t0)111112# THE LEBESGUE CONSTANT113114def lam_at(M, th):115 j = np.arange(M)116 d = np.minimum((j + th) / M, 1.0 - (j + th) / M)117 d = np.where(d <= 0.0, 1.0 / (2 * M), d)118 return float((abs(math.sin(math.pi * th)) / np.sin(np.pi * d)).sum())119120def lam_scan(M, grid=4000):121 ths = np.linspace(0.0, 0.5, grid + 1)[1:]122 vals = [lam_at(M, float(t)) for t in ths]123 k = int(np.argmax(vals))124 return vals[k], float(ths[k])125126def lebesgue(Ms):127 t0 = time.time()128 c1 = 2.0 / math.pi129 for M in Ms:130 v, arg = lam_scan(M)131 half = lam_at(M, 0.5)132 cap = c1 * math.log(M) + G0 + c1 / M133 print(f"M = {M}: L_M = {v:.6f} at offset {arg:.6f}, half-offset value {half:.6f}, L_M/M - (2/pi) log M = {v / M - c1 * math.log(M):.6f} under gamma' + (2/pi)/M = {G0 + c1 / M:.7f} {chk_le(v / M - cap, 0.0, f'lam {M}')} {chk_le(v - half, 1e-9 * M, f'half {M}')}")134 close(t0)135136def least_uniform(lo, hi, nd, m=1):137 t0 = time.time()138 rows = [(q, uniform_alpha(q, nd, m)[0]) for q in range(lo, hi)]139 first = next((q for q, e in rows if e < 0.25), None)140 if first is None:141 best = min(rows, key=lambda r: r[1])142 print(f"the digit-uniform window bound at {nd} digits and {m} excluded digits clears 1/4 nowhere in [{lo}, {hi}): the best row is q = {best[0]} at alpha_1 < {best[1]:.6f}, so the {m}-digit uniform window ceiling sits above {hi - 1}")143 close(t0)144 return145 bad = [q for q, e in rows if q > first and e >= 0.25]146 print(f"the digit-uniform window bound at {nd} digits and {m} excluded digits clears 1/4 first at q = {first}, and at every base above it in [{lo}, {hi}) {chk(float(len(bad)), 0.0, 0.0)}")147 print(" " + " ".join(f"q={q}:{e:.6f}" for q, e in rows if first - 4 <= q <= first + 3))148 close(t0)149150# THE PEELING BOUND151152def lam_max(M, grid=20000):153 return lam_scan(M, grid)[0]154155def a_seq(q, N, lams, m=1):156 a = {-1: 1.0, 0: 1.0}157 for n in range(1, N + 1):158 a[n] = m * a[n - 1] + m * sum(lams[l] * a[n - 1 - l] for l in range(1, n)) + lams[n]159 return a160161def peel(rows, m=1, grid=400):162 t0 = time.time()163 for q, N in rows:164 lams = {l: lam_max(q ** l) / q ** l for l in range(1, N + 1)}165 xs = (np.arange(grid) + 0.5) / (grid * q ** N)166 i = np.arange(q ** N)167 for msk in range(1, 1 << N):168 E = {j for j in range(N) if msk >> j & 1}169 rs = runs(E, N)170 best = 0.0171 for x in xs:172 t = x + i / q ** N173 pr = np.ones(len(i))174 for sr, l in rs:175 pr *= dirichlet(q ** l, (q ** sr) * t)176 best = max(best, float(pr.sum()))177 cap = q ** N * math.prod(lams[l] for _, l in rs)178 if best > cap * (1 + 1e-9):179 FAILS.append((q, N, sorted(E), best, cap))180 a = a_seq(q, N, lams, m)181 sig = max(sigma_u(q, N, float(x), m) for x in xs)182 cap = a[N] * (q / (q - m)) ** N183 if sig > cap * (1 + 1e-9):184 FAILS.append(("chain", q, N, sig, cap))185 print(f"base {q} levels {N} at {m} excluded digits: every one of {2 ** N - 1} run terms under q^N prod lambda {chk(0.0, 0.0, 0.0)}, and max Sigma_N^u = {sig:.6f} under a_N (q/(q-m))^N = {cap:.6f}")186 close(t0)187188# THE THRESHOLD189190C1 = 2.0 / math.pi191G0 = 0.9625229192193def root_z(q, m=1, g0=G0):194 L = math.log(q)195 f = lambda z: (z - m) * (z - 1) ** 2 - m * (C1 * L * z + g0 * (z - 1) + C1 * (z - 1) ** 2 / (q * z - 1))196 lo, hi = float(m), float(m) + 10.0197 while f(hi) < 0.0:198 hi *= 2.0199 for _ in range(300):200 mid = 0.5 * (lo + hi)201 if f(mid) < 0.0:202 lo = mid203 else:204 hi = mid205 return hi206207def clears(q, m=1, thr=0.25, g0=G0):208 return root_z(q, m, g0) * q / (q - m) < q ** thr209210def threshold(m=1, thr=0.25, name="1/4", g0=G0, hi=100000):211 t0 = time.time()212 qu = next(q for q in range(m + 2, hi) if all(clears(r, m, thr, g0) for r in range(q, min(q + 2000, hi))))213 bad = [q for q in range(qu, hi) if not clears(q, m, thr, g0)]214 z = root_z(qu, m, g0)215 print(f"{m} excluded digits against {name}, gamma' = {g0}: threshold q_u = {qu}, growth root z = {z:.6f} against q^({name})(1 - {m}/q) = {qu ** thr * (1 - float(m) / qu):.6f}, alpha_1 < {math.ceil(math.log(z * qu / (qu - m)) / math.log(qu) * 1e6) / 1e6:.6f}")216 print(f" the chain fails at q = {qu - 1} and at no base in [{qu}, {hi}) {chk(float(len(bad)), 0.0, 0.0)}")217 for q in (qu, 2 * qu, 10 * qu, 10 ** 6):218 z = root_z(q, m, g0)219 print(f" q = {q}: alpha_1 < {math.ceil(math.log(z * q / (q - m)) / math.log(q) * 1e6) / 1e6:.6f}")220 close(t0)221222# THE THRESHOLD, CERTIFIED223224def certify_threshold(qu, hi, m=1, thr=(1, 4), name="1/4", g0=G0):225 t0 = time.time()226 iv.prec = 120227 c1 = 2 / iv.pi228 e = iv.mpf(thr[0]) / thr[1]229 worst = None230 gam = c1 * (iv.euler + iv.log(8 / iv.pi))231 if not float(gam.b) <= g0:232 raise SystemExit("gamma' constant is below the true Lebesgue constant")233 def margin(q):234 w = iv.mpf(q) ** e * (1 - iv.mpf(m) / q)235 return (w - m) * (w - 1) ** 2 - iv.mpf(m) * (c1 * iv.log(q) * w + iv.mpf(g0) * (w - 1) + c1 * (w - 1) ** 2 / (q * w - 1))236 for q in range(qu, hi):237 g = float(margin(q).a)238 if g <= 0.0:239 FAILS.append((q, g))240 if worst is None or g < worst[1]:241 worst = (q, g)242 print(f"certified at 120 bits at {m} excluded digits, gamma' = {g0} above the true {float(gam.b):.10f}: (w - {m})(w - 1)^2 > {m}((2/pi) log q w + gamma'(w - 1) + (2/pi)(w-1)^2/(q w - 1)) at w = q^({name})(1 - {m}/q) for every q in [{qu}, {hi}), tightest margin {worst[1]:.6e} at q = {worst[0]} {chk(float(len(FAILS)), 0.0, 0.0)}")243 print(f" and the chain fails at q = {qu - 1}, margin {float(margin(qu - 1).b):.6e}, so {qu} is the least threshold this chain gives at {m} excluded digits")244 close(t0)245246# THE UNIFORM WINDOW BOUND247248def dn(x):249 return np.nextafter(np.asarray(x, float), -np.inf)250251def up(x):252 return np.nextafter(np.asarray(x, float), np.inf)253254def sup_sin(A, B):255 hit = np.ceil(A - 0.5) <= np.floor(B - 0.5)256 ends = np.maximum(np.abs(np.sin(np.pi * A)), np.abs(np.sin(np.pi * B)))257 return np.where(hit, 1.0, up(up(ends) * (1.0 + 2.0 ** -45)))258259def uniform_G(q, nd, m=1):260 W = q ** nd261 w = np.arange(W, dtype=np.float64)262 delta = np.minimum(w, W - 1 - w) / W263 sn = dn(dn(np.sin(np.pi * delta)) * (1.0 - 2.0 ** -45))264 with np.errstate(divide="ignore"):265 S2 = np.where(delta <= 0.0, np.inf, up(1.0 / np.maximum(sn, 0.0)))266 S1 = sup_sin(w / q ** (nd - 1), (w + 1) / q ** (nd - 1))267 D = np.minimum(float(q), up(S1 * S2))268 return np.minimum(1.0, up(up(D + float(m)) / (q - m)))269270def perron_up(G, q, nd, iters=20000):271 S = q ** (nd - 1)272 tgt = np.tile(np.arange(S), q)273 y = np.ones(S)274 lam = 0.0275 streak = 0276 for it in range(iters):277 z = (G * y[tgt]).reshape(S, q).sum(axis=1)278 nl = float(z.max())279 y = z / nl280 streak = streak + 1 if abs(nl - lam) <= 1e-13 * nl else 0281 lam = nl282 if streak >= 50 and it >= 300:283 break284 y = np.maximum(y, 1e-9 * float(y.max()))285 r = up(G * y[tgt]).reshape(S, q)286 acc = np.zeros(S)287 for c in range(q):288 acc = up(acc + r[:, c])289 return float(up(acc / y).max())290291def uniform_alpha(q, nd, m=1):292 mu = perron_up(uniform_G(q, nd, m), q, nd)293 e = math.ceil(math.log(mu) / math.log(q) * 1e6) / 1e6294 if not mu < q ** e:295 raise SystemExit("rounding unsafe")296 return e, mu297298def window(rows, m=1):299 t0 = time.time()300 for q, nd in rows:301 e, mu = uniform_alpha(q, nd, m)302 v = "CLEARS 1/4" if e < 0.25 else "does not clear 1/4"303 print(f"base {q} uniform windows {nd} digits at {m} excluded digits: lambda < {mu:.6f} -> alpha_1 < {e:.6f} at every excluded set [{v}]")304 close(t0)305306def compare(rows):307 t0 = time.time()308 for q, nd, a0, per in rows:309 e, mu = uniform_alpha(q, nd)310 ok = e >= per311 if not ok:312 FAILS.append((q, e, per))313 print(f"base {q} missing {a0}: uniform {e:.6f} against the per-digit ladder {per:.6f}, uniform weaker by {e - per:.6f} {chk(1.0 if ok else 0.0, 1.0, 0.0)}")314 close(t0)315316# VERBS317318def main():319 verb = sys.argv[1] if len(sys.argv) > 1 else "check"320 if verb == "check":321 domination([3, 4, 10, 20, 21, 34, 76, 99, 200])322 identity([(3, 4), (5, 3), (10, 3), (21, 2), (34, 2)])323 lebesgue([4, 21, 34, 76, 99, 200, 1000, 100000])324 peel([(3, 4), (5, 3), (10, 3)])325 threshold()326 certify_threshold(125, 3000)327 compare([(21, 4, 0, 0.250088), (34, 4, 16, 0.249371)])328 domination([4, 10, 21, 34], 2)329 identity([(4, 3), (10, 2), (21, 2)], 2)330 peel([(4, 3), (5, 3)], 2)331 threshold(2)332 certify_threshold(649, 3000, 2)333 if verb == "domination":334 domination([3, 4, 10, 20, 21, 34, 76, 99, 200])335 if verb == "peel":336 peel([(3, 4), (4, 4), (5, 3), (7, 3), (10, 3)])337 if verb == "least":338 least_uniform(35, 90, 3)339 if verb == "window":340 window([(21, 2), (21, 3), (21, 4), (34, 3), (34, 4), (76, 3), (99, 3), (126, 3), (200, 3)])341 if verb == "compare":342 compare([(9, 5, 0, 0.323432), (10, 5, 5, 0.350684), (20, 4, 6, 0.283414), (21, 4, 0, 0.250088), (33, 4, 15, 0.253000), (34, 4, 16, 0.249371)])343 if verb == "threshold":344 threshold()345 certify_threshold(125, 3000)346 if verb == "pairs":347 domination([4, 6, 10, 21, 34], 2)348 identity([(4, 3), (5, 3), (10, 2), (21, 2), (34, 2)], 2)349 peel([(4, 3), (5, 3), (7, 3)], 2)350 threshold(2)351 certify_threshold(649, 3000, 2)352 threshold(2, 1 / 3, "1/3")353 certify_threshold(105, 3000, 2, (1, 3), "1/3")354 threshold(1, 1 / 3, "1/3")355 certify_threshold(32, 3000, 1, (1, 3), "1/3")356 least_uniform(100, 300, 2, 2)357 if verb == "lebesgue":358 lebesgue([4, 9, 21, 34, 76, 99, 200, 1000, 5776, 100000])359 if verb == "identity":360 identity([(3, 4), (5, 3), (10, 3), (21, 2), (34, 2)])361362if __name__ == "__main__":363 main()