Babylonian spiral

Babylonian spiral is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The Babylonian spiral is a sequence of points in the plane that are created so as to continuously minimally increase in vector length and minimally bend in vector direction, while always moving from point to point on strictly integral coordinates. Of the two criteria of length and angle, the length has priority.

Examples

P(0) and P(1) are defined to be at (x = 0, y = 0) and (x = 0, y = 1). The first vector is from P(1) to P(2). It is vertical and of length 1. Note that the square of that length is 1.

Next in sequence is the vector from P(2) to P(3). This should be the smallest distance to a point with integral (x, y) which is longer than the last vector (that is, 1). It should also bend clockwise more than zero radians, but otherwise to the least degree.

The point chosen for P(2) that fits criteria is (x = 1, y = 1). Note the length of this vector is √2, which squared is 2. The lengths of the vectors thus determined can be given by a sorted list of sums of two integer squares, including 0 as a square.

Task

Find and show the first 40 (x, y) coordinates of the Babylonian spiral.

Stretch task

Show in your program how to calculate and plot the first 10000 points in the sequence. Your result should look similar to the page at https://oeis.org/plot2a?name1=A297346&name2=A297347&tform1=untransformed&tform2=untransformed&shift=0&radiop1=xy&drawlines=true".

See also

• OEIS:A256111 - squared distance to the origin of the n-th vertex on a Babylonian Spiral.
• OEIS:A297346 - List of successive x-coordinates in the Babylonian Spiral.
• OEIS:A297347 - List of successive y-coordinates in the Babylonian Spiral.

Python

<lang python>""" Rosetta Code task rosettacode.org/wiki/Babylonian_spiral """

from itertools import accumulate from math import isqrt, atan2, tau from matplotlib.pyplot import axis, plot, show

square_cache = []

def babylonian_spiral(nsteps):

   """
   Get the points for each step along a Babylonia spiral of `nsteps` steps.
   Origin is at (0, 0) with first step one unit in the positive direction along
   the vertical (y) axis. The other points are selected to have integer x and y
   coordinates, progressively concatenating the next longest vector with integer
   x and y coordinates on the grid. The direction of the new vector is chosen to
   be nonzero and clockwise in a direction that minimizes the change in direction
   from the previous vector.
   
   See also: oeis.org/A256111, oeis.org/A297346, oeis.org/A297347
   """
   if len(square_cache) <= nsteps:
       square_cache.extend([x * x for x in range(len(square_cache), nsteps)])
   xydeltas = [(0, 0), (0, 1)]
   δsquared = 1
   for _ in range(1, nsteps+1):
       x, y = xydeltas[-1]
       θ = atan2(y, x)
       candidates = []
       while not candidates:
           δsquared += 1
           for i, a in enumerate(square_cache):
               if a > (δsquared // 2):
                   break
               for j in range(isqrt(δsquared) + 1, 0, -1):
                   b = square_cache[j]
                   if a + b < δsquared:
                       break
                   if a + b == δsquared:
                       candidates.extend([(i, j), (-i, j), (i, -j), (-i, -j), (j, i), (-j, i),
                          (j, -i), (-j, -i)])

       p = min(candidates, key=lambda d: (θ - atan2(d[1], d[0])) % tau)
       xydeltas.append(p)

   return list(accumulate(xydeltas, lambda a, b: (a[0] + b[0], a[1] + b[1])))

points10000 = babylonian_spiral(10000) print("The first 40 Babylonian spiral points are:") for i, p in enumerate(points10000[:40]):

    print(str(p).ljust(10), end = '\n' if (i + 1) % 10 == 0 else )

1. stretch portion of task

plot(*zip(*points10000)) axis('scaled') show()

</lang>

Output:

The first 40 Babylonian spiral points are:
(0, 0)    (0, 1)    (1, 2)    (3, 2)    (5, 1)    (7, -1)   (7, -4)   (6, -7)   (4, -10)  (0, -10)  
(-4, -9)  (-7, -6)  (-9, -2)  (-9, 3)   (-8, 8)   (-6, 13)  (-2, 17)  (3, 20)   (9, 20)   (15, 19)
(21, 17)  (26, 13)  (29, 7)   (29, 0)   (28, -7)  (24, -13) (17, -15) (10, -12) (4, -7)   (4, 1)
(5, 9)    (7, 17)   (13, 23)  (21, 26)  (28, 21)  (32, 13)  (32, 4)   (31, -5)  (29, -14) (24, -22)