import functools from heapq import * import numpy as np """ General approach: generate all pairs (m, v) in order of increasing m with a priority queue, where v is an integer vector, m is the square magnitude of v, and m > 1. Define wurm(n) to be the n-th vector of the Babylonian Wurm and turn(n) to be A256175(n). Calculate wurm(n) by popping all vectors of minimal magnitude currently in the priority queue and adding the one most similar to the previous vector to wurm(n-1), breaking ties as necessary. """ # square magnitude of v def norm(v): return v[0]**2 + v[1]**2 # positive area means counterclockwise, negative means clockwise def area(v, u): return v[0]*u[1] - v[1]*u[0] # higher dot product means greater similarity def dot(v, u): return v[0]*u[0] + v[1]*u[1] visited, pq = set(), [] # add a vector to pq def add_v(v): global visited, pq v = tuple(v) if v not in visited: heappush(pq, (norm(v), v)) visited.add(tuple(v)) # initialize pq and visited for i in range(9): visited.add((i // 3 - 1, i % 3 - 1)) for i in range(4): add_v([0, 2, 0, -2, 0][i:i+2]) # pops all vectors of minimal magnitude def get_vs(): global pq m, vs = heappop(pq) vs = [np.array(vs)] while len(pq) > 0 and norm(pq[0][1]) == m: vs.append(np.array(heappop(pq)[1])) for v in vs: # add adjacent vectors to pq for i in range(4): add_v(v + [0, 1, 0, -1, 0][i:i+2]) return vs @functools.lru_cache(maxsize=None) def wurm(n): if n < 3: return np.array([(0, 0), (0, 1), (1, 2)][n]) # base case for i in range(n): wurm(i) # this keeps the recursion tree short vs = get_vs() u = wurm(n-1) - wurm(n-2) # previous vector dots = [dot(u, v) for v in vs] m = max(dots) ms = [i for i, val in enumerate(dots) if val == m] if len(ms) == 1 or turn(n-2) * -area(u, vs[ms[0]]) > 0: return wurm(n-1) + vs[ms[0]] return wurm(n-1) + vs[ms[1]] def turn(n): return np.sign(-area(wurm(n-1) - wurm(n), wurm(n) - wurm(n+1))) v = [turn(i) for i in range(1, 10001)] for i in range(len(v)): print(i+1, v[i])