shell.py

1from fractions import Fraction2from math import isqrt3import mpmath4import numpy as np5import sympy as sp67N2 = 2000008N_ROT = 60009DPS = 3010EXPONENTS = (2, 3, 4)1112def parity_classes(limit):13    bound = isqrt(limit)14    keys = []15    for a in range(-bound, bound + 1):16        square = a * a17        inner = isqrt(limit - square)18        b = np.arange(-inner, inner + 1, dtype=np.int64)19        norms = square + b * b20        slot = (a & 1) * 2 + (b & 1)21        keys.append(slot * (limit + 1) + norms)22    flat = np.bincount(np.concatenate(keys), minlength=4 * (limit + 1))23    return flat.reshape(4, limit + 1).astype(np.int64)2425def jacobi_r2(limit):26    out = np.zeros(limit + 1, dtype=np.int64)27    for d in range(1, limit + 1, 4):28        out[d::d] += 429    for d in range(3, limit + 1, 4):30        out[d::d] -= 431    return out3233def shifted(source, step, limit):34    out = np.zeros(limit + 1, dtype=np.int64)35    out[::step] = source[: limit // step + 1]36    return out3738def doubling_identity(cls, total, limit):39    predicted = shifted(total, 4, limit)40    return int(np.count_nonzero(cls[0][1:] != predicted[1:]))4142def rotation_identity(cls, limit):43    mixed = cls[1] + cls[2]44    predicted = shifted(mixed, 2, limit)45    return int(np.count_nonzero(cls[3][1:] != predicted[1:]))4647def swap_identity(cls):48    return int(np.count_nonzero(cls[1][1:] != cls[2][1:]))4950def odd_odd_points(norm):51    out = []52    bound = isqrt(norm)53    for i in range(-bound, bound + 1):54        if i % 2 == 0:55            continue56        rest = norm - i * i57        if rest < 0:58            continue59        j = isqrt(rest)60        if j * j != rest or j % 2 == 0:61            continue62        out.append((i, j))63        out.append((i, -j))64    return out6566def rotation_as_a_map(cls, limit):67    moved = 068    faults = 069    for norm in range(2, limit + 1, 2):70        points = odd_odd_points(norm)71        images = set()72        for i, j in points:73            u = (i + j) // 274            v = (i - j) // 275            if u * u + v * v != norm // 2 or (u & 1) == (v & 1):76                faults += 177            images.add((u, v))78        moved += len(points)79        if len(images) != len(points) or len(points) != int(cls[3][norm]):80            faults += 181        half = norm // 282        if len(images) != int(cls[1][half]) + int(cls[2][half]):83            faults += 184    return moved, faults8586def pinned_polynomial():87    t, whole = sp.symbols("t S")88    ee, eo, oe, oo = sp.symbols("S_ee S_eo S_oe S_oo")89    solution = sp.solve(90        [91            ee - t**2 * whole,92            oo - t * (eo + oe),93            eo - oe,94            ee + eo + oe + oo - whole,95        ],96        [ee, eo, oe, oo],97        dict=True,98    )[0]99    a_ee, a_eo, a_oe, a_oo = sp.symbols("a_ee a_eo a_oe a_oo")100    filled = (101        a_ee * solution[ee]102        + a_eo * solution[eo]103        + a_oe * solution[oe]104        + a_oo * solution[oo]105    )106    claimed = whole * (107        a_ee * t**2 + a_oo * t * (1 - t) + (a_eo + a_oe) * (1 - t) / 2108    )109    return sp.simplify(filled - claimed) == 0110111def dirichlet_beta(s):112    return (mpmath.zeta(s, mpmath.mpf(1) / 4) - mpmath.zeta(s, mpmath.mpf(3) / 4)) / 4**s113114def shell_sums(cls, total, limit):115    live = np.nonzero(total[: limit + 1])[0]116    live = live[live > 0]117    totals = {s: [mpmath.mpf(0)] * 4 for s in EXPONENTS}118    columns = [cls[slot] for slot in range(4)]119    for norm in live.tolist():120        base = mpmath.mpf(norm)121        powers = {2: base * base}122        powers[3] = powers[2] * base123        powers[4] = powers[3] * base124        for slot in range(4):125            weight = int(columns[slot][norm])126            if weight == 0:127                continue128            value = mpmath.mpf(weight)129            for s in EXPONENTS:130                totals[s][slot] += value / powers[s]131    return totals132133def indicators(code):134    return [(code >> k) & 1 for k in range(4)]135136def q_polynomial(a, s):137    t = Fraction(1, 2**s)138    return a[0] * t * t + a[3] * t * (1 - t) + Fraction(a[1] + a[2], 2) * (1 - t)139140def to_mpf(value):141    return mpmath.mpf(value.numerator) / mpmath.mpf(value.denominator)142143def closed_form_gaps(totals, limit):144    zeta = {s: mpmath.zeta(s) for s in EXPONENTS}145    beta = {s: dirichlet_beta(s) for s in EXPONENTS}146    worst = {}147    for s in EXPONENTS:148        biggest = mpmath.mpf(0)149        for code in range(1, 16):150            a = indicators(code)151            predicted = 4 * zeta[s] * beta[s] * to_mpf(q_polynomial(a, s))152            measured = mpmath.fsum(153                totals[s][slot] for slot in range(4) if a[slot]154            )155            tail = mpmath.mpf(0)156            if s == 2:157                tail = sum(a) * mpmath.pi / (4 * limit)158            gap = abs(measured + tail - predicted)159            biggest = max(biggest, gap)160        worst[s] = biggest161    return worst162163def sci(value):164    return mpmath.nstr(value, 2, strip_zeros=False)165166def main():167    mpmath.mp.dps = DPS168    cls = parity_classes(N2)169    total = cls.sum(axis=0)170    print(f"DOMAIN  parity classes of Z^2 to n = {N2}, rotation as a map to n = {N_ROT}")171    print()172    print("JACOBI")173    jacobi = jacobi_r2(N2)174    bad = int(np.count_nonzero(total[1:] != jacobi[1:]))175    print(f"  r2(n) = 4 (d1(n) - d3(n))            mismatches {bad} of {N2}")176    print()177    print("THE THREE SHELL RELATIONS AS INTEGER IDENTITIES ON THE COEFFICIENTS")178    bad_ee = doubling_identity(cls, total, N2)179    bad_oo = rotation_identity(cls, N2)180    bad_eo = swap_identity(cls)181    print(f"  S_ee = t^2 S    r2_ee(n) = r2(n/4)    mismatches {bad_ee} of {N2}")182    print(f"  S_oo = t S_mix  r2_oo(n) = r2_mix(n/2) mismatches {bad_oo} of {N2}")183    print(f"  S_eo = S_oe     r2_eo(n) = r2_oe(n)   mismatches {bad_eo} of {N2}")184    print()185    print("THE 45 DEGREE ROTATION AS A MAP")186    moved, faults = rotation_as_a_map(cls, N_ROT)187    print(f"  (i,j) -> ((i+j)/2, (i-j)/2) on every odd-odd point of even norm to n = {N_ROT}")188    print(f"  odd-odd points transported {moved}, faults {faults}")189    print()190    print("THE THREE RELATIONS PIN Q_c")191    print(f"  solving the four linear relations returns Q_c exactly: {pinned_polynomial()}")192    print()193    print("SELF TEST OF THE BETA VALUES")194    print(f"  beta(2) against Catalan   {sci(abs(dirichlet_beta(2) - mpmath.catalan))}")195    print(f"  beta(3) against pi^3/32   {sci(abs(dirichlet_beta(3) - mpmath.pi**3 / 32))}")196    print()197    print("THE CLOSED FORM AGAINST TRUNCATED LATTICE SUMS")198    totals = shell_sums(cls, total, N2)199    worst = closed_form_gaps(totals, N2)200    for s in EXPONENTS:201        note = " after the leading tail k pi / (4N)" if s == 2 else ""202        print(f"  s = {s}  worst gap over the 15 nonempty designs {sci(worst[s])}{note}")203204main()