design_zeta.py
12.2 kB · python · 329 lines
1import sys2import time34import mpmath as mp56# ENGINE78class Design:9 def __init__(self, q, F, P=None, dps=45, tol=mp.mpf(10) ** -22):10 mp.mp.dps = dps11 self.dps = dps12 self.q = q13 self.F = tuple(sorted(F))14 self.k = len(self.F)15 self.k1 = len([a for a in self.F if a != 0])16 self.amax = max(self.F)17 self.tol = tol18 if P is None:19 P = 220 while self.k > 1 and self.k1 * self.k ** P <= 100:21 P += 122 self.P = P23 self.alpha = mp.log(self.k) / mp.log(q)24 self.gam = [mp.mpf(sum(a ** l for a in self.F)) for l in range(0, 1200)]25 self.lo = sorted(n for j in range(1, P) for n in self.strings(j))26 self.mid = self.strings(P)27 self.loglo = [mp.log(n) for n in self.lo]28 self.logmid = [mp.log(n) for n in self.mid]2930 def strings(self, p):31 cur = [a for a in self.F if a != 0]32 for _ in range(p - 1):33 cur = sorted(self.q * m + a for m in cur for a in self.F)34 return sorted(cur)3536 def tail(self, p, sig):37 r = self.k * mp.power(self.q, -sig)38 if r >= 1:39 return mp.inf40 return self.k1 * r ** (p - 1) / (1 - r)4142 def poly_lo(self, w):43 return mp.fsum([mp.exp(-w * L) for L in self.loglo])4445 def poly_mid(self, w):46 return mp.fsum([mp.exp(-w * L) for L in self.logmid])4748 def cut(self, w, L):49 aw, sig = abs(w), mp.re(w)50 rho = max((aw + L + 1) / (L + 2), 1) * self.amax / mp.power(self.q, self.P)51 if rho >= mp.mpf("0.9"):52 return mp.inf53 maj = mp.binomial(aw + L, L + 1) * mp.power(self.q, -sig - (L + 1)) \54 * self.k * self.amax ** (L + 1) * self.tail(self.P, sig + L + 1)55 return maj / (1 - rho)5657 def ladder(self, s, W, L, top=True):58 q, k = self.q, self.k59 J = max(0, int(mp.ceil(W - mp.re(s))))60 g, e = {}, {}61 for j in range(J, J + L + 2):62 g[j] = mp.mpf(0)63 e[j] = self.tail(self.P, mp.re(s) + j)64 for j in range(J - 1, -1, -1):65 w = s + j66 acc = self.poly_mid(w)67 eacc = self.cut(w, L)68 c = mp.mpf(1)69 for l in range(1, L + 1):70 c = c * (-w - (l - 1)) / l71 co = mp.power(q, -w - l) * c * self.gam[l]72 acc += co * g[j + l]73 eacc += abs(co) * e[j + l]74 if not top and j == 0:75 return acc, eacc76 den = 1 - k * mp.power(q, -w)77 g[j] = acc / den78 e[j] = eacc / abs(den)79 return g[0], e[0]8081 def tune(self, s, top=True):82 W, L = mp.mpf(6), 883 while W <= 400:84 v, e = self.ladder(s, W, L, top)85 if e < self.tol:86 return v, e87 W += 688 L += 689 raise RuntimeError("FAIL tolerance " + mp.nstr(self.tol, 3) + " not reached at s = "90 + mp.nstr(s, 14) + ", best bound " + mp.nstr(e, 4))9192 def zeta(self, s):93 s = mp.mpmathify(s)94 old = mp.mp.dps95 mp.mp.dps = self.dps + 30 + int(abs(mp.im(s)) / 4)96 s = mp.mpmathify(s)97 g, e = self.tune(s)98 v = self.poly_lo(s) + g99 mp.mp.dps = old100 return +v, +e101102 def cofactor(self, s):103 s = mp.mpmathify(s)104 old = mp.mp.dps105 mp.mp.dps = self.dps + 30 + int(abs(mp.im(s)) / 4)106 s = mp.mpmathify(s)107 a, e = self.tune(s, top=False)108 v = (1 - self.k * mp.power(self.q, -s)) * self.poly_lo(s) + a109 mp.mp.dps = old110 return +v, +e111112 def residue(self, j):113 old = mp.mp.dps114 mp.mp.dps = self.dps + 30 + int(abs(j) * 2)115 sj = self.alpha + 2 * mp.pi * 1j * j / mp.log(self.q)116 a, e = self.tune(sj, top=False)117 r, er = a / mp.log(self.q), e / mp.log(self.q)118 mp.mp.dps = old119 return +r, +er120121# CONTOURS122123class Cache:124 def __init__(self, d, cof=False):125 self.d = d126 self.cof = cof127 self.m = {}128 self.n = 0129 self.emax = mp.mpf(0)130131 def __call__(self, z):132 key = (mp.nstr(mp.re(z), 22), mp.nstr(mp.im(z), 22))133 if key not in self.m:134 v, e = self.d.cofactor(z) if self.cof else self.d.zeta(z)135 self.m[key] = v136 self.emax = max(self.emax, e)137 self.n += 1138 return self.m[key]139140def phase(f, pt, n0, cap=1.0, budget=4000):141 us = [mp.mpf(i) / n0 for i in range(n0 + 1)]142 vs = [f(pt(u)) for u in us]143 for _ in range(30):144 bad = [i for i in range(len(us) - 1) if abs(mp.arg(vs[i + 1] / vs[i])) > cap]145 if not bad or len(us) > budget:146 break147 nu, nv = [], []148 for i in range(len(us) - 1):149 nu.append(us[i])150 nv.append(vs[i])151 if i in bad:152 um = (us[i] + us[i + 1]) / 2153 nu.append(um)154 nv.append(f(pt(um)))155 nu.append(us[-1])156 nv.append(vs[-1])157 us, vs = nu, nv158 tot = mp.mpf(0)159 mx = mp.mpf(0)160 for i in range(len(vs) - 1):161 dd = mp.arg(vs[i + 1] / vs[i])162 tot += dd163 mx = max(mx, abs(dd))164 return tot, mx165166def box_phase(f, x0, x1, y0, y1, n):167 a = phase(f, lambda u: mp.mpc(x0 + (x1 - x0) * u, y0), n)168 b = phase(f, lambda u: mp.mpc(x1, y0 + (y1 - y0) * u), n)169 c = phase(f, lambda u: mp.mpc(x1 - (x1 - x0) * u, y1), n)170 d = phase(f, lambda u: mp.mpc(x0, y1 - (y1 - y0) * u), n)171 tot = a[0] + b[0] + c[0] + d[0]172 return tot / (2 * mp.pi), max(a[1], b[1], c[1], d[1])173174def strip(f, x0, x1, y0, y1, nsub, n):175 cuts = [y0 + (y1 - y0) * i / nsub for i in range(nsub + 1)]176 hor = {}177 for y in cuts:178 hor[y] = phase(f, lambda u, y=y: mp.mpc(x0 + (x1 - x0) * u, y), n)179 out = []180 for i in range(nsub):181 a, b = cuts[i], cuts[i + 1]182 rt = phase(f, lambda u, a=a, b=b: mp.mpc(x1, a + (b - a) * u), n)183 lf = phase(f, lambda u, a=a, b=b: mp.mpc(x0, a + (b - a) * u), n)184 tot = hor[a][0] + rt[0] - hor[b][0] - lf[0]185 mx = max(hor[a][1], rt[1], hor[b][1], lf[1])186 out.append((a, b, tot / (2 * mp.pi), mx))187 return out188189def locate(f, box, steps=9):190 x0, x1, y0, y1 = [mp.mpf(v) for v in box]191 for _ in range(steps):192 xm, ym = (x0 + x1) / 2, (y0 + y1) / 2193 hit = False194 for (a, b, c, d) in [(x0, xm, y0, ym), (xm, x1, y0, ym), (x0, xm, ym, y1), (xm, x1, ym, y1)]:195 if abs(box_phase(f, a, b, c, d, 16)[0]) > mp.mpf("0.5"):196 x0, x1, y0, y1 = a, b, c, d197 hit = True198 break199 if not hit:200 break201 return mp.mpc((x0 + x1) / 2, (y0 + y1) / 2), max(x1 - x0, y1 - y0)202203def zero(d, box, cof=False):204 f = Cache(d, cof=cof)205 z0, w = locate(f, box)206 r = mp.findroot(f, [z0 - w, z0, z0 + w * 1j], solver="muller", tol=mp.mpf(10) ** -30, maxsteps=60)207 if abs(r - z0) > 4 * w:208 r = z0209 v = abs(f(r))210 return r, v, f.emax211212# STUDY213214def line(*a):215 print(" ".join(str(x) for x in a))216217def control():218 line("CONTROL base 2 full digit set, zeta_F = zeta")219 d = Design(2, (0, 1))220 for s in [mp.mpf(2), mp.mpf("0.5") + 14.134725141734693j, mp.mpf("0.3") + 40j,221 -1 + 2j, mp.mpf("0.5") + 100j]:222 v, e = d.zeta(s)223 mp.mp.dps = 60224 line(" s", mp.nstr(s, 8), "gap to mpmath zeta", mp.nstr(abs(v - mp.zeta(s)), 4),225 "proved bound", mp.nstr(e, 3))226 r, e = d.residue(0)227 line(" residue at alpha = 1", mp.nstr(r, 20), "bound", mp.nstr(e, 3))228 r, e = d.residue(1)229 line(" residue at alpha + 2 pi i/log 2", mp.nstr(abs(r), 4), "bound", mp.nstr(e, 3))230231def residues():232 line("RESIDUES against the certified enclosures of lab/py/burnol-residue")233 for q, F in [(3, (0, 1)), (3, (0, 2))]:234 d = Design(q, F)235 for j in (0, 1, 2):236 r, e = d.residue(j)237 line(" q", q, "F", F, "j", j, mp.nstr(r, 18), "bound", mp.nstr(e, 3))238239def scaling():240 line("SCALING zeta_(aF)(s) = a^(-s) zeta_F(s)")241 d1 = Design(3, (0, 1))242 d2 = Design(3, (0, 2))243 for s in [mp.mpf(2) + 3j, mp.mpf("0.8") + 23j, -mp.mpf("0.5") + 11j]:244 a = d2.zeta(s)[0]245 b = mp.power(2, -s) * d1.zeta(s)[0]246 line(" s", mp.nstr(s, 8), "gap", mp.nstr(abs(a - b), 4))247248def census(q, F, lo, hi, ymax, nsub, n=20, cof=True):249 d = Design(q, F, tol=mp.mpf(10) ** -10, dps=25)250 f = Cache(d, cof=cof)251 x0, x1 = d.alpha + mp.mpf(lo), d.alpha + mp.mpf(hi)252 t0 = time.time()253 res = strip(f, x0, x1, mp.mpf("0.02"), mp.mpf(ymax), nsub, n)254 tot = mp.fsum([w for a, b, w, mx in res])255 mxall = max([mx for a, b, w, mx in res])256 hits = [(a, b, w, mx) for a, b, w, mx in res if abs(w) > mp.mpf("0.2")]257 line(" q", q, "F", F, "alpha", mp.nstr(d.alpha, 12), "object",258 "Z = zeta_F (1 - k q^(-s))" if cof else "zeta_F",259 "strip Re in", mp.nstr(x0, 12), mp.nstr(x1, 12), "Im in 0.02", ymax)260 for a, b, w, mx in hits:261 line(" Im [", mp.nstr(a, 6), ",", mp.nstr(b, 6), "] winding", mp.nstr(w, 8),262 "max phase step", mp.nstr(mx, 4))263 line(" count", mp.nstr(tot, 10), "boxes", len(hits), "largest phase step",264 mp.nstr(mxall, 4), "largest bound on the contour", mp.nstr(f.emax, 4),265 "evaluations", f.n, "seconds", round(time.time() - t0, 1))266 return tot, hits267268def zeros(q, F, lo, hi, boxes, cof=False):269 d = Design(q, F, tol=mp.mpf(10) ** -22, dps=45)270 per = 2 * mp.pi / mp.log(q)271 out = []272 for (a, b) in boxes:273 r, v, e = zero(d, (d.alpha + mp.mpf(lo), d.alpha + mp.mpf(hi), a, b), cof=cof)274 out.append(r)275 line(" zero", mp.nstr(r, 16), "abs", "Z" if cof else "zeta_F", mp.nstr(v, 3),276 "largest bound", mp.nstr(e, 3),277 "Im/period", mp.nstr(mp.im(r) / per, 12))278 line(" period 2 pi/log q", mp.nstr(per, 14), "alpha", mp.nstr(d.alpha, 12),279 "alpha/2", mp.nstr(d.alpha / 2, 12))280 for i in range(len(out) - 1):281 g = mp.im(out[i + 1]) - mp.im(out[i])282 line(" gap", mp.nstr(g, 14), "gap minus period", mp.nstr(g - per, 8))283 return out284285def control_census(ymax):286 line("CONTROL base 2 full digit set, the census run on BOTH sides of the abscissa alpha = 1")287 a, _ = census(2, (0, 1), 0.02, 3.02, ymax, int(ymax / 2), cof=False)288 line(" object counted zeros of zeta_F = zeta right of the abscissa, count", mp.nstr(a, 6))289 b, _ = census(2, (0, 1), -0.98, -0.02, ymax, int(ymax / 2), cof=False)290 line(" object counted zeros of zeta_F = zeta left of the abscissa, count", mp.nstr(b, 6))291 c, _ = census(2, (0, 1), -0.92, 3.02, ymax, int(ymax / 2), cof=True)292 teeth = int(mp.floor(mp.mpf(ymax) * mp.log(2) / (2 * mp.pi)))293 line(" object counted zeros of the cofactor Z on one strip, count", mp.nstr(c, 6),294 "= zeros of zeta", mp.nstr(b, 6), "plus teeth of 1 - 2 q^(-s) on Re s = alpha",295 teeth, "predicted floor(T log q/2 pi)", teeth)296 line(" the Euler product forbids zeros of zeta in Re s >= 1, and the strip Re in",297 "[alpha + 0.02, alpha + 3.02] confirms it; the sliver alpha < Re s < alpha + 0.02",298 "carries the teeth of the cofactor and no zero of zeta")299300def main():301 full = "--full" in sys.argv302 ymax = 60 if full else 30303 y10 = ymax304 for a in sys.argv:305 if a.startswith("--y10="):306 y10 = float(a.split("=")[1])307 t0 = time.time()308 control()309 residues()310 scaling()311 control_census(ymax)312 line("CENSUS zeros of the cofactor Z = zeta_F (1 - k q^(-s)), analytic on Re s > alpha - 1")313 line(" the split is at Re s = alpha exactly, so no sliver is left unscanned")314 census(3, (0, 1), 0.0, 3.02, ymax, int(ymax / 2))315 census(3, (0, 1), -0.92, 0.0, ymax, int(ymax / 2))316 census(3, (0, 2), 0.0, 3.02, ymax, int(ymax / 2))317 if full:318 census(10, tuple(range(9)), 0.0, 3.02, y10, int(y10 / 2))319 census(10, tuple(range(9)), -0.92, 0.0, y10, int(y10 / 2))320 line("ZEROS polished, right of the abscissa")321 zeros(3, (0, 1), 0.02, 3.02, [(22, 24), (28, 30)])322 if full:323 zeros(10, tuple(range(9)), 0.02, 3.02, [(2, 4), (4, 6)])324 line("ZEROS polished on the cofactor, base 3 digits 0 1, left of the abscissa")325 zeros(3, (0, 1), -0.92, 0.0, [(6, 8), (10, 12)], cof=True)326 line("seconds", round(time.time() - t0, 1))327328if __name__ == "__main__":329 main()