eisenstein_stack.py

1from fractions import Fraction2from math import isqrt, lcm, log3from pathlib import Path45from sympy import N as N_6from sympy import Poly, Rational, Symbol, cos, cyclotomic_poly, minimal_polynomial, sin, sqrt7from sympy import pi as spi8from sympy import zeta as szeta910HERE = Path(__file__).resolve().parent1112X = Symbol("x")13CYC = 36014PHI360 = Poly(cyclotomic_poly(CYC, X), X, domain="QQ")15SPOT_DEGREES = [0, 30, 45, 60, 90, 120, 137, 180, 240, 270, 300, 359]16STACK_N = 5017COUNT_NS = [50, 200, 800, 1200, 2000, 3200]18CIRCLE_T = 40019CSL_M = 10020ROT_BOUND = 6021ROT_REACH = 2022FIG_N = 6023FIG_W = 90024FIG_H = 60025FIG_OUT = HERE / "eisenstein-stack.png"26ROOT3 = 3.0 ** 0.52728SURF = (255, 255, 255)29BLUE = (0, 140, 255)30GRAY = (186, 186, 191)31INK = (0, 0, 0)3233def enorm(z):34    a, b = z35    return a * a - a * b + b * b3637def emul(z, w):38    a, b = z39    c, d = w40    return (a * c - b * d, a * d + b * c - b * d)4142def econj(z):43    a, b = z44    return (a - b, -b)4546def edivides(w, z):47    p = emul(z, econj(w))48    n = enorm(w)49    return p[0] % n == 0 and p[1] % n == 05051def equo(z, w):52    p = emul(z, econj(w))53    n = enorm(w)54    return (p[0] // n, p[1] // n)5556def enearest(z, w):57    p = emul(z, econj(w))58    n = enorm(w)59    return ((2 * p[0] + n) // (2 * n), (2 * p[1] + n) // (2 * n))6061def egcd(z, w):62    while w != (0, 0):63        q = enearest(z, w)64        t = emul(q, w)65        z, w = w, (z[0] - t[0], z[1] - t[1])66    return z6768def associates(z):69    out = [z]70    for _ in range(2):71        a, b = out[-1]72        out.append((-b, a - b))73    return out + [(-a, -b) for a, b in out]7475def ecanon(z):76    if z == (0, 0):77        return z78    for c in associates(z):79        if c[0] > 0 and 0 <= c[1] < c[0]:80            return c81    return z8283def classes_up_to(n):84    out = set()85    r = isqrt(4 * n // 3) + 286    for a in range(-r, r + 1):87        for b in range(-r, r + 1):88            if (a, b) == (0, 0):89                continue90            if enorm((a, b)) <= n:91                out.add(ecanon((a, b)))92    return sorted(out, key=lambda z: (enorm(z), z))9394def hex_classes_closed(t):95    s = 096    j = 097    while 3 * j + 1 <= t:98        s += t // (3 * j + 1) - t // (3 * j + 2)99        j += 1100    return s101102def hex_classes_direct(t):103    return len(classes_up_to(t))104105def eisenstein_primes_dividing(z):106    n = enorm(z)107    rational = []108    m = n109    q = 2110    while q * q <= m:111        if m % q == 0:112            while m % q == 0:113                m //= q114            rational.append(q)115        q += 1116    if m > 1:117        rational.append(m)118    primes = []119    for q in rational:120        if q == 3:121            primes.append(ecanon((1, -1)))122        elif q % 3 == 2:123            primes.append((q, 0))124        else:125            r = 2126            root = None127            while root is None:128                c = pow(r, (q - 1) // 3, q)129                if (c * c + c + 1) % q == 0:130                    root = c131                r += 1132            p = egcd((q, 0), (root, -1))133            primes.append(ecanon(p))134            primes.append(ecanon(econj(p)))135    return [p for p in primes if edivides(p, z)]136137def eisenstein_totient(z):138    v = enorm(z)139    if v == 1:140        return 1141    for p in eisenstein_primes_dividing(z):142        v = v // enorm(p) * (enorm(p) - 1)143    return v144145def totient_sum(n_max):146    return sum(eisenstein_totient(z) for z in classes_up_to(n_max))147148def layer_nodes(z):149    n = enorm(z)150    c = econj(z)151    out = set()152    for m in range(n):153        for k in range(n):154            p, q = emul((m, k), c)155            out.add((Fraction(p % n, n), Fraction(q % n, n)))156    return out157158def literal_stack(n_max):159    hits = {}160    for z in classes_up_to(n_max):161        for p in layer_nodes(z):162            hits[p] = hits.get(p, 0) + 1163    return hits164165def reduced_denominator(node):166    x, y = node167    d = lcm(x.denominator, y.denominator)168    u = (x.numerator * (d // x.denominator), y.numerator * (d // y.denominator))169    g = egcd(u, (d, 0))170    return ecanon(equo((d, 0), g))171172def closed_brightness(n_max, d):173    return hex_classes_closed(n_max // enorm(d))174175def cyc_reduce(terms):176    acc = {}177    for e, c in terms:178        k = (e % CYC,)179        acc[k] = acc.get(k, 0) + c180    acc = {k: v for k, v in acc.items() if v}181    if not acc:182        return Poly(0, X, domain="QQ")183    return Poly(acc, X, domain="QQ").rem(PHI360)184185def field_cos_sin(d):186    c = cyc_reduce([(d, 1), (-d, 1)])187    s = cyc_reduce([(d + 120, 1), (d + 240, -1), (-d + 120, -1), (-d + 240, 1)])188    return c.degree() <= 0, s.degree() <= 0189190def field_rotation_degrees():191    return [d for d in range(CYC) if all(field_cos_sin(d))]192193def spot_check_degrees():194    out = []195    for d in SPOT_DEGREES:196        c = minimal_polynomial(cos(spi * Rational(d, 180)), X, polys=True).degree()197        s = minimal_polynomial(sin(spi * Rational(d, 180)) / sqrt(3), X, polys=True).degree()198        out.append((d, c == 1, s == 1, field_cos_sin(d)))199    return out200201def rotation_hits(bound, reach):202    seen = set()203    for a in range(-reach, reach + 1):204        for b in range(-reach, reach + 1):205            if (a, b) == (0, 0):206                continue207            n = enorm((a, b))208            seen.add((Fraction(2 * a * a - 2 * a * b - b * b, 2 * n), Fraction(2 * a * b - b * b, 2 * n)))209    direct = set()210    for r in range(1, bound + 1):211        for p in range(-r, r + 1):212            rest = r * r - p * p213            if rest < 0 or rest % 3:214                continue215            s = isqrt(rest // 3)216            if 3 * s * s == rest:217                for t in (s, -s):218                    direct.add((Fraction(p, r), Fraction(t, r)))219    return len(direct), direct <= seen220221def chi3(n):222    return (0, 1, -1)[n % 3]223224def dirichlet_mul(u, v, m):225    out = [0] * (m + 1)226    for i in range(1, m + 1):227        if u[i] == 0:228            continue229        for j in range(1, m // i + 1):230            out[i * j] += u[i] * v[j]231    return out232233def csl_zeta_ratio(m):234    a = [0] * (m + 1)235    for d in range(1, m + 1):236        c = chi3(d)237        if c:238            for n in range(d, m + 1, d):239                a[n] += c240    inv_square = [0] * (m + 1)241    for k in range(1, isqrt(m) + 1):242        f = k243        mu = 1244        q = 2245        while q * q <= f:246            if f % q == 0:247                f //= q248                if f % q == 0:249                    mu = 0250                    break251                mu = -mu252            q += 1253        if mu and f > 1:254            mu = -mu255        inv_square[k * k] = mu256    inv_three = [0] * (m + 1)257    e = 1258    j = 0259    while 3 ** j <= m:260        inv_three[3 ** j] = e261        e = -e262        j += 1263    return dirichlet_mul(dirichlet_mul(a, inv_square, m), inv_three, m)264265def csl_euler_product(m):266    out = [0] * (m + 1)267    out[1] = 1268    for q in range(2, m + 1):269        if q % 3 != 1 or any(q % r == 0 for r in range(2, isqrt(q) + 1)):270            continue271        local = [0] * (m + 1)272        local[1] = 1273        k = q274        while k <= m:275            local[k] = 2276            k *= q277        out = dirichlet_mul(out, local, m)278    return out279280def draw(n_max):281    from PIL import Image, ImageDraw282283    stack = literal_stack(n_max)284    top = max(stack.values())285    ss = 3286    w, h = FIG_W * ss, FIG_H * ss287    img = Image.new("RGB", (w, h), SURF)288    pen = ImageDraw.Draw(img)289    unit = 560.0 * ss290    ox = (w - 1.5 * unit) / 2 + unit / 2291    oy = (h + unit * ROOT3 / 2) / 2292    corners = [(0, 0), (1, 0), (1, 1), (0, 1)]293    place = lambda x, y: (ox + unit * (x - y / 2), oy - unit * y * ROOT3 / 2)294    pen.line([place(*c) for c in corners] + [place(0, 0)], fill=GRAY, width=2 * ss)295    for node, b in sorted(stack.items(), key=lambda kv: kv[1]):296        px, py = place(float(node[0]), float(node[1]))297        r = (1.4 + 7.0 * (b / top) ** 0.5) * ss298        pen.ellipse([px - r, py - r, px + r, py + r], fill=BLUE if b < top else INK)299    small = img.resize((FIG_W, FIG_H), Image.LANCZOS)300    small.convert("P", palette=Image.ADAPTIVE, colors=48).save(FIG_OUT, optimize=True)301    return len(stack), top302303def main():304    print("ring Z[omega], omega^2 = -1 - omega, norm a^2 - ab + b^2, units 6")305    print("coordinates in the basis 1, omega; fundamental domain the unit square of that chart")306307    rots = field_rotation_degrees()308    print("field rotation degrees", *rots)309    print("field rotation count", len(rots))310    spots = spot_check_degrees()311    print("spot degrees", len(spots))312    print("spot minpoly agrees with cyclotomic", all((c, s) == r for _, c, s, r in spots))313    hits, covered = rotation_hits(ROT_BOUND, ROT_REACH)314    print("rational rotations denominator <=", ROT_BOUND, hits)315    print("all are w/conj(w), Eisenstein w in box", ROT_REACH, covered)316317    print("hex circle count closed equals direct", all(hex_classes_closed(t) == hex_classes_direct(t) for t in range(CIRCLE_T + 1)))318    print("hex circle counts t=1..12", *[hex_classes_closed(t) for t in range(1, 13)])319320    stack = literal_stack(STACK_N)321    layers = classes_up_to(STACK_N)322    print("stack norm bound", STACK_N)323    print("stack layers", len(layers))324    print("stack nodes", len(stack))325    print("eisenstein totient sum", totient_sum(STACK_N))326    print("nodes equal totient sum", len(stack) == totient_sum(STACK_N))327    bad = 0328    for node, b in stack.items():329        if closed_brightness(STACK_N, reduced_denominator(node)) != b:330            bad += 1331    print("brightness comparisons", len(stack))332    print("brightness mismatches", bad)333    print("origin brightness", stack[(Fraction(0), Fraction(0))], "max", max(stack.values()))334335    lval = (szeta(2, Rational(1, 3)) - szeta(2, Rational(2, 3))) / 9336    z2 = spi ** 2 / 6337    rho = spi / (3 * sqrt(3))338    const = rho / (2 * z2 * lval)339    print("zeta(2)", "%.15f" % float(N_(z2, 30)))340    print("L(2, chi_-3)", "%.15f" % float(N_(lval, 30)))341    print("zeta_K(2) = zeta(2) L(2, chi_-3)", "%.15f" % float(N_(z2 * lval, 30)))342    print("residue of zeta_K at 1 = pi/(3 sqrt 3)", "%.15f" % float(N_(rho, 30)))343    c = float(N_(const, 30))344    print("constant pi/(6 sqrt3 zeta(2) L(2, chi_-3))", "%.15f" % c)345    for n in COUNT_NS:346        s = totient_sum(n)347        r = s / (n * n)348        print("node count N", n, "=", s, "ratio N^2", "%.6f" % r, "deviation", "%+.6f" % (r - c), "scaled by N/log N", "%+.3f" % ((r - c) * n / log(n)))349350    ratio = csl_zeta_ratio(CSL_M)351    euler = csl_euler_product(CSL_M)352    print("csl series bound", CSL_M)353    print("csl zeta ratio equals euler product", ratio == euler)354    print("csl nonzero coefficients", *[(n, ratio[n]) for n in range(1, CSL_M + 1) if ratio[n]])355356    nodes, top = draw(FIG_N)357    print("figure norm bound", FIG_N, "nodes", nodes, "brightest", top)358    print("figure bytes", FIG_OUT.stat().st_size)359360main()