smith_cascade.py

1import sys2import time3from fractions import Fraction4from math import comb56PREC = 2567TOP = 51189def jacobsthal(k):10    return (2 ** k - (-1) ** k) // 31112def polymul(a, b):13    out = [0] * (len(a) + len(b) - 1)14    for i, x in enumerate(a):15        if x:16            for j, y in enumerate(b):17                out[i + j] += x * y18    return out1920def digit_poly(D):21    base = [0] * (2 * D - 1)22    for k in range(D):23        base[2 * k] = comb(D - 1, k)24    return polymul(base, [1, D, 1])2526def m_even(D):27    P = digit_poly(D)28    R = (D - 1) // 229    n = R + 13031    def coefficient(x):32        return P[x] if 0 <= x < len(P) else 03334    M = [[0] * n for _ in range(n)]35    for cp in range(n):36        for c in range(n):37            v = coefficient(c + D - 3 * cp)38            if c:39                v += coefficient(-c + D - 3 * cp)40            M[cp][c] = v41    return M4243def fill(D):44    return 2 ** (D - 1) * (D + 2)4546def bareiss_det(A):47    n = len(A)48    M = [row[:] for row in A]49    sign = 150    prev = 151    for k in range(n - 1):52        if M[k][k] == 0:53            p = next((r for r in range(k + 1, n) if M[r][k]), None)54            if p is None:55                return 056            M[k], M[p] = M[p], M[k]57            sign = -sign58        for i in range(k + 1, n):59            for j in range(k + 1, n):60                M[i][j] = (M[i][j] * M[k][k] - M[i][k] * M[k][j]) // prev61        prev = M[k][k]62    return sign * M[n - 1][n - 1]6364def v_2(x):65    return (x & -x).bit_length() - 16667def smith_valuations(A, prec=PREC):68    n = len(A)69    cur = prec70    M = [[x % (1 << cur) for x in row] for row in A]71    out = []72    for k in range(n):73        bv = None74        best = None75        for i in range(k, n):76            for j in range(k, n):77                x = M[i][j]78                if x:79                    v = v_2(x)80                    if bv is None or v < bv:81                        bv, best = v, (i, j)82            if bv == 0:83                break84        if best is None:85            raise ValueError("singular modulo 2^%d at step %d" % (cur, k))86        if bv > cur - 64:87            raise ValueError("precision exhausted")88        i0, j0 = best89        M[k], M[i0] = M[i0], M[k]90        if j0 != k:91            for row in M:92                row[k], row[j0] = row[j0], row[k]93        out.append(bv)94        cur -= bv95        m2 = 1 << cur96        u = (M[k][k] >> bv) % m297        inv = pow(u, -1, m2)98        prow = M[k][k:]99        for i in range(k + 1, n):100            Mi = M[i]101            a = Mi[k] >> bv102            if a % m2:103                f = (a * inv) % m2104                M[i] = Mi[:k] + [(x - f * y) % m2 for x, y in zip(Mi[k:], prow)]105            elif bv:106                M[i] = [x % m2 for x in Mi]107            M[i][k] = 0108    return out109110TROUGHS = sorted({2 * jacobsthal(k) + s for k in range(2, 20) for s in (1, 3)})111112def tent(D):113    return min(abs(D - t) // 2 + 1 for t in TROUGHS)114115def layer(a, j):116    return sum(1 for x in a if x >= j)117118def profile(D, prec=PREC):119    a = smith_valuations(m_even(D), prec)120    n = len(a)121    return dict(D=D, n=n, a=a, v_2=sum(a), L1=layer(a, 1), L2=layer(a, 2), L3=layer(a, 3), L4=layer(a, 4), L5=layer(a, 5), amax=max(a), X=sum(a) - layer(a, 1), tent=tent(D))122123def octave(D):124    return D.bit_length() - 1125126def cone_checks(rows, key):127    seq = [r[key] for r in rows]128    lipschitz = sum(1 for i in range(len(seq) - 1) if abs(seq[i + 1] - seq[i]) > 1)129    minima = [rows[i]["D"] for i in range(1, len(seq) - 1) if seq[i] <= seq[i - 1] and seq[i] <= seq[i + 1] and (seq[i] < seq[i - 1] or seq[i] < seq[i + 1]) and seq[i] != 1]130    plateaus = [rows[i]["D"] for i in range(1, len(seq) - 1) if seq[i - 1] == seq[i] == seq[i + 1] and seq[i] != 1]131    return lipschitz, minima, plateaus132133def argmax(rows, key):134    m = max(r[key] for r in rows)135    return m, [r["D"] for r in rows if r[key] == m]136137def main():138    t0 = time.time()139    rows = [profile(D) for D in range(3, TOP + 1, 2)]140    print("odd D = 3..%d, rows %d, seconds %.0f" % (TOP, len(rows), time.time() - t0))141    print("tent law nullity = tent(D): %d/%d, misses %s" % (sum(r["L1"] == r["tent"] for r in rows), len(rows), [r["D"] for r in rows if r["L1"] != r["tent"]]))142    rows = [r for r in rows if r["D"] >= 5]143    print("rows at odd D = 5..%d: %d" % (TOP, len(rows)))144    for D in (5, 7):145        r = next(r for r in rows if r["D"] == D)146        print("D = %d profile %s v_2 = %d n = %d" % (D, r["a"], r["v_2"], r["n"]))147    print()148    print("octave  D range   maxL2 J(k-2)  maxL3 J(k-4)  maxL4  maxL5  max_amax k+4  maxX  max(v_2-ceil(n/3))  at D")149    for k in range(2, 9):150        W = [r for r in rows if octave(r["D"]) == k]151        l2, _ = argmax(W, "L2")152        l3, _ = argmax(W, "L3")153        l4, _ = argmax(W, "L4")154        l5, _ = argmax(W, "L5")155        am, _ = argmax(W, "amax")156        x, _ = argmax(W, "X")157        for r in W:158            r["slack"] = r["v_2"] - (-(-r["n"] // 3))159        sl, at = argmax(W, "slack")160        print("%-7d %3d..%-4d %5d %6d  %5d %6d  %5d  %5d  %8d %3d  %4d  %17d  %s" % (k, W[0]["D"], W[-1]["D"], l2, jacobsthal(k - 2), l3, jacobsthal(k - 4) if k >= 4 else 0, l4, l5, am, k + 4, x, sl, at))161    print()162    for j, key in ((1, "L1"), (2, "L2"), (3, "L3")):163        lip, mins, plat = cone_checks(rows, key)164        print("layer %d cones: Lipschitz breaks %d, interior local minima off 1 %s, plateaus off 1 %s" % (j, lip, mins, plat))165    print()166    print("v_2 <= ceil(n/3) + 9 at every row: %s, max slack %d at D = %s" % (all(r["slack"] <= 9 for r in rows), *argmax(rows, "slack")))167    print("v_2 > n at D = %s" % [r["D"] for r in rows if r["v_2"] > r["n"]])168    print("v_2 = n at D = %s" % [r["D"] for r in rows if r["v_2"] == r["n"]])169    cls = [r for r in rows if r["D"] % 6 == 1 and r["D"] >= 13]170    ratio = max(Fraction(r["v_2"], r["D"] - 1) for r in cls)171    print("class D = 1 mod 6, rows %d (%d..%d): max(v_2 - n) = %d, all v_2 < D - 1: %s, max v_2/(D-1) = %s at D = %s" % (len(cls), cls[0]["D"], cls[-1]["D"], max(r["v_2"] - r["n"] for r in cls), all(r["v_2"] < r["D"] - 1 for r in cls), ratio, [r["D"] for r in cls if Fraction(r["v_2"], r["D"] - 1) == ratio]))172    last = rows[-1]173    print("D = %d: v_2 = %d, ceil(n/3) + 13 = %d, ceil(n/3) + 9 = %d" % (last["D"], last["v_2"], -(-last["n"] // 3) + 13, -(-last["n"] // 3) + 9))174    print("max v_2/n over odd D >= 17: %.3f" % max(r["v_2"] / r["n"] for r in rows if r["D"] >= 17))175    print()176    print("first amax > 9 at D = %d, max amax %d at D = %s" % (next(r["D"] for r in rows if r["amax"] > 9), *argmax(rows, "amax")))177    print("first L3 != 1 at D = %d, max L3 %d at D = %s" % (next(r["D"] for r in rows if r["L3"] != 1), *argmax(rows, "L3")))178    print("first L2 > 5 at D = %d, max L2 %d at D = %s" % (next(r["D"] for r in rows if r["L2"] > 5), *argmax(rows, "L2")))179    print("max X %d at D = %s" % argmax(rows, "X"))180    print("X = v_2 - nullity at D = 255, 257: %s" % [r["X"] for r in rows if r["D"] in (255, 257)])181    print("second largest divisor valuation, max over rows %d; rows with L5 = 0: %s" % (max(sorted(r["a"])[-2] for r in rows), [r["D"] for r in rows if r["L5"] == 0]))182    print("v_2/(D-1) at D = 13, 19: %s" % [str(Fraction(r["v_2"], r["D"] - 1)) for r in rows if r["D"] in (13, 19)])183    print()184    bad = [D for D in range(5, 62, 2) if v_2(bareiss_det(m_even(D))) != sum(smith_valuations(m_even(D)))]185    print("exact determinant against Smith sum, odd D = 5..61: mismatches %s" % bad)186    t1 = time.time()187    stable = all(smith_valuations(m_even(D), 1024) == next(r["a"] for r in rows if r["D"] == D) for D in (255, 257, 511))188    print("profiles at D = 255, 257, 511 unchanged at precision 1024: %s, seconds %.0f" % (stable, time.time() - t1))189    E = m_even(7)190    pencil = [[fill(7) * (i == j) - 3 * E[i][j] for j in range(4)] for i in range(4)]191    d7 = bareiss_det(pencil)192    print("D = 7: det(fill I - 3 M_even) = %d, v_2 = %d, fill = %d" % (d7, v_2(d7), fill(7)))193    print("total seconds %.0f" % (time.time() - t0))194195if __name__ == "__main__":196    main()