ford_horocycle.py

1import argparse2import math3import sys4import time5from fractions import Fraction67import numpy as np89DELTA = {2: 0.5312805062772051416, 3: 0.705660908028}1011# TOTIENTS1213def totients(n):14    phi = np.arange(n + 1, dtype=np.int64)15    for p in range(2, n + 1):16        if phi[p] == p:17            phi[p::p] -= phi[p::p] // p18    return phi1920# QUOTIENTS2122def quotients(a, b):23    out = []24    while b:25        q, r = divmod(a, b)26        out.append(q)27        a, b = b, r28    return out[1:]2930def member(a, b, m):31    qs = quotients(a, b)32    if not qs:33        return True34    alt = qs[:-1] + [qs[-1] - 1, 1]35    return max(qs) <= m or max(alt) <= m3637def rule(a, b, m):38    qs = quotients(a, b)39    if not qs:40        return True41    return max(qs[:-1], default=1) <= m and qs[-1] <= m + 14243# CROSSINGS4445def egcd(a, b):46    if b == 0:47        return a, 1, 048    g, x, y = egcd(b, a % b)49    return g, y, x - (a // b) * y5051def mobius_image(a, ap, b, d, z):52    x, y = z53    num = (Fraction(a) * x + ap, Fraction(a) * y)54    den = (Fraction(b) * x + d, Fraction(b) * y)55    n2 = den[0] ** 2 + den[1] ** 256    re = (num[0] * den[0] + num[1] * den[1]) / n257    im = (num[1] * den[0] - num[0] * den[1]) / n258    return re, im5960def crossings(bound, heights):61    phi = totients(bound)62    checked = 063    for b in range(1, bound + 1):64        for a in range(b):65            if math.gcd(a, b) != 1:66                continue67            g, s, t = egcd(a, b)68            ap, d = -t, s69            assert a * d - ap * b == 170            centre = (Fraction(a, b), Fraction(1, 2 * b * b))71            r2 = Fraction(1, 4 * b ** 4)72            for x in (Fraction(0), Fraction(1, 3), Fraction(-7, 2), Fraction(5), Fraction(-d, b)):73                re, im = mobius_image(a, ap, b, d, (x, Fraction(1)))74                assert (re - centre[0]) ** 2 + (im - centre[1]) ** 2 == r275                assert im <= Fraction(1, b * b)76                if x == Fraction(-d, b):77                    assert im == Fraction(1, b * b) and re == Fraction(a, b)78                checked += 179    print(f"horocycle images checked exactly: {checked} points on {int(phi[1:bound + 1].sum())} Ford circles, b <= {bound}, every image on the Ford circle over a/b, top point 1/b^2 at x = -d/b")80    print("Q  height  closed-disk count  sum_{b<=Q} phi  open-disk count  sum_{b<Q} phi  count at 1/(2Q^2)  sum_{b<=floor(Q sqrt2)} phi")81    for Q in heights:82        h = Fraction(1, Q * Q)83        closed = sum(int(phi[b]) for b in range(1, bound + 1) if Fraction(1, b * b) >= h)84        opened = sum(int(phi[b]) for b in range(1, bound + 1) if Fraction(1, b * b) > h)85        half = sum(int(phi[b]) for b in range(1, bound + 1) if Fraction(1, b * b) >= h / 2)86        mQ = int(phi[1:Q + 1].sum())87        mQ1 = int(phi[1:Q].sum())88        m2 = int(phi[1:math.isqrt(2 * Q * Q) + 1].sum())89        assert closed == mQ and opened == mQ1 and half == m290        print(f"{Q}  1/{Q * Q}  {closed}  {mQ}  {opened}  {mQ1}  {half}  {m2}")9192# BRIDGE9394def bump_smooth(u):95    u = np.asarray(u, dtype=float)96    out = np.zeros_like(u)97    inside = (u > 1) & (u < 2)98    w = u[inside]99    out[inside] = np.exp(4 - 1 / ((w - 1) * (2 - w)))100    return out101102def bump_c2(u):103    u = np.asarray(u, dtype=float)104    out = np.zeros_like(u)105    inside = (u > 1) & (u < 2)106    w = u[inside]107    out[inside] = 64 * (w - 1) ** 3 * (2 - w) ** 3108    return out109110def gl(nodes):111    x, w = np.polynomial.legendre.leggauss(nodes)112    return x, w113114def integrate(f, lo, hi, x, w):115    mid, half = (lo + hi) / 2, (hi - lo) / 2116    return half * np.dot(w, f(mid + half * x))117118def g_transform(f, t, x, w):119    if t >= 1:120        return 0.0121    hi = math.sqrt(1 / (t * t) - 1)122    lo = math.sqrt(1 / (2 * t * t) - 1) if t * t <= 0.5 else 0.0123    return 2 * integrate(lambda v: f(1 / (t * t * (1 + v * v))), lo, hi, x, w)124125def horocycle_integral(f, h, x, w):126    Q = math.isqrt(int(1 / h))127    total = 0.0128    bumps = 0129    for c in range(1, Q + 1):130        s = h - c * c * h * h131        if s < 0:132            continue133        width = math.sqrt(s) / c134        inner = math.sqrt(max(h / 2 - c * c * h * h, 0.0)) / c135        dlo = math.ceil(-c * (1 + width))136        dhi = math.floor(c * width)137        for d in range(dlo, dhi + 1):138            if math.gcd(c, abs(d)) != 1:139                continue140            centre = -d / c141            pieces = ((centre - width, centre - inner), (centre + inner, centre + width))142            hit = False143            for lo, hi in pieces:144                lo, hi = max(0.0, lo), min(1.0, hi)145                if lo >= hi:146                    continue147                total += integrate(lambda t: f(h / ((c * t + d) ** 2 + c * c * h * h)), lo, hi, x, w)148                hit = True149            bumps += hit150    return total, bumps151152def bridge(kmax, nodes):153    x, w = gl(nodes)154    for name, f in (("smooth bump", bump_smooth), ("C2 bump", bump_c2)):155        main = 3 / math.pi * integrate(lambda u: f(u) / (u * u), 1, 2, x, w)156        print(f"{name}: main term (3/pi) int f(w) w^-2 dw = {main:.15f}")157        print("k  h=4^-k  Q  bumps  horocycle integral  h S_g(h^1/2)  abs gap  error E  E/h^(3/4)")158        Qmax = 2 ** kmax159        phi = totients(Qmax)160        for k in range(2, kmax + 1):161            h = 4.0 ** (-k)162            Q = 2 ** k163            t0 = time.time()164            lhs, bumps = horocycle_integral(f, h, x, w)165            rhs = h * sum(int(phi[c]) * g_transform(f, c * math.sqrt(h), x, w) for c in range(1, Q + 1))166            err = lhs - main167            print(f"{k}  {h:.3e}  {Q}  {bumps}  {lhs:.15f}  {rhs:.15f}  {abs(lhs - rhs):.1e}  {err:+.3e}  {err / h ** 0.75:+.4f}  ({time.time() - t0:.1f}s)")168169# CENSUS170171def census_walk(m, qmax):172    hist = np.zeros(qmax + 1, dtype=np.int64)173    hist[1] = 1174    stack = [(1, 2, 0, 1, 1, 1, 0, 1)]175    while stack:176        p, q, p1, q1, p2, q2, last, run = stack.pop()177        hist[q] += 1178        ql = q + q1179        if ql <= qmax:180            r = run + 1 if last == 0 else 1181            if r <= m:182                stack.append((p + p1, ql, p1, q1, p, q, 0, r))183        qr = q + q2184        if qr <= qmax:185            r = run + 1 if last == 1 else 1186            if r <= m:187                stack.append((p + p2, qr, p, q, p2, q2, 1, r))188    return np.cumsum(hist)189190def census_nodes(m, qmax):191    out = []192    stack = [(1, 2, 0, 1, 1, 1, 0, 1)]193    while stack:194        node = stack.pop()195        p, q, p1, q1, p2, q2, last, run = node196        out.append(node)197        ql, qr = q + q1, q + q2198        if ql <= qmax:199            r = run + 1 if last == 0 else 1200            if r <= m:201                stack.append((p + p1, ql, p1, q1, p, q, 0, r))202        if qr <= qmax:203            r = run + 1 if last == 1 else 1204            if r <= m:205                stack.append((p + p2, qr, p, q, p2, q2, 1, r))206    return out207208def census(jmax2, jmax3, jcontrol, jcheck):209    qc = 2 ** jcontrol210    t0 = time.time()211    cum = census_walk(10 ** 9, qc)212    phi = totients(qc)213    tot = np.cumsum(phi)214    assert all(int(cum[2 ** j]) == int(tot[2 ** j]) for j in range(0, jcontrol + 1))215    print(f"control A = N to Q = 2^{jcontrol}: the walk's count is sum_(b <= Q) phi(b) at every Q = 2^j, {int(cum[qc])} at the top ({time.time() - t0:.1f}s)")216    qk = 2 ** jcheck217    for m in (2, 3):218        walked = {(p, q) for p, q, *_ in census_nodes(m, qk)} | {(0, 1)}219        brute = {(a, b) for b in range(1, qk + 1) for a in range(b) if math.gcd(a, b) == 1 and member(a, b, m)}220        ruled = {(a, b) for b in range(1, qk + 1) for a in range(b) if math.gcd(a, b) == 1 and rule(a, b, m)}221        assert walked == brute == ruled222        both = one = 0223        for p, q, p1, q1, p2, q2, last, run in census_nodes(m, qk):224            qs = quotients(p, q)225            an = qs[-1]226            assert run == an - 1227            older, younger = ((p1, q1), (p2, q2)) if last == 0 else ((p2, q2), (p1, q1))228            assert older[1] < younger[1] or (p, q) == (1, 2)229            with_younger = member(p + younger[0], q + younger[1], m)230            with_older = member(p + older[0], q + older[1], m)231            assert with_younger232            assert with_older == (an <= m)233            if with_older:234                both += 1235            else:236                one += 1237        print(f"m = {m}, b <= 2^{jcheck}: walk, either-expansion test and the rule a_1..a_(n-1) <= m, a_n <= m+1 agree on {len(walked)} fractions; mediant with the younger neighbour stays in {both + one} of {both + one}, with the older in {both}, exactly the nodes with a_n <= m, and refused in {one}, exactly a_n = m + 1")238    for m, jmax in ((2, jmax2), (3, jmax3)):239        qmax = 2 ** jmax240        t0 = time.time()241        cum = census_walk(m, qmax)242        elapsed = time.time() - t0243        d = 2 * DELTA[m]244        print(f"m = {m}, A = {{1..{m}}}, 2 delta = {d:.13f}, walk to Q = 2^{jmax}: {int(cum[qmax])} fractions in [0, 1) ({elapsed:.1f}s)")245        print("j  Q  N(Q)  Q^(2 delta)  N/Q^(2 delta)  octave exponent log2(N(2Q)/N(Q))  octave exponent minus 2 delta  two-octave exponent log4(N(4Q)/N(Q))  min and max of N/Q^(2 delta) on the quarter octaves of [j, j+1)")246        for j in range(0, jmax + 1):247            Q = 2 ** j248            n = int(cum[Q])249            power = Q ** d250            one = f"{math.log2(int(cum[2 * Q]) / n):.4f}" if j + 1 <= jmax else "-"251            dev = f"{math.log2(int(cum[2 * Q]) / n) - d:+.4f}" if j + 1 <= jmax else "-"252            two = f"{math.log2(int(cum[4 * Q]) / n) / 2:.4f}" if j + 2 <= jmax else "-"253            band = [int(cum[int(round(2 ** (j + k / 4)))]) / 2 ** ((j + k / 4) * d) for k in range(4) if j + k / 4 <= jmax]254            print(f"{j}  {Q}  {n}  {power:.3f}  {n / power:.5f}  {one}  {dev}  {two}  {min(band):.5f} {max(band):.5f}")255256# MAIN257258def main():259    ap = argparse.ArgumentParser()260    ap.add_argument("verb", choices=["crossings", "bridge", "census", "all"])261    ap.add_argument("--bound", type=int, default=60)262    ap.add_argument("--kmax", type=int, default=10)263    ap.add_argument("--nodes", type=int, default=200)264    ap.add_argument("--jmax2", type=int, default=24)265    ap.add_argument("--jmax3", type=int, default=18)266    ap.add_argument("--jcontrol", type=int, default=10)267    ap.add_argument("--jcheck", type=int, default=9)268    args = ap.parse_args()269    t0 = time.time()270    if args.verb in ("crossings", "all"):271        crossings(args.bound, [1, 2, 3, 10, 32, 42])272    if args.verb in ("bridge", "all"):273        bridge(args.kmax, args.nodes)274    if args.verb in ("census", "all"):275        census(args.jmax2, args.jmax3, args.jcontrol, args.jcheck)276    print(f"total {time.time() - t0:.1f}s", file=sys.stderr)277278if __name__ == "__main__":279    main()