smoothed_novelty.py

1import math2import subprocess3import time4import numpy as np56N = 30_000_0007JMIN = 8.08JMAX = 23.59PER_OCTAVE = 1610ZERO_HEIGHT = 30011PI2 = math.pi ** 212T0 = time.time()1314def clock(label):15    print(f"[{time.time() - T0:6.1f}s] {label}")1617def totients(n):18    phi = np.arange(n + 1, dtype=np.int32)19    rem = np.arange(n + 1, dtype=np.int32)20    r = math.isqrt(n)21    small = np.ones(r + 1, dtype=bool)22    small[:2] = False23    for p in range(2, math.isqrt(r) + 1):24        if small[p]:25            small[p * p::p] = False26    for p in np.flatnonzero(small):27        p = int(p)28        phi[p::p] -= phi[p::p] // p29        pk = p30        while pk <= n:31            rem[pk::pk] //= p32            pk *= p33    big = rem > 134    phi[big] -= phi[big] // rem[big]35    return phi3637def totients_brute(n):38    return np.array([sum(1 for a in range(1, k + 1) if math.gcd(a, k) == 1) for k in range(n + 1)])3940def totient_sum_mobius(x):41    mu = np.ones(x + 1, dtype=np.int64)42    mu[0] = 043    prime = np.ones(x + 1, dtype=bool)44    prime[:2] = False45    for p in range(2, x + 1):46        if prime[p]:47            prime[p * p::p] = False48            mu[p::p] *= -149            mu[p * p::p * p] = 050    total = 051    for d in range(1, x + 1):52        if mu[d]:53            m = x // d54            total += int(mu[d]) * m * (m + 1) // 255    return total5657def bump_c2(u):58    v = (u - 1.0) * (2.0 - u)59    return np.where((u > 1.0) & (u < 2.0), 64.0 * v ** 3, 0.0)6061def bump_cinf(u):62    v = (u - 1.0) * (2.0 - u)63    inside = (u > 1.0) & (u < 2.0)64    safe = np.where(inside, v, 1.0)65    return np.where(inside, np.exp(4.0 - 1.0 / safe), 0.0)6667C2_POLY = 64.0 * np.polynomial.polynomial.polypow([-2.0, 3.0, -1.0], 3)6869def mellin_c2(s):70    s = np.asarray(s, dtype=complex)71    total = np.zeros_like(s)72    for k, a in enumerate(C2_POLY):73        total += a * (2.0 ** (s + k) - 1.0) / (s + k)74    return total7576GL_X, GL_W = np.polynomial.legendre.leggauss(2000)77GL_U = 1.5 + 0.5 * GL_X78GL_F = bump_cinf(GL_U) * 0.5 * GL_W7980def mellin_cinf(s):81    s = np.asarray(s, dtype=complex)82    return np.array([np.sum(GL_F * np.exp((z - 1.0) * np.log(GL_U))) for z in s.ravel()]).reshape(s.shape)8384def zeros_from_pari(height):85    script = f"""default(realprecision, 30);86z = lfunzeros(1, {height});87for(k = 1, #z, r = 1/2 + I*z[k]; a = zeta(r - 1); b = lfun(1, r, 1); print(z[k], " ", real(a), " ", imag(a), " ", real(b), " ", imag(b)))88"""89    out = subprocess.run(["gp", "-q"], input=script, capture_output=True, text=True, check=True).stdout90    rows = np.array([[float(t) for t in line.split()] for line in out.strip().splitlines()])91    gamma = rows[:, 0]92    zeta_left = rows[:, 1] + 1j * rows[:, 2]93    zeta_prime = rows[:, 3] + 1j * rows[:, 4]94    return gamma, zeta_left, zeta_prime9596def smoothed_sum(phi, y, bump):97    lo = math.ceil(1.0 / y)98    hi = math.floor(2.0 / y)99    n = np.arange(lo, hi + 1, dtype=np.float64)100    return float(np.sum(phi[lo:hi + 1].astype(np.float64) * bump(n * y)))101102def indicator_sum(prefix, y):103    lo = math.ceil(1.0 / y)104    hi = math.floor(2.0 / y)105    return int(prefix[hi] - prefix[lo - 1])106107def fit(logq, logerr):108    a, b = np.polyfit(logq, logerr, 1)109    resid = logerr - (a * logq + b)110    return a, float(np.sqrt(np.mean(resid ** 2)))111112def octave_rms(j, err):113    octaves = np.floor(j).astype(int)114    rows = []115    for k in sorted(set(octaves)):116        sel = octaves == k117        if sel.sum() >= PER_OCTAVE // 2:118            rows.append((k, float(np.sqrt(np.mean(err[sel] ** 2)))))119    return np.array(rows)120121def slopes(name, j, err):122    logq = -2.0 * j * math.log(2.0)123    rows = octave_rms(j, err)124    ok = np.abs(err) > 0125    a_all, r_all = fit(logq[ok], np.log(np.abs(err[ok])))126    lq = -2.0 * (rows[:, 0] + 0.5) * math.log(2.0)127    le = np.log(rows[:, 1])128    half = len(rows) // 2129    a_lo, r_lo = fit(lq[:half], le[:half])130    a_hi, r_hi = fit(lq[half:], le[half:])131    a_oct, r_oct = fit(lq, le)132    print(f"{name}: octave rms slope in q {a_oct:.4f} (residual {r_oct:.3f}) over j in [{rows[0, 0]:.0f}, {rows[-1, 0] + 1:.0f}]")133    print(f"{name}: lower window j in [{rows[0, 0]:.0f}, {rows[half - 1, 0] + 1:.0f}] slope {a_lo:.4f} (residual {r_lo:.3f}); upper window j in [{rows[half, 0]:.0f}, {rows[-1, 0] + 1:.0f}] slope {a_hi:.4f} (residual {r_hi:.3f})")134    print(f"{name}: all-sample slope of log|E| in q {a_all:.4f} (residual {r_all:.3f}) on {int(ok.sum())} samples")135    return a_oct136137def main():138    clock(f"sieve to {N}")139    phi = totients(N)140    prefix = np.cumsum(phi, dtype=np.int64)141    clock("sieve done")142    brute = totients_brute(2000)143    print(f"brute-force totients agree to 2000: {bool(np.array_equal(brute, phi[:2001]))}")144    x = 10 ** 6145    print(f"totient sum at {x}: sieve {int(prefix[x])}, mobius route {totient_sum_mobius(x)}, equal {int(prefix[x]) == totient_sum_mobius(x)}")146    j = np.arange(JMIN, JMAX + 1e-9, 1.0 / PER_OCTAVE)147    y = 2.0 ** (-j)148    c_ind = 6.0 / PI2 * 1.5149    c_c2 = 6.0 / PI2 * (24.0 / 35.0)150    f2_cinf = float(np.sum(GL_F * GL_U))151    c_cinf = 6.0 / PI2 * f2_cinf152    print(f"mellin transforms at 2: indicator 3/2, c2 bump 24/35 = {24 / 35:.12f} (closed form {mellin_c2(2.0).real:.12f}), cinf bump {f2_cinf:.12f}")153    e_ind = np.array([yy * yy * indicator_sum(prefix, yy) - c_ind for yy in y])154    clock("indicator sums done")155    e_c2 = np.array([yy * yy * smoothed_sum(phi, yy, bump_c2) - c_c2 for yy in y])156    clock("c2 bump sums done")157    e_cinf = np.array([yy * yy * smoothed_sum(phi, yy, bump_cinf) - c_cinf for yy in y])158    clock("cinf bump sums done")159    yy = y[-1]160    lo = math.ceil(1.0 / yy)161    hi = math.floor(2.0 / yy)162    n = np.arange(lo, hi + 1, dtype=np.float64)163    prod = phi[lo:hi + 1].astype(np.float64) * bump_cinf(n * yy)164    noise = abs(float(np.sum(prod)) - math.fsum(prod.tolist())) * yy * yy165    print(f"summation noise at j = {JMAX}: |pairwise - fsum| scaled by y^2 is {noise:.3e} against |E_cinf| = {abs(e_cinf[-1]):.3e}")166    print()167    print("SLOPES")168    a_ind = slopes("indicator", j, e_ind)169    a_c2 = slopes("c2 bump", j, e_c2)170    a_cinf = slopes("cinf bump", j, e_cinf)171    print(f"verdict: indicator {a_ind:.3f} against 1/2, c2 bump {a_c2:.3f} against 3/4, cinf bump {a_cinf:.3f} against 3/4")172    print()173    print("EXPLICIT FORMULA")174    gamma, zeta_left, zeta_prime = zeros_from_pari(ZERO_HEIGHT)175    rho = 0.5 + 1j * gamma176    print(f"{len(gamma)} zeros to height {ZERO_HEIGHT} from PARI, first {gamma[0]:.6f}, last {gamma[-1]:.6f}")177    for name, mellin, err in (("c2 bump", mellin_c2, e_c2), ("cinf bump", mellin_cinf, e_cinf)):178        coef = mellin(rho) * zeta_left / zeta_prime179        print(f"{name}: |c_rho| at zeros 1, 2, 10, 50, {len(gamma)}: " + ", ".join(f"{abs(coef[k]):.3e}" for k in (0, 1, 9, 49, len(gamma) - 1)))180        power = np.polyfit(np.log(gamma), np.log(np.abs(coef)), 1)[0]181        print(f"{name}: least-squares power of |c_rho| against gamma on {len(gamma)} zeros: {power:.2f}")182        for K in (1, 10, 30, len(gamma)):183            model = np.array([2.0 * np.sum((coef[:K] * yy ** (2.0 - rho[:K])).real) for yy in y])184            resid = np.max(np.abs(err - model)) / np.max(np.abs(err))185            print(f"{name}: first {K:3d} zeros, max |E - sum| / max |E| = {resid:.3e}")186        amp = np.abs(err) / y ** 1.5187        print(f"{name}: |E|/y^(3/2) over the grid: min {amp.min():.4e}, max {amp.max():.4e}, bound 2 sum |c_rho| = {2 * np.sum(np.abs(coef)):.4e}")188    print()189    print("SHARP CUTOFF")190    primes = [p for p in (10 ** 6 + 3, 10 ** 7 + 19) if phi[p] == p - 1]191    for p in primes:192        before = prefix[p - 1] - 3.0 / PI2 * (p - 1) ** 2193        after = prefix[p] - 3.0 / PI2 * p ** 2194        print(f"prime {p}: totient error jumps from {before:.1f} to {after:.1f}, jump {after - before:.1f} against phi(p) - 3(2p - 1)/pi^2 = {p - 1 - 3.0 / PI2 * (2 * p - 1):.1f}")195    clock("done")196197if __name__ == "__main__":198    main()