Babylonian spiral

Babylonian spiral
You are encouraged to solve this task according to the task description, using any language you may know.

f 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(1) and P(2) 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(3) that fits criteria is (x = 1, y = 2). Note the length of the vector from P(2) to P(3) is √2, which squared is 2. The lengths of the vectors thus determined can be given by a sorted list of possible 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 graph shown at the OEIS:

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.

C++

#include <cmath>
#include <iomanip>
#include <iostream>
#include <numbers>
#include <vector>

typedef std::pair<int32_t, int32_t> point;

std::vector<point> babylonian_spiral(int32_t step_count) {
	const double tau = 2 * std::numbers::pi;

	std::vector<int32_t> squares(step_count + 1);
	for ( int32_t i = 0; i < step_count; ++i ) {
		squares[i] = i * i;
	}
	std::vector<point> points { point(0, 0), point (0, 1) };
	int32_t norm = 1;

	for ( int32_t step = 0; step < step_count - 2; ++step ) {
		point previousPoint = points.back();
		const double theta = atan2(previousPoint.second, previousPoint.first);
		std::vector<point> candidates;
		while ( candidates.empty() ) {
			norm += 1;
			for ( int32_t i = 0; i < step_count; ++i ) {
				int32_t a = squares[i];
				if ( a > norm / 2 ) {
					break;
				}
				for ( int32_t j = sqrt(norm) + 1; j >= 0; --j ) {
					int32_t b = squares[j];
					if ( a + b < norm ) {
						break;
					}
					if ( a + b == norm ) {
						std::vector<point> next_points = { point(i, j), point(-i, j), point(i, -j), point(-i, -j),
								                           point(j, i), point(-j, i), point(j, -i), point(-j, -i) };
						candidates.reserve(next_points.size());
						std::move(next_points.begin(), next_points.end(), std::back_inserter(candidates));
					}
				}
			}
		}

		point minimum;
		double minimum_value = tau;
		for ( point candidate : candidates ) {
			double value = fmod(theta - atan2(candidate.second, candidate.first) + tau, tau);
			if ( value < minimum_value ) {
			    minimum_value = value;
			    minimum = candidate;
			}
		}

		points.push_back(minimum);
	}

	for ( uint64_t i = 0; i < points.size() - 1; ++i ) {
		points[i + 1] = point(points[i].first + points[i + 1].first, points[i].second + points[i + 1].second);
	}
	return points;
}

int main() {
	std::vector<point> points = babylonian_spiral(40);
	std::cout << "The first 40 points of the Babylonian spiral are:" << std::endl;
	for ( int32_t i = 0, column = 0; i < 40; ++i ) {
		std::string point_str = "(" + std::to_string(points[i].first) + ", " + std::to_string(points[i].second) + ")";
		std::cout << std::setw(10) << point_str;
		if ( ++column % 10 == 0 ) {
			std::cout << std::endl;
		}
	}
}

Output:

The first 40 points of the Babylonian spiral 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)

J

It's convenient, here, to use complex numbers to represent the coordinate pairs.

Implementation:

