smith_window.py

1import sys23# BITMASK POLYNOMIALS OVER F245def polymul2(a, b):6    r = 07    while b:8        lb = b & -b9        r ^= a * lb10        b ^= lb11    return r1213def pow2mask(m):14    r = 115    s = 116    while m:17        if m & 1:18            r ^= r << s19        s <<= 120        m >>= 121    return r2223def divmod2(a, b):24    db = b.bit_length() - 125    q = 026    while a and a.bit_length() - 1 >= db:27        sh = a.bit_length() - 1 - db28        q ^= 1 << sh29        a ^= b << sh30    return q, a3132def bits(v):33    out = []34    while v:35        lb = v & -v36        out.append(lb.bit_length() - 1)37        v ^= lb38    return out3940def echelon(vecs):41    B = []42    for v in vecs:43        for u in B:44            v = min(v, v ^ u)45        if v:46            B.append(v)47            B.sort(reverse=True)48    return B4950def reduce_by(v, B):51    for u in B:52        v = min(v, v ^ u)53    return v5455def pivots(vecs):56    piv = {}57    for v in vecs:58        while v:59            t = v.bit_length() - 160            if t in piv:61                v ^= piv[t]62            else:63                piv[t] = v64                break65    return piv6667def residue(v, piv, keys):68    for t in keys:69        if (v >> t) & 1:70            v ^= piv[t]71    return v7273def spans(v, piv):74    while v:75        t = v.bit_length() - 176        if t not in piv:77            return False78        v ^= piv[t]79    return True8081def nullspace_right(rows, n):82    piv = {}83    for r in rows:84        for p in sorted(piv, reverse=True):85            if (r >> p) & 1:86                r ^= piv[p]87        if r:88            q = r.bit_length() - 189            for pp in list(piv):90                if (piv[pp] >> q) & 1:91                    piv[pp] ^= r92            piv[q] = r93    out = []94    for f in [c for c in range(n) if c not in piv]:95        v = 1 << f96        for p, r in piv.items():97            if (r >> f) & 1:98                v |= 1 << p99        out.append(v)100    return out101102# JACOBSTHAL AND THE SLOT103104def jac(k):105    return 0 if k < 0 else (2 ** k - (-1) ** k) // 3106107def slot(R):108    b = 0109    while (1 << b) < 3 * R - 1:110        b += 1111    g = abs(2 * R - (1 << (b - 1)) - 1)112    e = 1113    while jac(e) < (g + 1) // 2:114        e += 1115    return b, g, e, b - 1 - e116117def window(R, b):118    A = 1 << b119    i0 = max(2, 4 * R - A)120    hi = (6 * R + 2 - A) // 3121    hi -= hi % 2122    return i0, (hi - i0) // 2123124def fibpoly(t):125    if t == 0:126        return 1127    a, b = 0b11, 1128    for _ in range(t - 1):129        a, b = polymul2(a, 2) ^ b, a130    return a131132def law_e(D):133    R = (D - 1) // 2134    b, gslot, e, k = slot(R)135    i0, K = window(R, b)136    tt = (jac(k) - 1) // 2137    C = K - (2 * jac(e - 1) if k % 2 == 0 else 0)138    w = C - tt * (1 << e)139    ce = 2 * jac(e - 2) - 1 if e >= 3 else 1140    m = max(0, 2 * w - ce)141    g = 0142    for a in bits(fibpoly(tt)):143        g ^= 1 << (a * (1 << e))144    return K, C, g << m, min(w + 1, ce + 1 - w)145146def tent(R):147    b, gslot, e, k = slot(R)148    N = jac(e) - jac(e - 1)149    u = (gslot + 1) // 2 - jac(e - 1) - 1150    p = u if R > (1 << (b - 2)) else N - 1 - u151    return p, min(p, N - 1 - p)152153def reach_law(R):154    b, gslot, e, k = slot(R)155    p, w = tent(R)156    if e == 1 and k % 2:157        return 1 if R == 1 << (b - 2) else (3 if R == 3 else 5)158    return 3 * w + (2 if e % 2 == 0 else 0) + ((1 + p % 2) if k % 2 else 0)159160# THE MOD-4 SYMBOL AND THE KERNEL FAMILY161162FW = 20163FM = (1 << FW) - 1164165def symbol4(D):166    out = [0] * (2 * D + 1)167    c = 1168    for k in range(D):169        v = c % 4170        if v:171            out[2 * k] = (out[2 * k] + v) % 4172            out[2 * k + 1] = (out[2 * k + 1] + v * D) % 4173            out[2 * k + 2] = (out[2 * k + 2] + v) % 4174        c = c * (D - 1 - k) // (k + 1)175    return out176177def symbol4pack(D):178    p = 0179    for i, v in enumerate(symbol4(D)):180        if v:181            p |= v << (FW * i)182    return p183184def family(D):185    R = (D - 1) // 2186    b, gslot, e, k = slot(R)187    i0, K = window(R, b)188    A = 1 << b189    G = polymul2(pow2mask(4 * R), 0b111)190    out = []191    for j in range(K + 1):192        i = i0 + 2 * j193        q = 1194        for a in bits(i):195            q = polymul2(q, 1 ^ (1 << (3 << a)))196        f = polymul2(1 << ((6 * R + 2 - 3 * i - A) // 2), q)197        f ^= f << A198        h, r = divmod2(f, G)199        assert r == 0, "family element off the ideal at D=%d" % D200        out.append(h)201    return out, K202203def colblock(R, G, wide):204    dec = [0, 0, 0]205    for a in bits(G):206        dec[a % 3] |= 1 << (a // 3)207    m = (1 << (R + 1)) - 1208209    def T(a):210        c = 1 - R + a211        r = c % 3212        q = (c - r) // 3213        v = dec[r]214        return ((v >> q) if q >= 0 else (v << -q)) & m215216    out = [T(0)]217    for j in range(1, wide):218        out.append(T(-j) ^ T(j))219    return out220221def columns(D):222    R = (D - 1) // 2223    return colblock(R, polymul2(pow2mask(4 * R), 0b111), R + 1)224225def obstruction(D, X, Ppack):226    R = (D - 1) // 2227    out = []228    for h in X:229        hp = 0230        for a in bits(h):231            hp |= 1 << (FW * a)232        pr = hp * Ppack233        o = 0234        for nu in range(R + 1):235            s = (pr >> (FW * (3 * nu + 1))) & FM236            assert s % 2 == 0, "family element not in the mod-2 kernel at D=%d" % D237            if s % 4 == 2:238                o |= 1 << nu239        out.append(o)240    return out241242# THE LAYER-2 WINDOW243244def layer2(D):245    X, K = family(D)246    cols = columns(D)247    piv = pivots(cols)248    keys = sorted(piv, reverse=True)249    raw = obstruction(D, X, symbol4pack(D))250    red = [residue(o, piv, keys) for o in raw]251    nb = max((o.bit_length() for o in red), default=0)252    rows = []253    for nu in range(nb):254        r = 0255        for j in range(K + 1):256            if (red[j] >> nu) & 1:257                r |= 1 << j258        if r:259            rows.append(r)260    V2 = echelon(nullspace_right(rows, K + 1))261    return K, V2, raw, cols, X, piv262263def shift_law(D, X, K, raw, piv):264    R = (D - 1) // 2265    G = polymul2(pow2mask(4 * R), 0b111)266    box = (1 << (2 * R + 1)) - 1267    for j in range(K):268        H = X[j]269        hi, lo = H << 3, H >> 3270        both = hi & lo271        assert X[j + 1] == (hi ^ lo), "family shift breaks at D=%d j=%d" % (D, j)272        Z = H ^ both273        assert (Z | box) == box, "Z leaves the box at D=%d j=%d" % (D, j)274        assert Z == rev_box(Z, R), "Z not palindromic at D=%d j=%d" % (D, j)275        ZG = polymul2(Z, G)276        az = 0277        for nu in range(R + 1):278            if (ZG >> (3 * nu + 1)) & 1:279                az |= 1 << nu280        assert spans(az, piv), "A(Z) off the image at D=%d j=%d" % (D, j)281        assert raw[j + 1] == fold_shift(raw[j], R) ^ az, "Lambda law breaks at D=%d j=%d" % (D, j)282    return True283284def rev_box(v, R):285    out = 0286    for a in bits(v):287        out |= 1 << (2 * R - a)288    return out289290def fold_shift(f, R):291    out = 0292    for nu in range(R + 1):293        a = f >> (nu - 1) & 1 if nu >= 1 else 0294        b = (f >> (nu + 1)) & 1 if nu + 1 <= R else ((f >> (R - 1)) & 1 if nu == R else 0)295        if a ^ b:296            out |= 1 << nu297    return out298299def corrector_index(cols, target):300    piv = {}301    for J, c in enumerate(cols):302        v = c303        while v:304            t = v.bit_length() - 1305            if t in piv:306                v ^= piv[t]307            else:308                piv[t] = v309                break310        if spans(target, piv):311            return J312    return None313314# THE SWEEP315316def sweep(lo, hi):317    n = bw = bg = bc = bm = br = bt = bf = ne = fr = 0318    cap, flo, tie = [], [], []319    r4 = [0] * 4320    r8 = [0] * 8321    for D in range(lo, hi + 1, 2):322        R = (D - 1) // 2323        b, gs, e, k = slot(R)324        K, V2, raw, cols, X, piv = layer2(D)325        assert V2, "empty layer-2 window at D=%d" % D326        assert X[K].bit_length() - 1 <= 2 * R, "family top leaves the box at D=%d" % D327        gen = min(V2)328        dg = gen.bit_length() - 1329        C = max(v.bit_length() - 1 for v in V2)330        n += 1331        r4[R % 4] += 1332        r8[R % 8] += 1333        if echelon([polymul2(gen, 1 << s) for s in range(C - dg + 1)]) != V2:334            bw += 1335            print("window fails at D=%d" % D)336        Kp, Cp, gp, L2p = law_e(D)337        if (Kp, Cp, gp, L2p) != (K, C, gen, len(V2)):338            print("law E fails at D=%d: %r against %r" % (D, (Kp, Cp, gp, L2p), (K, C, gen, len(V2))))339            if gp != gen:340                bg += 1341            if Cp != C or L2p != len(V2):342                bc += 1343        shift_law(D, X, K, raw, piv)344        om = 0345        for j in bits(gen):346            om ^= raw[j]347        J = corrector_index(cols, om)348        if k % 2 or e == 1:349            assert C == K, "upper half not free at D=%d" % D350            fr += 1351        else:352            assert K - C == 2 * jac(e - 1), "ceiling deficit is not 2J(e-1) at D=%d" % D353        if tent(R)[1] != C - dg:354            bt += 1355            print("tent identity fails at D=%d" % D)356        rl = reach_law(R)357        if rl != R - J:358            br += 1359            print("reach law fails at D=%d: %d against %d" % (D, rl, R - J))360        a, bq = K - dg, rl // 3361        if e == 1 and k % 2 and R > (1 << (b - 2)):362            ne += 1363            if (bq, C - dg, a) != (1, 0, 0):364                bf += 1365                print("escaping row misread at D=%d" % D)366        elif bq != C - dg + (1 if (k % 2 and e % 2 == 0) else 0):367            bf += 1368            print("floor identity fails at D=%d" % D)369        if C - dg != min(a, bq):370            bm += 1371            print("corrector law fails at D=%d" % D)372        elif a < bq:373            cap.append(D)374        elif bq < a:375            flo.append(D)376        else:377            tie.append(D)378    print("odd D = %d..%d: %d rows" % (lo, hi, n))379    print("R mod 4 classes %r, R mod 8 classes %r" % (r4, r8))380    print("window structure V2 = g F_2[z]_(<= C - deg g): %d/%d" % (n - bw, n))381    print("Law E generator z^m c_t(z^(2^e)): %d/%d" % (n - bg, n))382    print("ceiling C = K - 2J(e-1)[k even] and L_2: %d/%d" % (n - bc, n))383    print("the family shift law H^(j+1) = psi H^(j) - 2 Z_j and ob(X_(j+1)) = Lambda ob(X_j) + A(Z_j): %d/%d" % (n, n))384    print("the slot tent w = min(p, N - 1 - p) = C - deg g: %d/%d" % (n - bt, n))385    print("the reach law reach = 3w + 2[e even] + [k odd](1 + p mod 2), D = 4^m + 3 apart: %d/%d" % (n - br, n))386    print("floor(reach/3) = C - deg g + [k odd and e even] off the %d rows D = 4^m + 3, where it reads 1 against C - deg g = K - deg g = 0: %d/%d" % (ne, n - bf, n))387    print("corrector law C - deg g = min(K - deg g, floor(reach/3)) from the closed form: %d/%d" % (n - bm, n))388    print("upper half free where C = K, that is k odd or e = 1, else deficit K - C = 2J(e-1): free %d, open %d" % (fr, n - fr))389    print("branches: floor strict %d, K cap strict %d, tie %d" % (len(flo), len(cap), len(tie)))390    print("floor-strict rows are exactly the C < K rows: %s" % (sorted(flo) == sorted(d for d in range(lo, hi + 1, 2) if law_e(d)[1] < law_e(d)[0])))391392def unboxed(lo, hi):393    n = bad = co = 0394    for D in range(lo, hi + 1, 2):395        R = (D - 1) // 2396        G = polymul2(pow2mask(4 * R), 0b111)397        piv = pivots(colblock(R, G, 4 * R + 12))398        X, K = family(D)399        n += 1400        if len(piv) != R:401            co += 1402        if any(not spans(o, piv) for o in obstruction(D, X, symbol4pack(D))):403            bad += 1404    print("unboxed corrector image has corank exactly 1: %d/%d" % (n - co, n))405    print("every family obstruction meets that image, so the unboxed layer-2 window is all of the mod-2 kernel: %d/%d" % (n - bad, n))406407if __name__ == "__main__":408    a = sys.argv[1:]409    lo = int(a[0]) if a else 5410    hi = int(a[1]) if len(a) > 1 else 1601411    sweep(lo, hi)412    unboxed(5, min(hi, 601))