pairing.py
47.4 kB · python · 1121 lines
1import os2import sys34sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "mrly-euler"))56import euler78# FAMILIES910def full(q):11 return (1 << q) - 11213def miss(q, d):14 return full(q) & ~(1 << d)1516FAMILIES = [(3, 0b011, 14), (4, 0b0011, 11), (5, 0b00011, 9), (10, miss(10, 9), 6)]1718CONTROLS = [(2, 0b11, 20), (3, 0b111, 13)]1920THETAMAX = {(3, 0b011): 0.4869, (4, 0b0011): 0.4819, (5, 0b00011): 0.4633,21 (10, miss(10, 9)): 0.4871}2223def digs(q, F):24 return [d for d in range(q) if (F >> d) & 1]2526def alpha_of(q, F):27 from math import log28 return log(len(digs(q, F))) / log(q)2930def strings(q, F, L):31 import numpy32 out = numpy.zeros(1, dtype=numpy.int64)33 d = numpy.array(digs(q, F), dtype=numpy.int64)34 for i in range(L):35 out = (out[:, None] + d[None, :] * (q ** i)).ravel()36 out.sort()37 return out3839def indicator(q, F, L):40 import numpy41 v = numpy.zeros(q ** L, dtype=numpy.float64)42 v[strings(q, F, L)] = 1.043 return v4445# MOBIUS4647def mu_upto(N):48 import numpy49 mu = numpy.ones(N + 1, dtype=numpy.int8)50 mu[0] = 051 sieve = numpy.ones(N + 1, dtype=bool)52 sieve[:2] = False53 for p in range(2, int(N ** 0.5) + 1):54 if sieve[p]:55 sieve[p * p:: p] = False56 primes = numpy.nonzero(sieve)[0]57 for p in primes:58 mu[p:: p] = -mu[p:: p]59 for p in primes:60 if p * p > N:61 break62 mu[p * p:: p * p] = 063 return mu6465def mu_of(vals):66 import numpy67 rem = numpy.array(vals, dtype=numpy.int64)68 mu = numpy.ones(rem.shape, dtype=numpy.int8)69 top = int(rem.max())70 lim = int(top ** 0.5) + 171 sieve = numpy.ones(lim + 1, dtype=bool)72 sieve[:2] = False73 for p in range(2, int(lim ** 0.5) + 1):74 if sieve[p]:75 sieve[p * p:: p] = False76 for p in numpy.nonzero(sieve)[0]:77 d = rem % p == 078 if not d.any():79 continue80 rem[d] //= p81 mu[d] = -mu[d]82 again = d & (rem % p == 0)83 if again.any():84 mu[again] = 085 while again.any():86 rem[again] //= p87 again = again & (rem % p == 0)88 mu[rem > 1] = -mu[rem > 1]89 return mu9091VERBS = {}9293# SPLIT9495def levels(q, L):96 import numpy97 lev = numpy.full(q ** L, L, dtype=numpy.int16)98 for j in range(L - 1, -1, -1):99 lev[:: q ** (L - j)] = j100 return lev101102def split_one(q, F, L):103 import numpy104 from math import log105 N = q ** L106 k = len(digs(q, F))107 a = alpha_of(q, F)108 ind = indicator(q, F, L)109 mu = mu_upto(N - 1).astype(numpy.float64)[:N]110 direct = float(numpy.dot(ind, mu))111 Gh = numpy.fft.fft(ind)112 del ind113 Sh = numpy.fft.fft(mu)114 del mu115 paired = float((numpy.conj(Gh) * Sh).sum().real) / N116 zero = float((numpy.conj(Gh[0]) * Sh[0]).real) / N117 e2S = float((numpy.abs(Sh) ** 2).sum()) / N118 del Sh119 absG = numpy.abs(Gh)120 del Gh121 e2G = float((absG ** 2).sum()) / N122 lev = levels(q, L)123 w = numpy.bincount(lev, weights=absG, minlength=L + 1)124 del absG, lev125 C = float(w.sum())126 lq = log(q)127 ell1 = log(C / N) / (L * lq)128 cand = [((log(w[j]) / lq - L) / L + min(0.5 + j / (2.0 * L), 0.75), j)129 for j in range(L + 1) if w[j] > 0]130 best, arg = max(cand)131 return dict(q=q, F=F, L=L, k=k, alpha=a, direct=direct, paired=paired,132 e2G=e2G, e2S=e2S, zero=zero, w=w, C=C, ell1=ell1,133 unif=ell1 + 0.75, best=best, arg=arg,134 top=float(w[L]) / C, hi=float(w[(L + 1) // 2:].sum()) / C)135136def split():137 import numpy138 print("THE POSITION PAIRING ON THE GRID, EXACT")139 print("M_F(q^L) = q^(-L) sum_a G_L(a/q^L) S_L(a/q^L), both sides finite")140 print("q F L M_F(q^L) pairing int|G_L|^2 k^L int|S_L|^2 a=0 term")141 rows = []142 for q, F, L in FAMILIES + CONTROLS:143 r = split_one(q, F, L)144 rows.append(r)145 print("%2d %-8s %2d %10.1f %12.4f %10.0f %8d %11.0f %9.5f"146 % (q, euler.show(F, q), L, r["direct"], r["paired"], r["e2G"],147 r["k"] ** L, r["e2S"], r["zero"]))148 print()149 print("THE LEVEL PROFILE of the l^1 mass, C_L = sum_j k^(L-j) c_j")150 print("share of C_L at level j = L (primitive denominator q^L), at levels j >= L/2,")151 print("and the proved floor m/q for the top level")152 print("q F L C_L/q^L top share floor m/q share j >= L/2")153 for r in rows:154 m = r["q"] - r["k"]155 print("%2d %-8s %2d %12.4f %10.6f %10.6f %14.6f"156 % (r["q"], euler.show(r["F"], r["q"]), r["L"], r["C"] / r["q"] ** r["L"],157 r["top"], m / r["q"], r["hi"]))158 print()159 print("THE COST, exponents in log_q against the trivial alpha")160 print("cs = Cauchy-Schwarz (alpha+1)/2, unif = l^1 times the uniform GRH max x^(3/4),")161 print("den = the per-denominator split, GRH max x^(1/2) q^(j/2) at level j")162 print("mob = alpha times the thetamax of mobius.md, the measured exponent of M_F itself")163 print("q F L alpha alpha/2 mob cs unif den C_L/w_j q/m")164 for r in rows:165 m = r["q"] - r["k"]166 tm = THETAMAX.get((r["q"], r["F"]))167 print("%2d %-8s %2d %8.6f %8.6f %8s %8.6f %8.6f %8.6f %8.4f %8.4f"168 % (r["q"], euler.show(r["F"], r["q"]), r["L"], r["alpha"], r["alpha"] / 2,169 ("%.6f" % (r["alpha"] * tm)) if tm else "-",170 (r["alpha"] + 1) / 2, r["unif"], r["best"],171 r["C"] / r["w"][r["arg"]], r["q"] / m if m else float("inf")))172 print()173 print("THE LADDER IN L at base 3 F = {0,1}: the per-denominator gain is a constant")174 print("unif - den = log_q(C_L/w_L)/L and C_L/w_L <= q/m = 1.5 at every L")175 print(" L unif den unif-den L(unif-den) C_L/w_L top share")176 for L in range(6, 15):177 r = split_one(3, 0b011, L)178 d = r["unif"] - r["best"]179 print("%2d %8.6f %8.6f %9.6f %11.6f %8.4f %10.6f"180 % (L, r["unif"], r["best"], d, L * d, r["C"] / r["w"][r["arg"]], r["top"]))181 print()182 print("THE ONE-STEP CONSTANT read off the l^1 norms themselves")183 print("C_L = sum_(a mod q^L) |G_L(a/q^L)|, ratio C_L/C_(L-1) against the grid sup B_q(F)")184 print("q F L C_L C_L/C_(L-1) sup_t sum_r |g((t+r)/q)| top share")185 for q, F, L in FAMILIES:186 cs = []187 for l in range(1, min(L, 9) + 1):188 r = split_one(q, F, l)189 cs.append(r["C"])190 d = digs(q, F)191 t = numpy.linspace(0.0, 1.0, 200001)[:-1]192 h = numpy.zeros_like(t)193 for rr in range(q):194 z = numpy.zeros_like(t, dtype=numpy.complex128)195 for dd in d:196 z += numpy.exp(2j * numpy.pi * dd * (t + rr) / q)197 h += numpy.abs(z)198 for l in range(len(cs) - 2, len(cs)):199 print("%2d %-8s %2d %14.6f %15.9f %25.9f %10.6f"200 % (q, euler.show(F, q), l + 1, cs[l], cs[l] / cs[l - 1], h.max(),201 split_one(q, F, l + 1)["top"]))202 print()203 print("THE PRINCIPAL FIBRE carries almost none of the design meter")204 print("the a = 0 term is q^(-L) k^L M(q^L), exponent alpha - 1/2 under RH,")205 print("against the conjectured alpha/2 for M_F itself: alpha - 1/2 < alpha/2 iff alpha < 1")206 print("q F L M_F(q^L) a=0 term share alpha-1/2 alpha/2")207 for r in rows:208 sh = r["zero"] / r["direct"] if r["direct"] else float("nan")209 print("%2d %-8s %2d %10.1f %10.5f %10.6f %10.6f %8.6f"210 % (r["q"], euler.show(r["F"], r["q"]), r["L"], r["direct"], r["zero"],211 sh, r["alpha"] - 0.5, r["alpha"] / 2))212213VERBS["split"] = split214215# GLUE216217GLUE = [(3, 0b011, 18), (4, 0b0011, 14), (5, 0b00011, 12), (10, miss(10, 9), 7)]218219GLUE_CONTROLS = [(2, 0b11, 20), (3, 0b111, 12)]220221def coeffs(q, F, L):222 import numpy223 N = q ** L224 sf = strings(q, F, L)225 sf = sf[sf >= 1]226 mu = mu_of(sf).astype(numpy.int64)227 c = numpy.zeros(N + 1, dtype=numpy.int32)228 for e, m in zip(sf, mu):229 if m == 0:230 continue231 d = sf[: numpy.searchsorted(sf, N // e, "right")]232 c[d * e] += m233 return c234235PHASES = [0.0, 0.25, 0.5, 0.75]236237def partials(q, F, L):238 import numpy239 sf = strings(q, F, L)240 sf = sf[sf >= 1]241 mu = mu_of(sf).astype(numpy.int64)242 out = []243 for l in range(1, L):244 row = []245 for c in PHASES:246 x = int(float(q) ** (l + c))247 cut = numpy.searchsorted(sf, x, "right")248 a = numpy.searchsorted(sf, x // sf[:cut], "right")249 row.append((x, int((mu[:cut] * a).sum())))250 out.append(row)251 return out, sf, mu252253def glue():254 import numpy255 from math import log256 print("THE COEFFICIENTS OF ZETA_F M_F, c_F(n) = sum_(d e = n, d, e in S_F) mu(e)")257 print("q F L n<=q^L nonzero first n>1 c_F(n) max|c_F| at n")258 for q, F, L in [(3, 0b011, 12), (4, 0b0011, 9), (5, 0b00011, 8), (10, miss(10, 9), 6)]:259 c = coeffs(q, F, L)260 nz = numpy.nonzero(c[2:])[0] + 2261 am = int(numpy.abs(c).argmax())262 print("%2d %-8s %2d %10d %10d %10d %7d %9d %5d"263 % (q, euler.show(F, q), L, q ** L, len(nz) + (1 if c[1] else 0),264 int(nz[0]), int(c[nz[0]]), int(abs(c[am])), am))265 print()266 print("THE PARTIAL SUMS P(x) = sum_(n<=x) c_F(n) = sum_(e in S_F) mu(e) A_F(x/e)")267 print("the abscissa of D_F = zeta_F M_F - 1 is read off log|P|/log x against alpha")268 for q, F, L in GLUE + GLUE_CONTROLS:269 a = alpha_of(q, F)270 ps, sf, mu = partials(q, F, L)271 print("q = %d F = %s alpha = %.6f alpha/2 = %.6f" % (q, euler.show(F, q), a, a / 2))272 print(" L P(x) at log_q x = L, L+1/4, L+1/2, L+3/4 P(x)/x^alpha at the same four")273 for l in range(1, L):274 row = ps[l - 1]275 print(" %2d %10d %10d %10d %10d %9.6f %9.6f %9.6f %9.6f"276 % ((l,) + tuple(v for _, v in row)277 + tuple(v / float(x) ** a for x, v in row)))278 print()279 print("THE LIMIT TEST: sigma_c(D_F) < alpha forces M_F(sigma) -> 0 as sigma -> alpha+")280 print("M_F(sigma) = sum_(n in S_F) mu(n) n^(-sigma), summed to n <= q^L")281 print("q F L sigma-alpha M_F(sigma) tail bound q^(-L alpha/2)")282 for q, F, L in GLUE:283 a = alpha_of(q, F)284 sf = strings(q, F, L)285 sf = sf[sf >= 1]286 mu = mu_of(sf).astype(numpy.float64)287 x = sf.astype(numpy.float64)288 for eps in [0.2, 0.1, 0.05, 0.02, 0.0]:289 v = float((mu * x ** (-(a + eps))).sum())290 print("%2d %-8s %2d %12.4f %12.6f %18.2e"291 % (q, euler.show(F, q), L, eps, v, float(q) ** (-L * a / 2)))292293VERBS["glue"] = glue294295# INVERSE296297SEEDS = {(3, 0b011): ("0.720788", "28.60568"), (10, miss(10, 9)): ("1.001589", "2.7392")}298299def dirichlet_inverse(q, F, L):300 import numpy301 N = q ** L302 if F == full(q):303 sf = numpy.arange(2, N + 1, dtype=numpy.int64)304 else:305 sf = strings(q, F, L)306 sf = sf[sf >= 2]307 if (F >> 1) & 1 and (L == 0 or F & 1):308 sf = numpy.append(sf, numpy.int64(N))309 nu = numpy.zeros(N + 1, dtype=numpy.int32)310 nu[1] = 1311 lo = 1312 while lo <= N:313 hi = min(2 * lo, N + 1)314 for d in sf[: numpy.searchsorted(sf, hi - 1, "right")]:315 m0 = (lo + d - 1) // d316 m1 = (hi - 1) // d + 1317 if m1 <= m0:318 continue319 nu[m0 * d: (m1 - 1) * d + 1: d] -= nu[m0: m1]320 lo = hi321 return nu322323def refine(q, F):324 from mpmath import mp325 sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)),326 "..", "design-zeta"))327 import design_zeta328 d = design_zeta.Design(q, tuple(digs(q, F)))329 re0, im0 = SEEDS[(q, F)]330 s = mp.mpc(re0, im0)331 for _ in range(8):332 v, _ = d.zeta(s)333 h = mp.mpf(10) ** -8334 v2, _ = d.zeta(s + h)335 s = s - v * h / (v2 - v)336 v, e = d.zeta(s)337 return s, v, e, d.alpha338339def series_at(nu, q, L, sigmas):340 import numpy341 acc = [0.0] * len(sigmas)342 out = [[] for _ in sigmas]343 N = q ** L344 lo = 1345 marks = [q ** l for l in range(1, L + 1)]346 m = 0347 while lo <= N:348 hi = min(lo + (1 << 20), N + 1)349 if m < len(marks) and marks[m] < hi:350 hi = marks[m] + 1351 n = numpy.arange(lo, hi, dtype=numpy.float64)352 v = nu[lo:hi].astype(numpy.float64)353 for i, sg in enumerate(sigmas):354 acc[i] += float((v * n ** (-sg)).sum())355 if m < len(marks) and hi == marks[m] + 1:356 for i in range(len(sigmas)):357 out[i].append(acc[i])358 m += 1359 lo = hi360 return out361362def abscissa(q, F, L):363 from mpmath import mp364 sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)),365 "..", "design-zeta"))366 import design_zeta367 d = design_zeta.Design(q, tuple(digs(q, F)))368 nu = dirichlet_inverse(q, F, L)369 a = float(d.alpha)370 r = float(refine(q, F)[0].real)371 sigmas = [round(r + 0.08, 4), round(r + 0.20, 4)]372 out = series_at(nu, q, L, sigmas)373 tgt = [float(1 / d.zeta(mp.mpf(str(sg)))[0]) for sg in sigmas]374 return sigmas, out, tgt, a, r375376def ladder_stats(nu, q, L):377 import numpy378 N = q ** L379 marks = [q ** l for l in range(1, L + 1)]380 sig, mx = [], []381 run = 0382 best = 0383 lo, m = 1, 0384 while lo <= N:385 hi = min(lo + (1 << 20), N + 1)386 if m < len(marks) and marks[m] < hi:387 hi = marks[m] + 1388 c = numpy.cumsum(nu[lo:hi].astype(numpy.int64))389 c += run390 run = int(c[-1])391 numpy.abs(c, out=c)392 best = max(best, int(c.max()))393 if m < len(marks) and hi == marks[m] + 1:394 sig.append(run)395 mx.append(best)396 m += 1397 lo = hi398 return sig, mx399400def inverse():401 import numpy402 from math import log403 print("THE DIRICHLET INVERSE OF THE DESIGN INDICATOR, nu_F = 1_(S_F)^(-1)")404 print("zeta_F(s) N_F(s) = 1 exactly, the identity that replaces zeta M = 1 on a design")405 print("control: the full digit set gives nu_F = mu term for term")406 for q in [2, 3]:407 nu = dirichlet_inverse(q, full(q), 17 if q == 2 else 11)408 mu = mu_upto(len(nu) - 1)409 print(" q = %d full, n <= %d, nu_F = mu at every n: %s"410 % (q, len(nu) - 1, bool((nu[1:] == mu[1:]).all())))411 print()412 print("THE ZERO THAT FORCES THE ABSCISSA, refined on the sibling engine lab/py/design-zeta")413 print("sigma_c(N_F) >= Re rho for every zero rho of zeta_F with Re rho > alpha")414 print("q F zero rho |zeta_F(rho)| bound alpha")415 from mpmath import mp416 for q, F in [(3, 0b011), (10, miss(10, 9))]:417 r, v, e, a = refine(q, F)418 print("%2d %-8s %-42s %-9s %-9s %s"419 % (q, euler.show(F, q), mp.nstr(r, 16), mp.nstr(abs(v), 4),420 mp.nstr(e, 4), mp.nstr(a, 10)))421 print()422 print("THE DESIGN MERTENS OF nu_F against the design's own mass A_F(x) and against")423 print("the rightmost zero of zeta_F, whose real part is a lower bound for the abscissa")424 for q, F, L in [(3, 0b011, 16), (4, 0b0011, 12), (5, 0b00011, 10), (10, miss(10, 9), 7)]:425 a = alpha_of(q, F)426 k = len(digs(q, F))427 nu = dirichlet_inverse(q, F, L)428 raw, mxs = ladder_stats(nu, q, L)429 rz = float(refine(q, F)[0].real) if (q, F) in SEEDS else None430 print("q = %d F = %s alpha = %.6f rightmost censused zero Re = %s"431 % (q, euler.show(F, q), a, ("%.6f" % rz) if rz else "not censused"))432 print(" the level ratio of the running maximum against q^alpha = k = %d and"433 % k)434 print(" q^(Re rho) = %s"435 % ("%.6f" % (q ** rz) if rz else "not censused"))436 print(" L sum nu_F(n) max|sum| level ratio exponent max/A_F(q^L)")437 prev = 0438 for l in range(1, L + 1):439 x = q ** l440 mx = mxs[l - 1]441 ex = log(mx) / log(x) if mx > 1 else 0.0442 print(" %2d %13d %12d %12s %10.6f %14.4f"443 % (l, raw[l - 1], mx, ("%.4f" % (mx / prev)) if prev else "-", ex,444 mx / float(k ** l)))445 prev = mx446 del nu447448 print()449 print("THE ABSCISSA TEST: partial sums of N_F(sigma) = sum nu_F(n) n^(-sigma) against")450 print("1/zeta_F(sigma) from lab/py/design-zeta, at sigma above and below Re rho;")451 print("the transport theorem, not this table, places sigma_c(N_F) at or above Re rho")452 for q, F, L in [(3, 0b011, 16), (10, miss(10, 9), 7)]:453 sigmas, out, tgt, a, r = abscissa(q, F, L)454 print("q = %d F = %s alpha = %.7f Re rho = %.7f" % (q, euler.show(F, q), a, r))455 print(" sigma 1/zeta_F partial at q^(L-2) q^(L-1) q^L error")456 for i, sg in enumerate(sigmas):457 p3, p2, p1 = out[i][-3], out[i][-2], out[i][-1]458 print(" %-8.4f %12.6f %14.6f %12.6f %12.6f %11.2e"459 % (sg, tgt[i], p3, p2, p1, abs(p1 - tgt[i])))460461VERBS["inverse"] = inverse462463# BOX464465BOXES = [(10, miss(10, 9), "1.00150", "1.00168", "2.73915", "2.73925", 1),466 (3, 0b011, "0.72074", "0.72084", "28.60563", "28.60573", 1),467 (10, miss(10, 9), "0.99900", "1.00050", "2.73810", "2.74030", 0)]468469def box():470 from mpmath import mp471 sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)),472 "..", "design-zeta"))473 import design_zeta474 print("WINDING BOXES FOR THE ZEROS OF zeta_F, argument principle on lab/py/design-zeta")475 print("winding 1 on a rectangle certifies exactly one zero inside it, so Re rho is")476 print("pinned to the box edges; winding 0 certifies no zero there")477 for q, F, x0, x1, y0, y1, want in BOXES:478 d = design_zeta.Design(q, tuple(digs(q, F)))479 f = design_zeta.Cache(d)480 w, mx = design_zeta.box_phase(f, mp.mpf(x0), mp.mpf(x1), mp.mpf(y0), mp.mpf(y1), 16)481 lo = min(abs(v) for v in f.m.values())482 print("q = %d F = %s alpha = %s"483 % (q, euler.show(F, q), mp.nstr(d.alpha, 10)))484 print(" Re in [%s, %s] Im in [%s, %s]" % (x0, x1, y0, y1))485 print(" winding %s expected %d max phase step %s evaluations %d"486 % (mp.nstr(w, 8), want, mp.nstr(mx, 4), f.n))487 print(" min |zeta_F| on the contour %s engine bound %s ratio %s"488 % (mp.nstr(lo, 4), mp.nstr(f.emax, 4), mp.nstr(lo / f.emax, 4)))489490VERBS["box"] = box491492493# ONESTEP494495GAMMA = 0.5772156649015329496497ONESTEP_FAMILIES = [(3, 2), (10, 9), (11, 0), (13, 0), (100, 0), (100, 49), (1000, 0),498 (1000, 499), (2234, 0), (2234, 1116), (3690, 0), (3690, 1844)]499500def dn(x, d=6):501 from math import floor502 return floor(x * 10.0 ** d) / 10.0 ** d503504def up(x, d=6):505 from math import ceil506 return ceil(x * 10.0 ** d) / 10.0 ** d507508def psi_q(q):509 from math import log, pi510 n = -(-(q - 2) // 2)511 return (2 * q / pi) * (log(n) + GAMMA + 1 / (2 * n)) + (1 - 2 / pi) * q512513def phi_q(q):514 from math import log, pi515 n = -(-(q - 2) // 2)516 return (4 / pi) * q + (2 * q / pi) * (log(n) + GAMMA + 1 / (2 * n)) + (1 - 2 / pi) * (q - 2) + 0.727517518def sec_e(q, e0):519 from math import cos, pi520 return 1.0 / cos(pi * (e0 - (q - 1) / 2.0) / q)521522def bound_e(q, e0):523 from math import pi524 return (4 / pi) * q + psi_q(q) + q / 2.0 - sec_e(q, e0) / 2.0525526def harmonic(n):527 return sum(1.0 / j for j in range(1, int(n) + 1))528529def psi2_q(q, exact=False):530 from math import log, pi531 p = q // 2532 if exact:533 h = harmonic(p - 1)534 else:535 h = 0.0 if p < 3 else log(p - 1) + GAMMA + 1 / (2 * (p - 1))536 if q % 2 == 0:537 return (q / pi) * (2 * h - 1 + 1 / p) + (1 - 2 / pi) * q / 2538 return (q / pi) * (2 * h - 1 + 2 / p) + (1 - 2 / pi) * (q / 2 + 1 / (2 * q))539540def bound2_e(q, e0):541 from math import pi542 return (4 / pi) * q + psi2_q(q) + q / 2.0 - sec_e(q, e0) / 2.0543544def mono_floor(hi=400, exact=False):545 from math import pi546 ok = [psi2_q(q, exact) >= (1 + pi) * q / 2 for q in range(4, hi)]547 for i in range(len(ok)):548 if all(ok[i:]):549 return 4 + i550 return None551552def h_argmax_floor(hi=80, n=2001):553 import numpy554 from math import cos, pi555 tau = numpy.linspace(0.0, 0.5, n)556 c = numpy.cos(pi * tau)557 last = 3558 for q in range(4, hi):559 for e0 in range(q):560 b = pi * (e0 - (q - 1) / 2.0) / q561 h = c * (psi2_q(q) - q / 2.0 - numpy.cos(2 * b * tau) / (2 * cos(b)))562 if int(h.argmax()) != 0:563 last = q564 return last + 1, hi - 1565566def weight_cost(q):567 import numpy568 from math import pi569 r = numpy.arange(q, dtype=numpy.float64)570 a = 1.0 / numpy.sin(pi * (r + 0.5) / q)571 return 2.0 * float(numpy.sum(1.0 / (a + 1.0)))572573def terms(q, e0, t):574 import numpy575 from math import pi, sin576 r = numpy.arange(q, dtype=numpy.float64)577 ph = (t + r) / q578 a = ((-1.0) ** r) * sin(pi * t) / numpy.sin(pi * ph)579 c = e0 - (q - 1) / 2.0580 return numpy.sqrt(a * a + 1.0 - 2.0 * a * numpy.cos(2 * pi * c * ph))581582def phase_left(q, e0, t):583 import numpy584 from math import pi, sin585 r = numpy.arange(q, dtype=numpy.float64)586 ph = (t + r) / q587 sgn = numpy.sign(((-1.0) ** r) * sin(pi * t) / numpy.sin(pi * ph))588 c = e0 - (q - 1) / 2.0589 return float((1.0 + sgn * numpy.cos(2 * pi * c * ph)).sum())590591def phase_right(q, e0, t):592 from math import cos, pi593 c = e0 - (q - 1) / 2.0594 return q + cos(2 * pi * c * (t - 0.5) / q) / cos(pi * c / q)595596def sigma_t(q, e0, t):597 return float(terms(q, e0, t).sum())598599def kernel_t(q, t):600 import numpy601 from math import pi, sin602 r = numpy.arange(q, dtype=numpy.float64)603 return float(sin(pi * t) * (1.0 / numpy.sin(pi * (t + r) / q)).sum())604605def direct_t(q, e0, t):606 import numpy607 from math import pi608 f = numpy.array([d for d in range(q) if d != e0], dtype=numpy.float64)609 tot = 0.0610 for r in range(q):611 tot += abs(complex(numpy.exp(2j * pi * f * ((t + r) / q)).sum()))612 return tot613614ONESTEP_CUT = 2001615616def seat(q, e0, n=ONESTEP_CUT):617 import numpy618 best, at = -1.0, 0.0619 for t in numpy.linspace(1e-12, 0.5, n):620 v = sigma_t(q, e0, t)621 if v > best:622 best, at = v, t623 return best, at624625def ker_seat(q, n=ONESTEP_CUT):626 import numpy627 best, at = -1.0, 0.0628 for t in numpy.linspace(1e-12, 0.5, n):629 v = kernel_t(q, t)630 if v > best:631 best, at = v, t632 return best, at633def kernel_slack(q, n=ONESTEP_CUT):634 import numpy635 from math import pi, sin636 worst, at = 1e18, 0.0637 for t in numpy.linspace(1e-9, 0.5, n):638 v = (4 / pi) * q + sin(pi * t) * psi2_q(q) - kernel_t(q, t)639 if v < worst:640 worst, at = v, t641 return worst, at642643644BLOCK = 1 << 20645646def psi_grid(q, chord):647 import numpy648 from math import pi649 if not chord:650 n = numpy.ceil((q - 2) / 2)651 return (2 * q / pi) * (numpy.log(n) + GAMMA + 1 / (2 * n)) + (1 - 2 / pi) * q652 p = numpy.floor(q / 2)653 h = numpy.log(p - 1) + GAMMA + 1 / (2 * (p - 1))654 ex = numpy.where(q % 2 == 0, 1.0 / p, 2.0 / p)655 fl = numpy.where(q % 2 == 0, q / 2, q / 2 + 1 / (2 * q))656 return (q / pi) * (2 * h - 1 + ex) + (1 - 2 / pi) * fl657658def rhs_grid(q, worst, chord):659 import numpy660 from math import pi661 c = numpy.where(q % 2 == 1, 0.0, 0.5) if worst else (q - 1) / 2662 sec = 1.0 / numpy.cos(pi * c / q)663 return (4 / pi) + psi_grid(q, chord) / q + 0.5 - sec / (2 * q)664665def proved_wall(b, worst, hi, chord=False, lo=17):666 import numpy667 first, held = None, 0668 for a in range(lo, hi + 1, BLOCK):669 q = numpy.arange(a, min(a + BLOCK, hi + 1), dtype=numpy.float64)670 ok = numpy.nonzero((q - 1) * q ** (-b) - rhs_grid(q, worst, chord) > 0)[0]671 held += int(ok.size)672 if first is None and ok.size:673 first = int(q[ok[0]])674 return (first, held)675676def cost_cross(worst, chord, lo, hi, b=0.75):677 import numpy678 q = numpy.arange(lo, hi, dtype=numpy.float64)679 alpha = numpy.log(q - 1) / numpy.log(q)680 cq = numpy.log(rhs_grid(q, worst, chord)) / numpy.log(q)681 win = (alpha - b - cq) / alpha > 1.0 / (2 * (q - 1) * numpy.log(q))682 i = numpy.nonzero(win)[0]683 if not i.size:684 return None, None685 return int(q[i[0]]), bool(win[i[0]:].all())686687def half_max(q):688 import numpy689 from math import pi690 r = numpy.arange(q, dtype=numpy.float64)691 ph = (0.5 + r) / q692 a = ((-1.0) ** r) / numpy.sin(pi * ph)693 c = numpy.arange((q - 1) // 2 + 1, dtype=numpy.float64) - (q - 1) / 2.0694 s = numpy.sqrt(a * a + 1.0 - 2.0 * a * numpy.cos(2 * pi * c[:, None] * ph[None, :])).sum(axis=1)695 i = int(numpy.argmax(s))696 return float(s[i]), i697698def ceiling_wall(b, worst, hi):699 q = hi700 while q >= 17:701 v, e0 = half_max(q) if worst else (sigma_t(q, 0, 0.5), 0)702 if not ((q - 1) * q ** (-b) > v / q):703 return q + 1, q, e0704 q -= 1705 return 17, None, None706707CEIL_HI = 2000708709DRAWS = 4000710711DRAW_SEED = 1009712713DRAW_BASES = [17, 23, 60, 101, 333, 1000, 3690]714715MONO2 = 36716717RUNGS = [("1/2", 0.75, 3690, 4000000), ("13/25", 1417.0 / 1850, 8578, 8000000),718 ("11/20", 913.0 / 1160, 33547, 40000000)]719720def onestep():721 import numpy722 from math import pi723 print("THE EXACT ONE-STEP CONSTANT AT ONE EXCLUDED DIGIT")724 print("B_q(F) = sup_t sum_(r mod q) |g_F((t+r)/q)| with F = {0..q-1} less {e_0}")725 print("reduction: |g_F((t+r)/q)| = |A_r - e(c (t+r)/q)| with A_r = (-1)^r sin(pi t)/sin(pi (t+r)/q)")726 print("and c = e_0 - (q-1)/2, so the sum is exact in O(q) at every t")727 print(" q e_0 t reduction direct sum over F")728 for q, e0, t in [(7, 0, 0.13), (7, 3, 0.5), (11, 5, 0.77), (12, 4, 0.31)]:729 print("%2d %5d %6.2f %16.12f %16.12f" % (q, e0, t, sigma_t(q, e0, t), direct_t(q, e0, t)))730 print()731 print("THE PHASE IDENTITY the sharpening runs on, exact at every q, e_0 and t")732 print("sum_r (1 + sign(A_r) cos(2 pi c (t+r)/q)) = q + cos(2 pi c (t - 1/2)/q)/cos(pi c/q) >= q + 1")733 print(" q e_0 t left right q + 1")734 for q, e0, t in [(7, 0, 0.13), (11, 5, 0.77), (100, 0, 0.5), (100, 37, 0.29)]:735 print("%4d %5d %6.2f %16.9f %16.9f %10d" % (q, e0, t, phase_left(q, e0, t), phase_right(q, e0, t), q + 1))736 print()737 print("THE SEAT: the grid scan over t in [0, 1/2], cut 1/%d, Sigma symmetric about t = 1/2" % (2 * (ONESTEP_CUT - 1)))738 print("the t -> 0 value is 2(q-1) exactly, so B_q(F) >= 2(q-1) at every q and every e_0")739 print("the seat is t = 1/2 at every family printed here except q = 11 and q = 13, where it is interior,")740 print("so t = 1/2 is where the constant is read and not where it is proved to sit")741 print(" q e_0 sup on the cut at t Sigma(1/2) 2(q-1) seat")742 for q, e0 in ONESTEP_FAMILIES:743 v, at = seat(q, e0)744 h = sigma_t(q, e0, 0.5)745 tag = "t = 1/2" if abs(at - 0.5) < 1e-9 else ("t -> 0" if at < 0.4 else "interior")746 print("%5d %5d %18.10f %9.6f %18.10f %8d %s"747 % (q, e0, dn(v, 10), at, dn(h, 10), 2 * (q - 1), tag))748 print()749 print("THE SPLIT DEFECT: what the triangle |g_F| <= |D_q| + |g_E| throws away")750 print("K_q = sup_t sum_r |D_q((t+r)/q)| is the exact kernel sup and sum_r |g_E| = q exactly at m = 1,")751 print("so the triangle bound is K_q + q and the defect is (K_q + q) - B_q(F); Phi_q is the proved kernel bound")752 print(" q e_0 Phi_q/q K_q/q B_q(F)/q defect/q Phi_q slack/q proved bound/q")753 for q, e0 in ONESTEP_FAMILIES:754 if q < 50:755 continue756 k, _ = ker_seat(q)757 b, _ = seat(q, e0)758 print("%5d %5d %10.6f %10.6f %11.6f %10.6f %13.6f %14.6f"759 % (q, e0, up(phi_q(q) / q), dn(k / q), dn(b / q), up((k + q - b) / q),760 up((phi_q(q) - k) / q), up(bound_e(q, e0) / q)))761 print()762 print("THE CHORD KERNEL BOUND, what replaces Phi_q inside step 3")763 print("1/sin x <= 1/x + (2/pi)(1 - 2/pi) x on (0, pi/2] is the chord of the convex csc x - 1/x,")764 print("and pairing r with q-1-r puts every shifted-grid argument inside (0, pi/2] at t in (0, 1/2],")765 print("so K(t) = sin(pi t) sum_r 1/sin(pi (t+r)/q) <= (4/pi) q + sin(pi t) Psi'_q with")766 print("Psi'_q = (q/pi)(2 H(P-1) - 1 + 1/P) + (1 - 2/pi) q/2 at even q, P = floor(q/2),")767 print("and (q/pi)(2 H(P-1) - 1 + 2/P) + (1 - 2/pi)(q/2 + 1/(2q)) at odd q")768 print(" q Psi_q/q Psi'_q/q (K_q - (4/pi)q)/q old gap/q new gap/q lemma slack/q at t")769 for q, e0 in ONESTEP_FAMILIES:770 if e0 != 0 or q < 50:771 continue772 k, _ = ker_seat(q)773 v, at = kernel_slack(q)774 e = (k - (4 / pi) * q) / q775 print("%5d %10.6f %11.6f %18.6f %11.6f %11.6f %14.6f %5.3f"776 % (q, up(psi_q(q) / q), up(psi2_q(q) / q), dn(e), up(psi_q(q) / q - e),777 up(psi2_q(q) / q - e), dn(v / q), at))778 print("H is the harmonic upper bound ln n + gamma + 1/(2n) of mobius.md at every Psi and Psi' here")779 print("the monotone step of the sharpening needs Psi'_q >= (1 + pi) q/2, first true at q = %d with that H"780 % mono_floor())781 print("and at q = %d with the harmonic number itself, the two floors the convention separates"782 % mono_floor(exact=True))783 f, fhi = h_argmax_floor()784 print("the hypothesis is sufficient and not necessary: the max of h(tau) sits at tau = 0 at every e_0")785 print("from q = %d up on the exhaustive scan 4..%d, so either floor has room above the true one" % (f, fhi))786 print()787 print("THE WEIGHT ROUTE, priced: what dropping the singular terms would cost")788 print("the kept weight w_r = |A_r|/(|A_r| + 1) >= s/(1 + s) >= s/2 is worth q/2 at the seat,")789 print("and dropping it by (1 + sign(A_r) cos) <= 2 costs 2 sum_r 1/(|A_r| + 1) = (2 - 4/pi) q at t = 1/2")790 for q in (1000, 3690, 20000):791 print(" q = %6d cost/q %.6f limit 2(1 - 2/pi) = %.6f net q/2 - cost %.6f q"792 % (q, up(weight_cost(q) / q), up(2 * (1 - 2 / pi)), dn(0.5 - weight_cost(q) / q)))793 print(" the net is negative at every q printed, so the route loses more than it wins")794 print()795 print("THE PROVED SHARPENING, q >= 17 and m = 1")796 print("B_q(F) <= (4/pi) q + Psi_q + q/2 - sec(pi (e_0 - (q-1)/2)/q)/2 with Psi_q = Phi_q - (4/pi) q - 0.000239,")797 print("against the step 3 bound q PB_q(1) = q + Phi_q: a saving of q/2 at every e_0 and of q/2 + q/pi at e_0 in {0, q-1}")798 print("the wall q_0(a) is the least q with (q-1) q^(-b(a)) > B_q(F)/q, scanned exhaustively from q = 17")799 print("rung a b(a) q_0 at step 3 sharpened, every e_0 sharpened, e_0 in {0, q-1} scan")800 print("held is the count of q in the scan meeting the condition, so held = hi - w + 1 is an up-set")801 for name, b, old, hi in RUNGS:802 w1, n1 = proved_wall(b, True, hi)803 w2, n2 = proved_wall(b, False, hi)804 print("%-7s %.8f %14d %22s %28s %d..%d held %d %d"805 % (name, b, old, w1, w2, 17, hi, n1, n2))806 print()807 print("THE SAME WALLS WITH THE CHORD KERNEL BOUND, q >= %d and m = 1" % MONO2)808 print("B_q(F) <= (4/pi) q + Psi'_q + q/2 - sec(pi (e_0 - (q-1)/2)/q)/2, Psi'_q = Psi_q - q/2 + 2/pi at even q")809 print("rung a b(a) q_0 at step 3 chord, every e_0 chord, e_0 in {0, q-1} scan")810 for name, b, old, hi in RUNGS:811 w1, n1 = proved_wall(b, True, hi, True, MONO2)812 w2, n2 = proved_wall(b, False, hi, True, MONO2)813 print("%-7s %.8f %14d %18s %24s %d..%d held %d %d"814 % (name, b, old, w1, w2, MONO2, hi, n1, n2))815 print("the cost-out crossing of each chord wall is a kept number of lab/rs/mertens-numerology and is not printed here")816 print()817 print("THE CEILING OF THE LEVER: the same wall computed from the measured B_q(F) itself")818 print("this is a measurement of Sigma(1/2) and not an upper bound certificate,")819 print("so it names what an exact constant could ever buy and never enters a statement")820 print("every e_0 is scanned, not the middle one, and the scan runs downward from the top of")821 print("its range, so the printed wall is the last failure plus one and the tail above it is clear")822 print("rung a b(a) ceiling, every e_0 last failure ceiling, e_0 in {0, q-1} last failure scan")823 for name, b, old, hi in RUNGS[:1]:824 w1, f1, d1 = ceiling_wall(b, True, CEIL_HI)825 w2, f2, _ = ceiling_wall(b, False, 4 * CEIL_HI)826 print("%-7s %.8f %20s %8s at e_0 = %d %14s %13s %d..%d and %d..%d"827 % (name, b, w1, f1, d1, w2, f2, 17, CEIL_HI, 17, 4 * CEIL_HI))828 print()829 print("FALSIFICATION: a family with B_q(F) above the proved bound kills the lever")830 print("the ratio B_q(F)/bound over every e_0 at q = 17..60 on the cut and at the seats above")831 worst = None832 for q in range(17, 61):833 for e0 in range(q):834 v = seat(q, e0, 401)[0]835 if worst is None or v / bound_e(q, e0) > worst[0]:836 worst = (v / bound_e(q, e0), q, e0)837 big = max((seat(q, e0)[0] / bound_e(q, e0), q, e0) for q, e0 in ONESTEP_FAMILIES if q >= 100)838 rng = numpy.random.default_rng(DRAW_SEED)839 draws = None840 for _ in range(DRAWS):841 q = int(rng.choice(DRAW_BASES))842 e0 = int(rng.integers(0, q))843 t = float(rng.uniform(1e-9, 1.0 - 1e-9))844 v = sigma_t(q, e0, t) / bound_e(q, e0)845 if draws is None or v > draws[0]:846 draws = (v, q, e0, t)847 print(" q = 17..60, every e_0: worst ratio %.6f at q = %d, e_0 = %d" % worst)848 print(" the seats above: worst ratio %.6f at q = %d, e_0 = %d" % big)849 print(" %d draws, seed %d: worst ratio %.6f at q = %d, e_0 = %d, t = %.6f"850 % (DRAWS, DRAW_SEED, draws[0], draws[1], draws[2], draws[3]))851 print(" no family reaches 1, so nothing measured contradicts the sharpening")852 print()853 print("FALSIFICATION OF THE CHORD BOUND, the same three sweeps against (4/pi) q + Psi'_q + q/2 - sec/2")854 worst = None855 for q in range(MONO2, 61):856 for e0 in range(q):857 v = seat(q, e0, 401)[0]858 if worst is None or v / bound2_e(q, e0) > worst[0]:859 worst = (v / bound2_e(q, e0), q, e0)860 big = max((seat(q, e0)[0] / bound2_e(q, e0), q, e0) for q, e0 in ONESTEP_FAMILIES if q >= 100)861 rng = numpy.random.default_rng(DRAW_SEED)862 draws = None863 for _ in range(DRAWS):864 q = int(rng.choice(DRAW_BASES))865 e0 = int(rng.integers(0, q))866 t = float(rng.uniform(1e-9, 1.0 - 1e-9))867 v = sigma_t(q, e0, t) / bound2_e(q, e0)868 if draws is None or v > draws[0]:869 draws = (v, q, e0, t)870 print(" q = %d..60, every e_0: worst ratio %.6f at q = %d, e_0 = %d" % ((MONO2,) + worst))871 print(" the seats above: worst ratio %.6f at q = %d, e_0 = %d" % big)872 print(" %d draws, seed %d: worst ratio %.6f at q = %d, e_0 = %d, t = %.6f"873 % (DRAWS, DRAW_SEED, draws[0], draws[1], draws[2], draws[3]))874 print(" no family reaches 1, so nothing measured contradicts the chord bound either")875876VERBS["onestep"] = onestep877878# PERDEN879880PERDEN_FAMILIES = [(3, [2], 12), (10, [9], 6), (101, [0], 3), (1499, [749], 2)]881882PERDEN_POINTWISE = [(10, [9], 5), (10, [9], 6)]883884PERDEN_HONEST = [(3, [2], [6, 8, 10, 12]), (10, [9], [4, 5, 6])]885886PERDEN_BRUTE = [(3, [2], 4), (3, [2], 5), (10, [9], 2)]887888PERDEN_RUNGS = [((1, 2), (3, 4), "both"), ((13, 25), (1417, 1850), "Zhang"),889 ((11, 20), (913, 1160), "Zhang"), ((4, 7), (4, 5), "both"),890 ((3, 5), (4, 5), "BH"), ((2, 3), (5, 6), "BH")]891892def hat_abs(q, E, j):893 import numpy894 n = q ** j895 v = numpy.arange(n, dtype=numpy.int64)896 out = numpy.ones(n, dtype=numpy.complex128)897 for i in range(j):898 u = (v * pow(q, i, n)) % n899 t = u.astype(numpy.float64) / n900 z = numpy.exp(2j * numpy.pi * t)901 den = numpy.where(u == 0, 1.0 + 0j, 1.0 - z)902 g = numpy.where(u == 0, complex(q, 0), (1.0 - numpy.exp(2j * numpy.pi * q * t)) / den)903 for e in E:904 g = g - z ** e905 out *= g906 return numpy.abs(out)907908def levmass(q, E, J):909 import numpy910 C, c = [], []911 for j in range(J + 1):912 h = hat_abs(q, E, j)913 C.append(float(h.sum()))914 c.append(1.0 if j == 0 else float(h[numpy.arange(q ** j) % q != 0].sum()))915 del h916 return C, c917918def smalldens(q, E, L):919 import numpy920 n = q ** L921 w = hat_abs(q, E, L)922 a = numpy.arange(n, dtype=numpy.int64)923 d = n // numpy.gcd(a, n)924 bad = ((a % q) != 0) & (d.astype(numpy.float64) <= float(n) ** 0.5)925 return int(bad.sum()), float(w[bad].sum()) / float(w.sum())926927def nu(N):928 import numpy929 from math import sqrt930 a = numpy.arange(N, dtype=numpy.int64)931 best = numpy.full(N, float(N), dtype=numpy.float64)932 for Q in range(1, int(2.0 * sqrt(N)) + 2):933 u = (a * Q) % N934 numpy.minimum(best, numpy.minimum(u, N - u).astype(numpy.float64) + Q, out=best)935 return best936937def nu_slow(N):938 from math import gcd939 out = []940 for x in range(N):941 bst = float(N)942 for Q in range(1, N + 1):943 for r in range(Q + 1):944 if gcd(r, Q) == 1 or (Q == 1 and r == 0):945 bst = min(bst, Q + abs(x * Q - r * N))946 out.append(bst)947 return out948949def honest(q, E, L, a, b):950 import numpy951 N = q ** L952 w = hat_abs(q, E, L)953 con = numpy.minimum(float(N) ** b, float(N) ** a * numpy.sqrt(nu(N)))954 return float((w * con).sum()) / (float(w.sum()) * float(N) ** b)955956def brute(q, E, j):957 from cmath import exp958 from math import pi959 n, k = q ** j, q - len(E)960 D = [d for d in range(q) if d not in E]961 vals = [0]962 for i in range(j):963 vals = [v + d * q ** i for v in vals for d in D]964 assert len(vals) == k ** j965 tot, top = 0.0, 0.0966 for a in range(n):967 z = abs(sum(exp(2j * pi * a * v / n) for v in vals))968 tot += z969 if a % q:970 top += z971 return tot, top972973def shares(q, E, C, c, L):974 from math import log, exp975 k = q - len(E)976 return [exp((L - j) * log(k) + log(c[j]) - log(C[L])) for j in range(L + 1)]977978def perden_exp(q, sh, L, a, b):979 from math import log980 lq = log(q)981 gain = [log(sh[j]) / (L * lq) - max(0.0, b - a - j / (2.0 * L)) for j in range(L + 1)]982 ratio = sum(sh[j] * q ** (-L * max(0.0, b - a - j / (2.0 * L))) for j in range(L + 1))983 jm = max(range(L + 1), key=lambda j: gain[j])984 return jm, gain[jm], ratio985986def perden():987 import numpy988 from fractions import Fraction as Fr989 from math import log990 print("THE PER-DENOMINATOR MOBIUS INPUT IN STEP 5 (T-M10)")991 print("x = q^L, S_mu(theta) = sum_(n < x) mu(n) e(n theta), M_F(q^L) = q^(-L) sum_(a mod q^L) hatF_L(a/q^L) S_mu(-a/q^L)")992 print("group a by its q-power level j: a = a' q^(L-j) with q not dividing a', so a/q^L = a'/q^j")993 print("hatF_L(a'/q^j) = k^(L-j) hatF_j(a'/q^j) because g_F(integer) = k, so the level-j mass is k^(L-j) c_j")994 print("with c_j = sum over a' mod q^j, q not dividing a', of |hatF_j(a'/q^j)|, c_0 = 1, C_j = sum_(i <= j) k^(j-i) c_i")995 print("Baker-Harman PROPOSITION p.194 eq. 6 under hypothesis (4), L(s,chi) zero-free in sigma > a for EVERY")996 print("Dirichlet character: S_mu(theta) << x^(a+eps) Q^(1/2) (1 + x|theta - r/Q|)^(1/2) at (r,Q) = 1;")997 print("the reduced denominator of a'/q^j divides q^j, so the level-j constant is AT MOST")998 print("min(x^(b(a)), x^a q^(j/2)), taking (r,Q) = the frequency itself, where the second factor is 1,")999 print("and the first the uniform THEOREM p.193; that choice of (r,Q) is the exact-frequency COROLLARY,")1000 print("and the bracket [m/q, 1] below is proved for it and for it alone. The full PROPOSITION lets ANY")1001 print("reduced r/Q serve any frequency: per-frequency constant x^a nu(a)^(1/2) with nu(a) = min_Q (Q +")1002 print("||aQ||_(q^L)), since Q(1 + x|a/q^L - r/Q|) = Q + |aQ - r q^L|. That form is measured, not proved,")1003 print("in the honest block below, and it buys no exponent either")1004 print()1005 print("|M_F(q^L)| <<_(q,eps) x^eps q^(-L) sum_(j = 0..L) k^(L-j) c_j min(x^(b(a)), x^a q^(j/2))")1006 print()1007 print("THE TOP LEVEL PAYS THE UNIFORM PRICE AT EVERY RUNG")1008 print("the level-j charge is a + j/(2L), so the level beats b(a) only below j/L = 2(b - a);")1009 print("at j = L the charge reads a + 1/2, and a + 1/2 - b(a) > 0 at every rung, exact rationals")1010 print(" a b(a) source a + 1/2 - b(a) crossing 2(b - a) beats uniform at j = L")1011 for (an, ad), (bn, bd), src in PERDEN_RUNGS:1012 a, b = Fr(an, ad), Fr(bn, bd)1013 assert a + Fr(1, 2) > b and 2 * (b - a) <= Fr(1, 2)1014 print("%6s %9s %8s %15s %19s %s"1015 % (a, b, src, a + Fr(1, 2) - b, 2 * (b - a), "yes" if a + Fr(1, 2) < b else "no"))1016 print()1017 print("the charge is an upper bound, not the pointwise truth. At composite q the reduced denominator of")1018 print("a'/q^L with q not dividing a' can sit far below q^L: at q = 10, L = 6, a' = 5^6 = 15625 gives")1019 print("15625/10^6 = 1/64, denominator 2^6 = x^0.301 and charge x^(a + 0.1505), beneath x^(3/4) at a = 1/2.")1020 print("The displayed inequality is unharmed, being an upper bound, and the mass such frequencies carry is")1021 print(" q E L top-level a' with reduced denominator <= x^(1/2) their share of C_L")1022 for q, E, L in PERDEN_POINTWISE:1023 n, share = smalldens(q, E, L)1024 print("%5d %5d %3d %40d %24.6e" % (q, E[0], L, n, share))1025 print()1026 print("THE LEVEL PROFILE: the l^1 mass is NOT at the low denominators; it decays geometrically downward")1027 print("proved: sum_(s mod q) |g_F((t+s)/q)|^2 = qk exactly and |g_F| <= k, so sum_s |g_F| >= q,")1028 print("hence C_j >= q C_(j-1) at every j, hence c_j = C_j - k C_(j-1) >= (m/q) C_j,")1029 print("and sum_(i <= J) k^(j-i) c_i = k^(j-J) C_J <= (k/q)^(j-J) C_j")1030 print("the share column runs over j <= floor(L/2), the levels whose charge is at most x^(3/4), the tie at")1031 print("j = L/2 included, so it over-reports the levels the PROPOSITION strictly beats the uniform at;")1032 print("the q = 101 and q = 1499 families stop at j = 3 and j = 2, where C_j/C_(j-1) is still moving,")1033 print("244.658399 then 234.507307 at q = 101, so their top shares are short rows and not constants")1034 print(" q E j C_j C_j/C_(j-1) q c_j/C_j floor m/q den<=x^(1/2) cap (k/q)^(j-J)")1035 prof = {}1036 for q, E, J in PERDEN_FAMILIES:1037 C, c = levmass(q, E, J)1038 prof[(q, E[0])] = (C, c, J)1039 k, m = q - len(E), len(E)1040 for j in range(1, J + 1):1041 sh = shares(q, E, C, c, j)1042 Jc = j // 21043 low = sum(sh[:Jc + 1])1044 dec = sum(k ** (j - i) * c[i] for i in range(j + 1))1045 assert abs(dec - C[j]) <= 1e-9 * C[j]1046 assert C[j] >= q * C[j - 1] * (1 - 1e-12)1047 assert c[j] / C[j] >= m / q - 1e-121048 assert low <= (k / q) ** (j - Jc) + 1e-121049 print("%5d %5d %3d %13.4f %11.6f %6d %10.6f %11.6f %13.6f %14.6f"1050 % (q, E[0], j, C[j], C[j] / C[j - 1], q, c[j] / C[j], m / q, low, (k / q) ** (j - Jc)))1051 print()1052 print("THE EXPONENT IT BUYS, against the uniform step 5")1053 print("unif = log_q(C_L/q^L)/L + b(a), den = max_j [log_q(k^(L-j) c_j/q^L)/L + min(b(a), a + j/(2L))]")1054 print("ratio = the per-denominator sum divided by C_L x^(b(a)), proved to lie in [m/q, 1]")1055 print("proved: every level gain is at most 0 and the top level's is log_q(c_L/C_L)/L, so the saving is at")1056 print("MOST -log_q(top share)/L, itself at most log_q(q/m)/L; that the largest term is the top level is")1057 print("printed as argmax j and not proved. Read on the whole bound rather than on its largest term the")1058 print("saving is log_q(ratio)/L, smaller again, the two agreeing only through the x^eps that absorbs L+1")1059 print("terms. Either reading is a single factor at most q/m and never an exponent")1060 print("q E L a b(a) unif den den - unif L(den-unif) -log_q top cap log_q(q/m) argmax j ratio log_q(ratio)/L")1061 for q, E, J in PERDEN_FAMILIES:1062 C, c, _ = prof[(q, E[0])]1063 k, m = q - len(E), len(E)1064 a, b = 0.5, 0.751065 for L in ([J] if J < 6 else [J - 2, J]):1066 sh = shares(q, E, C, c, L)1067 jm, g, ratio = perden_exp(q, sh, L, a, b)1068 unif = log(C[L]) / (L * log(q)) - 1 + b1069 sp = -log(sh[L]) / log(q)1070 assert -log(q / m) / log(q) / L - 1e-12 <= -sp / L <= g <= 1e-121071 assert m / q - 1e-12 <= ratio <= 1 + 1e-121072 print("%5d %5d %3d %6.3f %8.5f %11.6f %11.6f %11.6f %12.6f %11.6f %15.6f %9d %9.6f %11.6f"1073 % (q, E[0], L, a, b, unif, unif + g, g, L * g, sp, log(q / m) / log(q), jm, ratio,1074 log(ratio) / (L * log(q))))1075 print()1076 print("THE FULL MINOR-ARC FORM, the same tool at full strength (measured, not proved)")1077 print("any reduced r/Q may serve any frequency, so the per-frequency constant is min(x^(b(a)), x^a nu^(1/2))")1078 print("with nu(a) = min_Q (Q + ||aQ||_(q^L)); dropping the coprimality changes nothing, a smaller Q being")1079 print("never worse, and the scan Q <= 2 q^(L/2) is exact because Dirichlet gives nu(a) <= 2 q^(L/2)")1080 assert float(numpy.abs(nu(81) - numpy.array(nu_slow(81))).max()) == 0.01081 print("nu checked against a full search over every reduced r/Q at q^L = 81: 0 mismatches")1082 print("ratio_nu = that weighted sum over C_L x^(b(a)); the exact-frequency ratio is the column beside it")1083 print("q E L ratio_nu log_q(ratio_nu)/L L times that exact-frequency ratio")1084 for q, E, Ls in PERDEN_HONEST:1085 C, c, _ = prof[(q, E[0])]1086 for L in Ls:1087 rn = honest(q, E, L, 0.5, 0.75)1088 gn = log(rn) / (L * log(q))1089 _, _, ratio = perden_exp(q, shares(q, E, C, c, L), L, 0.5, 0.75)1090 assert 0.0 < rn <= ratio + 1e-121091 print("%5d %5d %3d %10.6f %19.6f %14.6f %22.6f" % (q, E[0], L, rn, gn, L * gn, ratio))1092 print("ratio_nu rises with L while L times the gain falls, so the full form buys a bounded factor too and")1093 print("no exponent; that ratio_nu is bounded below in L is measured over these rows and is not proved")1094 print()1095 print("THE TRANSFORM, CHECKED WITHOUT ITSELF: brute force over the digit strings")1096 print("hatF_j(a/q^j) summed over the k^j allowed n directly, against the product form hat_abs")1097 print("q E j C_j brute C_j product c_j/C_j brute c_j/C_j product")1098 for q, E, j in PERDEN_BRUTE:1099 tot, top = brute(q, E, j)1100 C, c = levmass(q, E, j)1101 assert abs(tot - C[j]) <= 1e-9 * C[j] and abs(top - c[j]) <= 1e-9 * C[j]1102 print("%5d %5d %3d %13.6f %13.6f %14.6f %16.6f" % (q, E[0], j, tot, C[j], top / tot, c[j] / C[j]))1103 print()1104 print("WHERE THE PROPOSITION DOES BELONG: the d-form minor arc, not the q-power grid")1105 print("on a Dirichlet arc |theta - l/d| <= 1/d^2 with (l,d) = 1 the PROPOSITION reads")1106 print("S_mu(theta) << x^(a+eps) d^(1/2) (1 + x/d^2)^(1/2) = x^(a+eps) (d + x/d)^(1/2) <= x^(a+eps) (d^(1/2) + x^(1/2) d^(-1/2))")1107 print("and at a = 1/2 that is x^eps ((x d)^(1/2) + x d^(-1/2)), the first and third terms of the d-form")1108 print("minor-arc bound Theorem L5 owes, whose middle term x^(4/5) is the unconditional Vaughan piece;")1109 print("this row applies eq. 6 at an arbitrary arc denominator d, so it holds only within the printed")1110 print("range on Q that the source carries and the desk has not read, unlike the blocks above, where the")1111 print("min caps the PROPOSITION by the uniform THEOREM at the level that decides")1112 print()1113 print("THE WALL IT MOVES: none, and no wall row is printed here")1114 print("the exponent is b(a) + c_q with the SAME c_q, so the certificate (q-1) q^(-b) - PB_q(1, e_0) > 0")1115 print("takes no per-denominator quantity at all and the GRH walls stay exactly where verb onestep's chord")1116 print("put them; re-evaluating that certificate here would discriminate nothing")11171118VERBS["perden"] = perden11191120if __name__ == "__main__":1121 VERBS[sys.argv[1] if len(sys.argv) > 1 else "split"]()