require'stats'
bspir=: {{
  r=. 0
  e=. 2>.<.+:%:y
  for_qd.
    (y-1){.(</.~ *:@|) (/:|) (#~ (=<.)@:*:@:|) j./"1 (2 comb e),,.~1+i.e
  do.
    d=. ~.(,+)(,-)(,j.);qd
    ar=. 12 o. -~/_2{.r
    ad=. (- 2p1 * >:&ar) 12 o. d
    -r=. r, ({:r)+d{~ (i. >./) ad
  end.
}}

Task example:

   4 10$bspir 40 
    0   0j1   1j2   3j2   5j1   7j_1   7j_4   6j_7  4j_10  0j_10
_4j_9 _7j_6 _9j_2  _9j3  _8j8  _6j13  _2j17   3j20   9j20  15j19
21j17 26j13  29j7    29 28j_7 24j_13 17j_15 10j_12   4j_7    4j1
  5j9  7j17 13j23 21j26 28j21  32j13   32j4  31j_5 29j_14 24j_22

Also:

   require'plot'
   plot bspir 40
   plot bspir 1e4

There's an online example of these two cases here (follow the link, select the right pane, then hit control-R or select Edit Run > All Lines from the menu, this might take a couple seconds, depending on your machine, and this url may need to be updated in the future).

This could be made significantly faster with a better estimator (or with a better implementation of J compiled to javascript).

Java

import java.awt.Point;
import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Paths;
import java.util.Comparator;
import java.util.HashSet;
import java.util.List;
import java.util.Set;
import java.util.stream.Collectors;
import java.util.stream.IntStream;
import java.util.stream.Stream;

public final class BabylonianSpiral {

	public static void main(String[] aArgs) throws IOException {
		List<Point> points = babylonianSpiral(10_000);
		System.out.println("The first 40 points of the Babylonian spiral are:");
		for ( int i = 0, column = 0; i < 40; i++ ) {
			System.out.print(String.format("%9s%s",
				"(" + points.get(i).x + ", " + points.get(i).y + ")", ( ++column % 10 == 0 ) ? "\n" : " "));
		}
		System.out.println();
	
		String text = svgText(points, 800);		
		Files.write(Paths.get("C:/Users/psnow/Desktop/BabylonianSpiralJava.svg"), text.getBytes());
	}
	
	private static List<Point> babylonianSpiral(int aStepCount) {
		final double tau = 2 * Math.PI;
		List<Integer> squares = IntStream.rangeClosed(0, aStepCount).map( i -> i * i ).boxed().toList();
		List<Point> points = Stream.of( new Point(0, 0), new Point(0, 1) ).collect(Collectors.toList());
		int norm = 1;
		
		for ( int step = 0; step < aStepCount - 2; step++ ) {
		    Point previousPoint = points.get(points.size() - 1);
		    final double theta = Math.atan2(previousPoint.y, previousPoint.x);
		    Set<Point> candidates = new HashSet<Point>();
		    while ( candidates.isEmpty() ) {
		    	norm += 1;
			    for ( int i = 0; i < aStepCount; i++ ) {
			        int a = squares.get(i);
			        if ( a > norm / 2 ) {
			        	break;
			        }
			        for ( int j = (int) Math.sqrt(norm) + 1; j >= 0; j-- ) {
			        	int b = squares.get(j);	
			        	if ( a + b < norm ) {
			        		break;
			        	}
			        	if ( a + b == norm ) {
			        		candidates.addAll(
			        			List.of( new Point(i, j), new Point(-i, j), new Point(i, -j), new Point(-i, -j),
			        					 new Point(j, i), new Point(-j, i), new Point(j, -i), new Point(-j, -i) ));
			        	}
			        }
			    }
		    }
		    
		    Comparator<Point> comparatorPoint = (one, two) -> Double.compare(
		    	( theta - Math.atan2(one.y, one.x) + tau ) % tau, ( theta - Math.atan2(two.y, two.x) + tau ) % tau);
		    
		    Point minimum = candidates.stream().min(comparatorPoint).get();
		    points.add(minimum);
		}
		
		for ( int i = 0; i < points.size() - 1; i++ ) {
			points.set(i + 1, new Point(points.get(i).x + points.get(i + 1).x, points.get(i).y + points.get(i + 1).y));
		}
		return points;
	}

	private static String svgText(List<Point> aPoints, int aSize) {
    	StringBuilder text = new StringBuilder();
    	text.append("<svg xmlns='http://www.w3.org/2000/svg'");
        text.append(" width='" + aSize + "' height='" + aSize + "'>\n");
        text.append("<rect width='100%' height='100%' fill='cyan'/>\n");
        text.append("<path stroke-width='1' stroke='black' fill='cyan' d='");
        for ( int i = 0; i < aPoints.size(); i++ ) {
        	text.append(( i == 0 ? "M" : "L" ) +
        		( 150 + aPoints.get(i).x / 20 ) + ", " + ( 525 - aPoints.get(i).y / 20 ) + " ");
        }
        text.append("'/>\n</svg>\n");
        
        return text.toString();
    }
	
}

Output:

Media:BabylonianSpiralJava.svg

The first 40 points of the Babylonian spiral 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)

jq

Works with: jq

Also works with gojq, the Go implementation of jq

This entry is adapted from the Python and Wren entries but with an unabashedly stream-oriented approach, as this minimizes memory requirements.

In particular, the main function, `babylonianSpiral`, will generate points indefinitely, with the main limitation being the arithmetic precision imposed by the implementation of jq being used. Since gojq, the Go implementation of jq, supports unbounded integer arithmetic, this means that, from a practical perspective, the accuracy of the stream is assured when using gojq no matter how many points are generated.

Since gnuplot does not require the specification beforehand of the "xrange" or "yrange", the output of `babylonianSpiral` can be directly piped to gnuplot, as illustrated below.

# Emit an unbounded stream of points of the Babylonian spiral:
def babylonianSpiral:

  def heading($x):
    def tau: 6.283185307179586; 
    # Note: the Python modulo operator is not the same as Wren's or jq's
    fmod($x + tau; tau);

  [0,0],
  [0,1],
  ( { dxys: [0,1],   # the last point
       dsq : 1,
       sumx: 0,
       sumy: 1 }
    | foreach range(2; infinite) as $k (.;
         .dxys as [ $x, $y ]
        | atan2($y; $x) as $theta
        | .candidates = []
        | until(.candidates != [];
            .dsq += 1
            | .i = 0
            | until(.i >= $k;
                (.i*.i) as $a
                | if $a > ((.dsq/2)|floor) then .i = $k       # i.e., break
                  else .j = (.dsq|sqrt|floor) + 1
                  | until(.j <= 0;
                      (.j*.j) as $b
                      | if ($a + $b) < .dsq then .j = -1      # i.e., break
                        elif ($a + $b) == .dsq
                        then .candidates += [ [.i, .j], [-.i, .j], [.i, -.j], [-.i, -.j],
                                              [.j, .i], [-.j, .i], [.j, -.i], [-.j, -.i] ]
                        else .
                        end
                        | .j += -1 )
                  | .i += 1
                  end ) )
        # Python: lambda d: (θ - atan2(d[1], d[0])) % tau
        | .dxys = (.candidates | min_by( heading( ($theta - atan2(.[1]; .[0])) ) ))
        | .sumx += .dxys[0]
        | .sumy += .dxys[1];
        [.sumx, .sumy] )
  );

# Emit a stream of the first $n points:
def Points($n):
  limit($n; babylonianSpiral);

Print the first 40 in groups of 10 per line

"The first 40 Babylonian spiral points are:",
 ([Points(40)] | _nwise(10) | map(tostring) | join(" ") )

For plotting 10000 points

Points(10000) | join(" ")

Visualizing 10000 points with gnuplot

Assuming the .jq program generates points along the lines of:

   Points(10000) | join(" ") 

the following invocation of jq and gnuplot will produce a .png file showing the Babylonian spiral with the given number of points:

jq -rnf babylonian-spiral.jq | gnuplot -c <(cat <<EOF
reset
set terminal pngcairo
set output 'example-babylonian-spiral.png'
plot "/dev/stdin" with dots
EOF
)

Julia

Translation of: Python

""" Rosetta Code task rosettacode.org/wiki/Babylonian_spiral """

using GLMakie

const squarecache = Int[]

"""
Get the points for a Babylonian spiral of `nsteps` steps. Origin is at (0, 0)
with first step one unit in the positive direction along the vertical (y) axis.
See also: oeis.org/A256111, oeis.org/A297346, oeis.org/A297347
"""
function babylonianspiral(nsteps)
    if length(squarecache) <= nsteps
        append!(squarecache, map(x -> x * x, length(squarecache):nsteps))
    end
    # first line segment is 1 unit in vertical direction, with y vertical, x horizontal
    xydeltas, δ² = [(0, 0), (0, 1)], 1
    for _ in 1:nsteps-2
        x, y = xydeltas[end]
        θ = atan(y, x)
        candidates = Tuple{Int64,Int64}[]
        while isempty(candidates)
            δ² += 1
            for (k, a) in enumerate(squarecache)
                a > δ² ÷ 2 && break
                for j in isqrt(δ²)+1:-1:1
                    b = squarecache[j+1]
                    a + b < δ² && break
                    if a + b == δ²
                        i = k - 1
                        push!(
                            candidates,
                            (i, j),
                            (-i, j),
                            (i, -j),
                            (-i, -j),
                            (j, i),
                            (-j, i),
                            (j, -i),
                            (-j, -i),
                        )
                    end
                end
            end
        end
        _, idx = findmin(p -> mod1(θ - atan(p[2], p[1]), 2π), candidates)
        push!(xydeltas, candidates[idx])
    end
    return accumulate((a, b) -> (a[1] + b[1], a[2] + b[2]), xydeltas)
end

println("The first 40 Babylonian spiral points are:")
for (i, p) in enumerate(babylonianspiral(40))
    print(rpad(p, 10), i % 10 == 0 ? "\n" : "")
end

lines(babylonianspiral(10_000), linewidth = 1)

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)

iterator version

See Python generator version.

using DataStructures
using Plots

struct SquareSums end

function Base.iterate(ss::SquareSums, state = (1, PriorityQueue{Tuple{Int, Int, Int}, Int}()))
    i, sums = state
    while isempty(sums) || i^2 <= peek(sums)[2]
        enqueue!(sums, (i^2, i, 0) => i^2)
        i += 1
    end
    nextsum, xy = peek(sums)[2], Tuple{Int, Int}[]
    while !isempty(sums) && peek(sums)[2] == nextsum # pop all vectors with same length
        nextsum, a, b = dequeue!(sums)
        push!(xy, (a, b))
        if a > b
            hsq = a^2 + (b + 1)^2
            enqueue!(sums, (hsq, a, b + 1) => hsq)
        end
    end
    return xy, (i, sums)
end

function babylonian_spiral(N)
    dx, dy, points = 0, 1, [(0, 0)]
    for xys in SquareSums()
        for i in 1:length(xys)
            a, b = xys[i]
            a != b && push!(xys, (b, a))
            a != 0 && push!(xys, (-a, b), (b, -a))
            b != 0 && push!(xys, (a, -b), (-b, a))
            a * b != 0 && push!(xys, (-a, -b), (-b, -a))
        end
        filter!(p -> p[1] * dy - p[2] * dx >= 0, xys)
        _, idx = findmax(p -> p[1] * dx + p[2] * dy, xys)
        dx, dy = xys[idx]
        push!(points, (points[end][1] + dx, points[end][2] + dy))
        length(points) >= N && break
    end
    return @view points[begin:N]
end

println("The first 40 Babylonian spiral points are:")
for (i, p) in enumerate(babylonian_spiral(40))
    print(rpad(p, 10), i % 10 == 0 ? "\n" : "")
end

Plots.plot(babylonian_spiral(10_000))

MATLAB

Directly translate Julia code

Translation of: Julia

% Rosetta Code task rosettacode.org/wiki/Babylonian_spiral

clear all;close all;clc;
% Example usage
fprintf("The first 40 Babylonian spiral points are:\n");
spiral_points = babylonianspiral(40);
for i = 1:size(spiral_points, 1)
    fprintf('(%d, %d) ', spiral_points(i, 1), spiral_points(i, 2));
    if mod(i, 10) == 0
        fprintf('\n');
    end
end

% For plotting the spiral (requires MATLAB plotting functions)
spiral_points = babylonianspiral(10000);
plot(spiral_points(:, 1), spiral_points(:, 2), 'LineWidth', 1);


function points = babylonianspiral(nsteps)
% Get the points for a Babylonian spiral of `nsteps` steps. Origin is at (0, 0)
% with first step one unit in the positive direction along the vertical (y) axis.
% See also: oeis.org/A256111, oeis.org/A297346, oeis.org/A297347
    
    persistent squarecache;
    if isempty(squarecache)
        squarecache = [];
    end

    if length(squarecache) <= nsteps
        squarecache = [squarecache, arrayfun(@(x) x^2, length(squarecache):nsteps)];
    end

    xydeltas = [0, 0; 0, 1];
    deltaSq = 1;
    for i = 1:nsteps-2
        x = xydeltas(end, 1);
        y = xydeltas(end, 2);
        theta = atan2(y, x);
        candidates = [];
        while isempty(candidates)
            deltaSq = deltaSq + 1;
            for k = 1:length(squarecache)
                a = squarecache(k);
                if a > deltaSq / 2
                    break;
                end
                for j = floor(sqrt(deltaSq)):-1:1
                    b = squarecache(j+1);
                    if a + b < deltaSq
                        break;
                    end
                    if a + b == deltaSq
                        i = k - 1;
                        candidates = [candidates; i, j; -i, j; i, -j; -i, -j; ...
                                                     j, i; -j, i; j, -i; -j, -i];
                    end
                end
            end
        end
        [~, idx] = min(arrayfun(@(n) mod(theta - atan2(candidates(n, 2), candidates(n, 1)), 2*pi), 1:size(candidates, 1)));
        xydeltas = [xydeltas; candidates(idx, :)];
    end

    points = cumsum(xydeltas);
end

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) 

File:BabylonianSpiral.png

Nim

Translation of: Python

Library: gnuplotlib.nim

import std/[math, strformat, strutils]

type Vector = tuple[x, y: int]

func `+`(a, b: Vector): Vector =
  ## Return the sum of two vectors.
  (a.x + b.x, a.y + b.y)

func isqrt(n: int): int =
  ## Return the integer square root of "n".
  int(sqrt(n.toFloat))


proc babylonianSpiral(nsteps: int): seq[Vector] =
  ## 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 change of the  new vector is
  ## chosen to be nonzero and clockwise in a direction that minimizes the change
  ## in direction from the previous vector.

  var xyDeltas: seq[Vector] = @[(0, 0), (0, 1)]
  var δsquared = 1
  for _ in 0..nsteps - 3:
    let (x, y) = xyDeltas[^1]
    let θ = arctan2(y.toFloat, x.toFloat)
    var candidates: seq[Vector]
    while candidates.len == 0:
      inc δsquared, 1
      for i in 0..<nsteps:
        let a = i * i
        if a > δsquared div 2:
          break
        for j in countdown(isqrt(δsquared) + 1, 1):
          let b = j * j
          if a + b < δsquared:
            break
          if a + b == δsquared:
            candidates.add [(i, j), (-i, j), (i, -j), (-i, -j), (j, i), (-j, i), (j, -i), (-j, -i)]
      var p: Vector
      var minVal = TAU
      for candidate in candidates:
        let val = floorMod(θ - arctan2(candidate.y.toFloat, candidate.x.toFloat), TAU)
        if val < minVal:
          minVal = val
          p = candidate
      xyDeltas.add p

  result = cumsummed(xyDeltas)

let points10000 = babylonianSpiral(10_000)


### Task ###
echo "The first 40 Babylonian spiral points are:"
for i, p in points10000[0..39]:
  stdout.write alignLeft(&"({p.x}, {p.y})", 10)
  stdout.write if (i + 1) mod 10 == 0: "\n" else: ""


### Stretch task ###

import gnuplot

var x, y: seq[float]
for p in points10000:
  x.add p.x.toFloat
  y.add p.y.toFloat

withGnuPlot:
  cmd "set size ratio -1"
  plot(x, y, "Babylonian spiral", "with lines lc black lw 1")
  png("babylonian_spiral.png")

Output:

The first 40 Babylonian spiral points are:
(0, 0)    (0, 1)    (1, 2)    (1, 2)    (3, 2)    (5, 1)    (5, 1)    (5, 1)    (7, -1)   (7, -4)   
(6, -7)   (6, -7)   (6, -7)   (4, -10)  (4, -10)  (4, -10)  (0, -10)  (-4, -9)  (-7, -6)  (-7, -6)  
(-9, -2)  (-9, -2)  (-9, -2)  (-9, -2)  (-9, -2)  (-9, 3)   (-8, 8)   (-8, 8)   (-8, 8)   (-6, 13)  
(-6, 13)  (-6, 13)  (-2, 17)  (-2, 17)  (3, 20)   (3, 20)   (9, 20)   (15, 19)  (15, 19)  (15, 19)  

Perl

Translation of: Raku

use strict;
use warnings;
use feature <say state>;
use constant TAU => 2 * 2 * atan2(1, 0);

sub B_spiral {
    my($nsteps) = @_;
    my @squares = map $_**2, 0..$nsteps+1;
    my @dxys = ([0, 0], [0, 1]);
    my $dsq  = 1;

    for (1 .. $nsteps-2) {
        my ($x,$y) = @{$dxys[-1]};
        our $theta = atan2 $y, $x;
        my @candidates;

        until (@candidates) {
            $dsq++;
            for my $i (0..$#squares) {
                my $a = $squares[$i];
                next if $a > $dsq/2;
                for my $j ( reverse 0 .. 1 + int sqrt $dsq ) {
                    my $b = $squares[$j];
                    next if ($a + $b) < $dsq;
                    if ($dsq == $a + $b) {
                        push @candidates, ( [$i, $j], [-$i, $j], [$i, -$j], [-$i, -$j],
                                            [$j, $i], [-$j, $i], [$j, -$i], [-$j, -$i] );
                    }
                }
            }
        }

        sub comparer {
            my $i = ($theta - atan2 $_[1], $_[0]);
            my $z = $i - int($i / TAU) * TAU;
            $z < 0 ? TAU + $z : $z;
        }

        push @dxys, (sort { comparer(@$b) < comparer(@$a) } @candidates)[0];
    }

    map { state($x,$y); $x += $$_[0]; $y += $$_[1]; [$x,$y] } @dxys;
}

my @points = map { sprintf "(%3d,%4d)", @$_ } B_spiral(40);
say "The first 40 Babylonian spiral points are:\n" .
    join(' ', @points) =~ s/.{1,88}\K/\n/gr;

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)

Phix

Library: Phix/pGUI

Library: Phix/online

You can run this online here. Use left/right arrow keys to show less/more edges.

--
-- demo/rosetta/Babylonian_spiral.exw
-- ==================================
--
with javascript_semantics
function next_step(atom last_distance)
    // Find "the next longest vector with integral endpoints on a Cartesian grid"
    integer next_distance = 100*last_distance, // Set high so we can minimize
            nmax = floor(sqrt(last_distance)) + 2
            // ^ The farthest we could possibly go in one direction
    sequence next_steps = {}
    for n=0 to nmax do
        integer n2 = n*n,
              mmin = floor(sqrt(max(0,last_distance-n2)))
        for m=mmin to nmax do
            integer test_distance = n2 + m*m
            if test_distance>last_distance then
                if test_distance>next_distance then exit end if
                if test_distance<next_distance then
                    next_distance = test_distance
                    next_steps = {}
                end if
                next_steps &= {{m,n}}
                if m!=n then
                    next_steps &= {{n,m}}
                end if
            end if
        end for
    end for
    return {next_steps, next_distance}
end function

sequence x = {0,0},         -- first two points
         y = {0,1}          --  taken as given
integer distance = 1,
        px = 0, py = 1,     -- position
        pdx = 0, pdy = 1    -- previous delta

procedure make_spiral(integer npoints) // Make a Babylonian spiral of npoints.
    sequence deltas 
    atom t4 = time()+0.4
    for n=length(x)+1 to npoints do
        {deltas,distance} = next_step(distance)
        atom max_dot = 0, ldx = pdx, ldy = pdy
        for delta in deltas do
            integer {tx,ty} = delta
            for d in {{tx,ty},{-tx,ty},{tx,-ty},{-tx,-ty}} do
                integer {dx,dy} = d
                if ldx*dy-ldy*dx<0 then
                    atom dot = ldx*dx+ldy*dy
                    if dot>max_dot then
                        max_dot = dot
                        {pdx,pdy} = {dx,dy}
                    end if
                end if
            end for
        end for
        px += pdx
        py += pdy
        x &= px
        y &= py
        if time()>t4 then exit end if
    end for
end procedure

make_spiral(40)
printf(1,"The first 40 Babylonian spiral points are:\n%s\n",
         {join_by(columnize({x,y}),1,10," ",fmt:="(%3d,%3d)")})

include pGUI.e
include IupGraph.e
Ihandle dlg, graph, timer
constant mt = {{5,6,1}, -- {number, minmax, tick}
               {10,12,2},  -- add more/remove steps as desired
               {20,20,4},   
               {40,32,8},   
               {60,70,10},  
               {100,220,40},
               {200,350,50},
               {400,700,100},
               {800,1200,200},
               {1000,2000,400},
               {10000,12000,3000},
               {100000,150000,50000},
               {150000,150000,50000},
               {200000,400000,100000}}
                -- perfectly doable, but test yer patience:
--             {250000,400000,100000},
--             {500000,1200000,400000}}
integer mdx = 4, -- index to mt
        max_mdx = 4

