stack.py

1import hashlib2import time3from math import gcd45FAREY_N = 556Q = 99737M = 2 * Q8PROBES = 489PROBE_NS = [55, 5555, 10 ** 18]10SWEEP_NS = [1, 2, 55, 5555, 19945, 19946, 19947, 40001]11R = 51212D = 2 * R13IMAGE_N = 5514ODDS = list(range(1, IMAGE_N + 1, 2))15K = len(ODDS)1617def farey_nodes(n_max):18    out = []19    for b in range(1, n_max + 1):20        for a in range(1, b + 1):21            if gcd(a, b) == 1:22                out.append((a, b))23    return out2425def farey_by_form(n_max):26    return ["%d/%d:%d" % (a, b, n_max // b) for (a, b) in farey_nodes(n_max)]2728def farey_by_membership(n_max):29    scales = []30    for n in range(1, n_max + 1):31        scales.append(set((k // gcd(k, n), n // gcd(k, n)) for k in range(1, n + 1)))32    lines = []33    for (a, b) in farey_nodes(n_max):34        lines.append("%d/%d:%d" % (a, b, sum(1 for s in scales if (a, b) in s)))35    return lines3637def digest(lines):38    return hashlib.sha256("\n".join(lines).encode("utf-8")).hexdigest()3940def brightness_sum(lines):41    return sum(int(t.split(":")[1]) for t in lines)4243def probe_points():44    return [((137 * i + 61) % Q, (211 * i + 97) % Q) for i in range(PROBES)]4546def bad_residues(a, b):47    return [r for r in range(1, M, 2) if (r * a) % M >= Q and (r * b) % M >= Q]4849def brightness_closed(bad, n_max):50    dark = 051    for r in bad:52        if r <= n_max:53            dark += (n_max - r) // M + 154    return (n_max + 1) // 2 - dark5556def good_residues(a, b):57    return [r for r in range(1, M, 2) if (r * a) % M < Q or (r * b) % M < Q]5859def brightness_direct(good, n_max):60    lit = 061    for r in good:62        if r <= n_max:63            lit += (n_max - r) // M + 164    return lit6566def brightness_literal(a, b, n_max):67    c = 068    for n in range(1, n_max + 1, 2):69        if not (((n * a) // Q) % 2 == 1 and ((n * b) // Q) % 2 == 1):70            c += 171    return c7273def column_masks():74    masks = []75    for j in range(R):76        u = 2 * j + 177        m = 078        for k, n in enumerate(ODDS):79            if ((n * u) // D) % 2 == 1:80                m |= 1 << k81        masks.append(m)82    return masks8384def render_masks():85    masks = column_masks()86    buf = bytearray(R * R)87    for i in range(R):88        my = masks[i]89        row = i * R90        for j in range(R):91            buf[row + j] = K - (my & masks[j]).bit_count()92    return buf9394def render_layers():95    xs = [2 * j + 1 for j in range(R)]96    img = [[0] * R for _ in range(R)]97    for n in ODDS:98        odd_col = [((n * v) // D) % 2 == 1 for v in xs]99        for i in range(R):100            row = img[i]101            if ((n * xs[i]) // D) % 2 == 1:102                for j in range(R):103                    if not odd_col[j]:104                        row[j] += 1105            else:106                for j in range(R):107                    row[j] += 1108    buf = bytearray()109    for i in range(R):110        buf.extend(bytes(img[i]))111    return buf112113def render_pixels():114    buf = bytearray()115    for i in range(R):116        y = 2 * i + 1117        for j in range(R):118            x = 2 * j + 1119            c = 0120            for n in ODDS:121                if not (((n * x) // D) % 2 == 1 and ((n * y) // D) % 2 == 1):122                    c += 1123            buf.append(c)124    return buf125126def main():127    form = farey_by_form(FAREY_N)128    member = farey_by_membership(FAREY_N)129    print("farey N", FAREY_N)130    print("farey nodes", len(form))131    print("farey brightness sum", brightness_sum(form))132    print("farey N(N+1)/2", FAREY_N * (FAREY_N + 1) // 2)133    print("farey digest closed form", digest(form))134    print("farey digest membership", digest(member))135    print("farey digests equal", digest(form) == digest(member))136137    pts = probe_points()138    t0 = time.perf_counter()139    tables = [bad_residues(a, b) for (a, b) in pts]140    t_tables = time.perf_counter() - t0141    print("probe count", len(pts))142    print("probe modulus 2q", M)143    print("probe residue tables seconds", round(t_tables, 4))144145    values = {}146    for n_max in PROBE_NS:147        t0 = time.perf_counter()148        values[n_max] = [brightness_closed(t, n_max) for t in tables]149        dt = time.perf_counter() - t0150        print("closed form N", n_max, "seconds", round(dt, 4))151        print("closed form N", n_max, "first four", *values[n_max][:4])152153    t0 = time.perf_counter()154    twins = [good_residues(a, b) for (a, b) in pts]155    print("direct count residue tables seconds", round(time.perf_counter() - t0, 4))156    for n_max in PROBE_NS:157        t0 = time.perf_counter()158        other = [brightness_direct(t, n_max) for t in twins]159        dt = time.perf_counter() - t0160        print("direct count N", n_max, "seconds", round(dt, 4))161        print("direct count N", n_max, "equals closed form", other == values[n_max])162163    for n_max in (55, 5555):164        agree = sum(165            1166            for i, (a, b) in enumerate(pts)167            if values[n_max][i] == brightness_literal(a, b, n_max)168        )169        print("probes closed equal literal N", n_max, "%d/%d" % (agree, len(pts)))170171    bad = 0172    for i, (a, b) in enumerate(pts):173        for n_max in SWEEP_NS:174            if brightness_closed(tables[i], n_max) != brightness_literal(a, b, n_max):175                bad += 1176    print("sweep N values", *SWEEP_NS)177    print("sweep comparisons", len(pts) * len(SWEEP_NS))178    print("sweep mismatches", bad)179180    a = render_masks()181    b = render_layers()182    c = render_pixels()183    print("image side", R)184    print("image layers", K)185    print("image sha256 masks", hashlib.sha256(bytes(a)).hexdigest())186    print("image sha256 layers", hashlib.sha256(bytes(b)).hexdigest())187    print("image sha256 pixels", hashlib.sha256(bytes(c)).hexdigest())188    print("image three routes equal", a == b == c)189190main()