gaussian_zeta.py

1import math2import time3from fractions import Fraction45import numpy as np6from mpmath import mp, mpf, bernoulli, catalan, zeta78DIGITS = 909CORNERS = [(0, 0), (0, 1), (1, 0), (1, 1)]10DIAGONAL = [(0, 0), (1, 1)]11MATE = [(0, 1), (1, 0)]12EXACT_TO = 8013LIMITS = [49, 199, 999]14SHELL_MAX = 100000015RADIUS = 600016STATED = {17    1: "1/24", 2: "1/4", 3: "7/24", 4: "1/4", 5: "7/24", 6: "1/2", 7: "13/24",18    8: "1/8", 9: "1/6", 10: "3/8", 11: "5/12", 12: "3/8", 13: "5/12",19    14: "5/8", 15: "2/3",20}212223def code_of(filled):24    return sum(1 << k for k, cell in enumerate(CORNERS) if cell in filled)252627def fill(filled, n):28    table = np.zeros((2, 2), dtype=bool)29    for i, j in filled:30        table[i, j] = True31    parity = np.arange(n) % 232    return int(np.count_nonzero(table[parity[:, None], parity[None, :]]))333435def fluctuation(filled, n):36    return Fraction(fill(filled, n), n * n) - Fraction(len(filled), 4)373839def pi_machin():40    def arctan_inv(x):41        total = mpf(0)42        term = mpf(1) / x43        k = 044        while term > mpf(10) ** (-DIGITS - 10):45            total += term / (2 * k + 1) * (1 if k % 2 == 0 else -1)46            term /= x * x47            k += 148        return total4950    return 16 * arctan_inv(5) - 4 * arctan_inv(239)515253def odd_zeta(s, m=600, k=14):54    total = sum(mpf(1) / mpf(2 * j - 1) ** s for j in range(1, m))55    a = mpf(2 * m - 1)56    total += a ** (1 - s) / (2 * (s - 1)) + a ** (-s) / 257    for step in range(1, k + 1):58        poch = mpf(1)59        for i in range(2 * step - 1):60            poch *= s + i61        derivative = -(mpf(2) ** (2 * step - 1)) * poch * a ** (-(s + 2 * step - 1))62        total -= bernoulli(2 * step) / mp.factorial(2 * step) * derivative63    return total646566def lambda_part():67    print("DIAGONAL DESIGN, FILLED CORNERS (0,0) AND (1,1), CODE", code_of(DIAGONAL))68    even_ok = True69    odd_ok = True70    for n in range(1, EXACT_TO + 1):71        value = fluctuation(DIAGONAL, n)72        if n % 2 == 0:73            even_ok = even_ok and value == 074        else:75            odd_ok = odd_ok and value == Fraction(1, 2 * n * n)76        if n <= 11:77            print(f"  n={n:2d}  fill {fill(DIAGONAL, n):4d}  fluctuation {value}")78    print(f"  even n <= {EXACT_TO} fluctuation exactly 0: {even_ok}")79    print(f"  odd n <= {EXACT_TO} fluctuation exactly 1/(2n^2): {odd_ok}")80    forms_ok = all(81        fill(DIAGONAL, 2 * m) == 2 * m * m82        and fill(DIAGONAL, 2 * m - 1) == m * m + (m - 1) ** 283        for m in range(1, EXACT_TO // 2 + 1)84    )85    print(f"  fill(2m) = 2m^2 and fill(2m-1) = m^2 + (m-1)^2: {forms_ok}")86    mate_ok = all(87        fluctuation(MATE, n) == (0 if n % 2 == 0 else Fraction(-1, 2 * n * n))88        for n in range(1, 41)89    )90    print(f"  orbit mate code {code_of(MATE)} fluctuation -1/(2n^2) on odd n <= 40: {mate_ok}")91    print()9293    print("EXACT RATIONAL REDUCTION")94    z2 = Fraction(1, 2) * (1 - Fraction(1, 16)) / 9095    z4 = Fraction(1, 2) * (1 - Fraction(1, 64)) / 94596    print(f"  (1/2)(1-2^-4)/90  = {z2}   target 1/192   {z2 == Fraction(1, 192)}")97    print(f"  (1/2)(1-2^-6)/945 = {z4}   target 1/1920  {z4 == Fraction(1, 1920)}")98    print()99100    print(f"CONSTANTS AT {DIGITS} DIGITS")101    mp.dps = DIGITS102    machin = pi_machin()103    print(f"  pi machin          {mp.nstr(machin, 60)}")104    print(f"  pi library         {mp.nstr(mp.pi, 60)}")105    print(f"  pi difference      {mp.nstr(abs(machin - mp.pi), 3)}")106    for s in (4, 6):107        direct = odd_zeta(s)108        identity = (1 - mpf(2) ** (-s)) * zeta(s)109        print(f"  lambda({s}) two routes differ {mp.nstr(abs(direct - identity), 3)}")110    targets = {}111    cases = ((2, 4, 192, "0.50733901580"), (4, 6, 1920, "0.50072353832"))112    for s, power, denom, prefix in cases:113        left = odd_zeta(s + 2) / 2114        right = machin ** power / denom115        targets[s] = right116        text = mp.nstr(left, 50)117        print(f"  Z({s}) = lambda({s + 2})/2  {text}")118        print(f"  pi^{power}/{denom:<5}          {mp.nstr(right, 50)}")119        gap = mp.nstr(abs(left - right), 3)120        print(f"  difference {gap}  page prefix {prefix}: {text.startswith(prefix)}")121    print()122123    print("SERIES SUMMED FROM COUNTED FILLS, NO CLOSED FORM")124    counted = [fluctuation(DIAGONAL, n) for n in range(1, LIMITS[-1] + 1)]125    for s in (2, 4):126        total = mpf(0)127        for n, value in enumerate(counted, start=1):128            if value:129                total += mpf(value.numerator) / value.denominator / mpf(n) ** s130            if n in LIMITS:131                gap = abs(total - targets[s])132                print(f"  s={s}  n <= {n:4d}  partial {mp.nstr(total, 16)}  gap {mp.nstr(gap, 2)}")133    print()134135136def shell_counts(nmax):137    bound = math.isqrt(nmax)138    axis = np.arange(-bound, bound + 1)139    norm = axis[:, None] ** 2 + axis[None, :] ** 2140    cls = (np.abs(axis)[:, None] % 2) * 2 + (np.abs(axis)[None, :] % 2)141    keep = norm <= nmax142    return [np.bincount(norm[keep & (cls == c)], minlength=nmax + 1) for c in range(4)]143144145def chi4_excess(nmax):146    out = np.zeros(nmax + 1, dtype=np.int64)147    for d in range(1, nmax + 1, 2):148        out[d::d] += 1 if d % 4 == 1 else -1149    return out150151152def q_poly(a, t):153    return a[0] * t * t + a[3] * t * (1 - t) + (a[1] + a[2]) * (1 - t) / 2154155156def indicators(code):157    return [(code >> k) & 1 for k in range(4)]158159160def lattice_class_sums(radius):161    axis = np.arange(-radius, radius + 1, dtype=np.float64)162    limit = float(radius) ** 2163    out = [0.0] * 4164    for ri in range(2):165        for rj in range(2):166            rows = axis[np.abs(axis) % 2 == ri]167            cols = axis[np.abs(axis) % 2 == rj]168            squares = cols * cols169            pieces = []170            for i in rows:171                norm = i * i + squares172                norm = norm[(norm <= limit) & (norm > 0)]173                pieces.append(float(np.sum(1.0 / (norm * norm))))174            out[ri * 2 + rj] = math.fsum(pieces)175    return out176177178def gaussian_part():179    print(f"SHELL IDENTITIES ON Z^2, n <= {SHELL_MAX}")180    ee, eo, oe, oo = shell_counts(SHELL_MAX)181    r2 = ee + eo + oe + oo182    print(f"  r2(n) = 4(d1(n) - d3(n)): {bool(np.all(r2[1:] == 4 * chi4_excess(SHELL_MAX)[1:]))}")183    want = np.zeros(SHELL_MAX + 1, dtype=np.int64)184    want[4::4] = r2[1 : SHELL_MAX // 4 + 1]185    print(f"  S_ee(4m) = r2(m), zero elsewhere: {bool(np.all(ee[1:] == want[1:]))}")186    want = np.zeros(SHELL_MAX + 1, dtype=np.int64)187    want[2::2] = (eo + oe)[1 : SHELL_MAX // 2 + 1]188    print(f"  S_oo(2m) = S_mix(m), zero elsewhere: {bool(np.all(oo[1:] == want[1:]))}")189    print(f"  S_eo = S_oe: {bool(np.all(eo == oe))}")190    print("  so S_mix = (1-t) S and Q_c(t) = a_ee t^2 + a_oo t(1-t) + (a_eo + a_oe)(1-t)/2")191    print()192193    print("THE FIFTEEN NONEMPTY DESIGNS AT s = 2, Z_c(2) / (pi^2 G)")194    mp.dps = 50195    for code in range(1, 16):196        value = q_poly(indicators(code), Fraction(1, 4)) * Fraction(2, 3)197        ok = str(value) == STATED[code]198        print(f"  code {code:2d}  {str(value):>5}  page {STATED[code]:>5}  {ok}")199    z7 = mpf(13) / 24 * mp.pi ** 2 * catalan200    print(f"  code 7 (Sierpinski) 13 pi^2 G / 24 = {mp.nstr(z7, 20)}  page 4.8967847822")201    print()202203    print(f"TRUNCATED LATTICE SUMS, RADIUS {RADIUS}, TAIL pi/(4 R^2) PER CLASS")204    sums = lattice_class_sums(RADIUS)205    tail = math.pi / (4.0 * RADIUS * RADIUS)206    base = 4.0 * (math.pi ** 2 / 6.0) * float(catalan)207    worst = 0.0208    for code in range(1, 16):209        a = indicators(code)210        measured = sum(sums[k] for k in range(4) if a[k]) + tail * sum(a)211        predicted = base * float(q_poly(a, Fraction(1, 4)))212        gap = abs(measured - predicted)213        worst = max(worst, gap)214        print(f"  code {code:2d}  predicted {predicted:.12f}  measured {measured:.12f}  gap {gap:.1e}")215    print(f"  worst gap {worst:.1e}")216    print()217218219def main():220    start = time.time()221    lambda_part()222    gaussian_part()223    print(f"DOMAIN  fluctuation n <= {EXACT_TO}, series n <= {LIMITS[-1]}, "224          f"shells n <= {SHELL_MAX}, radius {RADIUS}")225    print(f"WALL {time.time() - start:.1f} s")226227228if __name__ == "__main__":229    main()