function get_data(Ihandle /*graph*/)
    integer {n,m,t} = mt[mdx]
    string title = sprintf("Babylonian spiral (%,d)", {n})
    if length(x)<mt[$][1] then
        title &= sprintf(" (calculating %,d/%,d)",{length(x),mt[$][1]})
    end if
    IupSetStrAttribute(graph, "GTITLE", title)
    sequence xn = x[1..n],
             yn = y[1..n]
    IupSetAttributes(graph,"XMIN=%d,XMAX=%d,XTICK=%d",{-m,m,t})
    IupSetAttributes(graph,"YMIN=%d,YMAX=%d,YTICK=%d",{-m,m,t})
    return {{xn,yn,CD_RED}}
end function

function key_cb(Ihandle /*ih*/, atom c)
    if c=K_ESC then return IUP_CLOSE end if -- (standard practice for me)
    if c=K_F5 then return IUP_DEFAULT end if -- (let browser reload work)
    if c=K_LEFT then mdx = max(mdx-1,1) end if
    if c=K_RIGHT then mdx = min(mdx+1,max_mdx) end if
    IupUpdate(graph)
    return IUP_CONTINUE
end function

function timer_cb(Ihandln /*ih*/)
    if max_mdx=length(mt) or not IupGetInt(dlg,"VISIBLE") then
        IupSetInt(timer,"RUN",false)
    else
        integer next_target = mt[max_mdx+1][1]
        make_spiral(next_target)
        if length(x)=next_target then
            max_mdx += 1
            if mdx=max_mdx-1 then
                mdx = max_mdx
            end if
        end if
    end if
    IupRedraw(graph)
    return IUP_IGNORE
