primes.py
18.8 kB · python · 469 lines
1import itertools2import json3import math4import os5import random6import sys7import time8from fractions import Fraction910import numpy as np11from mpmath import iv1213HERE = os.path.dirname(os.path.abspath(__file__))14DATA_DIR = os.path.join("data", os.path.relpath(HERE))15sys.path.insert(0, os.path.join(HERE, "..", "mobius-dissection"))16import dissection as MD17sys.path.insert(0, os.path.join(HERE, "..", "digit-uniform-bound"))18import ubound as UB1920C1 = 2.0 / math.pi21G0 = 0.962522922C0 = 0.972324# THE CHAIN AGAINST 1/52526def root_z(q, m=1):27 L = math.log(q)28 f = lambda z: (z - m) * (z - 1) ** 2 - m * (C1 * L * z + G0 * (z - 1) + C1 * (z - 1) ** 2 / (q * z - 1))29 lo, hi = float(m), float(m) + 10.030 while f(hi) < 0.0:31 hi *= 2.032 for _ in range(200):33 mid = 0.5 * (lo + hi)34 lo, hi = (mid, hi) if f(mid) < 0.0 else (lo, mid)35 return hi3637def clears(q, m=1, e=0.2):38 return root_z(q, m) < q ** e * (1 - m / q)3940def margin(q, m=1):41 c1 = 2 / iv.pi42 w = iv.mpf(q) ** (iv.mpf(1) / 5) * (1 - iv.mpf(m) / q)43 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))4445def cap_gap(q):46 c1 = 2 / iv.pi47 return iv.mpf(q) ** (iv.mpf(1) / 5) * (1 - iv.mpf(1) / q) - 1 - iv.sqrt(2 * c1 * iv.log(q) + iv.mpf(C0))4849def slope_gap(q):50 c1 = 2 / iv.pi51 return iv.mpf(q) ** (iv.mpf(1) / 5) * (1 - iv.mpf(1) / q) * iv.sqrt(2 * c1 * iv.log(q) + iv.mpf(C0)) - 5 * c15253def maynard_alpha(q):54 return math.log((q / (q - 1)) * math.log(q) + 3 * q / (q - 1)) / math.log(q)5556def maynard_s(q):57 s = 058 while math.log((1 + (2 + s + 1) / math.log(q)) * (q / (q - s - 1)) * math.log(q)) / math.log(q) < 0.2:59 s += 160 return s6162def w_budget(q):63 m = 064 while MD.pb_step3(q, m + 1) < (q - m - 1) * q ** -0.8:65 m += 166 return m6768def chain_budget(q):69 m = 070 while clears(q, m + 1):71 m += 172 return m7374def wall():75 t0 = time.time()76 iv.prec = 12077 gam = (2 / iv.pi) * (iv.euler + iv.log(8 / iv.pi))78 assert float(gam.b) <= G079 print("THE DIGIT-UNIFORM CHAIN AGAINST 1/5 at one excluded digit: z < base^(1/5)(1 - 1/base) gives alpha_1 < 1/5 at every digit")80 qu = next(q for q in range(3, 10 ** 5) if all(clears(r) for r in range(q, q + 3000)))81 qt = next(q for q in range(qu, 10 ** 5) if float(cap_gap(q).a) > 0.0)82 worst = min((float(margin(q).a), q) for q in range(qu, qt + 1))83 assert worst[0] > 0.0 and float(margin(qu - 1).b) < 0.084 assert float(cap_gap(qt).a) > 0.0 and float(slope_gap(100).a) > 0.085 z = root_z(qu)86 a1 = math.ceil(math.log(z * qu / (qu - 1)) / math.log(qu) * 1e6) / 1e687 print(" certified at 120 bits on [%d, %d]: tightest margin %.4e at base %d; the chain fails at base %d, margin %.4e"88 % (qu, qt, worst[0], worst[1], qu - 1, float(margin(qu - 1).b)))89 print(" from base %d the cap 1 + sqrt(2 (2/pi) log base + 0.97) sits below base^(1/5)(1 - 1/base), gap %.4e there,"90 " and the gap grows from base 100 on (slope test %.4f > 0)" % (qt, float(cap_gap(qt).a), float(slope_gap(100).a)))91 gap = 0.2 - math.log(z * qu / (qu - 1)) / math.log(qu)92 e = math.floor(math.log10(gap))93 gap = math.floor(gap / 10 ** e * 100) / 10094 print(" so the wall is %d: z = %.6f against %.6f and alpha_1 < %.6f at base %d, 1/5 - alpha_1 >= %.2fe%03d"95 % (qu, z, qu ** 0.2 * (1 - 1 / qu), a1, qu, gap, e))96 print()97 print("MAYNARD'S WRITTEN CONSTANT, alpha_q <= log((q/(q-1)) log q + 3q/(q-1))/log q, his Section 8")98 qm = next(q for q in range(2, 10 ** 7) if maynard_alpha(q) < 0.2 and all(maynard_alpha(r) < 0.2 for r in (q + 1, 2 * q, 10 * q)))99 print(" least q with alpha_q < 1/5: %d (alpha_q %.9f there, %.9f at %d); at q = 2000001 alpha_q = %.6f < 0.198"100 % (qm, maynard_alpha(qm), maynard_alpha(qm - 1), qm - 1, maynard_alpha(2000001)))101 assert maynard_alpha(2000001) < 0.198102 print()103 print("THE MISSING-DIGIT BUDGET at base 10^7 and 10^8: Maynard's C_(q,s) = 1 + (2+s)/log q, the wall condition (W), the chain")104 for q in (10 ** 7, 10 ** 8):105 print(" base %d: Maynard s <= %d, (W) m <= %d, the chain m <= %d" % (q, maynard_s(q), w_budget(q), chain_budget(q)))106 assert maynard_s(10 ** 8) >= 10107 print(" the chain's own wall at two and three excluded digits: %d and %d"108 % tuple(next(q for q in range(3, 10 ** 5) if all(clears(r, m) for r in range(q, q + 3000))) for m in (2, 3)))109 print("runtime %.1f s" % (time.time() - t0))110111def window():112 t0 = time.time()113 print("THE DIGIT-UNIFORM WINDOW AGAINST 1/5 at two window digits, one certificate per base for every excluded digit at once")114 rows = {}115 q = 583116 while True:117 rows[q] = UB.uniform_alpha(q, 2)[0]118 if rows[q] >= 0.2:119 break120 q -= 1121 lo = q + 1122 assert all(rows[r] < 0.2 for r in range(lo, 584))123 print(" alpha_1 < 1/5 certified at every base %d..583, the largest bound %.6f at base %d; base %d reads %.6f and fails"124 % (lo, max(rows[r] for r in range(lo, 584)), max(range(lo, 584), key=lambda r: rows[r]), q, rows[q]))125 print(" with the chain from 584 the one-missing-digit wall of the dissection is %d" % lo)126 print("runtime %.1f s" % (time.time() - t0))127128# THE WALL AT EVERY NUMBER OF EXCLUDED DIGITS129130def cap_gap_m(q, m):131 c1 = 2 / iv.pi132 return q ** (iv.mpf(1) / 5) * (1 - iv.mpf(m) / q) - m - iv.sqrt(m * (c1 * iv.log(q) + iv.mpf(C0)))133134def slope_m(q, m):135 c1 = 2 / iv.pi136 return q ** (iv.mpf(1) / 5) * iv.sqrt(c1 * iv.log(q) + iv.mpf(C0)) - iv.mpf(5) / 2 * iv.sqrt(iv.mpf(m)) * c1137138def w_exact(q, m):139 n = -(-(q - 2) // 2)140 x = iv.mpf(q)141 phq = 4 / iv.pi + (2 / iv.pi) * (iv.log(n) + iv.euler + iv.mpf(1) / (2 * n)) + (1 - 2 / iv.pi) * (x - 2) / x + iv.mpf("0.727") / x142 return (x - m) * x ** (-iv.mpf(4) / 5) - iv.sqrt(iv.mpf(m)) - phq143144def w_smooth(q, m):145 x = iv.mpf(q)146 s = 4 / iv.pi + (2 / iv.pi) * (iv.log(x / 2) + iv.euler + 1 / (x - 2)) + (1 - 2 / iv.pi) + iv.mpf("0.727") / x147 return (x - m) * x ** (-iv.mpf(4) / 5) - iv.sqrt(iv.mpf(m)) - s148149def certify(fn, lo, hi):150 stack, low = [(lo, hi)], None151 while stack:152 a, b = stack.pop()153 if b - a <= 4:154 for q in range(a, b + 1):155 v = float(fn(iv.mpf(q)).a)156 if v <= 0.0:157 return None158 if low is None or v < low[0]:159 low = (v, q)160 continue161 if float(fn(iv.mpf([a, b])).a) > 0.0:162 continue163 mid = (a + b) // 2164 stack += [(a, mid), (mid + 1, b)]165 return low166167def least(pred, lo):168 hi = lo169 while not pred(hi):170 hi *= 2171 while not all(pred(q) for q in range(hi, hi + 50)):172 hi *= 2173 lo = hi // 2174 while hi - lo > 1:175 mid = (lo + hi) // 2176 lo, hi = (lo, mid) if all(pred(q) for q in range(mid, mid + 50)) else (mid, hi)177 return hi178179def chain_wall(m):180 q0 = least(lambda q: clears(q, m), 86)181 while float(margin(q0 - 1, m).a) > 0.0:182 q0 -= 1183 assert float(margin(q0 - 1, m).b) < 0.0184 if m == 1:185 qt = next(q for q in range(q0, 10 ** 5) if float(cap_gap(q).a) > 0.0)186 assert float(slope_gap(100).a) > 0.0187 else:188 qt = least(lambda q: float(cap_gap_m(iv.mpf(q), m).a) > 0.0, q0)189 assert float(cap_gap_m(iv.mpf(qt), m).a) > 0.0 and float(slope_m(iv.mpf(qt), m).a) > 0.0190 low = certify(lambda q: margin(q, m), q0, qt)191 assert low is not None192 z = root_z(q0, m)193 gap = 0.2 - math.log(z * q0 / (q0 - m)) / math.log(q0)194 below = (q0 - 1) ** 0.2 * (1 - m / (q0 - 1)) - root_z(q0 - 1, m)195 return q0, qt, low, gap, below196197def w_wall(m):198 qs = least(lambda q: float(w_smooth(q, m).a) > 0.0, 327)199 assert qs >= 327 and float(w_smooth(qs, m).a) > 0.0200 q = qs - 1201 while float(w_exact(q, m).a) > 0.0:202 q -= 1203 assert float(w_exact(q, m).b) < 0.0204 at = w_exact(q + 1, m)205 wq = (q + 1 - m) * (q + 1) ** -0.8206 return q + 1, qs, float(at.a), math.log(wq / (wq - float(at.a))) / math.log(q + 1), float(w_exact(q, m).b)207208def walls():209 t0 = time.time()210 iv.prec = 120211 assert float(((2 / iv.pi) * (iv.euler + iv.log(8 / iv.pi))).b) <= G0212 print("THE WALL PER NUMBER m OF EXCLUDED DIGITS: the least base from which each certificate proves alpha_1 < 1/5 at every choice of m digits")213 print("chain: (z - m)(z - 1)^2 = m((2/pi) log(base) z + gamma'(z - 1) + (2/pi)(z - 1)^2/(base z - 1)) against z < base^(1/5)(1 - m/base)")214 print("(W): sqrt(m) + Phi_base/base < base^(1/5)(1 - m/base), Phi_base the step 3 constant of mobius")215 print(" m chain wall certified to root-equation margin 1/5 - alpha_1 w - z one below (W) wall smooth from w - PB at wall 1/5 - alpha_1 w - PB one below better")216 rows, thin, miss = [], [], []217 for m in range(1, 14):218 c, qt, low, gap, below = chain_wall(m)219 w, qs, wat, wgap, wbelow = w_wall(m)220 e = math.floor(math.log10(gap))221 g = math.floor(gap / 10 ** e * 100) / 100222 best = "chain" if c < w else "(W)"223 rows.append((m, c, w))224 ew = math.floor(math.log10(wgap))225 gw = math.floor(wgap / 10 ** ew * 100) / 100226 thin.append((gap, "chain", m, c))227 thin.append((wgap, "(W)", m, w))228 miss.append((-below, "chain", m, c - 1))229 miss.append((-wbelow, "(W)", m, w - 1))230 print("%2d %12d %14d %22.4e %12.2fe%03d %17.3e %12d %13d %16.4e %12.2fe%03d %18.3e %s"231 % (m, c, qt, low[0], g, e, below, w, qs, wat, gw, ew, wbelow, best))232 cross = next(m for m, c, w in rows if w < c)233 assert all(c < w for m, c, w in rows if m < cross) and all(w < c for m, c, w in rows if m >= cross)234 print("the chain is the better certificate at m <= %d and (W) at %d <= m <= 13" % (cross - 1, cross))235 t, u = min(thin), min(miss)236 et = math.floor(math.log10(t[0]))237 print("the thinnest pass over all 26 walls: 1/5 - alpha_1 >= %.2fe%03d, %s at m = %d, base %d; the closest failure one below: %.3e, %s at m = %d, base %d"238 % (math.floor(t[0] / 10 ** et * 100) / 100, et, t[1], t[2], t[3], -u[0], u[1], u[2], u[3]))239 h = w_smooth(iv.mpf(14) ** 5, 14)240 assert float(h.a) > 0.0241 print("from m = 14 on: the smooth (W) gap at base m^5 is at least %.4f at m = 14 and grows in m, so the (W) wall is below m^5,"242 " while the chain needs z < base^(1/5) with z > m, so its wall is above m^5" % (math.floor(float(h.a) * 1e4) / 1e4))243 print("runtime %.1f s" % (time.time() - t0))244245# THE SINGULAR SERIES246247def phi(n):248 r, m, p = n, n, 2249 while p * p <= m:250 if m % p == 0:251 r -= r // p252 while m % p == 0:253 m //= p254 p += 1255 return r - r // m if m > 1 else r256257def mob(n):258 k, p = 1, 2259 while p * p <= n:260 if n % p == 0:261 n //= p262 if n % p == 0:263 return 0264 k = -k265 p += 1266 return -k if n > 1 else k267268def ramanujan(d, f):269 return sum(mob(d // e) * e for e in range(1, d + 1) if d % e == 0 and f % e == 0)270271def kappa(q, F):272 return Fraction(q, phi(q)) * Fraction(sum(1 for f in F if math.gcd(f, q) == 1), len(F))273274def principal(q, F):275 return sum(Fraction(mob(d) * ramanujan(d, f), phi(d)) for d in range(1, q + 1) if q % d == 0 for f in F) / len(F)276277def series():278 t0 = time.time()279 print("THE MAIN TERM of C2: sum over squarefree d | base of mu(d)/phi(d) sum_((l,d)=1) hat F_k(l/d), divided by fill^k,")280 print("against kappa_F = (base/phi(base)) #{f in F : (f, base) = 1}/fill, in exact rationals")281 n = 0282 for q in range(3, 31):283 for m in (1, 2):284 for E in itertools.combinations(range(q), m):285 F = [v for v in range(q) if v not in E]286 assert principal(q, F) == kappa(q, F)287 n += 1288 print(" equal at all %d sets missing one or two digits of every base 3..30" % n)289 for q, E in ((10, (5,)), (10, (1,)), (10, (0, 5)), (10, (1, 3)), (12, (0, 6)), (30, (0, 15, 29))):290 F = [v for v in range(q) if v not in E]291 s1 = sum(1 for b in E if math.gcd(b, q) == 1)292 v1 = Fraction(q * (phi(q) - s1), (q - 1) * phi(q))293 print(" base %2d missing %-12s kappa_F = %-8s = %.6f; the printed q(phi(q) - s')/((q-1) phi(q)) = %.6f"294 % (q, str(E), str(kappa(q, F)), float(kappa(q, F)), float(v1)))295 print("runtime %.1f s" % (time.time() - t0))296297# THE REGIONS FOR LAMBDA298299def vonmangoldt(N):300 s = np.ones(N + 1, dtype=bool)301 s[:2] = False302 for p in range(2, math.isqrt(N) + 1):303 if s[p]:304 s[p * p:: p] = False305 ps = np.nonzero(s)[0]306 lam = np.zeros(N + 1)307 lam[ps] = np.log(ps)308 for p in ps[ps <= math.isqrt(N)]:309 pk = int(p) * int(p)310 while pk <= N:311 lam[pk] = math.log(p)312 pk *= int(p)313 return lam314315def regions_one(q, e0, k, Z):316 F = [v for v in range(q) if v != e0]317 fill = len(F)318 y = q ** k319 D = np.zeros(y)320 D[MD.strings(q, F, k)] = 1.0321 lam = vonmangoldt(y)[:y]322 hatF = np.conj(np.fft.fft(D))323 Sneg = np.fft.fft(lam)324 exact = float(lam[D > 0].sum())325 a = np.arange(y, dtype=np.int64)326 Q = int(y ** 0.6)327 l, d = MD.convergent(a, y, Q)328 h = np.abs(a * d - l * y)329 sm = np.array([MD.smooth(int(v), q) for v in range(Q + 1)])330 A = d >= y ** 0.4331 C = (~A) & (d < Z) & (h < Z)332 B = (~A) & (~C)333 C2 = C & sm[d]334 C1 = C & (~sm[d])335 term = hatF * Sneg / y336 parts = {nm: term[msk].sum().real for nm, msk in (("A", A), ("B", B), ("C1", C1), ("C2", C2))}337 tot = sum(parts.values())338 assert abs(tot - exact) < 1e-6 * y339 kap = float(kappa(q, F))340 h0 = C2 & (h == 0)341 mus = np.array([mob(int(v)) for v in d[h0]], dtype=float)342 phs = np.array([phi(int(v)) for v in d[h0]], dtype=float)343 pred = (hatF[h0] * mus / phs).sum().real344 assert abs(pred - kap * fill ** k) < 1e-6 * fill ** k345 print("base %d, excluded %d, level %d, y = %d, Z = %d, kappa_F = %.6f" % (q, e0, k, y, Z, kap))346 print(" exact sum of Lambda over the strings %.4f = %.6f kappa_F fill^k; the four regions sum to it, difference %.2e"347 % (exact, exact / (kap * fill ** k), abs(tot - exact)))348 print(" the principal characters at the C2 points with h = 0, y mu(d)/phi(d) in place of S: %.6f kappa_F fill^k" % (pred / (kap * fill ** k)))349 print(" C2 read on the grid at h = 0: %.6f kappa_F fill^k" % (term[h0].sum().real / (kap * fill ** k)))350 for nm in ("A", "B", "C1", "C2"):351 msk = {"A": A, "B": B, "C1": C1, "C2": C2}[nm]352 print(" region %-2s points %8d contribution %+.6f kappa_F fill^k" % (nm, msk.sum(), parts[nm] / (kap * fill ** k)))353354def regions():355 t0 = time.time()356 print("THE DISSECTION FOR LAMBDA at small base and level, every grid point a mod y, the regions of mobius")357 print()358 regions_one(10, 5, 6, 16)359 regions_one(10, 1, 6, 16)360 regions_one(5, 2, 9, 12)361 print("runtime %.1f s" % (time.time() - t0))362363# THE COUNT364365DESIGNS = [(10, (5,)), (10, (0,)), (10, (1,)), (10, (9,)), (10, (0, 5)), (7, (3,)), (3, (1,)), (5, (1, 3))]366367def count_upto(x, q, F):368 ds = []369 v = x370 while v:371 ds.append(v % q)372 v //= q373 ds = ds[::-1]374 L, fill, lead = len(ds), len(F), sum(1 for f in F if f > 0)375 tot = sum(lead * fill ** (j - 1) for j in range(1, L))376 for i, g in enumerate(ds):377 tot += sum(1 for f in F if f < g and (i > 0 or f > 0)) * fill ** (L - 1 - i)378 if g not in F:379 return tot380 return tot + 1381382def member(v, q, E):383 v = np.asarray(v, dtype=np.int64).copy()384 ok = np.ones(v.shape, dtype=bool)385 while (v > 0).any():386 ok &= ~((v > 0) & np.isin(v % q, E))387 v //= q388 return ok389390def census(N):391 s = np.ones(N + 1, dtype=bool)392 s[:2] = False393 for p in range(2, math.isqrt(N) + 1):394 if s[p]:395 s[p * p:: p] = False396 ps = np.nonzero(s)[0].astype(np.int64)397 del s398 pw, pl = [], []399 for p in ps[ps <= math.isqrt(N)]:400 pk = int(p) * int(p)401 while pk <= N:402 pw.append(pk)403 pl.append(math.log(p))404 pk *= int(p)405 pw = np.array(pw, dtype=np.int64)406 pl = np.array(pl)407 o = np.argsort(pw)408 pw, pl = pw[o], pl[o]409 rng = random.Random(2027)410 xs = sorted(rng.randrange(N // 100, N) for _ in range(12))411 out = []412 for q, E in DESIGNS:413 F = [v for v in range(q) if v not in E]414 mp = member(ps, q, E)415 P = ps[mp]416 cp = np.concatenate([[0.0], np.cumsum(np.log(P))])417 mw = member(pw, q, E)418 W = pw[mw]419 cw = np.concatenate([[0.0], np.cumsum(pl[mw])])420 pts = [q ** j for j in range(1, 40) if q ** j <= N] + xs421 rows = []422 for x in pts:423 psi = cp[np.searchsorted(P, x, side="right")] + cw[np.searchsorted(W, x, side="right")]424 rows.append([x, float(psi), count_upto(x, q, F)])425 out.append({"base": q, "excluded": list(E), "rows": rows})426 return out427428def count():429 t0 = time.time()430 N = 10 ** 8431 path = os.path.join(DATA_DIR, "count-%d.json" % N)432 if os.path.exists(path):433 with open(path) as fh:434 data = json.load(fh)435 src = "cached"436 else:437 data = census(N)438 os.makedirs(DATA_DIR, exist_ok=True)439 with open(path, "w") as fh:440 json.dump(data, fh)441 src = "computed"442 for q, E in DESIGNS:443 F = [v for v in range(q) if v not in E]444 cum = np.cumsum(member(np.arange(1, 10 ** 5 + 1), q, E))445 assert all(count_upto(x, q, F) == cum[x - 1] for x in (1, 9, 10, 11, 99, 100, 12345, 99999, 10 ** 5))446 print("PRIMES ON A DESIGN against the principal-character main term, sum_(n <= x, n in S_F) Lambda(n) / (kappa_F A_F(x)), x <= %d (%s)" % (N, src))447 top, seeded = [], []448 for rec in data:449 q, E = rec["base"], tuple(rec["excluded"])450 F = [v for v in range(q) if v not in E]451 kap = float(kappa(q, F))452 pw = [r for r in rec["rows"] if round(math.log(r[0], q), 9).is_integer()]453 rd = [r for r in rec["rows"] if r not in pw]454 tail = ", ".join("%.4f" % (r[1] / (kap * r[2])) for r in pw[-4:]) if kap else "kappa_F = 0"455 rr = [r[1] / (kap * r[2]) for r in rd] if kap else [0.0]456 print(" base %2d missing %-7s kappa_F %.6f at the last four powers of the base: %s; at 12 seeded x: %.4f..%.4f"457 % (q, str(E), kap, tail, min(rr), max(rr)))458 if kap and any(f + 1 in F for f in F):459 top.append(pw[-1][1] / (kap * pw[-1][2]))460 seeded.extend(rr)461 if not any(f + 1 in F for f in F):462 print(" outside (E): psi_F(%d) = %.4f against kappa_F A_F = %.1f" % (pw[-1][0], pw[-1][1], kap * pw[-1][2]))463 print(" the sets with two consecutive digits: at the largest power of the base %.4f..%.4f, at the seeded x %.4f..%.4f"464 % (min(top), max(top), min(seeded), max(seeded)))465 print("runtime %.1f s" % (time.time() - t0))466467if __name__ == "__main__":468 verb = sys.argv[1] if len(sys.argv) > 1 else "wall"469 {"wall": wall, "walls": walls, "window": window, "series": series, "regions": regions, "count": count}[verb]()