45 lines
1.1 KiB
Python
45 lines
1.1 KiB
Python
from itertools import islice, count
|
|
import matplotlib.pyplot as plt
|
|
import heapq
|
|
|
|
def twosquares():
|
|
q, n = [], 1
|
|
|
|
while True:
|
|
while not q or n*n <= q[0][0]:
|
|
heapq.heappush(q, (n*n, n, 0))
|
|
n += 1
|
|
|
|
s, xy = q[0][0], []
|
|
|
|
while q and q[0][0] == s: # pop all vectors with same length
|
|
s, a, b = heapq.heappop(q)
|
|
xy.append((a, b))
|
|
if a > b:
|
|
heapq.heappush(q, (a*a + (b+1)*(b+1), a, b + 1))
|
|
|
|
yield tuple(xy)
|
|
|
|
def gen_dirs():
|
|
d = (0, 1)
|
|
for v in twosquares():
|
|
# include symmetric vectors
|
|
v += tuple((b, a) for a, b in v if a != b)
|
|
v += tuple((a, -b) for a, b in v if b)
|
|
v += tuple((-a, b) for a, b in v if a)
|
|
|
|
# filter using dot and cross product
|
|
d = max((a*d[0] + b*d[1], a, b) for a, b in v if a*d[1] - b*d[0] >= 0)[1:]
|
|
yield d
|
|
|
|
def positions():
|
|
p = (0, 0)
|
|
for d in gen_dirs():
|
|
yield p
|
|
p = (p[0] + d[0], p[1] + d[1])
|
|
|
|
print(list(islice(positions(), 40)))
|
|
|
|
plt.plot(*zip(*list(islice(positions(), 100000))), lw=0.4)
|
|
plt.gca().set_aspect(1)
|
|
plt.show()
|