end function

IupOpen()
graph = IupGraph(get_data,"RASTERSIZE=640x480,XMARGIN=35")
dlg = IupDialog(graph,`TITLE="Babylonian spiral",MINSIZE=270x430`)
IupSetInt(graph,"GRIDCOLOR",CD_LIGHT_GREY)
IupShow(dlg)
IupSetCallback(dlg, "K_ANY", Icallback("key_cb"))
IupSetAttribute(graph,"RASTERSIZE",NULL)
timer = IupTimer(Icallback("timer_cb"), 30)
if platform()!=JS then
    IupMainLoop()
    IupClose()
end if

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)

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 change 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(nsteps - 2):
        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 '')

# stretch portion of task
plot(*zip(*points10000), color="navy", linewidth=0.2)
axis('scaled')
show()

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)

With priority queue

Use a priority queue to generate all x, y combinations. The advantage is that we don't need to do any real math, and it is much faster.

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()

Raku

Translation of: Wren

Translation

sub babylonianSpiral (\nsteps) {
    my @squareCache = (0..nsteps).hyper.map: *²;
    my @dxys = [0, 0], [0, 1];
    my $dsq  = 1;

    for ^(nsteps-2) {
        my \Θ = atan2 |@dxys[*-1][1,0];
        my @candidates;

        until @candidates.elems {
            $dsq++;
            for @squareCache.kv -> \i, \a {
                last if a > $dsq/2;
                for reverse 0 .. $dsq.sqrt.ceiling -> \j {
                    last if $dsq > (a + my \b = @squareCache[j]);
                    next if $dsq != a + b;
                    @candidates.append: [i, j], [-i, j], [i, -j], [-i, -j],
                                        [j, i], [-j, i], [j, -i], [-j, -i]
                }
            }
        }
        @dxys.push: @candidates.min: { ( Θ - atan2 |.[1,0] ) % τ };
    }

    [\»+«] @dxys
}

