census.py

1import subprocess2import sys3import time4from fractions import Fraction5from math import ceil, floor, log678def in_gasket(x, y, b):9    while x or y:10        u, v = x % b, y % b11        if u > 1 or v > 1 or (u and v):12            return False13        x //= b14        y //= b15    return True161718def direct(m):19    n = 3 ** m20    c = 021    for x in range(n):22        for y in range(n):23            if in_gasket(x, y, 2) and in_gasket(x, y, 3):24                c += 125    return c262728def count(m):29    pw = [3 ** i for i in range(m + 1)]30    bl = [p.bit_length() for p in pw]31    total = 03233    def walk(k, x, y):34        nonlocal total35        if k == 0:36            if x & y == 0:37                total += 138            return39        b = bl[k]40        hx = x >> b41        hy = y >> b42        if (hx & hy) and (hx & (hy + 1)) and ((hx + 1) & hy) and ((hx + 1) & (hy + 1)):43            return44        i = k - 145        p = pw[i]46        walk(i, x, y)47        walk(i, x + p, y)48        walk(i, x, y + p)4950    walk(m, 0, 0)51    return total525354def axis(m):55    return 2 ** (m + 1) - 1565758def terms(lo, hi):59    print("m  C(3^m)        2^(m+1)-1   C/prev    log_3(C)/m  sec")60    prev = None61    for m in range(lo, hi + 1):62        t0 = time.time()63        c = count(m)64        assert c >= axis(m)65        r = c / prev if prev else float("nan")66        e = log(c, 3) / m if m else float("nan")67        print("%-2d %-13d %-11d %-9.5f %-11.6f %.1f" % (m, c, axis(m), r, e, time.time() - t0))68        sys.stdout.flush()69        prev = c707172def control(hi):73    print("m  direct  pruned  axis  agree")74    for m in range(hi + 1):75        d = direct(m)76        c = count(m)77        print("%-2d %-7d %-7d %-5d %s" % (m, d, c, axis(m), d == c and c >= axis(m)))78        sys.stdout.flush()798081def down(x, k=6):82    return floor(x * 10 ** k) / 10 ** k838485def up(x, k=6):86    return ceil(x * 10 ** k) / 10 ** k878889def budget():90    da = log(3, 2)91    db = 1.092    ax = log(2, 3)93    print("dim A = log_2 3 = %.6f" % da)94    print("dim B = log_3 3 = %.6f" % db)95    print("budget dim A + dim B - 2 <= %.6f" % up(da + db - 2))96    print("axis exponent log_3 2 >= %.6f" % down(ax))97    print("excess >= %.6f" % down(ax - up(da + db - 2)))9899100def series(hi):101    return [count(m) for m in range(hi + 1)]102103104def bareiss(rows):105    n = len(rows)106    a = [row[:] for row in rows]107    sign = 1108    prev = 1109    for k in range(n - 1):110        if a[k][k] == 0:111            p = next((i for i in range(k + 1, n) if a[i][k]), None)112            if p is None:113                return 0114            a[k], a[p] = a[p], a[k]115            sign = -sign116        for i in range(k + 1, n):117            for j in range(k + 1, n):118                a[i][j] = (a[i][j] * a[k][k] - a[i][k] * a[k][j]) // prev119            a[i][k] = 0120        prev = a[k][k]121    return sign * a[n - 1][n - 1]122123124def solve(rows, rhs, r):125    aug = [[Fraction(v) for v in rows[i]] + [Fraction(rhs[i])] for i in range(len(rows))]126    piv = []127    row = 0128    for col in range(r):129        p = next((i for i in range(row, len(aug)) if aug[i][col]), None)130        if p is None:131            continue132        aug[row], aug[p] = aug[p], aug[row]133        d = aug[row][col]134        aug[row] = [v / d for v in aug[row]]135        for i in range(len(aug)):136            if i != row and aug[i][col]:137                f = aug[i][col]138                aug[i] = [aug[i][j] - f * aug[row][j] for j in range(r + 1)]139        piv.append(col)140        row += 1141    for i in range(row, len(aug)):142        if aug[i][r]:143            return "none", None144    c = [Fraction(0)] * r145    for i, col in enumerate(piv):146        c[col] = aug[i][r]147    return ("unique" if len(piv) == r else "free"), c148149150def charpoly_string(c):151    parts = ["x^%d" % len(c)]152    for i, v in enumerate(c):153        parts.append("- (%s)*x^%d" % (v, len(c) - 1 - i))154    return " ".join(parts)155156157def dominant_root(c):158    p = charpoly_string(c)159    out = subprocess.run(["gp", "-q"], input="v=polroots(%s); vecmax(abs(v))\n" % p,160                         capture_output=True, text=True)161    return out.stdout.strip()162163164def hankel(hi, rmax):165    a = series(hi)166    print("terms C(3^m), m = 0..%d" % hi)167    print(",".join(str(v) for v in a))168    print()169    print("k  det H_k (entries a[i+j], i,j < k)")170    for k in range(1, rmax + 2):171        if 2 * k - 1 > len(a):172            break173        print("%-2d %d" % (k, bareiss([[a[i + j] for j in range(k)] for i in range(k)])))174    print()175    print("r  eqs spare fit    coefficients c_1..c_r of a(n) = sum c_i a(n-i), n >= r")176    found = []177    for r in range(1, rmax + 1):178        rows = [[a[n - 1 - j] for j in range(r)] for n in range(r, len(a))]179        rhs = [a[n] for n in range(r, len(a))]180        st, c = solve(rows, rhs, r)181        print("%-2d %-3d %-5d %-6s %s" % (r, len(rows), len(rows) - r, st,182                                          "" if c is None else " ".join(str(v) for v in c)))183        if st != "none":184            found.append((r, c))185    print()186    if not found:187        print("no linear recurrence with constant coefficients of order r <= %d fits all %d terms"188              % (rmax, len(a)))189        print("last reading log_3 C(3^%d) / %d = %.6f" % (hi, hi, log(a[hi], 3) / hi))190        return191    for r, c in found:192        print("order %d characteristic polynomial %s" % (r, charpoly_string(c)))193        rho = dominant_root(c)194        print("order %d dominant root %s" % (r, rho))195        try:196            print("order %d log(root)/log(3) = %.6f" % (r, log(float(rho), 3)))197        except ValueError:198            print("order %d dominant root is not positive real" % r)199        print("order %d last reading log_3 C(3^%d) / %d = %.6f" % (r, hi, hi, log(a[hi], 3) / hi))200201202def first_order(a, rmax):203    for r in range(1, rmax + 1):204        rows = [[a[n - 1 - j] for j in range(r)] for n in range(r, len(a))]205        rhs = [a[n] for n in range(r, len(a))]206        st, c = solve(rows, rhs, r)207        if st != "none":208            return r, c209    return 0, None210211212def selftest():213    ok = True214    print("%-26s %-26s %-26s %s" % ("check", "expect", "got", "pass"))215216    def report(name, want, got):217        nonlocal ok218        hit = str(want) == str(got)219        ok = ok and hit220        print("%-26s %-26s %-26s %s" % (name, want, got, hit))221222    report("bareiss 1x1", 3, bareiss([[3]]))223    report("bareiss 2x2", -2, bareiss([[1, 2], [3, 4]]))224    report("bareiss 2x2 pivot swap", -6, bareiss([[0, 2], [3, 4]]))225    report("bareiss 2x2 singular", 0, bareiss([[1, 2], [2, 4]]))226    report("bareiss 3x3", 4, bareiss([[2, 1, 0], [1, 2, 1], [0, 1, 2]]))227    fib = [0, 1]228    while len(fib) < 20:229        fib.append(fib[-1] + fib[-2])230    r, c = first_order(fib, 6)231    report("fibonacci first order", 2, r)232    report("fibonacci charpoly", "x^2 - (1)*x^1 - (1)*x^0", charpoly_string(c))233    mix = [2 ** n + 3 ** n + 1 for n in range(20)]234    r, c = first_order(mix, 6)235    report("2^n + 3^n + 1 first order", 3, r)236    report("2^n + 3^n + 1 dominant root", "3.0000000", dominant_root(c)[:9])237    cube = [n ** 3 + 2 ** n for n in range(20)]238    r, c = first_order(cube, 8)239    report("n^3 + 2^n first order", 5, r)240    report("n^3 + 2^n charpoly", "x^5 - (6)*x^4 - (-14)*x^3 - (16)*x^2 - (-9)*x^1 - (2)*x^0",241           charpoly_string(c))242    print()243    print("all pass %s" % ok)244245246if __name__ == "__main__":247    v = sys.argv[1] if len(sys.argv) > 1 else "budget"248    if v == "terms":249        terms(int(sys.argv[2]), int(sys.argv[3]))250    elif v == "control":251        control(int(sys.argv[2]))252    elif v == "hankel":253        hankel(int(sys.argv[2]), int(sys.argv[3]) if len(sys.argv) > 3 else 12)254    elif v == "selftest":255        selftest()256    else:257        budget()