nyman_beurling.py
6.1 kB · python · 131 lines
1import os2import shutil3import subprocess4import sys5import time67HERE = os.path.dirname(os.path.abspath(__file__))8DATA_DIR = os.path.join("data", os.path.relpath(HERE))9PAIRS = [(1, 1), (1, 2), (2, 3), (3, 4), (5, 7), (1, 9), (4, 9), (13, 27), (10, 81), (40, 81), (6, 10)]10DESIGNS = [(3, (0, 1)), (3, (0, 2)), (3, (1, 2)), (4, (0, 1)), (4, (1, 2, 3)), (5, (0, 1)), (10, tuple(range(9))), (10, (0, 1, 2, 3, 4, 5, 6, 7, 9))]1112ENGINE = r"""13default(parisizemax, 6000000000);14default(realprecision, PREC);15LC = log(2*Pi) - Euler;16CT = Map();17cots(k) = my(c); if(mapisdefined(CT, k, &c), c, c = vector(k - 1, m, cotan(Pi*m/k)); mapput(CT, k, c); c);18VT = Map();19vv(h, k) = {20 my(r = h % k, key, val);21 if(k < 2, return(0));22 key = [r, k];23 if(mapisdefined(VT, key, &val), return(val));24 val = vector(k - 1, m, (m*r) % k) * cots(k)~ / k;25 mapput(VT, key, val);26 val;27}28nu(h, k) = LC/2 * (1/h + 1/k) + (k - h)/(2*h*k) * log(h/k) - Pi/(2*h*k) * (vv(h, k) + vv(k, h));29gram(m, n) = my(g = gcd(m, n)); nu(m/g, n/g)/g;30quad(h, k) = {31 my(pts = vecsort(concat(vector(h+1, j, (j-1)/h), vector(k+1, j, (j-1)/k)), , 8), tot = 0);32 for(i = 1, #pts - 1, tot += intnum(w = pts[i], pts[i+1], frac(k*w)*frac(h*w)*zetahurwitz(2, w)));33 tot/(h*k);34}35grammat(S) = my(n = #S, G = matrix(n, n)); for(i = 1, n, for(j = 1, i, G[i, j] = gram(S[i], S[j]); G[j, i] = G[i, j])); G;36target(S) = vector(#S, i, (log(S[i]) + 1 - Euler)/S[i])~;37cond1(G) = my(Gi = G^(-1), n = #G); vecmax(vector(n, j, vecsum(abs(G[, j])))) * vecmax(vector(n, j, vecsum(abs(Gi[, j]))));38dist(S) = {39 my(G = grammat(S), b = target(S), v = vector(#S, i, 1/S[i])~, x = matsolve(G, b), y = matsolve(G - v*v~, b));40 [1 - b~*x, 1 - b~*y, cond1(G), x, y];41}42design(base, F, N) = select(n -> setminus(Set(digits(n, base)), Set(F)) == [], [1..N]);43resid(S, c) = {44 my(f(u) = sum(i = 1, #S, c[i]*frac(u/S[i])), P = lcm(S), tot = intnum(u = 0, 1, f(u)^2/u^2));45 for(k = 1, P - 1, tot += intnum(u = k, k+1, (1 - f(u))^2/u^2));46 for(k = 0, P - 1, tot += intnum(v = k, k+1, (1 - f(v))^2*zetahurwitz(2, 1 + v/P))/P^2);47 tot;48}49floorq(S, K, tail) = {50 my(E = select(n -> n < K, S), d = #E + 1, L = vector(K - 1, k, log((k+1)/k)));51 my(Q(x) = my(s = x[1], t = if(tail, s^2, 0)); for(k = 1, K - 1, my(B = 1 + sum(i = 1, #E, x[i+1]*floor(k/E[i]))); t += s^2 - 2*s*B*L[k] + B^2/(k*(k+1))); t);52 my(e(i) = vector(d, j, j == i), c0 = Q(vector(d)), g = vector(d, i, (Q(e(i)) - Q(-e(i)))/4)~, H = matrix(d, d));53 for(i = 1, d, H[i, i] = (Q(e(i)) + Q(-e(i)))/2 - c0; for(j = 1, i - 1, H[i, j] = (Q(e(i) + e(j)) - Q(e(i)) - Q(e(j)) + c0)/2; H[j, i] = H[i, j]));54 c0 - g~*matsolve(H, g);55}56qzero(S, N) = my(r = 0); for(n = 1, N, if(issquarefree(n) && !setsearch(Set(S), n), r = n; break)); r;57"""585960def gp(lines, prec=60):61 if not shutil.which("gp"):62 raise SystemExit("PARI is needed: gp is not on PATH")63 os.makedirs(DATA_DIR, exist_ok=True)64 path = os.path.join(DATA_DIR, "run.gp")65 with open(path, "w") as fh:66 fh.write("PREC = {};".format(prec) + ENGINE + "\n".join(lines) + "\nquit;\n")67 done = subprocess.run(["gp", "-q", path], capture_output=True, text=True, stdin=subprocess.DEVNULL)68 if done.returncode != 0:69 raise RuntimeError(done.stderr)70 return [line for line in done.stdout.splitlines() if not line.startswith(" ***")]717273def verb_gram():74 lines = []75 for h, k in PAIRS:76 lines.append('my(a = gram({0}, {1}), b = quad({0}, {1})); print("{0} {1} ", a, " ", b, " ", abs(a - b))'.format(h, k))77 lines.append('print("norm ", gram(1, 1), " ", LC)')78 for S in ("[1, 2, 3]", "[1, 3]", "[1, 2, 3, 4]", "[1, 3, 4]"):79 lines.append('my(S = {0}, r = dist(S)); print("d2 {0} ", r[1], " ", resid(S, r[4]))'.format(S))80 lines.append('print("d2 1 ", dist([1])[1], " ", 1 - (1 - Euler)^2/LC)')81 for row in gp(lines, 40):82 print(row)838485def verb_table(levels="2-6", extra=7):86 lo, hi = (int(v) for v in levels.split("-"))87 for prec in (60, 90):88 lines = ['C = 2 + Euler - log(4*Pi); print("burnol C ", C)']89 for level in range(lo, hi + 1):90 lines.append('my(N = 3^{0}, r = dist([1..N])); print("full {0} ", N, " ", N, " ", r[1], " ", r[2], " ", r[3], " ", C/log(N), " ", r[1]*log(N))'.format(level))91 for level in range(1, max(hi, extra) + 1):92 lines.append('my(N = 3^{0}, D = design(3, [0, 1], N), r = dist(D)); print("design {0} ", N, " ", #D, " ", r[1], " ", r[2], " ", r[3], " ", C/log(N), " ", r[1]*log(N))'.format(level))93 start = time.time()94 rows = gp(lines, prec)95 print("precision", prec, "digits, {:.1f} s".format(time.time() - start))96 for row in rows:97 print(row)9899100def verb_zeros(height=1000):101 lines = ['my(z = lfunzeros(1, {0}), s = 2*sum(i = 1, #z, 1/(1/4 + z[i]^2)), tail = 2*(log({0}/(2*Pi)) + 1)/(2*Pi*{0})); print("zeros ", #z, " partial ", s, " tail ", tail, " sum ", s + tail, " C ", 2 + Euler - log(4*Pi))'.format(height)]102 for row in gp(lines, 30):103 print(row)104105106def verb_floor(kmax=40):107 lines = []108 for base, F in DESIGNS:109 lines.append('my(S = design({0}, {1}, 100000)); print("q0 {0} {1} ", qzero(S, 100000), " ", S[1..12])'.format(base, list(F)))110 lines.append('my(S = design(3, [0, 1], 1000)); for(K = 2, {0}, print("floor 3 [0,1] ", K, " ", floorq(S, K, 1), " ", floorq(S, K, 0)))'.format(kmax))111 lines.append('my(S = [1..1000]); for(K = 2, 12, print("floor full ", K, " ", floorq(S, K, 1), " ", floorq(S, K, 0)))')112 for base, F in DESIGNS[1:]:113 lines.append('my(S = design({0}, {1}, 1000)); foreach([{2}], K, print("floor {0} {1} ", K, " ", floorq(S, K, 1), " ", floorq(S, K, 0)))'.format(base, list(F), kmax))114 for row in gp(lines, 60):115 print(row)116117118if __name__ == "__main__":119 verb = sys.argv[1] if len(sys.argv) > 1 else "table"120 start = time.time()121 if verb == "gram":122 verb_gram()123 elif verb == "table":124 verb_table(*(sys.argv[2:3]))125 elif verb == "floor":126 verb_floor(*(int(v) for v in sys.argv[2:3]))127 elif verb == "zeros":128 verb_zeros(*(int(v) for v in sys.argv[2:3]))129 else:130 raise SystemExit("verbs: gram, table, floor, zeros")131 print("runtime {:.1f} s".format(time.time() - start))