euler.py
21.1 kB · python · 630 lines
1import sys2from math import gcd34# DIGITS56def dmask(n, q):7 m = 08 while n:9 m |= 1 << (n % q)10 n //= q11 return m1213def in_set(n, q, F):14 return n >= 1 and (dmask(n, q) & ~F) == 01516def repunit(q, L):17 return (q ** L - 1) // (q - 1)1819def show(F, q):20 return "".join(str(d) for d in range(q) if (F >> d) & 1)2122# WALL2324def theorem_witness(q, F):25 full = (1 << q) - 126 if not (F >> 1) & 1:27 return ("unit", None, None)28 if F == full:29 return ("full", None, None)30 missing = [c for c in range(q) if not (F >> c) & 1]31 high = [c for c in missing if c >= 2]32 if high:33 c = min(high)34 return ("repunit", repunit(q, c), repunit(q, c + 1))35 if q % 2 == 1:36 return ("odd", 2, (q * q + 1) // 2)37 return ("even", q * q - 1, q * q + 1)3839def breaks(q, F, m, n):40 if gcd(m, n) != 1 or m < 2 or n < 2:41 return False42 lhs = 1 if in_set(m * n, q, F) else 043 rhs = (1 if in_set(m, q, F) else 0) * (1 if in_set(n, q, F) else 0)44 return lhs != rhs4546def least_witness(q, F, pairs, masks):47 keep = ~F48 for m, n in pairs:49 a = 1 if (masks[m] & keep) == 0 else 050 b = 1 if (masks[n] & keep) == 0 else 051 c = 1 if (masks[m * n] & keep) == 0 else 052 if c != a * b:53 return (m, n)54 return None5556def coprime_pairs(bound):57 out = []58 for m in range(2, bound):59 for n in range(m + 1, bound // m + 1):60 if gcd(m, n) == 1:61 out.append((m, n))62 out.sort(key=lambda p: (p[0] * p[1], p[0]))63 return out6465def wall(qmax=12, bound=4000):66 pairs = coprime_pairs(bound)67 print("WALL: is 1_(S_F) multiplicative")68 print("q sets unit full witnessed kinds")69 total = {"unit": 0, "full": 0, "repunit": 0, "odd": 0, "even": 0}70 rows = []71 for q in range(2, qmax + 1):72 masks = [0] * (bound + 1)73 for n in range(1, bound + 1):74 masks[n] = dmask(n, q)75 kinds = {"unit": 0, "full": 0, "repunit": 0, "odd": 0, "even": 0}76 unfound = 077 for F in range(1, 1 << q):78 kind, m, n = theorem_witness(q, F)79 kinds[kind] += 180 total[kind] += 181 if kind in ("unit", "full"):82 continue83 assert breaks(q, F, m, n), (q, F, kind, m, n)84 lw = least_witness(q, F, pairs, masks)85 if lw is None:86 unfound += 187 print(" no witness below %d at q %d F {%s}" % (bound, q, show(F, q)))88 else:89 rows.append((q, F, kind, m, n, lw[0], lw[1]))90 print("%2d %5d %5d %5d %10d repunit %d odd %d even %d unfound %d"91 % (q, (1 << q) - 1, kinds["unit"], kinds["full"],92 kinds["repunit"] + kinds["odd"] + kinds["even"],93 kinds["repunit"], kinds["odd"], kinds["even"], unfound))94 hardest = max(rows, key=lambda r: r[5] * r[6])95 print("hardest least witness: q %d F {%s} pair (%d, %d) product %d"96 % (hardest[0], show(hardest[1], hardest[0]), hardest[5], hardest[6], hardest[5] * hardest[6]))97 print("totals unit %d full %d repunit %d odd %d even %d"98 % (total["unit"], total["full"], total["repunit"], total["odd"], total["even"]))99 print()100 print("least witness, every F with 1 in F and F not full, q <= 5")101 print("q F kind theorem pair least pair product")102 for (q, F, kind, m, n, a, b) in rows:103 if q > 5:104 continue105 if not (F >> 1) & 1:106 continue107 print("%d {%-8s} %-8s (%d, %d)%s(%d, %d) %s= %d"108 % (q, show(F, q), kind, m, n, " " * max(1, 16 - len("(%d, %d)" % (m, n))),109 a, b, " " * max(1, 12 - len("(%d, %d)" % (a, b))), a * b))110111VERBS = {"wall": wall}112113# LERCH114115def lerch_table(s, terms):116 from mpmath import mp, zeta, gamma, factorial117 return ([zeta(s - j) / factorial(j) for j in range(terms)], gamma(1 - s))118119def lerch(s, t, tab):120 from mpmath import mp, mpf, pi, j as I121 zs, g = tab122 t = mpf(t)123 if t <= mpf(1) / 2:124 mu = -2 * pi * I * t125 head = g * (2 * pi * I * t) ** (s - 1)126 else:127 u = 1 - t128 mu = 2 * pi * I * u129 head = g * (-2 * pi * I * u) ** (s - 1)130 acc = mp.mpc(0)131 p = mp.mpc(1)132 for c in zs:133 acc += c * p134 p *= mu135 return head + acc136137def gl_nodes(n):138 import numpy139 from mpmath import mpf140 x, w = numpy.polynomial.legendre.leggauss(n)141 return [mpf(float(a)) for a in x], [mpf(float(a)) for a in w]142143def design_sum(q, F, L, s):144 from mpmath import mp145 digs = [d for d in range(q) if (F >> d) & 1]146 vals = [0]147 for i in range(L):148 p = q ** i149 vals = [v + d * p for v in vals for d in digs]150 return mp.fsum([mp.mpf(n) ** (-s) for n in vals if n >= 1]), len(vals)151152def gseq(q, F, L, t):153 from mpmath import mp, pi, j as I154 digs = [d for d in range(q) if (F >> d) & 1]155 out = mp.mpc(1)156 for i in range(L):157 x = t * (q ** i)158 x = x - mp.floor(x)159 w = mp.exp(2 * pi * I * x)160 p = mp.mpc(1)161 acc = mp.mpc(0)162 for d in range(q):163 if (F >> d) & 1:164 acc += p165 p *= w166 out *= acc167 return out168169def cell_integral(q, F, L, s, tab, nodes, weights, lo, hi):170 from mpmath import mp171 half = (hi - lo) / 2172 mid = (hi + lo) / 2173 acc = mp.mpc(0)174 for x, w in zip(nodes, weights):175 t = mid + half * x176 acc += w * gseq(q, F, L, t) * lerch(s, t, tab)177 return acc * half178179def position_integral(q, F, L, s, npts=14, terms=80, grade=64):180 from mpmath import mp181 tab = lerch_table(s, terms)182 nodes, weights = gl_nodes(npts)183 N = q ** L184 h = mp.mpf(1) / N185 total = mp.mpc(0)186 for a in range(1, N - 1):187 total += cell_integral(q, F, L, s, tab, nodes, weights, a * h, (a + 1) * h)188 for j in range(grade):189 lo = h / (2 ** (j + 1))190 hi = h / (2 ** j)191 total += cell_integral(q, F, L, s, tab, nodes, weights, lo, hi)192 total += cell_integral(q, F, L, s, tab, nodes, weights, 1 - hi, 1 - lo)193 return total194195def position(qmax=None):196 from mpmath import mp, mpc, polylog, exp, pi, j as I197 mp.dps = 25198 print("LERCH EVALUATOR against mpmath polylog")199 for s in [mpc(3, 3) + mpc("0.3"), mpc("2.7", "1.9")]:200 tab = lerch_table(s, 80)201 for t in ["0.03125", "0.25", "0.4", "0.6", "0.9"]:202 a = lerch(s, mp.mpf(t), tab)203 b = polylog(s, exp(-2 * pi * I * mp.mpf(t)))204 print(" s %s t %s err %s" % (mp.nstr(s, 6), t, mp.nstr(abs(a - b), 3)))205 print()206 print("POSITION PRODUCT: sum over D_L of n^(-s) against the integral")207 print("q F L k^L s direct integral err")208 cases = [(10, 3), (10, 4), (3, 3), (3, 4), (3, 5)]209 for q, L in cases:210 F = ((1 << q) - 1) & ~(1 << (q - 1)) if q == 10 else 0b011211 for s in [mpc("3.3", 0), mpc("2.7", "1.9")]:212 direct, kL = design_sum(q, F, L, s)213 integ = position_integral(q, F, L, s)214 print("%2d %-8s %d %6d %-14s %-26s %-26s %s"215 % (q, show(F, q), L, kL, mp.nstr(s, 5), mp.nstr(direct, 14),216 mp.nstr(integ, 14), mp.nstr(abs(direct - integ), 3)))217218VERBS["position"] = position219220# PAIR221222def mobius_sieve(N):223 mu = [1] * (N + 1)224 primes = []225 comp = [False] * (N + 1)226 for i in range(2, N + 1):227 if not comp[i]:228 primes.append(i)229 mu[i] = -1230 for p in primes:231 if i * p > N:232 break233 comp[i * p] = True234 if i % p == 0:235 mu[i * p] = 0236 break237 mu[i * p] = -mu[i]238 return mu239240def pair_coeff(n, q, F, mu, divs, masks):241 keep = ~F242 tot = 0243 for d in divs[n]:244 e = n // d245 if (masks[d] & keep) == 0 and (masks[e] & keep) == 0:246 tot += mu[e]247 return tot248249def pair(N=4000):250 mu = mobius_sieve(N)251 divs = [[] for _ in range(N + 1)]252 for d in range(1, N + 1):253 for m in range(d, N + 1, d):254 divs[m].append(d)255 print("ZETA_F TIMES M_F IS NOT 1")256 print("coefficient c(n) = sum over a b = n with a, b in S_F of mu(b)")257 print("q F k first n > 1 with c(n) nonzero c(n)")258 cases = []259 for q in range(2, 9):260 for F in range(1, 1 << q):261 if (F >> 1) & 1:262 cases.append((q, F))263 for f in [((1 << 10) - 1) & ~(1 << 9), ((1 << 10) - 1) & ~1, (1 << 10) - 1]:264 cases.append((10, f))265 shown = {(3, 0b011), (3, 0b101), (3, 0b110), (4, 0b0011), (5, 0b00011),266 (10, ((1 << 10) - 1) & ~(1 << 9)), (10, ((1 << 10) - 1) & ~1),267 (3, 0b111), (10, (1 << 10) - 1)}268 masks_by_q = {}269 firsts = []270 for q, F in cases:271 if q not in masks_by_q:272 masks_by_q[q] = [0] + [dmask(n, q) for n in range(1, N + 1)]273 masks = masks_by_q[q]274 assert pair_coeff(1, q, F, mu, divs, masks) == 1275 hit = None276 for n in range(2, N + 1):277 c = pair_coeff(n, q, F, mu, divs, masks)278 if c != 0:279 hit = (n, c)280 break281 firsts.append((q, F, hit))282 if (q, F) in shown:283 print("%2d %-10s %2d %-30s %s"284 % (q, show(F, q), bin(F).count("1"),285 "none below %d" % N if hit is None else str(hit[0]),286 "-" if hit is None else str(hit[1])))287 full = [(q, F, h) for (q, F, h) in firsts if F == (1 << q) - 1]288 part = [(q, F, h) for (q, F, h) in firsts if F != (1 << q) - 1]289 print("full digit sets tested %d, all with c(n) = 0 for 1 < n <= %d: %s"290 % (len(full), N, all(h is None for (_, _, h) in full)))291 print("proper sets tested %d, all with a nonzero c(n): %s, largest first witness %d"292 % (len(part), all(h is not None for (_, _, h) in part),293 max(h[0] for (_, _, h) in part)))294295VERBS["pair"] = pair296297# WORD298299def mu_int(n):300 r = 1301 d = 2302 m = n303 while d * d <= m:304 if m % d == 0:305 m //= d306 if m % d == 0:307 return 0308 r = -r309 d += 1310 if m > 1:311 r = -r312 return r313314def lyndon(k, L):315 t = 0316 for d in range(1, L + 1):317 if L % d == 0:318 t += mu_int(d) * k ** (L // d)319 return t // L320321def series_mul(a, b, n):322 out = [0] * n323 for i, x in enumerate(a):324 if x == 0:325 continue326 for j, y in enumerate(b):327 if i + j >= n:328 break329 out[i + j] += x * y330 return out331332def binom(n, r):333 t = 1334 for i in range(r):335 t = t * (n - i) // (i + 1)336 return t337338def word(order=17):339 print("WORD EULER PRODUCT: 1/(1 - k u) = prod over L of (1 - u^L)^(-c_k(L))")340 print("c_k(L) is the Lyndon count (1/L) sum over d dividing L of mu(d) k^(L/d)")341 print("k c_k(1..10)")342 for k in [2, 3, 4, 9, 10]:343 print("%2d %s" % (k, " ".join(str(lyndon(k, L)) for L in range(1, 11))))344 print()345 print("k product to u^%d equals 1/(1 - k u)" % (order - 1))346 for k in [2, 3, 4, 9, 10]:347 acc = [0] * order348 acc[0] = 1349 for L in range(1, order):350 c = lyndon(k, L)351 f = [0] * order352 j = 0353 while L * j < order:354 f[L * j] = binom(c + j - 1, j)355 j += 1356 acc = series_mul(acc, f, order)357 print("%2d %s" % (k, acc == [k ** i for i in range(order)]))358 print()359 print("word Mobius: 1 - k u has coefficients 1, -k and zero after, so the word")360 print("Mertens is 1 at norm 1 and 1 - k at every norm above, bounded in the norm")361 print("word zeta 1/(1 - k q^(-s)) has no zero and simple poles exactly at")362 print("s = alpha + 2 pi i m / log q with alpha = log_q k, the design pole lattice")363364VERBS["word"] = word365366# CHARACTERS367368def prime_factors(n):369 out = {}370 d = 2371 while d * d <= n:372 while n % d == 0:373 out[d] = out.get(d, 0) + 1374 n //= d375 d += 1376 if n > 1:377 out[n] = out.get(n, 0) + 1378 return out379380def primitive_root(p):381 if p == 2:382 return 1383 fs = list(prime_factors(p - 1))384 g = 2385 while True:386 if all(pow(g, (p - 1) // f, p) != 1 for f in fs):387 return g388 g += 1389390def crt_lift(x, m, Q):391 o = Q // m392 if o == 1:393 return x % Q394 inv = pow(o % m, -1, m)395 return (1 + o * ((x - 1) * inv % m)) % Q396397def cyclic_parts(Q):398 parts = []399 for p, e in prime_factors(Q).items():400 m = p ** e401 if p == 2:402 if e == 1:403 continue404 if e == 2:405 parts.append((crt_lift(3, 4, Q), 2))406 else:407 parts.append((crt_lift(m - 1, m, Q), 2))408 parts.append((crt_lift(5, m, Q), 1 << (e - 2)))409 else:410 g = primitive_root(p)411 if e > 1 and pow(g, p - 1, p * p) == 1:412 g += p413 parts.append((crt_lift(g, m, Q), (p - 1) * p ** (e - 1)))414 return parts415416def char_table(Q):417 import cmath418 parts = cyclic_parts(Q)419 orders = [n for _, n in parts]420 logs = {}421 def walk(i, cur, exps):422 if i == len(parts):423 logs[cur] = tuple(exps)424 return425 g, n = parts[i]426 x = cur427 for a in range(n):428 walk(i + 1, x, exps + [a])429 x = x * g % Q430 walk(0, 1 % Q, [])431 idx = []432 def tuples(i, cur):433 if i == len(orders):434 idx.append(tuple(cur))435 return436 for a in range(orders[i]):437 tuples(i + 1, cur + [a])438 tuples(0, [])439 chars = []440 for a in idx:441 tab = [0j] * Q442 for r, l in logs.items():443 ph = sum(a[i] * l[i] / orders[i] for i in range(len(orders)))444 tab[r] = cmath.exp(2j * cmath.pi * ph)445 chars.append(tab)446 return chars447448# FIBRE449450def fibre(N=3000):451 import cmath452 from math import gcd as G453 mu = mobius_sieve(N)454 print("LERCH-MOBIUS AT A RATIONAL: M(s, a/Q) as inverse L-functions")455 print("M(s, a/Q) = sum over g dividing Q of mu(g) g^(-s) (1/phi(Q/g))")456 print(" times sum over chi mod Q/g of Gauss(chi, a) A_chi(s)")457 print("A_chi(s) = 1/L(s, chi) times prod over p dividing Q not Q/g of (1 - chi(p) p^(-s))^(-1)")458 print("Q a divisors used characters max coefficient error to n = %d Euler step" % N)459 for Q, a in [(3, 1), (9, 1), (9, 2), (27, 5), (5, 2), (10, 3), (10, 7), (100, 21), (7, 1), (8, 3), (12, 5)]:460 lhs = [0j] * (N + 1)461 for n in range(1, N + 1):462 lhs[n] = mu[n] * cmath.exp(-2j * cmath.pi * n * a / Q)463 rhs = [0j] * (N + 1)464 gs = [g for g in range(1, Q + 1) if Q % g == 0 and mu[g] != 0]465 nchars = 0466 euler_err = 0.0467 for g in gs:468 Qp = Q // g469 chars = char_table(Qp) if Qp > 1 else [[1j * 0 + 1]]470 nchars += len(chars)471 phi = len([r for r in range(Qp) if G(r, Qp) == 1]) if Qp > 1 else 1472 extra = [p for p in prime_factors(Q) if Qp % p != 0]473 for chi in chars:474 cv = (lambda m: chi[m % Qp]) if Qp > 1 else (lambda m: 1.0 + 0j)475 gauss = sum((cv(r).conjugate()) * cmath.exp(-2j * cmath.pi * r * a / Qp)476 for r in range(Qp) if G(r, Qp) == 1) if Qp > 1 else 1.0 + 0j477 A = [0j] * (N // g + 1)478 for m in range(1, N // g + 1):479 if G(m, Q) == 1:480 A[m] = mu[m] * cv(m)481 for m in range(1, N // g + 1):482 if A[m] != 0j:483 rhs[g * m] += mu[g] * gauss * A[m] / phi484 B = [0j] * (N // g + 1)485 for m in range(1, N // g + 1):486 B[m] = mu[m] * cv(m)487 S = [0j] * (N // g + 1)488 S[1] = 1.0 + 0j489 for p in extra:490 T = list(S)491 pk = p492 while pk <= N // g:493 for m in range(1, N // g // pk + 1):494 T[m * pk] += S[m] * cv(pk)495 pk *= p496 S = T497 C = [0j] * (N // g + 1)498 for m in range(1, N // g + 1):499 if B[m] != 0j:500 for l in range(1, N // g // m + 1):501 if S[l] != 0j:502 C[m * l] += B[m] * S[l]503 euler_err = max(euler_err, max(abs(C[m] - A[m]) for m in range(1, N // g + 1)))504 err = max(abs(lhs[n] - rhs[n]) for n in range(1, N + 1))505 print("%-4d %-2d %-14s %-11d %-33s %s"506 % (Q, a, ",".join(str(g) for g in gs), nchars, "%.3e" % err, "%.3e" % euler_err))507 print()508 print("t = 0 fibre: M(s, 0) = 1/zeta(s), so the classical RH is one fibre of the family")509510VERBS["fibre"] = fibre511512# BEURLING513514def primes_to(X):515 sieve = bytearray([1]) * (X + 1)516 sieve[0] = sieve[1] = 0517 i = 2518 while i * i <= X:519 if sieve[i]:520 sieve[i * i::i] = bytearray(len(sieve[i * i::i]))521 i += 1522 return [i for i in range(2, X + 1) if sieve[i]]523524def generated(P, X):525 out = [(1, 1)]526 def go(i, val, sign, sqf):527 for j in range(i, len(P)):528 p = P[j]529 if val * p > X:530 break531 v = val * p532 out.append((v, -sign if sqf else 0))533 go(j + 1, v, -sign, sqf)534 w = v * p535 while w <= X:536 out.append((w, 0))537 go(j + 1, w, 0, False)538 w *= p539 go(0, 1, 1, True)540 out.sort()541 return out542543def beurling(X=10 ** 6):544 from math import log545 print("BEURLING SYSTEM ON A DESIGN: N_F is the free semigroup on the primes of S_F")546 print("N_F is not S_F. mu_B is the restriction of mu, M_B(x) = sum of mu over N_F below x")547 allp = primes_to(X)548 cases = [(10, ((1 << 10) - 1) & ~(1 << 9), "base 10 missing 9"),549 (3, 0b011, "base 3 {0,1}"),550 (3, 0b101, "base 3 {0,2}"),551 (10, (1 << 10) - 1, "base 10 full, control")]552 for q, F, name in cases:553 k = bin(F).count("1")554 alpha = log(k) / log(q)555 P = [p for p in allp if in_set(p, q, F)]556 gen = generated(P, X)557 print()558 print("%s k %d alpha %.6f alpha/2 %.6f primes in S_F below %d: %d"559 % (name, k, alpha, alpha / 2, X, len(P)))560 print(" x pi_F(x) N_F(x) M_B(x) max abs M_B log max / log x log max / log N_F log_q x")561 idx = 0562 run = 0563 mx = 0564 rows = []565 checks = []566 j = 2567 while True:568 for r in range(4):569 x = int(round(q ** (j + r / 4.0)))570 if x > X:571 break572 if x >= 10:573 checks.append(x)574 if q ** j > X:575 break576 j += 1577 checks = sorted(set(c for c in checks if c <= X))578 pi_at = 0579 pidx = 0580 for x in checks:581 while idx < len(gen) and gen[idx][0] <= x:582 run += gen[idx][1]583 if abs(run) > mx:584 mx = abs(run)585 idx += 1586 while pidx < len(P) and P[pidx] <= x:587 pidx += 1588 pi_at = pidx589 nf = idx590 e = log(mx) / log(x) if mx > 0 else float("nan")591 e2 = log(mx) / log(nf) if mx > 0 and nf > 1 else float("nan")592 rows.append((log(x) / log(q), e, e2))593 print(" %-12d %-9d %-10d %-9d %-13d %-16.6f %-18.6f %.4f"594 % (x, pi_at, nf, run, mx, e, e2, log(x) / log(q)))595 print(" exponent log max / log x by residue class of log_q x, last four checkpoints each")596 for r in [0.0, 0.25, 0.5, 0.75]:597 cls = [row for row in rows if abs((row[0] % 1.0) - r) < 0.02 or abs((row[0] % 1.0) - r - 1) < 0.02]598 if cls:599 print(" r = %.2f %s" % (r, " ".join("%.4f" % c[1] for c in cls[-4:])))600601VERBS["beurling"] = beurling602603# DUAL604605def dual():606 from mpmath import mp, mpc, mpf, pi, gamma, zeta, exp, j as I607 mp.dps = 30608 print("REFLECTION: Z(s, t) from the Hurwitz formula, DLMF 25.13.3 solved for F(-t, s)")609 print("Z(s, t) = (2 pi)^s Gamma(1-s) / (2 pi i) times")610 print(" [ e^(pi i s/2) zeta(1-s, t) - e^(-pi i s/2) zeta(1-s, 1-t) ]")611 print("s t series value reflection value err")612 for s in [mpc("3.3", 0), mpc("2.7", "1.9"), mpc("0.6", "4.1")]:613 tab = lerch_table(s, 120)614 for t in ["0.125", "0.3", "0.5", "0.77"]:615 t = mpf(t)616 a = lerch(s, t, tab)617 b = ((2 * pi) ** s * gamma(1 - s) / (2 * pi * I)) * (618 exp(pi * I * s / 2) * zeta(1 - s, t) - exp(-pi * I * s / 2) * zeta(1 - s, 1 - t))619 print("%-14s %-9s %-26s %-26s %s"620 % (mp.nstr(s, 5), mp.nstr(t, 4), mp.nstr(a, 14), mp.nstr(b, 14), mp.nstr(abs(a - b), 3)))621 print()622 print("the reflection moves the kernel and not the design measure, so the position")623 print("identity returns a dual integral of the same G_L against zeta(1-s, t) and")624 print("zeta(1-s, 1-t), never a relation between zeta_F(s) and zeta_F(1-s)")625626VERBS["dual"] = dual627628if __name__ == "__main__":629 v = sys.argv[1] if len(sys.argv) > 1 else "wall"630 VERBS[v]()