eisenstein_visibility.py

1import math2import random3from decimal import Decimal, getcontext4from fractions import Fraction56import numpy as np78PREC = 609EULER_TERMS = 10010BERNOULLI = [11    Fraction(1, 6),12    Fraction(-1, 30),13    Fraction(1, 42),14    Fraction(-1, 30),15    Fraction(5, 66),16    Fraction(-691, 2730),17    Fraction(7, 6),18    Fraction(-3617, 510),19    Fraction(43867, 798),20    Fraction(-174611, 330),21]2223PAIRS = 20000024CUTOFF = 400000025SIDE = 300026TRIALS = 40000027RADIUS = 50028SEED = 20260725293031def atan_inv(x):32    base = Decimal(x)33    term = Decimal(1) / base34    total = term35    square = base * base36    eps = Decimal(10) ** -5837    k = 138    sign = -139    while True:40        term = term / square41        add = Decimal(sign) * term / Decimal(2 * k + 1)42        if abs(add) < eps:43            break44        total += add45        sign = -sign46        k += 147    return total484950def machin_pi():51    return 16 * atan_inv(5) - 4 * atan_inv(239)525354def trigamma(x):55    total = Fraction(0)56    for n in range(EULER_TERMS):57        shifted = n + x58        total += 1 / (shifted * shifted)59    m = EULER_TERMS + x60    total += 1 / m61    total += Fraction(1, 2) / (m * m)62    for index, bern in enumerate(BERNOULLI):63        total += bern / m ** (2 * (index + 1) + 1)64    return total656667def dec(value):68    return Decimal(value.numerator) / Decimal(value.denominator)697071def contfrac(value, terms):72    out = []73    x = value74    for _ in range(terms):75        a = int(x)76        out.append(a)77        frac = x - a78        if frac == 0:79            break80        x = Decimal(1) / frac81    return out828384def denominators(quotients):85    out = []86    previous, current = 0, 187    for index, a in enumerate(quotients):88        if index == 0:89            out.append(current)90            continue91        previous, current = current, a * current + previous92        out.append(current)93    return out949596def cut(value, width):97    return str(value)[:width]9899100def precision_block():101    getcontext().prec = PREC102    pi = machin_pi()103    pi2 = pi * pi104    psi13 = trigamma(Fraction(1, 3))105    psi23 = trigamma(Fraction(2, 3))106    psi14 = trigamma(Fraction(1, 4))107    psi34 = trigamma(Fraction(3, 4))108    l3 = dec(psi13 - psi23) / 9109    catalan = dec(psi14 - psi34) / 16110    z2 = pi2 / 6111    zk3 = z2 * l3112    zki = z2 * catalan113    sqrt3 = Decimal(3).sqrt()114    reflection = dec(psi13 + psi23) - 4 * pi2 / 3115    cl2 = l3 * 3 * sqrt3 / 4116    r1 = l3 * sqrt3 / pi2117    r2 = l3 / pi2118    print("pi           =", cut(pi, 45))119    print("zeta(2)      =", cut(z2, 45))120    print("6/pi^2       =", cut(6 / pi2, 45))121    print("L(2,chi_-3)  =", cut(l3, 45))122    print("Catalan G    =", cut(catalan, 45))123    print("zeta_K3(2)   =", cut(zk3, 45))124    print("1/zeta_K3(2) =", cut(1 / zk3, 45))125    print("zeta_Qi(2)   =", cut(zki, 45))126    print("1/zeta_Qi(2) =", cut(1 / zki, 45))127    print("psi1(1/3)+psi1(2/3)-4pi^2/3 =", reflection)128    print("Cl2(pi/3) implied =", cut(cl2, 30))129    for label, value, terms in (130        ("L3*sqrt3/pi^2", r1, 22),131        ("L3/pi^2      ", r2, 22),132        ("G/pi^2       ", catalan / pi2, 22),133        ("zeta(2)/pi^2 ", z2 / pi2, 8),134    ):135        quotients = contfrac(value, terms)136        print(f"{label} = {cut(value, 30)} CF: {quotients}")137        print(f"{label}   last convergent denominator: {denominators(quotients)[-1]}")138    return float(zk3), float(zki), float(1 / zk3)139140141def dirichlet(pairs):142    m = np.arange(pairs, dtype=np.float64)143    odd = 3 * m + 1144    even = 3 * m + 2145    return math.fsum(1.0 / (odd * odd) - 1.0 / (even * even))146147148def lattice_radius(cutoff):149    return int(2.2 * math.sqrt(cutoff))150151152def lattice_sums(cutoff):153    r = lattice_radius(cutoff)154    b = np.arange(-r, r + 1, dtype=np.int64)155    hexagonal = 0.0156    square = 0.0157    for a in range(-r, r + 1):158        hexnorm = a * a - a * b + b * b159        sqnorm = a * a + b * b160        take = hexnorm[(hexnorm > 0) & (hexnorm <= cutoff)].astype(np.float64)161        hexagonal += float(np.sum(1.0 / (take * take)))162        take = sqnorm[(sqnorm > 0) & (sqnorm <= cutoff)].astype(np.float64)163        square += float(np.sum(1.0 / (take * take)))164    return hexagonal / 6.0, square / 4.0165166167def visible_pairs(side):168    axis = np.arange(1, side + 1, dtype=np.int64)169    count = 0170    for a in range(1, side + 1):171        count += int(np.count_nonzero(np.gcd(a, axis) == 1))172    return count173174175def emul(z, w):176    a, b = z177    c, d = w178    return (a * c - b * d, a * d + b * c - b * d)179180181def enorm(z):182    a, b = z183    return a * a - a * b + b * b184185186def edivmod(z, w):187    c, d = w188    n = enorm(w)189    p = emul(z, (c - d, -d))190    q = (round(p[0] / n), round(p[1] / n))191    s = emul(q, w)192    return q, (z[0] - s[0], z[1] - s[1])193194195def egcd(z, w):196    while enorm(w) != 0:197        z, w = w, edivmod(z, w)[1]198    return z199200201def seeded_sieve(trials, radius):202    rng = random.Random(SEED)203    coprime = 0204    for _ in range(trials):205        z = (rng.randint(-radius, radius), rng.randint(-radius, radius))206        w = (rng.randint(-radius, radius), rng.randint(-radius, radius))207        if enorm(z) == 0 or enorm(w) == 0:208            continue209        if enorm(egcd(z, w)) == 1:210            coprime += 1211    return coprime212213214def phi(n):215    r, m, p = n, n, 2216    while p * p <= m:217        if m % p == 0:218            r -= r // p219            while m % p == 0:220                m //= p221        p += 1222    if m > 1:223        r -= r // m224    return r225226227def fresh(base, level):228    return sum(1 for k in range(1, base**level + 1) if k % base)229230231def sieve_block(zk3, zki, inv_zk3):232    print(f"L(2,chi_-3) direct Dirichlet sum ({PAIRS} pairs) = {dirichlet(PAIRS):.12f}")233    hexagonal, square = lattice_sums(CUTOFF)234    hextail = math.pi / (3 * math.sqrt(3)) / CUTOFF235    sqtail = math.pi / 4 / CUTOFF236    print(237        f"zeta_K3(2) lattice sum (norm<={CUTOFF}) = {hexagonal:.10f}"238        f"  (+tail~{hextail:.2e} -> {hexagonal + hextail:.10f})"239    )240    print(f"zeta_K3(2) expected zeta(2)*L       = {zk3:.10f}")241    print(242        f"zeta_Qi(2) lattice sum (norm<={CUTOFF}) = {square:.10f}"243        f"  (+tail~{sqtail:.2e} -> {square + sqtail:.10f})"244    )245    print(f"zeta_Qi(2) expected zeta(2)*G       = {zki:.10f}")246    seen = visible_pairs(SIDE)247    print(248        f"square-lattice visible density N={SIDE}: {seen / SIDE**2:.6f}"249        f"  vs 6/pi^2 = {6 / math.pi**2:.6f}"250    )251    coprime = seeded_sieve(TRIALS, RADIUS)252    print(253        f"Eisenstein sieve: {coprime}/{TRIALS} coprime = {coprime / TRIALS:.5f}"254        f"  vs 1/zeta_K3(2) = {inv_zk3:.5f}"255    )256    assert enorm(emul((3, 2), (1, 5))) == enorm((3, 2)) * enorm((1, 5))257    print(f"gcd(6,9) in Z[omega] has norm {enorm(egcd((6, 0), (9, 0)))} (expect 9 = N(3))")258    for base in (2, 3, 5):259        for level in range(1, 5):260            new = fresh(base, level)261            assert new == base**level - base ** (level - 1) == phi(base**level)262    print(263        "fresh nodes per level = b^L - b^(L-1) = phi(b^L)"264        " for b=2,3,5, L=1..4: OK; fraction (b-1)/b"265    )266    print(267        f"b=6,L=2: new={fresh(6, 2)}, b^L-b^(L-1)={6**2 - 6}, phi(b^L)={phi(36)}"268        " (phi identity needs prime b)"269    )270271272def main():273    print(f"domain: decimal prec {PREC}, Euler-Maclaurin {EULER_TERMS} terms through B_20")274    print(f"domain: Dirichlet {PAIRS} pairs, lattice norm <= {CUTOFF}"275          f" (radius {lattice_radius(CUTOFF)}), census [1,{SIDE}]^2,"276          f" sieve {TRIALS} pairs in [-{RADIUS},{RADIUS}]^2")277    zk3, zki, inv_zk3 = precision_block()278    sieve_block(zk3, zki, inv_zk3)279280281main()