peel.py

1import json2import os3import sys4from time import perf_counter56from mpmath import mp, mpc, log, sqrt78mp.dps = 409P = 1410SHIFT = 6011CUT = 3012T = [[1, 1], [1, 0]]13GAM = [[0, 0], [1, 0]]14LN2 = log(2)1516low, mid = [], [[], []]17for n in range(1, 1 << P):18    if n & (n << 1):19        continue20    if n.bit_length() < P:21        low.append(n)22    else:23        mid[n & 1].append(n)242526def mv(m, v):27    return [sum(mpc(m[i][j]) * v[j] for j in range(2)) for i in range(2)]282930def solve(x, rhs):31    a = [[1 - x * T[0][0], -x * T[0][1]], [-x * T[1][0], 1 - x * T[1][1]]]32    det = a[0][0] * a[1][1] - a[0][1] * a[1][0]33    return [(a[1][1] * rhs[0] - a[0][1] * rhs[1]) / det,34            (a[0][0] * rhs[1] - a[1][0] * rhs[0]) / det]353637def epoly(w):38    return [sum(mpc(n) ** (-w) for n in mid[u]) for u in range(2)]394041def ladder(s):42    levels = int(mp.ceil(SHIFT - s.real))43    rows = levels + CUT + 244    g = [[mpc(0), mpc(0)] for _ in range(rows)]45    num = None46    for j in range(levels - 1, -1, -1):47        w = s + j48        acc = epoly(w)49        c = mpc(1)50        for l in range(1, CUT + 1):51            c = c * (-w - (l - 1)) / l52            weight = c * mpc(2) ** (-w - l)53            lift = mv(GAM, g[j + l])54            acc = [acc[i] + weight * lift[i] for i in range(2)]55        if j == 0:56            num = acc57        g[j] = solve(mpc(2) ** (-w), acc)58    return g[0], num596061def zeta(s):62    g, _ = ladder(s)63    return sum(mpc(n) ** (-s) for n in low) + g[0] + g[1]646566def cofactor(s):67    _, num = ladder(s)68    x = mpc(2) ** (-s)69    det = 1 - x - x * x70    adj = [[mpc(1), x], [x, 1 - x]]71    return det * sum(mpc(n) ** (-s) for n in low) + sum(mv(adj, num))727374CONTROL = os.path.join(os.path.dirname(os.path.abspath(__file__)), "control.json")757677def control():78    start = perf_counter()79    phi = (1 + sqrt(5)) / 280    alpha = log(phi) / LN281    period = 2 * mp.pi / LN282    half = mp.pi / LN283    rows = []84    for kind, s in [("zeta", mpc(3)), ("zeta", mpc("0.8")), ("Z", mpc(3)), ("Z", mpc(2)),85                    ("Z", mpc("0.8")), ("Z", mpc("-0.95", 20)), ("Z", mpc(0, 30)),86                    ("Z", mpc("-0.693919320380", "9.492981163426")),87                    ("Z", mpc("-0.737611737911", "14.343440171295")),88                    ("Z", mpc("-0.708500086380", "27.443172842651")),89                    ("Z", mpc("0.665353220372", "18.164471712166")),90                    ("Z", mpc("-0.315485550754", "23.209873760385")),91                    ("residue", mpc(alpha, period)), ("residue", mpc(-alpha, half * 5))]:92        v = zeta(s) if kind == "zeta" else cofactor(s) if kind == "Z" else residue(s)93        rows.append({"kind": kind, "re": mp.nstr(s.real, 25), "im": mp.nstr(s.imag, 25),94                     "value": [mp.nstr(v.real, 25), mp.nstr(v.imag, 25)]})95    json.dump({"source": "lab/py/memory-zeta peel.py control", "dps": mp.dps, "peel": P,96               "shift": SHIFT, "cut": CUT, "rows": rows}, open(CONTROL, "w"), indent=1)97    print("control.json", len(rows), "rows,", round(perf_counter() - start, 2), "s")9899100def residue(w0):101    _, num = ladder(w0)102    x = mpc(2) ** (-w0)103    adj = [[mpc(1), x], [x, 1 - x]]104    top = sum(mv(adj, num))105    den = -x * LN2 * (-1 - 2 * x)106    return top / den107108109def main():110    start = perf_counter()111    phi = (1 + sqrt(5)) / 2112    alpha = log(phi) / LN2113    period = 2 * mp.pi / LN2114    half = mp.pi / LN2115    print("code 7, k = 2, dps", mp.dps, "peel", P, "shift", SHIFT, "cut", CUT)116    print("alpha", mp.nstr(alpha, 25))117    for s in [mpc(3), mpc(2), mpc("0.8"), mpc("1.2", 9)]:118        print("zeta", s, mp.nstr(zeta(s), 25))119    for j in range(3):120        print("comb one j", j, mp.nstr(residue(mpc(alpha, period * j)), 25))121    for j in range(3):122        print("comb two j", j, mp.nstr(residue(mpc(-alpha, half * (2 * j + 1))), 25))123    print("runtime", round(perf_counter() - start, 2), "s")124125126if len(sys.argv) > 1 and sys.argv[1] == "control":127    control()128else:129    main()