zeta_locus.py
20.9 kB · python · 545 lines
1import os2import sys3import time45import mpmath as mp67HERE = os.path.dirname(os.path.abspath(__file__))8sys.path.insert(0, os.path.join(HERE, "..", "design-zeta"))910from design_zeta import Design, phase, strip1112# COFACTOR1314def build(q, F, dps=45, tol=None, P=None):15 if tol is None:16 tol = mp.mpf(10) ** -2217 if P is None and len(F) == 1:18 P = 1219 return Design(q, F, P=P, dps=dps, tol=tol)2021def ztune(d, s):22 key = int(abs(mp.im(s)) / 8)23 W, L = getattr(d, "wl", {}).get(key, (mp.mpf(6), 8))24 v, e = d.ladder(s, W, L, False)25 while e >= d.tol and W <= 54:26 W += 627 L += 628 v, e = d.ladder(s, W, L, False)29 if e < d.tol:30 if not hasattr(d, "wl"):31 d.wl = {}32 d.wl[key] = (max(mp.mpf(6), W - 6), max(8, L - 6))33 return v, e3435def zwork(d, s):36 a, e = ztune(d, s)37 return (1 - d.k * mp.power(d.q, -s)) * d.poly_lo(s) + a, e3839def cofactor(d, s):40 old = mp.mp.dps41 mp.mp.dps = d.dps + 30 + int(abs(mp.im(s)) / 4)42 v, e = zwork(d, mp.mpmathify(s))43 mp.mp.dps = old44 return +v, +e4546class ZCache:47 def __init__(self, d):48 self.d = d49 self.m = {}50 self.n = 05152 def __call__(self, z):53 key = (mp.nstr(mp.re(z), 22), mp.nstr(mp.im(z), 22))54 if key not in self.m:55 self.m[key] = cofactor(self.d, z)[0]56 self.n += 157 return self.m[key]5859def laurent(d, j):60 old = mp.mp.dps61 mp.mp.dps = d.dps + 40 + int(abs(j) * 2)62 L = mp.log(d.q)63 s0 = mp.log(d.k) / L + 2 * mp.pi * 1j * j / L64 h = mp.mpf(10) ** -565 z0, e0 = zwork(d, s0)66 zp = zwork(d, s0 + h)[0]67 zm = zwork(d, s0 - h)[0]68 z1 = (zp - zm) / (2 * h)69 z2 = (zp + zm - 2 * z0) / (2 * h * h)70 r = z0 / L71 c0 = z1 / L + z0 / 272 c1 = z2 / L + z1 / 2 + z0 * L / 1273 mp.mp.dps = old74 return +s0, +r, +c0, +c1, +e07576def predict(r, c0, c1):77 u1 = -r / c078 if c1 == 0:79 return u1, u180 disc = mp.sqrt(c0 * c0 - 4 * c1 * r)81 roots = [(-c0 + disc) / (2 * c1), (-c0 - disc) / (2 * c1)]82 u2 = min(roots, key=lambda z: abs(z - u1))83 return u1, u28485def disc(f, s0, rho, n=32):86 t, mx = phase(f, lambda u: s0 + rho * mp.exp(2 * mp.pi * 1j * u), n, budget=280)87 return t / (2 * mp.pi), mx8889def ring(f, s0, rho, nr=5, na=12):90 pts = [(abs(f(s0)), s0)]91 for i in range(1, nr + 1):92 rr = rho * mp.mpf(i) / (nr + 1)93 for m in range(na):94 pts.append((abs(f(s0 + rr * mp.exp(2 * mp.pi * 1j * mp.mpf(m) / na))), i, m))95 pts.sort(key=lambda t: t[0])96 out = []97 for t in pts:98 if len(t) == 2:99 out.append(t[1])100 else:101 rr = rho * mp.mpf(t[1]) / (nr + 1)102 out.append(s0 + rr * mp.exp(2 * mp.pi * 1j * mp.mpf(t[2]) / na))103 return out104105def teeth(f, dp, s0, rho, want, skip):106 out = []107 for z0 in ring(f, s0, rho)[:want + 6]:108 if len(out) >= want:109 break110 z, v = polish(f.d, z0, rho / 40, rho, dp)111 if z is None or abs(z - s0) > rho:112 continue113 if skip and abs(z - s0) < mp.mpf("1e-9"):114 continue115 if any(abs(z - y) < mp.mpf("1e-9") for y in out):116 continue117 out.append(z)118 return out119120def seek(f, x0, x1, y0, y1, take, nx=11, ny=7):121 pts = [(abs(f(mp.mpc(x0 + (x1 - x0) * i / (nx - 1), y0 + (y1 - y0) * m / (ny - 1)))),122 i, m) for i in range(nx) for m in range(ny)]123 pts.sort()124 return [mp.mpc(x0 + (x1 - x0) * i / (nx - 1), y0 + (y1 - y0) * m / (ny - 1))125 for v, i, m in pts[:take]]126127def muller(d, z0, w, span, tol, steps):128 def f(z):129 if abs(z - z0) > span:130 raise ValueError131 return cofactor(d, z)[0]132 try:133 r = mp.findroot(f, [z0 - w, z0, z0 + w * 1j], solver="muller", tol=tol, maxsteps=steps)134 except Exception:135 return None136 return r if abs(r - z0) <= span else None137138def polish(d, z0, w, span, dp=None):139 r = muller(d, z0, w, span, mp.mpf(10) ** -18, 40)140 if r is None:141 return None, None142 if dp is not None:143 r2 = muller(dp, r, mp.mpf(10) ** -9, mp.mpf("1e-4"), mp.mpf(10) ** -20, 15)144 if r2 is not None:145 r = r2146 return r, abs(cofactor(dp, r)[0])147 return r, abs(cofactor(d, r)[0])148149# SWEEP150151def line(*a):152 print(" ".join(str(x) for x in a))153154def tag(q, F):155 return "q" + str(q) + "F" + "".join(str(a) for a in F)156157def alpha_of(q, F):158 return mp.log(len(F)) / mp.log(q)159160def row(q, F, alpha, per, z, kind, pred):161 t = mp.im(z) / per162 line(" ROW", tag(q, F), "Re", mp.nstr(mp.re(z), 12), "Im", mp.nstr(mp.im(z), 12),163 "Im/per", mp.nstr(t, 10), "frac", mp.nstr(t - mp.floor(t), 8),164 "alpha", mp.nstr(alpha, 10), "k/q", mp.nstr(mp.mpf(len(F)) / q, 8),165 "Re-alpha", mp.nstr(mp.re(z) - alpha, 10), kind, pred)166167# STUDY168169Q3 = [(3, (1,)), (3, (0, 1)), (3, (1, 2)), (3, (0, 1, 2))]170Q4 = [(4, (1,)), (4, (0, 1)), (4, (1, 2)), (4, (1, 3)), (4, (2, 3)),171 (4, (0, 1, 2)), (4, (0, 1, 3)), (4, (0, 2, 3)), (4, (1, 2, 3)), (4, (0, 1, 2, 3))]172Q5 = [(5, (0, 1)), (5, (1, 2))]173HALF = [(9, (0, 1, 2)), (16, (0, 1, 2, 3))]174WIDE = [(10, tuple(range(9)))]175CONTROL = [(2, (0, 1))]176177def comb(q, F, ymax, dps=40, tol=mp.mpf(10) ** -24, verbose=True):178 d = build(q, F, dps=dps, tol=tol)179 dp = build(q, F, dps=32, tol=mp.mpf(10) ** -18)180 fw = ZCache(build(q, F, dps=25, tol=mp.mpf(10) ** -10))181 rho = mp.mpf("0.45")182 L = mp.log(q)183 per = 2 * mp.pi / L184 alpha = alpha_of(q, F)185 if verbose:186 line(" ", tag(q, F), "alpha", mp.nstr(alpha, 12), "period", mp.nstr(per, 12),187 "k/q", mp.nstr(mp.mpf(len(F)) / q, 8), "disc radius", mp.nstr(rho, 4))188 out = []189 for j in range(0, int(mp.ceil(mp.mpf(ymax) / per)) + 1):190 s0, r, c0, c1, e = laurent(d, j)191 u1, u2 = predict(r, c0, c1)192 w, wmx = disc(fw, s0, rho)193 null = abs(r) < 1000 * e194 nz = int(mp.nint(w)) - (1 if null else 0)195 zl = teeth(fw, dp, s0, rho, nz, null) if nz > 0 else []196 z = min(zl, key=lambda y: abs(y - s0)) if zl else None197 rec = {"j": j, "s0": s0, "r": r, "R": c0, "u1": u1, "u2": u2, "z": z,198 "w": w, "wmx": wmx, "null": null, "nz": nz, "all": zl,199 "v": None if z is None else abs(fw(z))}200 out.append(rec)201 if not verbose:202 continue203 msg = [" j", j, "r", mp.nstr(r, 12), "R", mp.nstr(c0, 12),204 "u1", mp.nstr(u1, 10), "u2-u1", mp.nstr(abs(u2 - u1), 8),205 "ladder bound", mp.nstr(e, 3), "Z zeros in disc", mp.nstr(w, 8),206 "step", mp.nstr(wmx, 4), "residue null", null,207 "zeta_F zeros in disc", nz, "located", len(zl)]208 if z is None:209 msg += ["tooth", "none"]210 else:211 m1 = abs(z - s0 - u1)212 m2 = abs(z - s0 - u2)213 msg += ["tooth", mp.nstr(z, 14), "|u|", mp.nstr(abs(z - s0), 8),214 "|Z|", mp.nstr(rec["v"], 3), "miss1", mp.nstr(m1, 6),215 "miss2", mp.nstr(m2, 6), "miss2/miss1",216 mp.nstr(m2 / max(m1, mp.mpf(10) ** -40), 5)]217 line(*msg)218 return out219220def census(q, F, ymax, nsub, lo=-0.92, hi=3.02, n=20, cmb=None, deep=False):221 d = build(q, F, dps=25, tol=mp.mpf(10) ** -10)222 dp = build(q, F, dps=32, tol=mp.mpf(10) ** -18)223 L = mp.log(q)224 per = 2 * mp.pi / L225 alpha = alpha_of(q, F)226 x0, x1 = alpha + mp.mpf(lo), alpha + mp.mpf(hi)227 f = ZCache(d)228 t0 = time.time()229 res = strip(f, x0, x1, mp.mpf("0.02"), mp.mpf(ymax), nsub, n)230 tot = sum(r[2] for r in res)231 mx = max(r[3] for r in res)232 bx = [r for r in res if abs(r[2]) > mp.mpf("0.2")]233 nulls = [mp.im(c["s0"]) for c in (cmb or []) if c["null"] and c["j"] > 0234 and mp.im(c["s0"]) <= ymax]235 line(" ", tag(q, F), "alpha", mp.nstr(alpha, 10), "strip Re", mp.nstr(x0, 8),236 mp.nstr(x1, 8), "Im to", ymax, "Z winding", mp.nstr(tot, 8), "null teeth",237 len(nulls), "zeta_F zeros", mp.nstr(tot - len(nulls), 8), "boxes", len(bx),238 "max phase step", mp.nstr(mx, 4), "evals", f.n, "sec", round(time.time() - t0, 1))239 cum, marks = mp.mpf(0), []240 for a, b, w, mxx in res:241 cum += w242 marks.append((b, cum - len([y for y in nulls if y <= b])))243 line(" N_F(T)", " ".join("T " + mp.nstr(T, 6) + " N " + mp.nstr(c, 6)244 for T, c in marks[nsub // 4 - 1::max(1, nsub // 4)]),245 "per period", mp.nstr((tot - len(nulls)) * per / mp.mpf(ymax), 8))246 zs = [(c["z"], "comb j " + str(c["j"])) for c in (cmb or [])247 if c["z"] is not None and mp.im(c["z"]) <= ymax + 1]248 for c in (cmb or []):249 for y in c["all"]:250 if y is c["z"] or mp.im(y) > ymax + 1:251 continue252 if not any(abs(y - t[0]) < mp.mpf("1e-9") for t in zs):253 zs.append((y, "comb j " + str(c["j"]) + " second tooth"))254 for a, b, w, mxx in (bx if deep else []):255 want = int(mp.nint(abs(w)))256 hit = 0257 for z0 in seek(f, x0, x1, a, b, want + 3):258 if hit >= want:259 break260 z, v = polish(d, z0, (b - a) / 40, x1 - x0, dp)261 if z is None or mp.re(z) < x0 or mp.re(z) > x1:262 continue263 if mp.im(z) < a - 1 or mp.im(z) > b + 1:264 continue265 hit += 1266 if any(abs(z - y[0]) < mp.mpf("1e-9") for y in zs):267 continue268 if any(c["null"] and abs(z - c["s0"]) < mp.mpf("1e-9") for c in (cmb or [])):269 continue270 zs.append((z, "second"))271 line(" located", len(zs), "of zeta_F zeros", mp.nstr(tot - len(nulls), 8),272 "complete", len(zs) == int(mp.nint(tot)) - len(nulls), "deep", deep,273 "teeth", sum(1 for z, k in zs if k.startswith("comb")))274 for z, kind in sorted(zs, key=lambda y: mp.im(y[0])):275 pred = "-"276 if cmb and kind.startswith("comb"):277 best = min(cmb, key=lambda c: abs(z - c["s0"]))278 pred = mp.nstr(best["s0"] + best["u1"], 12) + " miss " + \279 mp.nstr(abs(z - best["s0"] - best["u1"]), 6)280 row(q, F, alpha, per, z, kind, pred)281 return zs, tot, res282283def laws(got):284 line("LAWS")285 groups = {}286 for q, F in got:287 groups.setdefault(mp.nstr(alpha_of(q, F), 10), []).append((q, F))288 for a in sorted(groups):289 fam = groups[a]290 if len(fam) < 2:291 continue292 line(" EQUAL ALPHA", a, "designs", " | ".join(tag(q, F) for q, F in fam))293 for q, F in fam:294 line(" ", tag(q, F), "k/q", mp.nstr(mp.mpf(len(F)) / q, 8),295 "zeros", len(got[(q, F)]), "complete", got[(q, F)][0][2] if got[(q, F)] else "-")296 for i in range(len(fam)):297 for m in range(i + 1, len(fam)):298 A = [z for z, k, c in got[fam[i]]]299 B = [z for z, k, c in got[fam[m]]]300 if not A or not B:301 continue302 pairs = [(abs(mp.im(z) - mp.im(y)), abs(mp.re(z) - mp.re(y)), z, y)303 for z in A for y in B]304 shared = min(abs(z - y) for z in A for y in B)305 pairs.sort()306 di, dr, z, y = pairs[0]307 near = [p for p in pairs if p[0] < mp.mpf("0.02")]308 worst = max(near, key=lambda p: p[1]) if near else None309 line(" ", tag(*fam[i]), "vs", tag(*fam[m]), "closest in Im: dIm",310 mp.nstr(di, 6), "dRe", mp.nstr(dr, 6), "at Im", mp.nstr(mp.im(z), 10),311 "| shared zero", mp.nstr(shared, 6),312 "| worst dRe at dIm<0.02", "none" if worst is None else313 mp.nstr(worst[1], 6) + " at Im " + mp.nstr(mp.im(worst[2]), 10))314 line(" SECOND FAMILY, the zeros left when every tooth is stripped")315 for q, F in got:316 sec = [mp.re(z) for z, k, c in got[(q, F)] if k == "second"]317 a2 = alpha_of(q, F) / 2318 line(" ", tag(q, F), "alpha/2", mp.nstr(a2, 10), "n", len(sec),319 "Re range", "-" if not sec else mp.nstr(min(sec), 8),320 "-" if not sec else mp.nstr(max(sec), 8),321 "complete", got[(q, F)][0][2] if got[(q, F)] else "-")322 line(" COMB in frac(Im s log q/2 pi): worst Re gap between zeros of equal frac")323 for q, F in got:324 per = 2 * mp.pi / mp.log(q)325 fr = sorted(((mp.im(z) / per) % 1, mp.re(z)) for z, k, c in got[(q, F)])326 if len(fr) < 2:327 continue328 pairs = [abs(fr[i][1] - fr[i + 1][1]) for i in range(len(fr) - 1)329 if fr[i + 1][0] - fr[i][0] < mp.mpf("0.03")]330 line(" ", tag(q, F), "n", len(fr), "Re range",331 mp.nstr(min(x for f, x in fr), 8), mp.nstr(max(x for f, x in fr), 8),332 "worst gap at equal frac", "none" if not pairs else mp.nstr(max(pairs), 8))333334# ROUCHE335336def explicit(d, s0, M=60):337 L = mp.log(d.q)338 dm = [mp.mpf(0)] * (M + 1)339 em = [mp.mpf(0)] * (M + 1)340 for arr, tgt in ((d.loglo, dm), (d.logmid, em)):341 for t in arr:342 v = mp.exp(-s0 * t)343 for m in range(M + 1):344 tgt[m] += v345 v = v * (-t) / (m + 1)346 c = [mp.mpf(0)] + [-(-L) ** m / mp.factorial(m) for m in range(1, M + 1)]347 P = [em[m] + mp.fsum([c[i] * dm[m - i] for i in range(1, m + 1)]) for m in range(M + 1)]348 return P, c349350def ptail(d, s0, rho, P, M=60):351 sig = mp.re(s0)352 L = mp.log(d.q)353 h = M // 2 + 1354 sd = mp.fsum([mp.exp((rho - sig) * t) for t in d.loglo])355 cd = mp.fsum([mp.exp((rho - sig) * t) * (rho * t) ** h / mp.factorial(h) for t in d.loglo])356 ce = mp.fsum([mp.exp((rho - sig) * t) * (rho * t) ** h / mp.factorial(h) for t in d.logmid])357 C = mp.expm1(L * rho)358 Ch = (L * rho) ** h * mp.exp(L * rho) / mp.factorial(h)359 rem = ce + Ch * sd + C * cd360 return mp.fsum([abs(P[m]) * rho ** m for m in range(2, M + 1)]) + rem361362def gbound(d, sig, extra=2):363 if not hasattr(d, "lv"):364 cur = list(d.mid)365 d.lv = []366 for i in range(extra):367 d.lv.append([mp.log(n) for n in cur])368 cur = sorted(d.q * m + a for m in cur for a in d.F)369 t = d.tail(d.P + extra, sig)370 if t == mp.inf:371 return mp.inf372 for lg in d.lv:373 t += mp.fsum([mp.exp(-sig * x) for x in lg])374 return t375376def tbound(d, s0, R2):377 sig = mp.re(s0) - R2378 aw = abs(s0) + R2379 rat = d.amax * mp.power(d.q, -d.P)380 tot = mp.mpf(0)381 for l in range(1, 600):382 t = gbound(d, sig + l)383 if t == mp.inf:384 return mp.inf385 term = mp.binomial(aw + l - 1, l) * mp.power(d.q, -sig - l) * d.gam[l] * t386 tot += term387 rl = (aw + l) / (l + 1) * rat388 if rl < 1:389 return tot + term * rl / (1 - rl)390 return mp.inf391392def certify(d, s0, Z0c, e0, samp, R, N, rhos, r2s):393 P, c = explicit(d, s0)394 p1 = mp.log(d.q) * d.poly_lo(s0) - mp.fsum([mp.exp(-s0 * t) * t for t in d.logmid])395 maxe = max(e for v, e in samp)396 w = mp.exp(2 * mp.pi * 1j / N)397 T = []398 for m in range(N):399 u = R * w ** m400 Pu = (1 - mp.exp(-mp.log(d.q) * u)) * d.poly_lo(s0 + u) + d.poly_mid(s0 + u)401 T.append(samp[m][0] - Pu)402 c1 = mp.fsum([T[m] * w ** (-m) for m in range(N)]) / (N * R)403 U0 = abs(Z0c) + e0404 best = None405 for R2 in r2s:406 if R2 <= R:407 continue408 B = tbound(d, s0, mp.mpf(R2))409 if B == mp.inf:410 continue411 al = (mp.mpf(R) / R2) ** N412 e1 = maxe / R + (B / R2) * al / (1 - al)413 Z1 = p1 + c1414 L1 = abs(Z1) - e1415 if L1 <= 0:416 continue417 for rho in rhos:418 rho = mp.mpf(rho)419 if rho >= R2:420 continue421 minM = L1 * rho - U0422 if minM <= 0:423 continue424 tau = rho / R2425 marg = minM - ptail(d, s0, rho, P) - B * tau * tau / (1 - tau)426 if best is None or marg > best[0]:427 best = (marg, rho, mp.mpf(R2), minM, B, L1, U0)428 return best, U0429430def rouche(q, F, ymax, R=0.2, N=24, verbose=True, pool=1200, secs=None):431 base = build(q, F, dps=40, tol=mp.mpf(10) ** -24)432 fw = ZCache(build(q, F, dps=25, tol=mp.mpf(10) ** -10))433 L = mp.log(q)434 per = 2 * mp.pi / L435 depths, P = [], base.P436 while len(F) ** P <= pool and len(depths) < 3:437 depths.append(P)438 P += 1439 ds = {base.P: base}440 rhos = [mp.mpf(x) / 400 for x in range(6, 181)]441 r2s = [0.22, 0.25, 0.28, 0.31, 0.34, 0.37, 0.4, 0.45, 0.5, 0.55, 0.6, 0.7, 0.8, 0.9]442 if verbose:443 line(" ", tag(q, F), "alpha", mp.nstr(alpha_of(q, F), 10), "sample circle R",444 R, "samples", N, "peel depths tried", depths)445 t0 = time.time()446 out = []447 for j in range(0, int(mp.ceil(mp.mpf(ymax) / per)) + 1):448 old = mp.mp.dps449 mp.mp.dps = base.dps + 40 + int(abs(j) * 2)450 s0 = mp.log(base.k) / L + 2 * mp.pi * 1j * j / L451 Z0c, e0 = zwork(base, s0)452 mp.mp.dps = old453 if abs(Z0c) < 1000 * e0 or abs(Z0c) < mp.mpf(10) ** -15:454 if verbose:455 line(" j", j, "residue null, no tooth to certify")456 out.append({"j": j, "null": True, "cert": False, "skip": False})457 continue458 if secs is not None and time.time() - t0 > secs:459 if verbose:460 line(" j", j, "SKIPPED, design time budget spent")461 out.append({"j": j, "null": False, "cert": False, "skip": True})462 continue463 best, used = None, None464 for P in depths:465 if P not in ds:466 ds[P] = build(q, F, dps=40, tol=mp.mpf(10) ** -24, P=P)467 d = ds[P]468 mp.mp.dps = d.dps + 40 + int(abs(j) * 2)469 z0, ee = zwork(d, s0)470 samp = [zwork(d, s0 + mp.mpf(R) * mp.exp(2 * mp.pi * 1j * mp.mpf(m) / N))471 for m in range(N)]472 b, U0 = certify(d, s0, z0, ee, samp, mp.mpf(R), N, rhos, r2s)473 mp.mp.dps = old474 if b is not None and (best is None or b[0] > best[0]):475 best, used = b, P476 if b is not None and b[0] > 0:477 break478 n45 = disc(fw, s0, mp.mpf("0.45"))[0]479 rec = {"j": j, "null": False, "skip": False, "n45": n45, "best": best,480 "P": used, "cert": best is not None and best[0] > 0}481 if rec["cert"]:482 rec["nrho"] = disc(fw, s0, best[1])[0]483 out.append(rec)484 if not verbose:485 continue486 if best is None:487 line(" j", j, "no admissible radius at any depth", depths,488 "zeros in 0.45", mp.nstr(n45, 6))489 else:490 marg, rho, R2, minM, B, L1, U0b = best491 line(" j", j, "CERTIFIED" if marg > 0 else "failed", "peel P", used,492 "margin", mp.nstr(marg, 8), "rho", mp.nstr(rho, 5), "R2", mp.nstr(R2, 4),493 "min|model|", mp.nstr(minM, 6), "B_T", mp.nstr(B, 6),494 "|Z_1| low", mp.nstr(L1, 6), "|Z_0| up", mp.nstr(U0b, 6),495 "|Z_0/Z_1|", mp.nstr(U0b / L1, 6), "samples", N, "zeros in rho",496 "-" if marg <= 0 else mp.nstr(rec["nrho"], 6),497 "zeros in 0.45", mp.nstr(n45, 6))498 return out499500501def main():502 argv = sys.argv[1:]503 verb = argv[0] if argv else "shadow"504 ymax = float(argv[1]) if len(argv) > 1 else 30.0505 sets = {"q3": Q3, "q4": Q4, "q5": Q5, "half": HALF, "wide": WIDE, "ctl": CONTROL,506 "all": Q3 + Q4 + Q5 + HALF + CONTROL + WIDE,507 "rou": Q3 + Q5 + HALF + WIDE}508 which = sets[argv[2]] if len(argv) > 2 else sets["all"]509 t0 = time.time()510 if verb == "rouche":511 line("ROUCHE certificate: exactly one zero of Z, hence of zeta_F, in abs(u) < rho")512 tot = {"cert": 0, "fail": 0, "one": 0, "not one": 0, "null": 0, "skip": 0}513 for q, F in which:514 for rec in rouche(q, F, ymax, secs=300 if q == 10 else 105):515 if rec.get("null"):516 tot["null"] += 1517 continue518 if rec.get("skip"):519 tot["skip"] += 1520 continue521 tot["cert" if rec["cert"] else "fail"] += 1522 if rec["cert"]:523 tot["one" if abs(rec["nrho"] - 1) < mp.mpf("0.2") else "not one"] += 1524 line("CERTIFIED", tot["cert"], "FAILED", tot["fail"], "null poles", tot["null"],525 "skipped for budget", tot["skip"],526 "certified discs whose argument principle count is one", tot["one"],527 "certified discs disagreeing", tot["not one"])528 elif verb == "shadow":529 line("SHADOW derived from the residue and the regular part, no fit")530 for q, F in which:531 comb(q, F, ymax)532 elif verb == "census":533 line("CENSUS zeros of the cofactor Z, one strip of width 4, every design deep")534 got = {}535 for q, F in which:536 c = comb(q, F, ymax, verbose=True)537 zs, tot, res = census(q, F, ymax, int(ymax / 2), cmb=c, deep=True)538 nul = sum(1 for x in c if x["null"] and x["j"] > 0 and mp.im(x["s0"]) <= ymax)539 done = len(zs) == int(mp.nint(tot)) - nul540 got[(q, F)] = [(z, k, done) for z, k in zs]541 laws(got)542 line("seconds", round(time.time() - t0, 1))543544if __name__ == "__main__":545 main()