pascal_shear.py

1import sys2import urllib.error3import urllib.request4from pathlib import Path56HERE = Path(__file__).resolve().parent7GOULD_URL = "https://oeis.org/A001316/b001316.txt"8GASKET_URL = "https://oeis.org/A047999/b047999.txt"9GASKET_CACHE = "a047999.txt"10CODE = 711TOP_LEVEL = 912DIGIT_LEVEL = 1413EXACT = 12814ROWS = 102415SHEAR_LEVEL = 616TRUNCATION_LEVEL = 8171819def popcount(value):20    return bin(value).count("1")212223def tile_of(code):24    grid = [[0, 0], [0, 0]]25    for index, (x, y) in enumerate([(0, 0), (0, 1), (1, 0), (1, 1)]):26        grid[x][y] = (code >> index) & 127    return grid282930def kron(a, b):31    ra, ca = len(a), len(a[0])32    rb, cb = len(b), len(b[0])33    out = [[0] * (ca * cb) for _ in range(ra * rb)]34    for i in range(ra):35        for j in range(ca):36            if a[i][j]:37                for u in range(rb):38                    for v in range(cb):39                        out[i * rb + u][j * cb + v] = b[u][v]40    return out414243def fractal(tile, level):44    out = [row[:] for row in tile]45    for _ in range(1, level):46        out = kron(out, tile)47    return out484950def pascal_mod2(rows):51    out = [[1]]52    for n in range(1, rows):53        previous = out[n - 1]54        row = [1]55        for k in range(1, n):56            row.append((previous[k - 1] + previous[k]) & 1)57        row.append(1)58        out.append(row)59    return out606162def carry_count(i, j):63    count = 064    carry = 065    while i or j or carry:66        total = (i & 1) + (j & 1) + carry67        carry = 1 if total >= 2 else 068        count += carry69        i >>= 170        j >>= 171    return count727374def binomial(n, k):75    k = min(k, n - k)76    out = 177    for i in range(k):78        out = out * (n - i) // (i + 1)79    return out808182def valuation2(value):83    return (value & -value).bit_length() - 1848586def parse_terms(text):87    terms = []88    for line in text.splitlines():89        line = line.strip()90        if not line or line.startswith("#"):91            continue92        parts = line.split()93        terms.append((int(parts[0]), int(parts[1])))94    return terms959697def read_cache(name):98    return parse_terms((HERE / name).read_text())99100101def read_live(url):102    request = urllib.request.Request(url, headers={"User-Agent": "curl/8"})103    try:104        with urllib.request.urlopen(request, timeout=30) as handle:105            return parse_terms(handle.read().decode())106    except (urllib.error.URLError, OSError, TimeoutError):107        return None108109110def triangular_rows(count):111    rows = 0112    seen = 0113    while seen < count:114        rows += 1115        seen += rows116    return rows117118119def level_sets():120    print("KRONECKER LEVEL SET OF CODE 7")121    tile = tile_of(CODE)122    print(f"  tile [[{tile[0][0]}, {tile[0][1]}], [{tile[1][0]}, {tile[1][1]}]]")123    for level in range(1, TOP_LEVEL + 1):124        grid = fractal(tile, level)125        side = 1 << level126        cells = {(i, j) for i in range(side) for j in range(side) if grid[i][j]}127        anded = {(i, j) for i in range(side) for j in range(side) if i & j == 0}128        agree = cells == anded and len(cells) == 3**level129        print(130            f"  L={level:2d}  side {side:4d}  cells {len(cells):6d}"131            f"  3^L {3 ** level:6d}  set == (i AND j == 0) {agree}"132        )133    same = all(134        sum(1 << (level - popcount(i)) for i in range(1 << level)) == 3**level135        for level in range(1, DIGIT_LEVEL + 1)136    )137    print(f"  digit sum gives 3^L for L = 1..{DIGIT_LEVEL}: {same}")138139140def kummer():141    print("KUMMER ON EXACT BINOMIALS")142    faults = 0143    for i in range(EXACT):144        for j in range(EXACT):145            value = binomial(i + j, i)146            if valuation2(value) != carry_count(i, j):147                faults += 1148            if (value & 1 == 1) != (i & j == 0):149                faults += 1150    print(f"  0 <= i, j < {EXACT}: {EXACT * EXACT} binomials, {faults} faults")151152153def recurrence(triangle):154    print("PASCAL MOD 2 FROM THE ADDITIVE RECURRENCE")155    faults = 0156    entries = 0157    for n, row in enumerate(triangle):158        for k, value in enumerate(row):159            entries += 1160            if (value == 1) != (k & (n - k) == 0):161                faults += 1162    print(f"  rows 0..{ROWS - 1}: {entries} entries, {faults} mismatched cells")163164165def shear(triangle):166    print("THE SHEAR (i, j) -> (i, i + j)")167    side = 1 << SHEAR_LEVEL168    cells = {(i, j) for i in range(side) for j in range(side) if i & j == 0}169    image = {(i, i + j) for (i, j) in cells}170    odd = {171        (k, n)172        for n, row in enumerate(triangle[: 2 * side - 1])173        for k, value in enumerate(row)174        if value and k < side and n - k < side175    }176    print(177        f"  L={SHEAR_LEVEL}: {len(cells)} cells map onto {len(image)} points,"178        f" odd entries in range {len(odd)}, bijection {image == odd}"179    )180181182def gould(triangle):183    print("ANTIDIAGONAL POPULATION")184    faults = sum(1 for n, row in enumerate(triangle) if sum(row) != 1 << popcount(n))185    print(f"  row sums == 2^popcount(n) for n = 0..{ROWS - 1}, {faults} faults")186    side = 1 << TRUNCATION_LEVEL187    counts = [0] * (2 * side - 1)188    for i in range(side):189        for j in range(side):190            if i & j == 0:191                counts[i + j] += 1192    inside = all(counts[n] == 1 << popcount(n) for n in range(side))193    outside = all(counts[n] < 1 << popcount(n) for n in range(side, 2 * side - 1))194    print(195        f"  inside [0, 2^{TRUNCATION_LEVEL})^2 the count is 2^popcount(n)"196        f" for n < 2^L {inside}, strictly smaller above {outside}"197    )198199200def oeis(triangle):201    print("OEIS")202    flat = [value for row in triangle for value in row]203    sums = [sum(row) for row in triangle]204    terms = read_live(GOULD_URL)205    if terms is None:206        print("  A001316 live read unavailable, its b-file is too large to keep here")207    else:208        wrong = sum(209            1210            for n, (index, value) in enumerate(terms)211            if index != n or value != 1 << popcount(n)212        )213        print(214            f"  A001316 live: {len(terms)} terms, n = 0..{len(terms) - 1},"215            f" {wrong} differences"216        )217        head = [value for _, value in terms[:ROWS]]218        print(f"  A001316 first {ROWS} terms == our row sums: {head == sums}")219    cache = read_cache(GASKET_CACHE)220    gasket = [("cache", cache)]221    live = read_live(GASKET_URL)222    if live is None:223        print("  A047999 live read unavailable, cache only")224    else:225        gasket.append(("live", live))226    for label, terms in gasket:227        wrong = sum(228            1229            for n, (index, value) in enumerate(terms)230            if index != n or value != flat[n]231        )232        rows = triangular_rows(len(terms))233        print(234            f"  A047999 {label}: {len(terms)} terms, rows 0..{rows - 1},"235            f" {wrong} differences"236        )237    if live is not None:238        print(f"  A047999 cache == live term for term: {cache == live}")239240241def main():242    print(f"DOMAIN levels 1..{TOP_LEVEL}, binomials < {EXACT}, Pascal rows < {ROWS}")243    level_sets()244    kummer()245    triangle = pascal_mod2(ROWS)246    recurrence(triangle)247    shear(triangle)248    gould(triangle)249    oeis(triangle)250    return 0251252253if __name__ == "__main__":254    sys.exit(main())