toolpath.py

1import sys2import time34U = [(1, 0), (0, 1), (-1, 1), (-1, 0), (0, -1), (1, -1)]5F, R, M, MR = 0, 1, 2, 36REV = [R, F, MR, M]7MIR = [M, MR, F, R]8NAMES = "F R M MR".split()910def norm(c):11    a, b = c12    return a * a + a * b + b * b1314def mul(p, q):15    x, y = p16    s, t = q17    return (x * s - y * t, x * t + y * s + y * t)1819def mirror_index(c):20    c2 = mul(c, c)21    n = norm(c)22    for m in range(6):23        if (U[m][0] * n, U[m][1] * n) == c2:24            return m25    return None2627def turn(d1, d2):28    t = (d2 - d1) % 629    return t - 6 if t > 3 else t3031def walks(n, c, gentle):32    out = []33    ta, tb = c34    def dfs(path, seen, a, b, d):35        k = len(path)36        if k == n:37            if (a, b) == (ta, tb):38                out.append(tuple(path))39            return40        rem = n - k - 141        for nd in range(6):42            if k:43                t = abs(turn(d, nd))44                if t == 3 or (gentle and t > 1):45                    continue46            x, y = a + U[nd][0], b + U[nd][1]47            dx, dy = ta - x, tb - y48            if max(abs(dx), abs(dy), abs(dx + dy)) > rem:49                continue50            if (x, y) in seen:51                continue52            seen.add((x, y))53            path.append(nd)54            dfs(path, seen, x, y, nd)55            path.pop()56            seen.discard((x, y))57    dfs([], {(0, 0)}, 0, 0, None)58    return out5960def copy(d, f, g, fl, m):61    n = len(g)62    if f == F:63        return [((d + g[j]) % 6, fl[j]) for j in range(n)]64    if f == R:65        return [((d + g[n - 1 - j]) % 6, REV[fl[n - 1 - j]]) for j in range(n)]66    if f == M:67        return [((d + m - g[j]) % 6, MIR[fl[j]]) for j in range(n)]68    return [((d + m - g[n - 1 - j]) % 6, REV[MIR[fl[n - 1 - j]]]) for j in range(n)]6970def expand(seq, g, fl, m):71    out = []72    for d, f in seq:73        out.extend(copy(d, f, g, fl, m))74    return out7576def points(seq):77    a = b = 078    P = [(0, 0)]79    for d, _ in seq:80        a, b = a + U[d][0], b + U[d][1]81        P.append((a, b))82    return P8384def avoids(seq):85    a = b = 086    seen = {(0, 0)}87    for d, _ in seq:88        a, b = a + U[d][0], b + U[d][1]89        if (a, b) in seen:90            return False91        seen.add((a, b))92    return True9394def counts(seq):95    c = [0, 0, 0, 0]96    for (d1, _), (d2, _) in zip(seq, seq[1:]):97        c[abs(turn(d1, d2))] += 198    return c99100def ends(g, fl, m):101    out = {}102    for f in (F, R, M, MR):103        if m is None and f >= M:104            continue105        s = copy(0, f, g, fl, m)106        out[f] = (s[0][0], s[0][1], s[-1][0], s[-1][1])107    return out108109def orbit(t, f, h, E):110    seen = set()111    x = (t % 6, f, h)112    worst = 0113    while x not in seen:114        seen.add(x)115        worst = max(worst, abs(turn(0, x[0])))116        t, f, h = x117        x = ((t + E[h][0] - E[f][2]) % 6, E[f][3], E[h][1])118    return worst, seen119120def closure(g, fl, m):121    E = ends(g, fl, m)122    types = set()123    flags = set(fl)124    todo = [(g[j + 1] - g[j], fl[j], fl[j + 1]) for j in range(len(g) - 1)]125    fdone = set()126    while todo or flags - fdone:127        while todo:128            x = todo.pop()129            x = (x[0] % 6, x[1], x[2])130            if x in types:131                continue132            types.add(x)133            t, f, h = x134            todo.append(((t + E[h][0] - E[f][2]) % 6, E[f][3], E[h][1]))135        for f in list(flags - fdone):136            fdone.add(f)137            s = copy(0, f, g, fl, m)138            flags.update(x[1] for x in s)139            todo.extend((b[0] - a[0], a[1], b[1]) for a, b in zip(s, s[1:]))140    return max(abs(turn(0, t)) for t, _, _ in types)141142def flag_search(g, flags, m, gentle=True):143    n = len(g)144    shapes = {}145    for d in range(6):146        for f in flags:147            steps = copy(d, f, g, [F] * n, m)148            shapes[(d, f)] = (points(steps)[1:], steps[0][0], steps[-1][0])149    tj = [(g[j + 1] - g[j]) % 6 for j in range(n - 1)]150    res = []151    for f0 in flags:152        for fz in flags:153            proto = [f0] + [F] * (n - 2) + [fz]154            E = ends(g, proto, m)155            good = {(t, f, h): not gentle or orbit(t, f, h, E)[0] <= 1 for t in range(6) for f in flags for h in flags}156            def dfs(j, fl, used, o, last):157                if j == n:158                    res.append(tuple(fl))159                    return160                cand = [f0] if j == 0 else ([fz] if j == n - 1 else flags)161                for f in cand:162                    if j and not good[(tj[j - 1], fl[-1], f)]:163                        continue164                    K, first, end = shapes[(g[j], f)]165                    if last is not None and abs(turn(last, first)) > (1 if gentle else 2):166                        continue167                    P = [(o[0] + p[0], o[1] + p[1]) for p in K]168                    if any(p in used for p in P):169                        continue170                    used.update(P)171                    fl.append(f)172                    dfs(j + 1, fl, used, P[-1], end)173                    fl.pop()174                    used.difference_update(P)175            dfs(0, [], {(0, 0)}, (0, 0), None)176    return res177178def canon(g, fl, mm):179    orb = [(g, fl), (g[::-1], fl[::-1])]180    if mm is not None:181        h = tuple((mm - d) % 6 for d in g)182        orb += [(h, fl), (h[::-1], fl[::-1])]183    return min(orb)184185A229214 = (1, 2, -1, 3, 1, 1, -3, 1, 2, 2, -1, -2, 3, 2, 3, -1, -1, -3, 1, -2, -1, 3, -1, -3, -2, 3, 3, 2, 1, 2, -1, 3, 1, 1, -3, 1, 2, -1, 3, 1, 1, -3, -2, -3, -3, 2, 3, 1, -3)186187CHORDS = [(3, (1, 1)), (4, (2, 0)), (7, (2, 1)), (9, (3, 0)), (12, (2, 2)), (13, (3, 1)), (16, (4, 0)), (19, (3, 2))]188189def census(orders, top, stop, mirrors=True):190    for n, c in CHORDS:191        if n not in orders:192            continue193        t0 = time.time()194        mm = mirror_index(c)195        m = mm if mirrors else None196        flags = [F, R] if m is None else [F, R, M, MR]197        G = walks(n, c, True)198        W = G[:stop] if stop else G199        lv2 = sharp = 0200        cross = {}201        alive = []202        kept = set()203        for g in W:204            for fl in flag_search(g, flags, m):205                lv2 += 1206                if closure(g, fl, m) > 1:207                    sharp += 1208                    continue209                kept.add(canon(g, fl, mm))210                seq = expand([(d, f) for d, f in zip(g, fl)], g, fl, m)211                L = 2212                while L < top:213                    seq = expand(seq, g, fl, m)214                    L += 1215                    if not avoids(seq):216                        cross[L] = cross.get(L, 0) + 1217                        break218                else:219                    alive.append((g, fl))220        mir = "F R" if m is None else "F R M MR"221        sym = "reversal" if mm is None else "reversal and mirror"222        if len(W) == len(G):223            tail = f"classes up to {sym} {len(kept)}, exhausted"224        else:225            tail = f"CUT after walk {len(W)} of {len(G)}"226        print(f"order {n} chord {c} flags {mir}: gentle walks {len(G)}, (generator, flags) pairs gentle and point-avoiding at level 2 by the typed search {lv2}, sharp at a later level by the junction closure {sharp}, gentle at every level {lv2 - sharp}, first crossing level {sorted(cross.items())}, alive at level {top} {len(alive)}, {tail}, {time.time() - t0:.1f}s")227        for g, fl in alive[:4]:228            print("  alive", g, [NAMES[f] for f in fl])229        sys.stdout.flush()230231def curve(g, fl, m, k):232    seq = [(0, F)]233    for _ in range(k):234        seq = expand(seq, g, fl, m)235    return seq236237def harmonic(seq):238    c = [0, 0, 0]239    for d, _ in seq:240        c[d % 3] += 1241    return c, c[0] ** 2 + c[1] ** 2 + c[2] ** 2 - c[0] * c[1] - c[1] * c[2] - c[0] * c[2]242243def gosper():244    g, fl = (0, 5, 3, 4, 0, 0, 1), (F, R, R, F, F, F, R)245    for k in range(1, 7):246        seq = curve(g, fl, None, k)247        tc = counts(seq)248        want = [7 ** (k - 1), 4 * 7 ** (k - 1) - 1, 2 * 7 ** (k - 1), 0]249        assert tc == want, (k, tc, want)250        assert k > 5 or avoids(seq)251        ax, h2 = harmonic(seq)252        assert h2 == 7 ** k253        if k == 2:254            assert [(-d) % 6 for d, _ in seq] == [x - 1 if x > 0 else 2 - x for x in A229214]255        print(f"gosper level {k}: steps {len(seq)}, turns 0/60/120 {tc[:3]}, axis counts {ax}, |h2|^2 {h2}")256257def lsys(rules, axiom, k):258    s = axiom259    for _ in range(k):260        s = "".join(rules.get(ch, ch) for ch in s)261    return s262263def square(rules, axiom, k):264    d = 0265    h = v = 0266    tc = [0, 0, 0]267    last = None268    for ch in lsys(rules, axiom, k):269        if ch == "+":270            d = (d + 1) % 4271        elif ch == "-":272            d = (d - 1) % 4273        elif ch == "F":274            if last is not None:275                tc[min((d - last) % 4, (last - d) % 4)] += 1276            last = d277            if d % 2:278                v += 1279            else:280                h += 1281    return h, v, tc282283def rate(orders):284    for n, c in CHORDS:285        if n not in orders:286            continue287        t0 = time.time()288        m = mirror_index(c)289        flags = [F, R] if m is None else [F, R, M, MR]290        pairs = 0291        rows = {}292        for g in walks(n, c, False):293            for fl in flag_search(g, flags, m, False):294                seq = curve(g, fl, m, 4)295                if not avoids(seq) or counts(seq)[3]:296                    continue297                pairs += 1298                s = [counts(curve(g, fl, m, k))[2] for k in range(1, 6)]299                rows[canon(g, fl, m)] = s300        sym = "reversal" if m is None else "reversal and mirror"301        print(f"order {n} chord {c}: (generator, flags) pairs avoiding points through level 4 {pairs}, classes up to {sym} {len(rows)}, {time.time() - t0:.1f}s")302        for (g, fl), s in sorted(rows.items(), key=lambda x: (x[1], x[0])):303            print(f"  {g} {[NAMES[f] for f in fl]}: 120-degree turns at levels 1..5 {s}, per unit step at level 5 {s[-1]}/{n ** 5}")304        sys.stdout.flush()305306def bead():307    hil = {"A": "+BF-AFA-FB+", "B": "-AF+BFB+FA-"}308    pea = {"X": "XFYFX+F+YFXFY-F-XFYFX", "Y": "YFXFY-F-XFYFX+F+YFXFY"}309    for k in range(1, 7):310        h, v, tc = square(hil, "A", k)311        assert h + v == 4 ** k - 1312        print(f"hilbert level {k}: steps {h + v}, horizontal {h}, vertical {v}, h2 {h - v}, turns 0/90 {tc[:2]}")313    for k in range(1, 5):314        h, v, tc = square(pea, "X", k)315        assert h + v == 9 ** k - 1316        print(f"peano level {k}: steps {h + v}, horizontal {h}, vertical {v}, h2 {h - v}, turns 0/90 {tc[:2]}")317318if __name__ == "__main__":319    verb = sys.argv[1] if len(sys.argv) > 1 else "census"320    orders = [int(x) for x in sys.argv[2].split(",")] if len(sys.argv) > 2 else None321    stop = int(sys.argv[3]) if len(sys.argv) > 3 else None322    if verb == "census":323        census(orders or [3, 4, 7, 9, 13], 4, stop)324    if verb == "plain":325        census(orders or [12], 4, stop, False)326    if verb == "gosper":327        gosper()328    if verb == "bead":329        bead()330    if verb == "rate":331        rate(orders or [7])