shifted_sums.py

1import subprocess2import sys3import time4from concurrent.futures import ThreadPoolExecutor56TOP = 267LOW = 168WORKERS = 89SHOW = 4010CONTROL_TOP = 1611CHECK = 1812PAGE_T = [0, -1, -3, -4, 0, -6, -5, -27, -21, -25, -51, -118, -168, -178, -257, -169]13PAGE_ENDPOINT = [0, -1, -1, -1, 3, -6, 0, -22, 6, -4, -27, -67, -49, -10, -78, 88, 82, 209, 543, 858, 335, -407, 3030]14LOWS = f"L=vector(2^{LOW},r,6*fromdigits(binary(r-1),3));"1516# PARI171819def gp(script):20    out = subprocess.run(["gp", "-q"], input=script + "\n", capture_output=True, text=True, check=True).stdout21    return [[int(x) for x in line.split()] for line in out.split("\n") if line.strip()]222324def fan(scripts):25    with ThreadPoolExecutor(len(scripts)) as pool:26        return [row for part in pool.map(gp, scripts) for row in part]272829def head(j):30    return f"h=6*fromdigits(binary({j}),3)*3^{LOW}+1;"313233def sums_script(js):34    body = [LOWS]35    for j in js:36        body.append(37            f"{head(j)} p=0; mn=0; mx=0; "38            f"for(r=1,2^{LOW}, p+=moebius(h+L[r]); if(p<mn,mn=p); if(p>mx,mx=p)); "39            f'print({j}," ",p," ",mn," ",mx);'40        )41    return "\n".join(body)424344def cross_script(jobs):45    body = [LOWS]46    for j, start, ls in jobs:47        body.append(48            f"{head(j)} p={start}; ls={ls}; "49            f'for(r=1,2^{LOW}, p+=moebius(h+L[r]); if(p!=0&&sign(p)!=ls, if(ls!=0, print({j}," ",r," ",sign(p)," ",h+L[r])); ls=sign(p)));'50            f'print({j}," ",0," ",ls," ",p);'51        )52    return "\n".join(body)535455def low_levels():56    return gp(LOWS + f" p=0; s=0; for(r=1,2^{LOW}, p+=moebius(1+L[r]); if(r==2^s, print(s,\" \",p); s++));")575859def sequential(top):60    return gp(f"p=0; ls=0; k=0; at=0; for(i=0,2^{top}-1, p+=moebius(6*fromdigits(binary(i),3)+1); if(p!=0&&sign(p)!=ls, if(ls!=0, k++; at=i+1); ls=sign(p))); print(p,\" \",k,\" \",at);")[0]616263def spread(items):64    return [items[k::WORKERS] for k in range(WORKERS) if items[k::WORKERS]]6566# THE WALK676869def chunk_sums():70    got = fan([sums_script(js) for js in spread(list(range(2 ** (TOP - LOW))))])71    return {j: (p, mn, mx) for j, p, mn, mx in got}727374def candidates(sums):75    todo, starts, ls, start = [], {}, 0, 076    for j in range(len(sums)):77        p, mn, mx = sums[j]78        starts[j] = (start, ls)79        if ls == 0 or (ls > 0 and start + mn < 0) or (ls < 0 and start + mx > 0):80            todo.append(j)81            end = start + p82            if end != 0:83                ls = 1 if end > 0 else -184            else:85                ls = gp(cross_script([(j, start, ls)]))[-1][2]86        start += p87    return todo, starts888990def crossings(sums):91    todo, starts = candidates(sums)92    got = fan([cross_script([(j, *starts[j]) for j in js]) for js in spread(todo)])93    flips = sorted((j * 2**LOW + r, sg, n) for j, r, sg, n in got if r > 0)94    return flips, len(todo)959697def positive_share(flips, n):98    pos, sign, at = 0, 1, 199    for m, sg, _ in flips:100        if m > n:101            break102        if sign > 0:103            pos += m - at104        sign, at = sg, m105    if sign > 0:106        pos += n + 1 - at107    return pos / n108109110def verb_walk():111    t0 = time.time()112    low = dict(low_levels())113    sums = chunk_sums()114    t1 = time.time()115    flips, todo = crossings(sums)116    t2 = time.time()117    level = dict(low)118    run = 0119    for j in range(len(sums)):120        run += sums[j][0]121        if (j + 1) & j == 0:122            level[LOW + (j + 1).bit_length() - 1] = run123    assert [level[s] for s in range(1, 17)] == PAGE_T124    mu = dict(gp("\n".join(f'print({t}," ",moebius(3^{t}+1));' for t in range(2, TOP + 2))))125    ends = [level[t - 1] - level[t - 2] + mu[t] for t in range(2, 25)]126    assert ends == PAGE_ENDPOINT127    early = [m for m, _, _ in flips if m <= 2**CHECK]128    assert sequential(CHECK) == [level[CHECK], len(early), early[-1]]129    print(f"T_s = sum_(w in B_s) mu(6w + 1), B_s the 2^s base 3 {{0,1}}-strings of length s, walked in ascending w to s = {TOP}")130    print(f"checks: T_1..T_16 equal the page's list; T_(t-1) - T_(t-2) + mu(3^t + 1) equals the page's M_F(x_t; R_t) at t = 2..24; a sequential walk to 2^{CHECK} repeats T_{CHECK}, the flip count and the last flip")131    print("  s |        2^s |       top element |       T_s | T_s/2^(s/2) | flips n <= 2^s | share of n <= 2^s with last nonzero sign +")132    for s in range(TOP + 1):133        n = 2**s134        k = sum(1 for m, _, _ in flips if m <= n)135        print(f"{s:3d} | {n:10d} | {3 ** (s + 1) - 2:17d} | {level[s]:9d} | {level[s] / 2 ** (s / 2):11.4f} | {k:14d} | {positive_share(flips, n):.4f}")136    neg = [s for s in range(1, TOP + 1) if level[s] < 0]137    pos = [s for s in range(1, TOP + 1) if level[s] > 0]138    print(f"T_s < 0 at s = {neg}")139    print(f"T_s > 0 at s = {pos}")140    print(f"T_s = 0 at s = {[s for s in range(1, TOP + 1) if level[s] == 0]}")141    print(f"cut: n = 2^{TOP} = {2**TOP} terms, largest element 3^{TOP + 1} - 2 = {3 ** (TOP + 1) - 2}, P = {level[TOP]}")142    print(f"sign changes of the prefix walk P_n (last nonzero sign flips): {len(flips)} in all, {todo} chunks of 2^{LOW} rewalked")143    print(f"the last {min(SHOW, len(flips))}: n, new sign, element 6w + 1")144    for m, sg, n in flips[-SHOW:]:145        print(f"  n = {m:10d}  {'+' if sg > 0 else '-'}  6w + 1 = {n}")146    print(f"walk {t1 - t0:.1f} s, rewalk {t2 - t1:.1f} s, {WORKERS} gp processes")147148# THE CONTROL149150151def smooth23(x):152    count, a = 0, 1153    while a <= x:154        b = a155        while b <= x:156            count, b = count + 1, 3 * b157        a *= 2158    return count159160161def verb_control():162    t0 = time.time()163    got = gp(164        "p=0; ls=0; ch=0; ps=0; c=0; s=0; nx=1; "165        f"forstep(n=1,3^{CONTROL_TOP}-2,6, p+=moebius(n); c++; if(p!=0&&sign(p)!=ls, if(ls!=0,ch++); ls=sign(p)); if(ls>0,ps++); "166        'if(n==nx, print(s," ",p," ",ch," ",ps," ",c); s++; nx=3^(s+1)-2));'167    )168    print("control: the whole class, M(x; 6, 1) = sum_(n <= x, n = 1 mod 6) mu(n), at x = 3^(s+1) - 2")169    print("drift = 3/4 - Psi_(2,3)(x), the s = 0 residue of M(x; 6, 1) with M's constant -2 and 1/(L(0, chi_-3)(1 + 2^0)) = 3/2")170    print("  s |                x | M(x; 6, 1) |   drift | M/sqrt(x/6) | flips to x | share with last nonzero sign +")171    for s, p, ch, ps, c in got:172        x = 3 ** (s + 1) - 2173        print(f"{s:3d} | {x:16d} | {p:10d} | {0.75 - smooth23(x):7.2f} | {p / (x / 6) ** 0.5:11.4f} | {ch:10d} | {ps / c:.4f}")174    print(f"control {time.time() - t0:.1f} s")175176177def main():178    verbs = {"walk": verb_walk, "control": verb_control}179    for v in sys.argv[1:] or list(verbs):180        verbs[v]()181182183if __name__ == "__main__":184    main()