I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

# 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(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.

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

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".

## 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.}} `

`   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  15j1921j17 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 the 1e4 case here (follow the link, then hit run, this might take a couple seconds, depending on your machine).

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

## 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, p), 2π), candidates)        push!(xydeltas, candidates[idx])    end    return accumulate((a, b) -> (a + b, a + b), 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)
```

## 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 \$_, \$_);            my \$z = \$i - int(\$i / TAU) * TAU;            \$z < 0 ? TAU + \$z : \$z;        }         push @dxys, (sort { comparer(@\$b) < comparer(@\$a) } @candidates);    }     map { state(\$x,\$y); \$x += \$\$_; \$y += \$\$_; [\$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[\$] then
title &= sprintf(" (calculating %,d/%,d)",{length(x),mt[\$]})
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]
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 accumulatefrom math import isqrt, atan2, taufrom 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, d)) % tau)        xydeltas.append(p)     return list(accumulate(xydeltas, lambda a, b: (a + b, a + b)))  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 taskplot(*zip(*points10000))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)
```

## 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 ( a + my \b = @squareCache[j] ) < \$dsq;                    if ((a + b) == \$dsq) {                        @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 tasksay "The first \$_ Babylonian spiral points are:\n",(babylonianSpiral(\$_).map: { sprintf '(%3d,%4d)', @\$_ }).batch(10).join("\n") given 40; # Stretchuse 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<hyp> == @this<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 tasksay "The first \$_ Babylonian spiral points are:\n",@spiral[^\$_].map({ sprintf '(%3d,%4d)', |\$_ }).batch(10).join: "\n" given 40; # Stretchuse 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: Wren-trait
Library: Wren-seq
Library: Wren-fmt
`import "./trait" for Indexedimport "./seq" for Lstimport "./fmt" for Fmtimport "io" for File // 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]        var y = dxys[-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.atan(d), Num.tau) }        candidates.sort { |a, b| comparer.call(a) < comparer.call(b) }        dxys.add(candidates)    }     var accs = []    var sumx = 0    var sumy = 0    for (dxy in dxys) {        sumx = sumx + dxy        sumy = sumy + dxy        accs.add([sumx, sumy])    }    return accs} // find first 10,000 pointsvar points10000 = babylonianSpiral.call(9998) // first two added automatically // print first 40 to terminalSystem.print("The first 40 Babylonian spiral points are:")for (chunk in Lst.chunks(points10000[0..39], 10)) Fmt.print("\$-9s", chunk) // create .csv file for all 10,000 points for display by an external plotterFile.create("babylonian_spiral.csv") { |file|    for (p in points10000) {        file.writeBytes("%(p), %(p)\n")    }}`
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]
```