transport_census.py

1import os2import sys3import time45import mpmath as mp67HERE = os.path.dirname(os.path.abspath(__file__))8sys.path.insert(0, os.path.join(HERE, "..", "design-zeta"))9sys.path.insert(0, os.path.join(HERE, "..", "zeta-locus"))1011import design_zeta as dz12import zeta_locus as zl1314# DESIGNS1516LOCUS = zl.Q3 + zl.Q4 + zl.Q5 + zl.HALF + zl.CONTROL + zl.WIDE17SHADOW = [(20, tuple(range(19))), (50, tuple(range(49)))]18RUNGS = [(5, (0, 1, 2, 3)), (10, tuple(range(8)))]19FULL = [(2, (0, 1)), (3, (0, 1, 2)), (4, (0, 1, 2, 3))]20SETS = {"locus": LOCUS, "shadow": SHADOW, "ctl": FULL, "rungs": RUNGS,21        "b20": SHADOW[:1], "b50": SHADOW[1:], "all": LOCUS + RUNGS}2223HALF = mp.mpf("5e-5")24EDGE = mp.mpf(10) ** -625YMAX = mp.mpf(40)2627def line(*a):28    print(" ".join(str(x) for x in a))2930def down(x, n=7):31    return mp.nstr(mp.floor(x * mp.mpf(10) ** n) / mp.mpf(10) ** n, n + 2)3233def build(q, F):34    return zl.build(q, F, dps=25, tol=mp.mpf(10) ** -10)3536# THE ZERO FREE EDGE3738def edge(d):39    amin = min(a for a in d.F if a != 0)40    s = d.alpha + mp.mpf("0.05")41    while s < mp.mpf(14):42        v, e = d.zeta(s)43        if mp.power(amin, s) * (v + e) < 2:44            return s, v, e45        s += mp.mpf("0.05")46    return None, None, None4748# THE POLES4950def poledist(d, x0, x1, y0, y1):51    per = 2 * mp.pi / mp.log(d.q)52    best = mp.inf53    for m in range(0, 3):54        xp = d.alpha - m55        for j in range(int(mp.floor(y0 / per)) - 1, int(mp.ceil(y1 / per)) + 2):56            yp = j * per57            dx = max(x0 - xp, xp - x1, mp.mpf(0))58            dy = max(y0 - yp, yp - y1, mp.mpf(0))59            best = min(best, mp.sqrt(dx * dx + dy * dy))60    return best6162# THE SWEEP6364def sweep(q, F, ymax, nsub):65    d = build(q, F)66    dp = zl.build(q, F, dps=32, tol=mp.mpf(10) ** -18)67    s1, v1, e1 = edge(d)68    x0, x1 = d.alpha + EDGE, s169    f = dz.Cache(d, cof=True)70    t0 = time.time()71    res = dz.strip(f, x0, x1, mp.mpf("0.02"), ymax, nsub, 20)72    tot = sum(r[2] for r in res)73    mx = max(r[3] for r in res)74    zs = []75    for a, b, w, m in res:76        want = int(mp.nint(abs(w)))77        if want == 0:78            continue79        hit = 080        for z0 in zl.seek(f, x0, x1, a, b, 2 * want + 8, 15, 11):81            if hit >= want:82                break83            z, val = zl.polish(d, z0, (b - a) / 40, x1 - x0, dp)84            if z is None or mp.re(z) < x0 or mp.re(z) > x1:85                continue86            if mp.im(z) < a - mp.mpf("0.5") or mp.im(z) > b + mp.mpf("0.5"):87                continue88            if any(abs(z - y) < mp.mpf("1e-9") for y in zs):89                continue90            hit += 191            zs.append(z)92    return {"d": d, "s1": s1, "v1": v1, "e1": e1, "x0": x0, "x1": x1, "res": res,93            "n": int(mp.nint(tot)), "wind": tot, "mx": mx, "zs": zs,94            "evals": f.n, "sec": time.time() - t0}9596# THE BOX9798def certify(q, F, z, half=HALF, n=16):99    d = build(q, F)100    f = dz.Cache(d)101    x0, x1 = mp.re(z) - half, mp.re(z) + half102    y0, y1 = mp.im(z) - half, mp.im(z) + half103    w, mx = dz.box_phase(f, x0, x1, y0, y1, n)104    lo = min(abs(v) for v in f.m.values())105    return {"x0": x0, "x1": x1, "y0": y0, "y1": y1, "w": w, "mx": mx, "lo": lo,106            "e": f.emax, "pole": poledist(d, x0, x1, y0, y1), "alpha": d.alpha}107108VERBS = {}109110# CENSUS111112def census(which=None, ymax=YMAX):113    which = which or LOCUS114    line("TRANSPORT CENSUS: zeros of zeta_F right of alpha, and the Mertens exponent")115    line("of the design's own Mobius nu_F they force through the transport theorem")116    line("sigma1 is a proved zero free edge: a_min^sigma zeta_F(sigma) < 2 there, so the")117    line("census over alpha < Re s < sigma1 is complete in Re and cut only in height;")118    line("the count is the winding and is exact, the rightmost is the rightmost LOCATED")119    got = {}120    for q, F in which:121        r = sweep(q, F, ymax, int(ymax / 2))122        d = r["d"]123        zs = sorted(r["zs"], key=lambda z: -mp.re(z))124        line(" ", zl.tag(q, F), "alpha", mp.nstr(d.alpha, 10), "k/q",125             mp.nstr(mp.mpf(d.k) / q, 8), "sigma1", mp.nstr(r["s1"], 6),126             "a_min^s zeta_F", mp.nstr(mp.power(min(a for a in d.F if a != 0), r["s1"])127                                       * r["v1"], 8), "bound", mp.nstr(r["e1"], 3))128        line("    strip Re", mp.nstr(r["x0"], 10), mp.nstr(r["x1"], 8), "Im to",129             mp.nstr(ymax, 6), "Z winding", mp.nstr(r["wind"], 8), "zeros",130             r["n"], "located", len(zs), "complete", len(zs) == r["n"],131             "max phase step", mp.nstr(r["mx"], 4), "evals", r["evals"],132             "sec", round(r["sec"], 1))133        if not zs:134            line("    NO ZERO RIGHT OF alpha to height", mp.nstr(ymax, 6),135                 "- no transport bound on theta(nu_F) from this census")136            got[(q, F)] = None137            continue138        for z in zs:139            line("    ZERO", zl.tag(q, F), "Re", mp.nstr(mp.re(z), 12), "Im",140                 mp.nstr(mp.im(z), 12), "Re-alpha", mp.nstr(mp.re(z) - d.alpha, 10))141        z = zs[0]142        b = certify(q, F, z)143        line("    BOX Re [", mp.nstr(b["x0"], 12), ",", mp.nstr(b["x1"], 12), "] Im [",144             mp.nstr(b["y0"], 12), ",", mp.nstr(b["y1"], 12), "] winding",145             mp.nstr(b["w"], 8), "step", mp.nstr(b["mx"], 4), "min |zeta_F|",146             mp.nstr(b["lo"], 4), "engine bound", mp.nstr(b["e"], 4),147             "pole distance", mp.nstr(b["pole"], 6))148        ok = abs(b["w"] - 1) < mp.mpf("0.05") and b["lo"] > b["e"]149        line("    THETA", zl.tag(q, F), "theta(nu_F) >=", down(b["x0"]), "alpha",150             mp.nstr(d.alpha, 10), "gain", down(b["x0"] - d.alpha), "above 1",151             b["x0"] > 1, "certified", ok)152        got[(q, F)] = (z, b, len(zs), r["n"])153    return got154155VERBS["census"] = census156157# LAW158159def nullbox(q, F, j, half=HALF):160    d = build(q, F)161    s0 = d.alpha + 2 * mp.pi * 1j * j / mp.log(q)162    out = []163    for cof in (True, False):164        f = dz.Cache(d, cof=cof)165        w, mx = dz.box_phase(f, mp.re(s0) - half, mp.re(s0) + half,166                             mp.im(s0) - half, mp.im(s0) + half, 16)167        out.append((w, mx, min(abs(v) for v in f.m.values()), f.emax))168    return s0, out169170LANDED = [(3, (0, 1), "0.72074", "0.72084", "28.60563", "28.60573", 1),171          (10, tuple(range(9)), "1.00150", "1.00168", "2.73915", "2.73925", 1),172          (10, tuple(range(9)), "0.99900", "1.00050", "2.73810", "2.74030", 0)]173174def repro():175    line("REPRODUCTION of the two landed boxes of lab/py/mrly-pairing and its control,")176    line("same edges, same argument principle, evaluated here")177    for q, F, x0, x1, y0, y1, want in LANDED:178        d = build(q, F)179        f = dz.Cache(d)180        w, mx = dz.box_phase(f, mp.mpf(x0), mp.mpf(x1), mp.mpf(y0), mp.mpf(y1), 16)181        lo = min(abs(v) for v in f.m.values())182        line(" ", zl.tag(q, F), "Re [", x0, ",", x1, "] Im [", y0, ",", y1,183             "] winding", mp.nstr(w, 8), "expected", want, "step", mp.nstr(mx, 4),184             "min |zeta_F|", mp.nstr(lo, 4), "engine bound", mp.nstr(f.emax, 4),185             "pole distance", mp.nstr(poledist(d, mp.mpf(x0), mp.mpf(x1),186                                               mp.mpf(y0), mp.mpf(y1)), 6),187             "reproduced", abs(w - want) < mp.mpf("0.05"))188189def law(which=None, ymax=YMAX):190    which = which or LOCUS191    got = census(which, ymax)192    line("LAW: the columns, no fit")193    line("  DESIGN alpha k/q sigma1 zeros-right-of-alpha rightmost-Re gain=Re-alpha above-1 1-in-F")194    rows = []195    for q, F in which:196        g = got[(q, F)]197        d = build(q, F)198        one = 1 in F199        if g is None:200            rows.append((q, F, d.alpha, mp.mpf(d.k) / q, 0, None, None, one))201            continue202        z, b, nloc, ntot = g203        rows.append((q, F, d.alpha, mp.mpf(d.k) / q, ntot, mp.re(z),204                     mp.re(z) - d.alpha, one))205    for r in sorted(rows, key=lambda t: (t[3], t[2])):206        q, F, a, kq, n, re, gain, one = r207        line("  ROW", zl.tag(q, F), "alpha", mp.nstr(a, 10), "k/q", mp.nstr(kq, 8),208             "zeros", n, "rightmost located", "-" if re is None else mp.nstr(re, 12),209             "gain", "-" if gain is None else mp.nstr(gain, 10),210             "above 1", "-" if re is None else (re > 1), "1 in F", one)211    line("  BY k/q, then by alpha: the gain column above is printed in k/q order and")212    line("  again in alpha order below; equal-key pairs with unequal gain kill a law")213    for r in sorted(rows, key=lambda t: (t[2], t[3])):214        q, F, a, kq, n, re, gain, one = r215        line("  ALPHAORDER", zl.tag(q, F), "alpha", mp.nstr(a, 10), "k/q",216             mp.nstr(kq, 8), "gain", "-" if gain is None else mp.nstr(gain, 10))217    keys = {}218    for r in rows:219        keys.setdefault((mp.nstr(r[2], 10), mp.nstr(r[3], 8)), []).append(r)220    for key, fam in sorted(keys.items()):221        if len(fam) < 2:222            continue223        gs = [r[6] for r in fam if r[6] is not None]224        line("  TIE alpha", key[0], "k/q", key[1], "designs",225             " ".join(zl.tag(r[0], r[1]) for r in fam), "gains",226             " ".join("-" if r[6] is None else mp.nstr(r[6], 8) for r in fam),227             "spread", "-" if len(gs) < 2 else mp.nstr(max(gs) - min(gs), 8))228    ab = [r for r in rows if r[5] is not None and r[5] > 1]229    line("  ABOVE ONE: designs whose rightmost censused zero has Re rho > 1, so the")230    line("  Mertens of nu_F outruns x itself:", " ".join(zl.tag(r[0], r[1]) for r in ab),231         "| their k/q", " ".join(mp.nstr(r[3], 8) for r in ab),232         "| their alpha", " ".join(mp.nstr(r[2], 10) for r in ab))233    line("  COROLLARY: 1 in F and a certified zero rho with Re rho > alpha give")234    line("  sigma_c(N_F) >= Re rho > alpha >= alpha/2, so sum_(n <= x) nu_F(n) is not")235    line("  O(x^(alpha/2 + eps)) and not O(x^(alpha - eps)): no square-root-shaped bound")236    for r in sorted(rows, key=lambda t: t[2]):237        q, F, a, kq, n, re, gain, one = r238        if re is None:239            line("  NOBOUND", zl.tag(q, F), "no zero right of alpha to height",240                 mp.nstr(ymax, 6))241        elif not one:242            line("  NOINVERSE", zl.tag(q, F), "1 not in F, so 1 not in S_F and the",243                 "Dirichlet inverse nu_F does not exist; the zero stands, the",244                 "transport does not")245        else:246            b = got[(q, F)][1]247            line("  BOUND", zl.tag(q, F), "theta(nu_F) >=", down(b["x0"]), ">",248                 "alpha", mp.nstr(a, 10), ">= alpha/2", mp.nstr(a / 2, 10),249                 "square-root shape refuted for nu_F", True)250    line("  CONTROLS: the cofactor only teeth of the full sets at Re s = alpha = 1")251    for q, F in FULL:252        s0, out = nullbox(q, F, 1)253        (wz, mz, lz, ez), (wf, mf, lf, ef) = out254        line("  NULLBOX", zl.tag(q, F), "pole", mp.nstr(s0, 12), "winding of Z",255             mp.nstr(wz, 8), "winding of zeta_F", mp.nstr(wf, 8), "min |Z|",256             mp.nstr(lz, 4), "min |zeta_F|", mp.nstr(lf, 4), "steps",257             mp.nstr(mz, 4), mp.nstr(mf, 4))258    repro()259260VERBS["law"] = law261262def main():263    argv = sys.argv[1:]264    verb = argv[0] if argv else "census"265    ymax = mp.mpf(argv[1]) if len(argv) > 1 else YMAX266    which = SETS[argv[2]] if len(argv) > 2 else LOCUS267    t0 = time.time()268    VERBS[verb](which, ymax)269    line("seconds", round(time.time() - t0, 1))270271if __name__ == "__main__":272    main()