ratio.py

1import argparse2import time3from math import ceil, floor, gcd, log45import numpy as np67def up(x, d):8    return ceil(x * 10 ** d) / 10 ** d910def down(x, d):11    return floor(x * 10 ** d) / 10 ** d1213# INCREMENT AUTOMATON1415def level(z1, z2):16    start = z1 // 317    if start == 0:18        return 0, 019    dist = {start: 1}20    frontier = [start]21    d = 122    while frontier:23        d += 124        nxt = []25        for j in frontier:26            if j % 3 == 0:27                moves = (j // 3, (j + z1) // 3)28            elif (j - z2) % 3 == 0:29                moves = ((j - z2) // 3,)30            else:31                continue32            for k in moves:33                if k == 0:34                    return d, len(dist)35                if k not in dist:36                    dist[k] = d37                    nxt.append(k)38        frontier = nxt39    return 0, len(dist)4041def reach_size(z1, z2):42    start = z1 // 343    if start == 0:44        return 045    seen = {0, start}46    stack = [start]47    while stack:48        j = stack.pop()49        if j % 3 == 0:50            moves = (j // 3, (j + z1) // 3)51        elif (j - z2) % 3 == 0:52            moves = ((j - z2) // 3,)53        else:54            continue55        for k in moves:56            if k not in seen:57                seen.add(k)58                stack.append(k)59    return len(seen)6061def band_cap(z1, z2):62    return (z1 - 1) // 2 + (z2 - 1) // 2 + 16364def succ(j, z1, z2):65    if j % 3 == 0:66        return (j // 3, (j + z1) // 3)67    if (j - z2) % 3 == 0:68        return ((j - z2) // 3,)69    return ()7071def reach(z1, z2):72    start = z1 // 373    seen = {start}74    stack = [start]75    hit = 076    while stack:77        for k in succ(stack.pop(), z1, z2):78            if k == 0:79                hit = 180            if k not in seen:81                seen.add(k)82                stack.append(k)83    return hit, len(seen)8485def core(z1, z2):86    lo = -((z2 - 1) // 2)87    hi = (z1 - 1) // 288    deg = {}89    for j in range(lo, hi + 1):90        deg[j] = sum(1 for k in succ(j, z1, z2) if lo <= k <= hi)91    dead = [j for j in deg if j and not deg[j]]92    while dead:93        j = dead.pop()94        del deg[j]95        for p in (3 * j, 3 * j - z1, 3 * j + z2):96            if lo <= p <= hi and p in deg and p != 0 and j in succ(p, z1, z2):97                deg[p] -= 198                if not deg[p]:99                    dead.append(p)100    live = set(deg)101    start = z1 // 3102    seen = {start}103    stack = [start]104    ok = start in live105    while stack:106        for k in succ(stack.pop(), z1, z2):107            if k in live:108                ok = True109            if k not in seen:110                seen.add(k)111                stack.append(k)112    return ok, len(live)113114def base3(n):115    s = ""116    while n:117        s = str(n % 3) + s118        n //= 3119    return s or "0"120121BETA = log(2) / log(3)122123def binaries(k):124    out = [0]125    for i in range(k):126        out = out + [a + 3 ** i for a in out]127    return out128129def pairs(cap):130    for z1 in range(3, cap + 1, 3):131        for z2 in range(1, cap + 1):132            if z2 % 3 and gcd(z1, z2) == 1:133                yield z1, z2134135# PAIR CARRY AUTOMATON136137OUT = ((0, 0), (0, 1), (1, 0))138139def carry_occupied(z1, z2):140    seen = set()141    stack = []142    for d in (1, 2):143        s1 = d * z1144        s2 = d * z2145        if (s1 % 3, s2 % 3) in OUT:146            key = (s1 // 3) * z2 + (s2 // 3)147            if key not in seen:148                seen.add(key)149                stack.append((s1 // 3, s2 // 3))150    while stack:151        c1, c2 = stack.pop()152        e1 = c1 % 3153        if e1 == 2:154            continue155        for e2 in ((0, 1) if e1 == 0 else (0,)):156            d = (e2 - c2) * pow(z2 % 3, -1, 3) % 3157            n1 = (d * z1 + c1) // 3158            n2 = (d * z2 + c2) // 3159            if n1 == 0 and n2 == 0:160                return True, len(seen)161            key = n1 * z2 + n2162            if key not in seen:163                seen.add(key)164                stack.append((n1, n2))165    return False, len(seen)166167# BRUTE GASKET168169def brute(n, cap):170    out = {}171    for lab in range(3 ** n):172        u = v = 0173        w = lab174        for i in range(n):175            d = w % 3176            w //= 3177            if d == 1:178                u += 3 ** i179            elif d == 2:180                v += 3 ** i181        if u and v:182            g = gcd(u, v)183            z = (u // g, v // g)184            if max(z) <= cap:185                k = 1186                while 3 ** k <= max(u, v):187                    k += 1188                if z not in out or k < out[z]:189                    out[z] = k190    return out191192# SWEEP193194def cmd_sweep(args):195    cap = args.cap196    t = time.time()197    hs = []198    ws = []199    ls = []200    for z1, z2 in pairs(cap):201        d, _ = level(z1, z2)202        if d:203            hs.append(max(z1, z2))204            ws.append(z1 + z2)205            ls.append(d)206    secs = time.time() - t207    h = np.array(hs, dtype=np.int64)208    w = np.array(ws, dtype=np.int64)209    lv = np.array(ls, dtype=np.int64)210    print("convention height max(z1,z2), ordered directions, occupied at some level, no level cap")211    print("x A(x) theta local")212    ladder = []213    x = 32214    while x <= cap:215        ladder.append(x)216        x *= 2217    loc = []218    prev = None219    for x in ladder:220        a = 2 * int((h <= x).sum())221        th = log(a) / log(x)222        if prev:223            loc.append(log(a / prev[1]) / log(x / prev[0]))224        print(x, a, round(th, 4), round(loc[-1], 4) if prev else "-")225        prev = (x, a)226    print("theta band", down(min(log(2 * (h <= x).sum()) / log(x) for x in ladder[-4:]), 4),227          up(max(log(2 * (h <= x).sum()) / log(x) for x in ladder[-4:]), 4),228          "local band", down(min(loc[-3:]), 4), up(max(loc[-3:]), 4),229          "over x", ladder[-4], ladder[-1])230    print("W Zsum Zmean Zmax argmax logZmax/logW")231    cnt = np.bincount(w, minlength=cap + 1)[: cap + 1]232    zexp = []233    for x in ladder:234        seg = cnt[: x + 1]235        top = 2 * int(seg.max())236        zexp.append(log(top) / log(x))237        print(x, int(2 * seg.sum()), round(2 * seg.sum() / x, 4), top, int(seg.argmax()),238              round(zexp[-1], 4))239    print("Zmax exponent band", down(min(zexp), 4), up(max(zexp), 4),240          "over W", ladder[0], ladder[-1], "needed below", 0.8073)241    print("x meanlev maxlev meanc maxc share_lev_below_1.8073_log3_x")242    for x in ladder:243        sel = h <= x244        c = lv[sel] / (np.log(h[sel]) / log(3))245        print(x, round(float(lv[sel].mean()), 3), int(lv[sel].max()),246              round(float(c.mean()), 3), round(float(c.max()), 3),247              round(float((lv[sel] <= 1.8073 * log(x) / log(3)).mean()), 4))248    print("secs", round(secs, 1))249250# LEVELS251252def cmd_levels(args):253    cap = args.cap254    t = time.time()255    tab = np.zeros(200, dtype=np.int64)256    for z1, z2 in pairs(cap):257        d, _ = level(z1, z2)258        if d:259            tab[d] += 2260    cum = np.cumsum(tab)261    print("cap", cap, "A(inf)", int(cum[-1]))262    print("n A(n,cap) new")263    for n in range(args.lo, args.hi + 1):264        print(n, int(cum[n]), int(tab[n]))265    print("secs", round(time.time() - t, 1))266267# MISSING DIGIT MULTIPLES268269def dcount(n, q):270    v = np.zeros(q, dtype=np.int64)271    v[0] = 1272    idx = np.arange(q)273    a = (3 * idx) % q274    b = (3 * idx + 1) % q275    for _ in range(n):276        nv = np.zeros(q, dtype=np.int64)277        nv[a] = v278        nv[b] += v279        v = nv280    return int(v[0])281282def cmd_multiples(args):283    n = args.n284    print("uniform test D_n(q) q / 2^n over q <= Q coprime to 3, n =", n)285    worst = (0.0, 0)286    for q in range(2, args.qmax + 1):287        if q % 3 == 0:288            continue289        r = dcount(n, q) * q / 2 ** n290        if r > worst[0]:291            worst = (r, q)292    print("Q", args.qmax, "max ratio", round(worst[0], 4), "at q", worst[1])293    print("h q=1+3^h D_2h(q) 2^h 2^(2h)/q ratio")294    for h in range(1, args.hmax + 1):295        q = 1 + 3 ** h296        d = dcount(2 * h, q)297        print(h, q, d, 2 ** h, round(4 ** h / q, 2), round(d * q / 4 ** h, 3))298299# BAND300301def cmd_band(args):302    t = time.time()303    big = (0, None)304    frac = (0, 1, None)305    bad = 0306    for z1, z2 in pairs(args.cap):307        r = reach_size(z1, z2)308        cap = band_cap(z1, z2)309        if r > cap:310            bad += 1311        if r > big[0]:312            big = (r, (z1, z2))313        if r * frac[1] > frac[0] * cap:314            frac = (r, cap, (z1, z2))315    print("cap", args.cap, "violations of the halved cap", bad)316    print("largest reachable set", big[0], "at", big[1], "halved cap", band_cap(*big[1]),317          "fraction", up(big[0] / band_cap(*big[1]), 4))318    print("fullest reachable set", frac[0], "of", frac[1], "at", frac[2],319          "fraction", up(frac[0] / frac[1], 4))320    print("secs", round(time.time() - t, 1))321322# WEIGHT LAYERS323324def cmd_layers(args):325    cap = args.cap326    t = time.time()327    Z = [0] * (cap + 1)328    bad = 0329    for z1 in range(3, cap + 1, 3):330        for z2 in range(1, cap + 1):331            if z2 % 3 == 0 or gcd(z1, z2) != 1:332                continue333            w = z1 + z2334            if w > cap:335                break336            if level(z1, z2)[0]:337                Z[w] += 2338                if w <= 3 * z1 <= 2 * w:339                    bad += 1340    secs = time.time() - t341    print("cap", cap, "occupied directions with z1/w inside [1/3,2/3]", bad)342    print("octave argmax_Z Zmax base3 argmax_exp exp base3 argmax_const const base3")343    x = 32344    binary = True345    while 2 * x - 1 <= cap:346        seg = range(x, 2 * x)347        mz = max(seg, key=lambda w: Z[w])348        me = max(seg, key=lambda w: log(Z[w]) / log(w) if Z[w] else 0.0)349        mc = max(seg, key=lambda w: Z[w] / w ** BETA)350        print(x, mz, Z[mz], base3(mz), me, up(log(Z[me]) / log(me), 4), base3(me),351              mc, up(Z[mc] / mc ** BETA, 4), base3(mc))352        binary = binary and all(set(base3(a)) <= {"0", "1"} for a in (mz, me, mc))353        x *= 2354    print("every argmax of all three columns binary base 3", binary)355    e = max(range(4, cap + 1), key=lambda w: log(Z[w]) / log(w) if Z[w] else 0.0)356    c = max(range(4, cap + 1), key=lambda w: Z[w] / w ** BETA)357    print("max exp", up(log(Z[e]) / log(e), 4), "at", e,358          "max const", up(Z[c] / c ** BETA, 4), "at", c,359          "needed exp below", 0.8073, "secs", round(secs, 1))360361def cmd_weights(args):362    print("w Z(w) phi3 logZ/logw Z/w^(log2/log3) meanreach/sqrt(w) secs")363    for w in args.w:364        t = time.time()365        z = n = r = 0366        for z1 in range(3, w, 3):367            if gcd(z1, w) != 1:368                continue369            n += 1370            hit, s = reach(z1, w - z1)371            z += 2 * hit372            r += s373        print(w, z, 2 * n, up(log(z) / log(w), 4), up(z / w ** BETA, 4),374              up(r / n / w ** 0.5, 4), round(time.time() - t, 1))375376def cmd_core(args):377    print("w Z(w) Zinf(w) Zinf/Z coremean band secs")378    for w in args.w:379        t = time.time()380        z = h = n = c = 0381        for z1 in range(3, w, 3):382            if gcd(z1, w) != 1:383                continue384            n += 1385            z += 2 * reach(z1, w - z1)[0]386            ok, sz = core(z1, w - z1)387            h += 2 * ok388            c += sz389        print(w, z, h, up(h / z, 4), round(c / n, 1), (w - 1) // 2,390              round(time.time() - t, 1))391392# SLOPE COVER393394def tri(k, base):395    a = np.zeros(1, dtype=np.int64)396    b = np.zeros(1, dtype=np.int64)397    for i in range(k):398        p = 3 ** (base + i)399        a = np.concatenate([a, a + p, a])400        b = np.concatenate([b, b, b + p])401    return a, b402403def cover(n, e, low=13):404    m = n + e405    if 3 * 3 ** (n + m) // 2 + 3 ** n >= 2 ** 63:406        raise SystemExit("int64 overflow: reduce --extra or n")407    low = min(low, m)408    al, bl = tri(low, 0)409    ah, bh = tri(m - low, low)410    N = 3 ** n411    M = 3 ** m412    hit = np.zeros(N + 1, dtype=bool)413    u = np.empty(len(al), dtype=np.int64)414    v = np.empty(len(al), dtype=np.int64)415    c = np.empty(len(al), dtype=np.int64)416    d = np.empty(len(al), dtype=np.int64)417    for i in range(len(ah)):418        np.add(al, ah[i], out=u)419        np.add(bl, bh[i], out=v)420        np.add(u, M, out=c)421        np.add(c, v, out=d)422        c *= N423        c //= d424        hit[c] = True425        np.add(u, v, out=d)426        d += M427        np.multiply(u, N, out=c)428        c //= d429        hit[c] = True430    return int(hit[:N].sum())431432def cmd_box(args):433    print("cover counts both swap halves and is a lower estimate of the saturated cover")434    print("saturation at n =", args.lo, "over extra digits 0 ..", args.extra)435    for e in range(args.extra + 1):436        print(args.lo, e, cover(args.lo, e))437    print("n cover(3^-n) 3^n cover*n/3^n log_3 cover / n local secs")438    prev = None439    for n in range(args.lo, args.hi + 1):440        t = time.time()441        c = cover(n, args.extra)442        print(n, c, 3 ** n, down(c * n / 3 ** n, 4), down(log(c) / (n * log(3)), 4),443              "-" if prev is None else down(log(c / prev) / log(3), 4),444              round(time.time() - t, 1))445        prev = c446447# REPUNIT448449def witness(z1, z2):450    start = z1 // 3451    par = {start: (None, z1)}452    frontier = [start]453    while frontier:454        nxt = []455        for j in frontier:456            for k in succ(j, z1, z2):457                a = k * 3 - j458                if k == 0:459                    path = [a]460                    cur = j461                    while cur is not None:462                        path.append(par[cur][1])463                        cur = par[cur][0]464                    e = sum(3 ** i for i, b in enumerate(reversed(path)) if b)465                    return e // (z1 + z2)466                if k not in par:467                    par[k] = (j, a)468                    nxt.append(k)469        frontier = nxt470    return 0471472def factor(n):473    out = []474    p = 2475    while p * p <= n:476        if n % p == 0:477            out.append(p)478            while n % p == 0:479                n //= p480        p += 1481    if n > 1:482        out.append(n)483    return out484485def order3(q):486    d = 1487    x = 3 % q488    while x != 1:489        x = x * 3 % q490        d += 1491    return d492493def count_multiples(k, q):494    v = np.zeros(q, dtype=np.int64)495    v[0] = 1496    idx = np.arange(q)497    for i in range(k):498        nv = v.copy()499        np.add.at(nv, (idx + pow(3, i, q)) % q, v)500        v = nv501    return int(v[0])502503def fourier_multiples(k, q):504    d = order3(q)505    t = np.arange(q)506    P = np.ones(q, dtype=complex)507    for r in range(d):508        P *= 1 + np.exp(2j * np.pi * t * pow(3, r, q) / q)509    return (P ** (k // d)).sum() / q510511def floor_formula(k):512    ps = factor((3 ** k - 1) // 2)513    exact = 0514    four = 0515    for mask in range(2 ** len(ps)):516        q = 1517        for i, p in enumerate(ps):518            if mask >> i & 1:519                q *= p520        sign = (-1) ** bin(mask).count("1")521        exact += sign * count_multiples(k, q)522        four += sign * (2 ** k if q == 1 else fourier_multiples(k, q))523    return ps, exact, four524525def lift_union(k):526    w = (3 ** k - 1) // 2527    pw = np.array([3 ** i for i in range(k)], dtype=np.int64)528    S = np.arange(2 ** k)529    bits = ((S[:, None] >> np.arange(k)) & 1).astype(np.int64)530    occ = set()531    total = 0532    heur = 0.0533    for T in range(0, 2 ** k, 2):534        tb = np.array([(T >> i) & 1 for i in range(k)], dtype=np.int64)535        m = 1 + 2 * int(tb.dot(pw))536        heur += 2 ** k / m537        A = bits.dot(pw * (1 - tb) + pw * tb * 3 ** k)538        for a in A[A % m == 0]:539            z = int(a) // m540            if 0 < z < w and gcd(z, w) == 1:541                occ.add(z)542                total += 1543    return occ, total, heur544545def cmd_repunit(args):546    print("floor Phi_k = #{submask a of R_k : 0 < a < R_k, gcd(a, R_k) = 1}; Occ_T = directions with a witness K_T, "547          "T inside [1, k-1]; lift = |union of Occ_T|, liftsum = Sum |Occ_T|, heur = Sum 2^k / m_T; "548          "deep = Z - lift; delta = Phi / 2^k; X = (Z - Phi) / Phi; maxdigits = longest minimal witness m R_k, maxcol = most uses of one column mod k by a minimal witness")549    print("k R_k primes Phi_k formula fourier Z(R_k) excess lift liftsum heur deep Z/R_k^beta delta X maxdigits maxcol secs")550    lists = {}551    for k in range(2, args.kmax + 1):552        t = time.time()553        w = (3 ** k - 1) // 2554        subs = set(a for a in binaries(k) if 0 < a < w and gcd(a, w) == 1)555        ps, exact, four = floor_formula(k)556        occ = []557        for z1 in range(3, w, 3):558            if gcd(z1, w) != 1:559                continue560            m = witness(z1, w - z1)561            if m:562                occ.append((z1, m))563        Z = 2 * len(occ)564        if k <= args.liftmax:565            lu, liftsum, heur = lift_union(k)566            lift, deep, heur = len(lu), Z - len(lu), round(heur, 1)567        else:568            lift = deep = liftsum = heur = "-"569        digits = maxcol = 0570        for z1, m in occ:571            if z1 in subs:572                continue573            cols = [0] * k574            for i, ch in enumerate(reversed(base3(m * w))):575                cols[i % k] += ch == "1"576            digits = max(digits, len(base3(m * w)))577            maxcol = max(maxcol, max(cols))578        print(k, w, ps, len(subs), exact, round(four.real, 3) if abs(four.imag) < 1e-6 else four,579              Z, Z - len(subs), lift, liftsum, heur, deep, up(Z / w ** BETA, 4), down(len(subs) / 2 ** k, 4),580              down((Z - len(subs)) / len(subs), 4), digits, maxcol, round(time.time() - t, 1))581        if k in (7, 8, 9):582            lists[k] = [(z1, w - z1, m) for z1, m in occ if z1 not in subs]583    for k, rows in lists.items():584        print("non-submask occupied directions of weight R_" + str(k), "=", (3 ** k - 1) // 2,585              "one of each swap pair, 3 | z1; columns z1 z2 m base3(z1) base3(z2) base3(m) base3(m z1) base3(m z2) base3(m w)")586        for z1, z2, m in rows:587            print(z1, z2, m, base3(z1), base3(z2), base3(m), base3(m * z1), base3(m * z2), base3(m * (z1 + z2)))588589# LIFT UNION590591def subset_sums(ws):592    out = np.zeros(1, dtype=np.int64)593    for x in ws:594        out = np.concatenate((out, out + x))595    return out596597def submask_floor(k):598    w = (3 ** k - 1) // 2599    b = subset_sums([3 ** i for i in range(k)])600    return int(((b > 0) & (b < w) & (np.gcd(b, w) == 1)).sum())601602def lift_scan(k):603    w = (3 ** k - 1) // 2604    h = max(k // 2, 1)605    v1 = {}606    for tl in range(0, 1 << h, 2):607        v1[tl] = subset_sums([3 ** (i + k) if tl >> i & 1 else 3 ** i for i in range(h)])608    v2 = {}609    for th in range(1 << (k - h)):610        v2[th] = subset_sums([3 ** (i + k) if th >> (i - h) & 1 else 3 ** i for i in range(h, k)])611    occ = set()612    total = 0613    model = 0.0614    best = (0.0, 0, 0)615    for th in range(1 << (k - h)):616        ah = sum(3 ** i for i in range(h, k) if th >> (i - h) & 1)617        V2 = v2[th]618        for tl in range(0, 1 << h, 2):619            al = sum(3 ** i for i in range(h) if tl >> i & 1)620            T = tl | (th << h)621            m = 1 + 2 * (al + ah)622            model += 2 ** k / m623            z = mod_match(v1[tl], V2, m, w)624            total += int(z.size)625            occ.update(z.tolist())626            ratio = z.size * m / 2 ** k627            if ratio > best[0]:628                best = (ratio, m, T)629    return occ, total, model, best630631def occ_of(k, T):632    w = (3 ** k - 1) // 2633    m = 1 + 2 * sum(3 ** i for i in range(k) if T >> i & 1)634    A = subset_sums([3 ** (i + k) if T >> i & 1 else 3 ** i for i in range(k)])635    z = A[A % m == 0] // m636    z = z[(z > 0) & (z < w) & (np.gcd(z, w) == 1)]637    return m, set(z.tolist())638639def cyclotomic_occ(t):640    k = 2 * t + 1641    w = (3 ** k - 1) // 2642    m = 3 ** (2 * t) - 3 ** t + 1643    out = set()644    for S in range(1, 1 << (t - 1)):645        c = sum(3 ** (1 + i) for i in range(t - 1) if S >> i & 1)646        z = c * (3 ** t + 1)647        if gcd(z, w) == 1:648            out.add(z)649            out.add(w - z)650    return m, out651652def cmd_lifts(args):653    print("Occ_T = {A / m_T : A submask of K_T = m_T R_k, m_T | A, 0 < A < K_T, gcd(A / m_T, R_k) = 1}, T inside [1, k-1], m_T = 1 + 2 a_T; "654          "U = |union of Occ_T|, sum = Sum_T |Occ_T|, model = Sum_T 2^k / m_T, peak = max_T |Occ_T| m_T / 2^k at multiplier mstar; "655          "Phi = submask floor, agg = sum 2^k / (Phi model) is the aggregate against the model after the coprime cut; "656          "cyc = |Occ_T| at T = [t, 2t-1] when k = 2t+1, pred its cyclotomic lower bound; deep = Z(R_k) - U where the band automaton reaches")657    print("k R_k Phi U sum L U/2^k U/Phi U/sum agg peak mstar cyc pred Z deep secs")658    band = []659    for k in range(2, args.kmax + 1):660        t = time.time()661        w = (3 ** k - 1) // 2662        phi = submask_floor(k)663        occ, total, model, best = lift_scan(k)664        if k % 2 and k >= 5:665            j = k // 2666            mstar, low = cyclotomic_occ(j)667            cyc = len(occ_of(k, ((1 << j) - 1) << j)[1])668            pred = len(low)669        else:670            cyc = pred = mstar = "-"671        if k <= args.zmax:672            z = 2 * sum(1 for z1 in range(3, w, 3) if gcd(z1, w) == 1 and witness(z1, w - z1))673            deep = z - len(occ)674        else:675            z = deep = "-"676        agg = total * 2 ** k / (phi * model)677        band.append((len(occ) / phi, agg, model / 2 ** k, best[0], len(occ) / total))678        print(k, w, phi, len(occ), total, down(model / 2 ** k, 5), down(len(occ) / 2 ** k, 4),679              down(len(occ) / phi, 4), down(len(occ) / total, 4), down(agg, 4),680              down(best[0], 3), best[1], cyc, pred, z, deep, round(time.time() - t, 1))681    tail = band[args.lo - 2:]682    names = ("U/Phi", "agg", "L", "peak", "U/sum")683    print("bands over k =", args.lo, "..", args.kmax, "".join(684        " " + n + " [" + str(down(min(r[i] for r in tail), 5)) + ", " + str(up(max(r[i] for r in tail), 5)) + "]"685        for i, n in enumerate(names)))686687def raw_scan(k):688    h = max(k // 2, 1)689    v1 = {}690    for tl in range(0, 1 << h, 2):691        v1[tl] = subset_sums([3 ** (i + k) if tl >> i & 1 else 3 ** i for i in range(h)])692    v2 = {}693    for th in range(1 << (k - h)):694        v2[th] = subset_sums([3 ** (i + k) if th >> (i - h) & 1 else 3 ** i for i in range(h, k)])695    ms = []696    ns = []697    tops = []698    for th in range(1 << (k - h)):699        ah = sum(3 ** i for i in range(h, k) if th >> (i - h) & 1)700        V2 = v2[th]701        for tl in range(0, 1 << h, 2):702            al = sum(3 ** i for i in range(h) if tl >> i & 1)703            T = tl | (th << h)704            m = 1 + 2 * (al + ah)705            if m == 1:706                n = 1 << k707            else:708                r2s = np.sort(V2 % m)709                key = (m - v1[tl] % m) % m710                n = int((np.searchsorted(r2s, key, side="right") - np.searchsorted(r2s, key, side="left")).sum())711            ms.append(m)712            ns.append(n)713            tops.append(T.bit_length() - 1)714    return np.array(ms, dtype=np.int64), np.array(ns, dtype=np.int64), np.array(tops)715716def residue_dist(k, T):717    m = 1 + 2 * sum(3 ** i for i in range(k) if T >> i & 1)718    v = np.zeros(m, dtype=np.int64)719    v[0] = 1720    for i in range(k):721        p = i + k if T >> i & 1 else i722        v = v + np.roll(v, pow(3, p, m))723    return m, v724725def fourier_direct(k, T, m):726    u = np.arange(m)727    F = np.ones(m, dtype=complex)728    for i in range(k):729        p = i + k if T >> i & 1 else i730        F *= 1 + np.exp(2j * np.pi * (u * pow(3, p, m) % m) / m)731    return F732733def phi2p_T(p, t):734    k = (p - 1) * t + 1735    T = 0736    for i in range(1, p - 1, 2):737        for j in range(t):738            T |= 1 << (t * i + j)739    return k, T740741def cyclotomic_poly(n, x):742    num = 1743    den = 1744    for d in range(1, n + 1):745        if n % d == 0:746            mu = mobius(n // d)747            if mu == 1:748                num *= x ** d - 1749            elif mu == -1:750                den *= x ** d - 1751    return num // den752753def mobius(n):754    ps = factor(n)755    m = n756    for p in ps:757        if m % (p * p) == 0:758            return 0759    return (-1) ** len(ps)760761def raw_count(k, T):762    m = 1 + 2 * sum(3 ** i for i in range(k) if T >> i & 1)763    A = subset_sums([3 ** (i + k) if T >> i & 1 else 3 ** i for i in range(k)])764    return m, int((A % m == 0).sum())765766def no_one_digits(u, m, lo, hi):767    ok = np.ones(u.size, dtype=bool)768    for j in range(lo, hi + 1):769        ok &= (u * 3 ** j // m) % 3 != 1770    return ok771772def cmd_agg(args):773    print("N_T = #{A submask of K_T : m_T | A} with A = 0 and K_T kept, T over every subset of [1, k-1] with {} included; "774          "F_T(u) = Prod_{p in supp K_T} (1 + e(u 3^p / m_T)); identity N_T = (1/m_T) Sum_u F_T(u), the u = 0 term 2^k / m_T; "775          "M = Sum_T N_T, L = Sum_T 1/m_T, share = (M - 2^k L) / 2^k is the aggregate u != 0 part, "776          "agg' = Sum_T (N_T - 2) / (2^k L) the cut-free aggregate against the model, "777          "Abs = Sum_T (1/m_T) Sum_{u != 0} |F_T(u)| / 2^k, top = max_T max_{u != 0} F_T(u) / 2^k at T*, "778          "iderr = max over T and u of |F_T from the product - F_T from the FFT|, imag = max |Im F_T|")779    print("k M M/2^k L share agg' Abs Abs/prev top T* iderr imag secs")780    prev = None781    for k in range(args.ulo, args.umax + 1):782        t0 = time.time()783        M = 0784        L = 0.0785        absum = 0.0786        top = (0.0, 0)787        iderr = 0.0788        imag = 0.0789        for T in range(0, 1 << k, 2):790            m, v = residue_dist(k, T)791            F = np.fft.fft(v)792            imag = max(imag, float(np.abs(F.imag).max()))793            Fr = F.real794            n = int(v[0])795            if abs(Fr.mean() - n) > 1e-6 * max(n, 1):796                raise SystemExit("identity fails at k %d T %d" % (k, T))797            if k <= args.directmax:798                D = fourier_direct(k, T, m)799                iderr = max(iderr, float(np.abs(D - np.conj(F)).max()))800            M += n801            L += 1 / m802            if m > 1:803                absum += (np.abs(Fr).sum() - 2 ** k) / m804                mx = float(Fr[1:].max()) / 2 ** k805                if mx > top[0]:806                    top = (mx, T)807        share = (M - 2 ** k * L) / 2 ** k808        aggp = (M - 2 * (1 << (k - 1))) / (2 ** k * L)809        absk = absum / 2 ** k810        Tset = [i for i in range(k) if top[1] >> i & 1]811        print(k, M, down(M / 2 ** k, 5), down(L, 5), down(share, 5), down(aggp, 5), down(absk, 4),812              "-" if prev is None else down(absk / prev, 4), down(top[0], 5), Tset,813              "%.1e" % iderr if k <= args.directmax else "-", "%.1e" % imag, round(time.time() - t0, 1))814        prev = absk815    print("mass at the cyclotomic T = [t, 2t-1], k = 2t+1, m = Phi_6(3^t): N_T against 2^t, F_T(1)/2^k against its bound 1 - 13 9^(-t) for t >= 2, "816          "share_no1 = the share of Sum_{u != 0} F_T(u) carried by the u whose digits d_2..d_t of u/m avoid 1, frac_no1 = how many such u per m")817    print("t k m N_T 2^t F(1)/2^k bound share_no1 frac_no1")818    for t in range(2, args.cycmax + 1):819        k = 2 * t + 1820        T = ((1 << t) - 1) << t821        m, v = residue_dist(k, T)822        F = np.fft.fft(v).real823        u = np.arange(m)824        ok = no_one_digits(u, m, 2, t)825        ok[0] = False826        rest = F[1:].sum()827        print(t, k, m, int(v[0]), 2 ** t, down(F[1] / 2 ** k, 6), down(1 - 13 * 9.0 ** (-t), 6),828              down(F[ok].sum() / rest, 4), down(ok.sum() / m, 5))829    print("the antipodal family: p odd, t >= 1, 0 <= s <= t, k = (p-1) t + s, T = Union_{i odd <= p-2} [ti, ti+t-1], "830          "m_T = (3^(pt) + 1) / (3^t + 1), N_T = 2^((p-1)(t-s)/2 + s) exactly, 2^t at s = t; the verb raises on any mismatch; "831          "printed per (p, t): the k range, m_T, whether m_T is that quotient, and N_T against pred at every s")832    print("p t k_lo..k_hi m_T quotient N_T(s=0..t) pred(s=0..t) equal")833    for p in range(3, 20, 2):834        for t in range(1, 20):835            if (p - 1) * t > args.kmax:836                break837            T = phi2p_T(p, t)[1]838            m = 1 + 2 * sum(3 ** i for i in range(64) if T >> i & 1)839            ns = []840            preds = []841            ok = True842            for s in range(0, t + 1):843                k = (p - 1) * t + s844                if k > args.kmax:845                    break846                n = raw_count(k, T)[1]847                pred = 2 ** ((p - 1) * (t - s) // 2 + s)848                ns.append(n)849                preds.append(pred)850                ok &= n == pred851            print(p, t, "%d..%d" % ((p - 1) * t, (p - 1) * t + len(ns) - 1), m, m == (3 ** (p * t) + 1) // (3 ** t + 1), ns, preds, ok)852            if not ok or m != (3 ** (p * t) + 1) // (3 ** t + 1):853                raise SystemExit("antipodal family fails at p %d t %d" % (p, t))854    print("the cut-free aggregate to k = kmax by meet in the middle: M, M/2^k, L, share, agg', "855          "and the excess by top digit t = max T, X_t = Sum_{max T = t} (N_T - 2 - 2^k / m_T) / 2^k at its argmax")856    print("k M M/2^k L share agg' t* X_t* secs")857    band = []858    for k in range(2, args.kmax + 1):859        t0 = time.time()860        ms, ns, tops = raw_scan(k)861        M = int(ns.sum())862        L = float((1.0 / ms).sum())863        share = (M - 2 ** k * L) / 2 ** k864        aggp = (M - 2 * (1 << (k - 1))) / (2 ** k * L)865        ex = (ns - 2 - 2.0 ** k / ms) / 2 ** k866        ex[ms == 1] = 0867        xt = np.array([ex[tops == t].sum() for t in range(k)])868        ts = int(xt.argmax())869        band.append((M / 2 ** k, share, aggp))870        print(k, M, down(M / 2 ** k, 5), down(L, 5), down(share, 5), down(aggp, 5), ts, down(xt[ts], 5), round(time.time() - t0, 1))871    tail = band[args.lo - 2:]872    names = ("M/2^k", "share", "agg'")873    print("bands over k =", args.lo, "..", args.kmax, "".join(874        " " + n + " [" + str(down(min(r[i] for r in tail), 5)) + ", " + str(up(max(r[i] for r in tail), 5)) + "]"875        for i, n in enumerate(names)))876877# DEEP TAIL878879def col_sums(opts):880    out = np.zeros(1, dtype=np.int64)881    for o in opts:882        out = np.concatenate([out + v for v in o])883    return out884885def block_half(k, cols, slots, key):886    opts = []887    g = 0888    for i, r in enumerate(cols):889        e = key // 3 ** i % 3890        if slots == 1:891            opts.append([0, 3 ** (r + k * e)])892        else:893            ps = [3 ** (r + k * j) for j in range(3) if j != e]894            opts.append([0, ps[0], ps[1], ps[0] + ps[1]])895        g += 3 ** r * (0 if e == 0 else 1 if e == 1 else 3 ** k + 1)896    return col_sums(opts), g897898def mod_match(V1, V2, m, w):899    r2 = V2 % m900    order = np.argsort(r2, kind="stable")901    r2s = r2[order]902    key = (m - V1 % m) % m903    lo = np.searchsorted(r2s, key, side="left")904    hi = np.searchsorted(r2s, key, side="right")905    cnt = hi - lo906    nz = np.nonzero(cnt)[0]907    if nz.size == 0:908        return np.zeros(0, dtype=np.int64)909    reps = cnt[nz]910    left = np.repeat(nz, reps)911    base = np.repeat(np.cumsum(reps) - reps, reps)912    pos = np.repeat(lo[nz], reps) + np.arange(reps.sum()) - base913    z = (V1[left] + V2[order[pos]]) // m914    return z[(z > 0) & (z < w) & (np.gcd(z, w) == 1)]915916def block3_scan(k, slots):917    if 3 * k * log(3) / log(2) > 62:918        raise ValueError("the depth-3 census is int64 bound to k <= 13")919    w = (3 ** k - 1) // 2920    X = 3 ** k921    h = max(k // 2, 1)922    low = [block_half(k, list(range(h)), slots, key) for key in range(3 ** h)]923    occ = set()924    for kh in range(3 ** (k - h)):925        V2, g2 = block_half(k, list(range(h, k)), slots, kh)926        for kl, (V1, g1) in enumerate(low):927            if (kl % 3 == 0) != (slots == 1):928                continue929            g = g1 + g2930            m = 1 + 2 * g if slots == 1 else X * (X + 1) - 2 * g931            occ.update(mod_match(V1, V2, m, w).tolist())932    return occ933934def column_vectors(k, b):935    n = b * k936    S = np.arange(1, 1 << n)937    bits = ((S[:, None] >> np.arange(n)) & 1).astype(np.int64)938    val = bits.dot(np.array([3 ** i for i in range(n)], dtype=np.int64))939    keep = val % ((3 ** k - 1) // 2) == 0940    return bits[keep].reshape(-1, b, k).sum(axis=1)941942def lift_lengths(k):943    w = (3 ** k - 1) // 2944    dep = {}945    for z1 in range(3, w, 3):946        if gcd(z1, w) != 1:947            continue948        m = witness(z1, w - z1)949        if m:950            d = len(base3(m * w))951            dep[z1] = d952            dep[w - z1] = d953    return dep954955def cmd_tail(args):956    print("d(z) is the base-3 length of m(z) R_k for the minimal witness m(z) of an occupied direction of weight R_k, "957          "and the depth b(z) = ceil(d(z) / k) is the number of k-blocks that lift fills; "958          "U = |Union_T Occ_T| is depth at most 2, V = #{b(z) <= 3} the depth-3 census, Z = Z(R_k), tail = Z - U the deep tail; "959          "cap = V - U is the tail captured at depth 3 and share = cap / tail; "960          "fam = |U union the depth-3 lift families|, new1 and new2 what the all-1 and the all-2 family add beyond U; "961          "the depth histogram of the tail is printed under each row")962    print("k R_k U V Z tail cap share fam new1 new2 maxdepth secs")963    rows = []964    for k in range(2, args.kmax + 1):965        t = time.time()966        w = (3 ** k - 1) // 2967        occ, _, _, _ = lift_scan(k)968        dep = lift_lengths(k)969        deep = sorted(set(dep) - occ)970        hist = {}971        for z in deep:972            b = -(-dep[z] // k)973            hist[b] = hist.get(b, 0) + 1974        V = len(occ) + sum(1 for z in deep if dep[z] <= 3 * k)975        fam = new1 = new2 = "-"976        o1 = set()977        if k <= min(max(args.onemax, args.famax), 13):978            o1 = block3_scan(k, 1)979            new1 = len(o1 - occ)980        if k <= min(args.famax, 13):981            o2 = block3_scan(k, 2)982            new2 = len(o2 - occ)983            fam = len(occ | o1 | o2)984            if fam != V:985                raise ValueError("depth-3 family census " + str(fam) + " against the automaton " + str(V) + " at k = " + str(k))986        share = down((V - len(occ)) / len(deep), 4) if deep else "-"987        rows.append((k, len(deep), V - len(occ)))988        print(k, w, len(occ), V, len(dep), len(deep), V - len(occ), share, fam, new1, new2,989              max(hist) if hist else 0, round(time.time() - t, 1))990        if hist:991            print("  depth histogram of the tail at k =", k, sorted(hist.items()))992    live = [r for r in rows if r[1]]993    print("tail", [r[1] for r in live], "captured", [r[2] for r in live], "at k =", [r[0] for r in live])994    print("tail growth over two steps", [down(live[i + 2][1] / live[i][1], 4) for i in range(len(live) - 2)],995          "against 4 for 2^k, captured growth", [down(live[i + 2][2] / live[i][2], 4) for i in range(len(live) - 2) if live[i][2]])996997# CHECKS998999def cmd_check(args):1000    cap = 601001    for n in (6, 9, 12):1002        b = brute(n, cap)1003        miss = [z for z in b if level(*(z if z[0] % 3 == 0 else z[::-1]))[0] == 0]1004        levbad = [z for z in b1005                  if level(*(z if z[0] % 3 == 0 else z[::-1]))[0] != b[z]]1006        print("brute n", n, "rays", len(b), "not occupied by automaton", len(miss),1007              "level mismatch", len(levbad))1008    b = brute(12, cap)1009    extra = [(a, c) for a in range(1, cap + 1) for c in range(1, cap + 1)1010             if gcd(a, c) == 1 and (a % 3 == 0) != (c % 3 == 0)1011             and level(*((a, c) if a % 3 == 0 else (c, a)))[0]1012             and (a, c) not in b]1013    print("automaton occupied but absent from brute n=12", len(extra))1014    bad = sym = gap = core_bad = 01015    for z1, z2 in pairs(120):1016        d, size = level(z1, z2)1017        o, _ = carry_occupied(z1, z2)1018        if (d > 0) != o:1019            bad += 11020        if reach_size(z1, z2) > band_cap(z1, z2):1021            sym += 11022        if d and z1 + z2 <= 3 * z1 <= 2 * (z1 + z2):1023            gap += 11024        if d and not core(z1, z2)[0]:1025            core_bad += 11026    print("j vs pair-carry disagreements to 120", bad, "band violations", sym,1027          "occupied with z1/w inside [1/3,2/3]", gap,1028          "occupied not reaching the core", core_bad)1029    print("gasket rays with u/(u+v) inside [1/3,2/3] at n=12",1030          sum(1 for a, c in brute(12, 3 ** 12) if a + c <= 3 * a <= 2 * (a + c)))1031    b = brute(12, 60)1032    print("gasket rays symmetric under swap", all((c, a) in b for a, c in b),1033          "every ray has exactly one coordinate divisible by 3",1034          all((a % 3 == 0) != (c % 3 == 0) for a, c in b),1035          "so no direction with 3 dividing neither coordinate is occupied")1036    print("k w submask_floor Z(w) w^(log2/log3) max_lev floor_beats_w^0.6309")1037    for k in range(2, 9):1038        wt = (3 ** k - 1) // 21039        sub = [a for a in binaries(k) if 0 < a < wt and gcd(a, wt) == 1]1040        got = [a for a in range(1, wt) if gcd(a, wt) == 11041               and (a % 3 == 0) != ((wt - a) % 3 == 0)1042               and level(*((a, wt - a) if a % 3 == 0 else (wt - a, a)))[0]]1043        lv = max(level(*((a, wt - a) if a % 3 == 0 else (wt - a, a)))[0] for a in got)1044        print(k, wt, len(sub), len(got), round(wt ** (log(2) / log(3)), 1), lv,1045              len(sub) > wt ** (log(2) / log(3)))1046    R = {}1047    for k in range(2, 7):1048        s = set()1049        for lab in range(3 ** k):1050            u = v = 01051            m = lab1052            for i in range(k):1053                d = m % 31054                m //= 31055                if d == 1:1056                    u += 3 ** i1057                elif d == 2:1058                    v += 3 ** i1059            if u and v % 3:1060                s.add(u * pow(v, -1, 3 ** k) % 3 ** k)1061        R[k] = s1062    Z = [0] * 7291063    recov = miss = 01064    for z1 in range(3, 729, 3):1065        for z2 in range(1, 729):1066            if z2 % 3 == 0 or gcd(z1, z2) != 1 or z1 + z2 > 728:1067                continue1068            if not level(z1, z2)[0]:1069                continue1070            w = z1 + z21071            Z[w] += 21072            k = 11073            while 3 ** k <= w:1074                k += 11075            r = z1 * pow(z2, -1, 3 ** k) % 3 ** k1076            if r * w * pow(1 + r, -1, 3 ** k) % 3 ** k != z1:1077                recov += 11078            if 2 <= k <= 6 and r not in R[k]:1079                miss += 11080    over = 01081    for w in range(4, 729):1082        k = 11083        while 3 ** k <= w:1084            k += 11085        if Z[w] > 2 * len(R[k]):1086            over += 11087    print("|R_k| for k = 2..6", [len(R[k]) for k in range(2, 7)],1088          "recovery failures", recov, "residues outside R_k", miss,1089          "weights breaking Z(w) <= 2|R_k|", over)1090    for k in (7, 8):1091        occ, total, model, best = lift_scan(k)1092        direct = set()1093        dtot = 01094        for T in range(0, 1 << k, 2):1095            m, o = occ_of(k, T)1096            direct |= o1097            dtot += len(o)1098        print("lift meet-in-the-middle against direct submask enumeration at k", k,1099              "union", len(occ), len(direct), "sum", total, dtot, "agree", occ == direct and total == dtot)1100    for t in range(2, 7):1101        k = 2 * t + 11102        m, low = cyclotomic_occ(t)1103        mt, full = occ_of(k, ((1 << t) - 1) << t)1104        print("cyclotomic lift t", t, "k", k, "m_T", m, "= 3^(2t)-3^t+1", m == mt,1105              "divides 3^(3t)+1", (3 ** (3 * t) + 1) % m == 0, "cofactor 3^t+1 coprime to R_k",1106              gcd(3 ** t + 1, (3 ** k - 1) // 2) == 1, "predicted", len(low), "actual", len(full),1107              "equal", low == full)1108    for k in (7, 8):1109        ms, ns, tops = raw_scan(k)1110        bad = 01111        for j, T in enumerate(range(0, 1 << k, 2)):1112            m, v = residue_dist(k, T)1113            F = fourier_direct(k, T, m)1114            if m != ms[j] or int(v[0]) != ns[j] or abs(F.sum() / m - ns[j]) > 1e-6:1115                bad += 11116        fam = all(raw_count((p - 1) * t + s, phi2p_T(p, t)[1])[1] == 2 ** ((p - 1) * (t - s) // 2 + s)1117                  for p, t in ((3, 2), (3, 3), (5, 2), (7, 1), (9, 1)) for s in range(t) if (p - 1) * t + s <= 8)1118        print("Fourier identity N_T = (1/m_T) Sum_u Prod_p (1 + e(u 3^p / m_T)) against the residue DP and the meet-in-the-middle count at k",1119              k, "over", len(ms), "sets T, failures", bad, "antipodal family N_T = 2^((p-1)(t-s)/2 + s) at k <= 8", fam)1120    rigid = sum(1 for k in range(2, 6) for c in column_vectors(k, 3) if len(set(c.tolist())) > 1)1121    branch = [c.tolist() for c in column_vectors(3, 4) if len(set(c.tolist())) > 1]1122    cen = []1123    for k in (7, 8, 9):1124        fam = len(lift_scan(k)[0] | block3_scan(k, 1) | block3_scan(k, 2))1125        Z = len(lift_lengths(k))1126        assert fam == Z, (k, fam, Z)1127        cen.append((k, fam))1128    print("binary multiples of R_k below 3^(3k) carry a constant column vector at k = 2..5, exceptions", rigid,1129          "and at four blocks the rigidity breaks, non-constant vectors at k = 3", len(branch), "first", branch[0],1130          "; the depth-3 family census equals Z(R_k) at", cen)1131    worst = max(sum(1 / (1 + 2 * sum(3 ** i for i in range(k) if T >> i & 1))1132                    for T in range(0, 1 << k, 2)) for k in range(2, 13))1133    print("Sum_T 1/m_T over T inside [1, k-1] at k = 2..12 stays below 3/2, largest", down(worst, 5))1134    print("regression A(inf, 3^5) >= 474 and A(9, 3^7) = 2818 are checked by levels")11351136def main():1137    p = argparse.ArgumentParser()1138    s = p.add_subparsers(dest="cmd", required=True)1139    a = s.add_parser("sweep")1140    a.add_argument("cap", type=int)1141    a.set_defaults(fn=cmd_sweep)1142    b = s.add_parser("levels")1143    b.add_argument("cap", type=int)1144    b.add_argument("--lo", type=int, default=2)1145    b.add_argument("--hi", type=int, default=24)1146    b.set_defaults(fn=cmd_levels)1147    c = s.add_parser("multiples")1148    c.add_argument("n", type=int)1149    c.add_argument("--qmax", type=int, default=500)1150    c.add_argument("--hmax", type=int, default=8)1151    c.set_defaults(fn=cmd_multiples)1152    e = s.add_parser("band")1153    e.add_argument("cap", type=int)1154    e.set_defaults(fn=cmd_band)1155    f = s.add_parser("layers")1156    f.add_argument("cap", type=int)1157    f.set_defaults(fn=cmd_layers)1158    g = s.add_parser("weights")1159    g.add_argument("w", type=int, nargs="+")1160    g.set_defaults(fn=cmd_weights)1161    h = s.add_parser("core")1162    h.add_argument("w", type=int, nargs="+")1163    h.set_defaults(fn=cmd_core)1164    i = s.add_parser("box")1165    i.add_argument("lo", type=int)1166    i.add_argument("hi", type=int)1167    i.add_argument("--extra", type=int, default=3)1168    i.set_defaults(fn=cmd_box)1169    j = s.add_parser("repunit")1170    j.add_argument("--kmax", type=int, default=13)1171    j.add_argument("--liftmax", type=int, default=15)1172    j.set_defaults(fn=cmd_repunit)1173    l = s.add_parser("lifts")1174    l.add_argument("--kmax", type=int, default=15)1175    l.add_argument("--zmax", type=int, default=13)1176    l.add_argument("--lo", type=int, default=11)1177    l.set_defaults(fn=cmd_lifts)1178    q = s.add_parser("tail")1179    q.add_argument("--kmax", type=int, default=13)1180    q.add_argument("--famax", type=int, default=12)1181    q.add_argument("--onemax", type=int, default=13)1182    q.set_defaults(fn=cmd_tail)1183    n = s.add_parser("agg")1184    n.add_argument("--kmax", type=int, default=19)1185    n.add_argument("--umax", type=int, default=11)1186    n.add_argument("--ulo", type=int, default=5)1187    n.add_argument("--directmax", type=int, default=9)1188    n.add_argument("--cycmax", type=int, default=5)1189    n.add_argument("--lo", type=int, default=11)1190    n.set_defaults(fn=cmd_agg)1191    d = s.add_parser("check")1192    d.set_defaults(fn=cmd_check)1193    args = p.parse_args()1194    args.fn(args)11951196main()