interval.py

1import math2import sys3import time45import numpy as np6from mpmath import iv, mp78GAMMA = 0.57721566490153299ZETA3 = 1.202056903159594210C_P = math.sqrt(2) - 4 / math.pi11K_2 = 2 * (4 / 3 + (36 / 35) * (7 * ZETA3 / 8 - 1))12TAIL = 3.41314# CLOSED FORMS1516def x0_float(N):17    return (N / math.pi) * (math.log(N + 3) + GAMMA + math.log(math.tan(3 * math.pi / 8 + math.pi / (4 * N)))) + C_P * (N + 1) ** 2 / (8 * N)1819def lam_float(N):20    return (N + 1) / 2 + x0_float(N) + 0.5 / math.sin(math.pi / (2 * N))2122def rate_iv(N):23    iv.prec = 12024    N = iv.mpf(N)25    pi = iv.pi26    cp = iv.sqrt(2) - 4 / pi27    x0 = (N / pi) * (iv.log(N + 3) + iv.euler + iv.log(iv.tan(3 * pi / 8 + pi / (4 * N)))) + cp * (N + 1) ** 2 / (8 * N)28    n = (N + 1) / 229    return (n + x0 + 1 / (2 * iv.sin(pi / (2 * N)))) / n3031def c_inf_iv():32    iv.prec = 12033    pi = iv.pi34    return 1 + 2 / pi + (2 / pi) * (iv.euler + iv.log(1 + iv.sqrt(2))) + (iv.sqrt(2) - 4 / pi) / 43536def tail_gap_iv(N, d, c):37    iv.prec = 12038    N = iv.mpf(N)39    return N ** (iv.mpf(1) / d) - (2 / iv.pi) * iv.log(N) - c - TAIL / N4041def exact_gap_iv(N, d):42    iv.prec = 12043    return iv.mpf(N) ** (iv.mpf(1) / d) - rate_iv(N)4445def first_tail(e, c_hi, lo, hi):46    f = lambda N: N ** e - (2 / math.pi) * math.log(N) - c_hi - TAIL / N47    while hi - lo > 1:48        mid = (lo + hi) // 249        if f(mid) > 0:50            hi = mid51        else:52            lo = mid53    return hi5455def down(x, d=5):56    v = mp.mpf(x.a.a)57    k = d - 1 - int(mp.floor(mp.log10(abs(v))))58    return mp.nstr(mp.floor(v * 10 ** k) / mp.mpf(10) ** k, d)5960def up(x, d=5):61    v = mp.mpf(x.b.b)62    k = d - 1 - int(mp.floor(mp.log10(abs(v))))63    return mp.nstr(mp.ceil(v * 10 ** k) / mp.mpf(10) ** k, d)6465def fup(x, d=6):66    return f"{math.ceil(x * 10 ** d) / 10 ** d:.{d}f}"6768def fdown(x, d=6):69    return f"{math.floor(x * 10 ** d) / 10 ** d:.{d}f}"7071def gdown(x, d=5):72    k = d - 1 - math.floor(math.log10(x))73    return f"{math.floor(x * 10 ** k) / 10 ** k:.{d - 1}e}"7475def gup(x, d=3):76    k = d - 1 - math.floor(math.log10(x))77    return f"{math.ceil(x * 10 ** k) / 10 ** k:.{d - 1}e}"7879def odd_up(N):80    return N if N % 2 else N + 18182# WALL8384def wall_at(d, label, lo, hi, mono):85    e = 1 / d86    c = c_inf_iv()87    c_hi = float(c.b)88    na = first_tail(e, c_hi, lo, hi)89    while tail_gap_iv(na, d, c).a <= 0:90        na += 191    assert tail_gap_iv(na - 1, d, c).a <= 0 or na - 1 < mono92    assert na >= max(mono, 101)93    N = odd_up(na)94    checked = 095    while True:96        g = exact_gap_iv(N - 2, d)97        if g.a <= 0:98            break99        N -= 2100        checked += 1101    below = exact_gap_iv(N - 2, d)102    top = odd_up(na)103    n0 = N104    for M in range(n0, top + 1, 2):105        assert exact_gap_iv(M, d).a > 0106    margin = iv.mpf(1) / d - iv.log(rate_iv(n0)) / iv.log(n0)107    print(f"  {label}: tail bound holds for every base >= {na} (monotone from {mono}), closed form certified at every odd base {n0}..{top} ({(top - n0) // 2 + 1} bases), fails at {n0 - 2} with gap <= {up(below)}")108    print(f"  {label}: wall {n0}, bar - alpha_1 >= {down(margin)} there, gap >= {down(exact_gap_iv(n0, d))}")109    return n0110111def maynard_alpha(q, consecutive_half):112    L = math.log(q)113    if consecutive_half:114        return math.log((2 + 2 / L) * (2 * q / (q + 1)) * L) / L115    return math.log((1 + 3 / L) * (q / (q - 1)) * L) / L116117def maynard_cross(half):118    lo, hi = 10, 10 ** 14119    while hi - lo > 1:120        mid = (lo + hi) // 2121        if maynard_alpha(mid, half) < 0.2:122            hi = mid123        else:124            lo = mid125    return hi126127def wall():128    t0 = time.time()129    c = c_inf_iv()130    print(f"c_inf = 1 + 2/pi + (2/pi)(gamma + log(1 + sqrt 2)) + (sqrt 2 - 4/pi)/4 <= {up(c, 9)}, tail {TAIL}/base from base 101")131    worst = None132    for N in list(range(101, 3002, 2)) + [94939, 200001, 10 ** 6 + 1, 10 ** 8 + 1]:133        r = rate_iv(N)134        tl = (2 / iv.pi) * iv.log(N) + c + TAIL / N135        assert r.b < tl.a136        gap = tl - r137        worst = gap if worst is None or gap.a < worst.a else worst138    print(f"  direct: lambda/fill <= (2/pi) log base + c_inf + {TAIL}/base at every odd base 101..3001 and at 94939, 200001, 10^6 + 1, 10^8 + 1, smallest gap >= {down(worst)}")139    n5 = wall_at(5, "bar 1/5", 1000, 10 ** 7, 327)140    n4 = wall_at(4, "bar 1/4", 100, 10 ** 7, 43)141    iv.prec = 120142    low = lambda N: (2 * iv.mpf(N) / (iv.pi * (N + 1))) * iv.log(iv.mpf(N + 2) / 5) - iv.mpf(1) / 4 - ((2 / iv.pi) * iv.log(N) - iv.mpf(131) / 100)143    assert all(low(N).a > 0 for N in range(9, 10 ** 4, 2))144    assert ((2 / iv.pi) * iv.log(9) - iv.mpf(131) / 100).a > 0 and ((2 / iv.pi) * iv.log(7) - iv.mpf(131) / 100).b < 0145    tail = iv.mpf(106) / 100 - (2 / iv.pi) * iv.log(5) - (2 / iv.pi) * iv.log(iv.mpf(103) / 5) / 102146    assert tail.a > 0147    print(f"  lower bound: (2 base/(pi (base+1))) log((base+2)/5) - 1/4 >= (2/pi) log base - 1.31 > 0 at every odd base 9..9999, and from 101 by 1.06 - (2/pi) log 5 - (2/pi) log((base+2)/5)/(base+1) >= {down(tail)}; the factor is negative at 7")148    floor = lambda N: (iv.mpf(N) ** (iv.mpf(1) / 5) - 2 * iv.mpf(N) / (N + 1))149    f27 = min(N for N in range(3, 200, 2) if all(floor(M).a > 0 for M in range(N, 200, 2)))150    assert floor(f27 - 2).b < 0151    print(f"  density floor: 1 - alpha_base < 1/5 at every odd base >= {f27}, fails at {f27 - 2}")152    for N in (100003, 10 ** 6 + 3, 10 ** 9 + 7):153        r = rate_iv(N)154        a1 = iv.log(r) / iv.log(N)155        print(f"  base {N}: base^alpha_1 = lambda/fill <= {up(r, 8)}, (2/pi) log base = {2 / math.pi * math.log(N):.6f}, alpha_1 <= {up(a1, 7)}")156    one = maynard_cross(False)157    half = maynard_cross(True)158    assert one == 1520573159    print(f"  Maynard 2022 alpha_q below 1/5: one missing digit from {one} (calibration), consecutive half interval from {half}")160    print(f"wall in {time.time() - t0:.1f}s")161162# CHECK163164def grid_sums(N, ts):165    n = (N + 1) // 2166    r = np.arange(N, dtype=np.float64)167    G = np.empty(len(ts))168    S = np.empty(len(ts))169    step = max(1, 4_000_000 // N)170    for i in range(0, len(ts), step):171        t = ts[i:i + step, None]172        u = (t + r[None, :]) / N173        num = np.abs(np.sin(np.pi * n * u))174        den = np.sin(np.pi * u)175        with np.errstate(divide="ignore", invalid="ignore"):176            h = np.where(den > 0, num / np.where(den > 0, den, 1), n)177        G[i:i + step] = h.sum(1)178        S[i:i + step] = num.sum(1)179    return G, S180181def check_base(N, T):182    n = (N + 1) // 2183    ts = np.linspace(0.0, 0.5, T)184    G, S = grid_sums(N, ts)185    sig = np.mod(n * ts, 1.0)186    Sc = np.cos(np.pi * (sig - 0.5) / N) / np.sin(np.pi / (2 * N))187    assert np.max(np.abs(S - Sc) / Sc) < 1e-9188    x0 = x0_float(N)189    lam = lam_float(N)190    sn = np.sin(np.pi * ts)191    gb = n + x0 + sn * (x0 / 2 + K_2 * N / (4 * math.pi))192    low = (N / math.pi) * math.log((N + 2) / 5) - n / 4193    r1 = np.max(G / gb)194    r2 = np.max((G + S / 2) / (lam * (1 + sn / 2)))195    r3 = np.min(G) / low if low > 0 else float("inf")196    assert r1 < 1 and r2 < 1 and r3 > 1197    assert K_2 * N / (2 * math.pi) <= n + np.min(S) / 2198    L = math.log(N)199    return r1, r2, r3, G[0] / n - 2 / math.pi * L, np.max(G) / n - 2 * math.sqrt(2) / math.pi * L, ts[np.argmax(G)], np.min(G) / n - 2 / math.pi * L200201def levels_mp(N, i, s):202    mp.dps = 40203    n = (N + 1) // 2204    tot = mp.mpf(0)205    Y = N ** i206    for a in range(Y):207        v = mp.mpf(s) + mp.mpf(a) / Y208        p = mp.mpf(1)209        for j in range(i):210            w = v * N ** j211            w = w - mp.floor(w)212            if w == 0:213                p *= n214            else:215                p *= abs(mp.sin(mp.pi * n * w) / mp.sin(mp.pi * w))216        tot += p217    return tot218219def pieces(N, ts):220    n = (N + 1) // 2221    J = (n - 2) // 2222    K = (n - 1) // 2223    worst = [-1e9, -1e9, -1e9]224    for t in ts:225        m = 2 * np.arange(J + 1) + 1.0226        bp = np.sum(1 / (m + t) + 1 / (m - t)) - (math.log(N + 3) + GAMMA + K_2 * t * t)227        k = np.arange(1, K + 1, dtype=np.float64)228        le = np.sum(1 / (t + 2 * k)) if K else 0.0229        ho = 2 * k - t230        sg = np.sum(np.where(ho / N <= t, 1 / ho, -1 / ho)) if K else 0.0231        ap = le + sg - (math.log(N + 4) - 0.0757)232        dl = (t + 2 * k) / N233        dh = (2 * k - t) / N234        q = np.sum(1 / (2 * np.cos(np.pi * dl / 2))) + np.sum(1 / (2 * np.cos(np.pi * dh / 2))) if K else 0.0235        bq = q - (N / math.pi) * math.log(math.tan(3 * math.pi / 8 + math.pi / (4 * N)))236        worst = [max(worst[0], bp), max(worst[1], ap), max(worst[2], bq)]237    assert max(worst) < 0238    return worst239240def check():241    t0 = time.time()242    w = [-1e9] * 3243    for N in list(range(3, 402, 2)) + [1001, 10001, 100001]:244        p = pieces(N, np.linspace(0.0, 0.5, 201))245        w = [max(a, b) for a, b in zip(w, p)]246    print(f"pieces at every odd base 3..401 and 1001, 10001, 100001 on 201 shifts, smallest margin under the bound: bP >= {gdown(-w[0])}, aP >= {gdown(-w[1])}, bQ >= {gdown(-w[2])}")247    worst = [0, 0, 1e9]248    for N in range(3, 402, 2):249        r1, r2, r3, *_ = check_base(N, 801)250        worst = [max(worst[0], r1), max(worst[1], r2), min(worst[2], r3)]251    print(f"every odd base 3..401 at 801 shifts in [0, 1/2]: max G/bound <= {fup(worst[0])}, max T phi/(lambda phi) <= {fup(worst[1])}, min G/lower >= {fdown(worst[2])}")252    print("base, max G/bound, max T phi/(lambda phi), min G/lower, G(0)/fill - (2/pi) log base, max G/fill - (2 sqrt2/pi) log base, argmax t, min G/fill - (2/pi) log base")253    for N in (1001, 4001, 10001, 30001, 100001):254        r1, r2, r3, g0, gm, tm, gmin = check_base(N, 401)255        print(f"  {N}: <= {fup(r1)} <= {fup(r2)} >= {fdown(r3)}, readings {g0:.6f} {gm:.6f} {tm:.4f} {gmin:.6f}")256    print(f"  gamma' + 1 = {(2 / math.pi) * (GAMMA + math.log(8 / math.pi)) + 1:.6f}")257    print("levels at 40 digits, every level 1..top: base, top level, sum at the last shift, (3/2) lambda^top, its ratio")258    worst = 0.0259    for N, top in ((3, 8), (5, 5), (7, 4), (9, 4), (11, 3), (13, 3), (15, 3), (17, 3), (21, 3)):260        lam = lam_float(N)261        for i in range(1, top + 1):262            for s in (0, 0.5, 1 / (2 * N), (math.sqrt(5) - 1) / 2):263                v = levels_mp(N, i, s)264                bound = 1.5 * lam ** i265                worst = max(worst, float(v) / bound)266                assert v < bound267            if i == top:268                print(f"  {N} {i} {mp.nstr(v, 10)} {fup(bound, 2)} {fup(float(v) / bound)}")269    print(f"  largest sum/bound over all 4 shifts and levels 1..top: <= {fup(worst)}")270    err = 0.0271    for N in (5, 7, 9):272        n = (N + 1) // 2273        for s in (0, 0.5, 1 / (2 * N), (math.sqrt(5) - 1) / 2):274            t = math.fmod(N * N * s, 1.0)275            u = (t + np.arange(N)) / N276            h = np.array([abs(math.sin(math.pi * n * w) / math.sin(math.pi * w)) if w > 0 else n for w in u])277            g, _ = grid_sums(N, u)278            two = float(np.sum(h * g))279            err = max(err, abs(two - float(levels_mp(N, 2, s))) / two)280    assert err < 1e-9281    print(f"  transfer identity sum_2(s) = (T^2 1)(base^2 s) at bases 5, 7, 9 and 4 shifts, largest relative error <= {gup(err)}")282    print(f"check in {time.time() - t0:.1f}s")283284# RATE285286def power(N, M, it):287    n = (N + 1) // 2288    t = (np.arange(M) + 0.5) / M289    r = np.arange(N)290    phi = np.ones(M)291    lo = hi = 0.0292    for _ in range(it):293        new = np.zeros(M)294        for k in range(0, M, max(1, 2_000_000 // N)):295            u = (t[k:k + max(1, 2_000_000 // N), None] + r[None, :]) / N296            h = np.abs(np.sin(np.pi * n * u) / np.sin(np.pi * u))297            new[k:k + u.shape[0]] = (h * np.interp(u.ravel(), t, phi, period=1.0).reshape(u.shape)).sum(1)298        q = new / phi299        lo, hi = q.min(), q.max()300        phi = new / new.max()301    return lo / n, hi / n, phi.min()302303def rate():304    t0 = time.time()305    print("power iteration on 1000 cells, readings: base, lambda/fill - (2/pi) log base low/high, min phi/max phi, base^(1/5) - lambda/fill")306    for N in (101, 1001, 10001):307        lo, hi, pm = power(N, 1000, 30)308        print(f"  {N}: {lo - 2 / math.pi * math.log(N):.5f} {hi - 2 / math.pi * math.log(N):.5f} {pm:.4f} {N ** 0.2 - hi:.4f}")309    print("power iteration on 300 cells near the route's own crossing")310    for N in (60001, 70001, 80001):311        lo, hi, pm = power(N, 300, 12)312        print(f"  {N}: {lo - 2 / math.pi * math.log(N):.5f} {hi - 2 / math.pi * math.log(N):.5f} {pm:.4f} {N ** 0.2 - hi:.4f}")313    one = lambda N: grid_sums(N, np.array([0.5]))[0][0] / ((N + 1) // 2)314    for e in (0.25, 0.2):315        lo, hi = 1001, 2_000_001316        while hi - lo > 2:317            mid = odd_up((lo + hi) // 2)318            if one(mid) < mid ** e:319                hi = mid320            else:321                lo = mid322        print(f"  one-step reading G(1/2)/fill < base^{e} first at odd base {hi}, G(1/2)/fill - (2 sqrt2/pi) log base = {one(hi) - 2 * math.sqrt(2) / math.pi * math.log(hi):.5f} there")323    print(f"rate in {time.time() - t0:.1f}s")324325# METER326327def meter():328    t0 = time.time()329    X = 10 ** 7330    mu = np.ones(X + 1, dtype=np.int8)331    mu[0] = 0332    isp = np.ones(X + 1, dtype=bool)333    isp[:2] = False334    for p in range(2, int(X ** 0.5) + 1):335        if isp[p]:336            isp[p * p::p] = False337    for p in np.nonzero(isp)[0]:338        mu[p::p] *= -1339        if p * p <= X:340            mu[p * p::p * p] = 0341    k = np.arange(X + 1, dtype=np.int64)342    logk = np.log(np.maximum(k, 1))343    print("sanity only: base, x, A(x), M(x), max |M|/A, max |M|/sqrt(A), primes: sum log p/(kappa A) at x, kappa = p/(p+1)")344    for N in (101, 1009, 10007):345        h = (N - 1) // 2346        ok = np.ones(X + 1, dtype=bool)347        y = k.copy()348        while y.any():349            ok &= (y % N) <= h350            y //= N351        ok[0] = False352        A = np.cumsum(ok)353        M = np.cumsum(np.where(ok, mu, 0).astype(np.int64))354        sel = A >= 100355        th = np.sum(logk[ok & isp])356        kap = N / (N + 1)357        print(f"  {N}: {X} {A[-1]} {M[-1]} {np.max(np.abs(M[sel]) / A[sel]):.4f} {np.max(np.abs(M[sel]) / np.sqrt(A[sel])):.3f} {th / (kap * A[-1]):.4f}")358    print(f"meter in {time.time() - t0:.1f}s")359360if __name__ == "__main__":361    verb = sys.argv[1] if len(sys.argv) > 1 else "wall"362    {"wall": wall, "check": check, "rate": rate, "meter": meter}[verb]()