# The task
say "The first $_ Babylonian spiral points are:\n",
(babylonianSpiral($_).map: { sprintf '(%3d,%4d)', @$_ }).batch(10).join("\n") given 40;

# Stretch
use SVG;

'babylonean-spiral-raku.svg'.IO.spurt: SVG.serialize(
    svg => [
        :width<100%>, :height<100%>,
        :rect[:width<100%>, :height<100%>, :style<fill:white;>],
        :polyline[ :points(flat babylonianSpiral(10000)),
          :style("stroke:red; stroke-width:6; fill:white;"),
          :transform("scale (.05, -.05) translate (1000,-10000)")
        ],
    ],
);

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)

Stretch goal:

(offsite SVG image - 96 Kb) - babylonean-spiral-raku.svg

Independent implementation

Exact same output; about one tenth the execution time.

my @next = { :x(1), :y(1), :2hyp },;

sub next-interval (Int $int) {
     @next.append: (0..$int).map: { %( :x($int), :y($_), :hyp($int² + .²) ) };
     @next = |@next.sort: *.<hyp>;
}

my @spiral = [\»+«] lazy gather {
    my $interval = 1;
    take [0,0];
    take my @tail = 0,1;
    loop {
        my \Θ = atan2 |@tail[1,0];
        my @this = @next.shift;
        @this.push: @next.shift while @next and @next[0]<hyp> == @this[0]<hyp>;
        my @candidates = @this.map: {
            my (\i, \j) = .<x y>;
            next-interval(++$interval) if $interval == i;
            |((i,j),(-i,j),(i,-j),(-i,-j),(j,i),(-j,i),(j,-i),(-j,-i))
        }
        take @tail = |@candidates.min: { ( Θ - atan2 |.[1,0] ) % τ };
    }
}

