totient_constant.py

1from decimal import Decimal, getcontext2from fractions import Fraction3from math import comb, gcd, isqrt45getcontext().prec = 5067BOUNDS = [50, 200, 800, 1200, 2000, 3200, 12800, 51200, 102400]8SET_BOUNDS = [50, 200, 800]9PUBLISHED = {10    "Q(i)": {50: 672, 200: 10608, 800: 168088},11    "Q(sqrt -3)": {50: 630, 200: 9606, 800: 151020, 1200: 337026, 2000: 945486, 3200: 2419950},12}13FIELDS = [14    {"name": "Q(i)", "ram": 2, "mod": 4, "w": 4, "disc": 4, "D": -4},15    {"name": "Q(sqrt -3)", "ram": 3, "mod": 3, "w": 6, "disc": 3, "D": -3},16]17DISCS = [-4, -3, -20, -23]18EULER_M = 6019EULER_J = 102021def bernoulli(m):22    b = [Fraction(0)] * (m + 1)23    b[0] = Fraction(1)24    for n in range(1, m + 1):25        s = sum(comb(n + 1, k) * b[k] for k in range(n))26        b[n] = -s / (n + 1)27    return b2829def dec(f):30    return Decimal(f.numerator) / Decimal(f.denominator)3132def hurwitz2(a, bern):33    base = dec(a)34    s = sum(1 / ((base + k) * (base + k)) for k in range(EULER_M))35    x = base + EULER_M36    s += 1 / x + 1 / (2 * x * x)37    for j in range(1, EULER_J + 1):38        s += dec(bern[2 * j]) / x ** (2 * j + 1)39    return s4041def arctan_inv(n):42    total = Decimal(0)43    term = Decimal(1) / Decimal(n)44    sq = Decimal(n) * Decimal(n)45    k = 046    while True:47        t = term / (2 * k + 1)48        if t < Decimal(10) ** -46:49            return total50        total += t if k % 2 == 0 else -t51        term /= sq52        k += 15354def dpi():55    return 16 * arctan_inv(5) - 4 * arctan_inv(239)5657def jacobi(a, n):58    a %= n59    r = 160    while a:61        while a % 2 == 0:62            a //= 263            if n % 8 in (3, 5):64                r = -r65        a, n = n, a66        if a % 4 == 3 and n % 4 == 3:67            r = -r68        a %= n69    return r if n == 1 else 07071def kron(d, n):72    if n == 0:73        return 074    r = 175    while n % 2 == 0:76        if d % 2 == 0:77            return 078        n //= 279        r *= 1 if d % 8 in (1, 7) else -180    return r * jacobi(d % n, n) if n > 1 else r8182def units(d):83    return 4 if d == -4 else (6 if d == -3 else 2)8485def class_number(d):86    q = -d87    return units(d) * -sum(kron(d, a) * a for a in range(1, q)) // (2 * q)8889def lvalue(d, bern):90    q = -d91    s = sum(kron(d, a) * hurwitz2(Fraction(a, q), bern) for a in range(1, q))92    return s / (q * q)9394def spf_sieve(n):95    spf = list(range(n + 1))96    for p in range(2, isqrt(n) + 1):97        if spf[p] == p:98            for m in range(p * p, n + 1, p):99                if spf[m] == m:100                    spf[m] = p101    return spf102103def factor(n, spf):104    out = []105    while n > 1:106        p = spf[n]107        e = 0108        while n % p == 0:109            n //= p110            e += 1111        out.append((p, e))112    return out113114def kind(fld, p):115    if p == fld["ram"]:116        return "r"117    return "s" if p % fld["mod"] == 1 else "i"118119def mulc(fld, u, v):120    x, y = u121    c, d = v122    if fld["ram"] == 2:123        return (x * c - y * d, x * d + y * c)124    return (x * c - y * d, x * d + y * c - y * d)125126def conjc(fld, z):127    a, b = z128    return (a, -b) if fld["ram"] == 2 else (a - b, -b)129130def nrm(fld, z):131    a, b = z132    return a * a + b * b if fld["ram"] == 2 else a * a - a * b + b * b133134def classes(fld, m):135    if fld["ram"] == 2:136        return [(a, b) for a in range(1, isqrt(m) + 1) for b in range(isqrt(m - a * a) + 1)]137    r = isqrt(4 * m // 3) + 1138    out = []139    for a in range(-r, r + 1):140        for b in range(-r, r + 1):141            n = a * a - a * b + b * b142            if 0 < n <= m:143                z = (a, b)144                best = z145                for _ in range(3):146                    z = (-z[1], z[0] - z[1])147                    best = min(best, z, (-z[0], -z[1]))148                if best == (a, b):149                    out.append((a, b))150    return out151152def totient(fld, z, n, spf):153    v = n154    a, b = z155    for p, e in factor(n, spf):156        k = kind(fld, p)157        if k == "r":158            v = v // p * (p - 1)159        elif k == "i":160            v = v // (p * p) * (p * p - 1)161        else:162            v = v // p * (p - 1)163            if a % p == 0 and b % p == 0:164                v = v // p * (p - 1)165    return v166167def mobius(fld, z, n, spf):168    v = 1169    a, b = z170    for p, e in factor(n, spf):171        k = kind(fld, p)172        if k == "r":173            if e != 1:174                return 0175            v = -v176        elif k == "i":177            if e != 2:178                return 0179            v = -v180        else:181            if e == 1:182                v = -v183            elif e == 2 and a % p == 0 and b % p == 0:184                pass185            else:186                return 0187    return v188189def tables(fld, m, spf):190    cls = classes(fld, m)191    phi = [0] * (m + 1)192    nsum = [0] * (m + 1)193    mob = [0] * (m + 1)194    for z in cls:195        n = nrm(fld, z)196        phi[n] += totient(fld, z, n, spf)197        nsum[n] += n198        mob[n] += mobius(fld, z, n, spf)199    for n in range(1, m + 1):200        phi[n] += phi[n - 1]201        nsum[n] += nsum[n - 1]202    return cls, phi, nsum, mob203204def norm_arrays(d, m):205    chi = [kron(d, n) for n in range(m + 1)]206    ideals = [0] * (m + 1)207    for k in range(1, m + 1):208        if chi[k]:209            for n in range(k, m + 1, k):210                ideals[n] += chi[k]211    prime = [True] * (m + 1)212    prime[0] = prime[1] = False213    for k in range(2, isqrt(m) + 1):214        if prime[k]:215            for n in range(k * k, m + 1, k):216                prime[n] = False217    mu = [1] * (m + 1)218    sq = [True] * (m + 1)219    for k in range(2, m + 1):220        if prime[k]:221            for n in range(k, m + 1, k):222                mu[n] = -mu[n]223            for n in range(k * k, m + 1, k * k):224                sq[n] = False225    mk = [0] * (m + 1)226    for k in range(1, m + 1):227        if sq[k] and mu[k]:228            for j in range(1, m // k + 1):229                if sq[j] and chi[j]:230                    mk[k * j] += mu[k] * mu[j] * chi[j]231    tt = [0] * (m + 1)232    for n in range(1, m + 1):233        tt[n] = tt[n - 1] + n * ideals[n]234    return mk, tt235236def norm_totient_sum(n, mk, tt):237    return sum(mk[k] * tt[n // k] for k in range(1, n + 1) if mk[k])238239def totient_sum(phi, m):240    return phi[m]241242def convolution_sum(m, nsum, mob):243    return sum(mob[n] * nsum[m // n] for n in range(1, m + 1) if mob[n])244245def hnf(r1, r2):246    a1, b1 = r1247    a2, b2 = r2248    x, y, d = 1, 0, a1249    u, v, e = 0, 1, a2250    while e:251        q = d // e252        d, e = e, d - q * e253        x, u = u, x - q * u254        y, v = v, y - q * v255    if d < 0:256        d, x, y = -d, -x, -y257    return d, x * b1 + y * b2, abs((a2 // d) * b1 - (a1 // d) * b2)258259def farey_set(fld, m):260    seen = set()261    for z in classes(fld, m):262        n = nrm(fld, z)263        r1 = mulc(fld, (1, 0), z)264        r2 = mulc(fld, (0, 1), z)265        h11, h12, h22 = hnf(r1, r2)266        cj = conjc(fld, z)267        for x in range(h11):268            for y in range(h22):269                u, v = mulc(fld, (x, y), cj)270                u %= n271                v %= n272                g = gcd(gcd(u, v), n)273                seen.add((u // g, v // g, n // g))274    return len(seen)275276def d12(x):277    return str(+x.quantize(Decimal("1.000000000000")))278279def main():280    bern = bernoulli(2 * EULER_J)281    pi = dpi()282    z2 = pi * pi / 6283    lv = {d: lvalue(d, bern) for d in DISCS}284    print("pi", d12(pi))285    print("zeta(2)", d12(z2))286    print("L(2, chi_-4) Catalan", d12(lv[-4]))287    print("L(2, chi_-3)", d12(lv[-3]))288    const = {}289    for fld in FIELDS:290        dk = fld["disc"]291        zk = z2 * lv[fld["D"]]292        rho = 2 * pi / (fld["w"] * Decimal(dk).sqrt())293        c = rho / (2 * zk)294        const[fld["name"]] = c295        print("field", fld["name"], "w", fld["w"], "|D|", dk)296        print("  zeta_K(2)", d12(zk))297        print("  rho_K", d12(rho))298        print("  c = rho_K/(2 zeta_K(2))", d12(c))299        print("  Sayous c_K = pi/(sqrt|D| zeta_K(2))", d12(pi / (Decimal(dk).sqrt() * zk)), "= w c", d12(fld["w"] * c))300    spf = spf_sieve(max(BOUNDS))301    byclass = {}302    for fld in FIELDS:303        name = fld["name"]304        c = const[name]305        cls, phi, nsum, mob = tables(fld, max(BOUNDS), spf)306        byclass[fld["D"]] = phi307        print("field", name, "classes to", max(BOUNDS), len(cls))308        for m in BOUNDS:309            s = totient_sum(phi, m)310            dn = Decimal(m)311            main_term = c * dn * dn312            ratio = Decimal(s) / main_term313            dev = Decimal(s) - main_term314            mark = PUBLISHED[name].get(m)315            tag = "published " + str(mark) + " match " + str(mark == s) if mark else "new"316            print("  N", m, "count", s, "count/(c N^2)", d12(ratio), "dev/N^1.5", d12(dev / (dn * dn.sqrt())), "dev/(N ln N)", d12(dev / (dn * dn.ln())), tag)317        bad = [m for m in BOUNDS if convolution_sum(m, nsum, mob) != totient_sum(phi, m)]318        print("  convolution identity mismatches", len(bad))319        for m in SET_BOUNDS:320            k = farey_set(fld, m)321            print("  Farey set size at N", m, k, "equals totient sum", k == totient_sum(phi, m))322    mx = max(BOUNDS)323    for d in DISCS:324        h = class_number(d)325        w = units(d)326        zk = z2 * lv[d]327        rho = 2 * pi * h / (w * Decimal(-d).sqrt())328        c = rho / (2 * zk)329        mk, tt = norm_arrays(d, mx)330        print("discriminant", d, "h", h, "w", w, "L(2, chi_D)", d12(lv[d]), "rho_K", d12(rho), "c", d12(c))331        for m in BOUNDS:332            s = norm_totient_sum(m, mk, tt)333            dn = Decimal(m)334            main_term = c * dn * dn335            tail = " class route " + str(s == totient_sum(byclass[d], m)) if d in byclass else ""336            print("  N", m, "count", s, "count/(c N^2)", d12(s / main_term), "dev/N^1.5", d12((s - main_term) / (dn * dn.sqrt())) + tail)337338main()