# The task
say "The first $_ Babylonian spiral points are:\n",
@spiral[^$_].map({ sprintf '(%3d,%4d)', |$_ }).batch(10).join: "\n" given 40;

# Stretch
use SVG;

'babylonean-spiral-raku.svg'.IO.spurt: SVG.serialize(
    svg => [
        :width<100%>, :height<100%>,
        :rect[:width<100%>, :height<100%>, :style<fill:white;>],
        :polyline[ :points(flat @spiral[^10000]),
          :style("stroke:red; stroke-width:6; fill:white;"),
          :transform("scale (.05, -.05) translate (1000,-10000)")
        ],
    ],
);

Same output:

Wren

Translation of: Python

Library: DOME

Library: Wren-iterate

Library: Wren-seq

Library: Wren-fmt

Library: Wren-plot

Generates an image similar to the OEIS one.

import "dome" for Window
import "graphics" for Canvas, Color
import "./iterate" for Indexed, Stepped
import "./seq" for Lst
import "./fmt" for Fmt
import "./plot" for Axes

// Python modulo operator (not same as Wren's)
var pmod = Fn.new { |x, y| ((x % y) + y) % y }

var squareCache = []

"""
    Get the points for each step along a Babylonian 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 change 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
"""
var babylonianSpiral = Fn.new { |nsteps|
    for (x in 0...nsteps) squareCache.add(x*x)
    var dxys = [[0, 0], [0, 1]]
    var dsq = 1
    for (i in 0...nsteps) {
        var x = dxys[-1][0]
        var y = dxys[-1][1]
        var theta = y.atan(x)
        var candidates = []
        while (candidates.isEmpty) {
            dsq = dsq + 1
            for (se in Indexed.new(squareCache)) {
                var i = se.index
                var a = se.value
                if (a > (dsq/2).floor) break
                for (j in dsq.sqrt.floor + 1...0) {
                    var b = squareCache[j]
                    if ((a + b) < dsq) break
                    if ((a + b) == dsq) {
                        candidates.addAll([ [i, j], [-i, j], [i, -j], [-i, -j],
                                            [j, i], [-j, i], [j, -i], [-j, -i] ])
                    }
                }
            }
        }
        var comparer = Fn.new { |d| pmod.call(theta - d[1].atan(d[0]), Num.tau) }
        candidates.sort { |a, b| comparer.call(a) < comparer.call(b) }
        dxys.add(candidates[0])
    }

    var accs = []
    var sumx = 0
    var sumy = 0
    for (dxy in dxys) {
        sumx = sumx + dxy[0]
        sumy = sumy + dxy[1]
        accs.add([sumx, sumy])
    }
    return accs
}

// find first 10,000 points
var Points10000 = babylonianSpiral.call(9998) // first two added automatically

// print first 40 to terminal
System.print("The first 40 Babylonian spiral points are:")
for (chunk in Lst.chunks(Points10000[0..39], 10)) Fmt.print("$-9s", chunk)

class Main {
    construct new() {
        Window.title = "Babylonian spiral"
        Canvas.resize(1000, 1000)
        Window.resize(1000, 1000)
        Canvas.cls(Color.white)
        var axes = Axes.new(100, 900, 800, 800, -1000..11000, -5000..10000)
        axes.draw(Color.black, 2)
        var xMarks = Stepped.new(0..10000, 2000)
        var yMarks = Stepped.new(-5000..10000, 5000)
        axes.label(xMarks, yMarks, Color.black, 2, Color.black)
        axes.line(-1000, 10000, 11000, 10000, Color.black, 2)
        axes.line(11000, -5000, 11000, 10000, Color.black, 2)
        axes.lineGraph(Points10000, Color.black, 2)
    }

    init() {}

    update() {}

    draw(alpha) {}
}

var Game = Main.new()

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]

File:Wren-Babylonian spiral.png