Babylonian spiral
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 page at https://oeis.org/plot2a?name1=A297346&name2=A297347&tform1=untransformed&tform2=untransformed&shift=0&radiop1=xy&drawlines=true".
- See also
Phix
You can run this online here.
-- -- 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 nmax = floor(sqrt(last_distance)) + 2 // ^^ The farthest we could possibly go in one direction atom next_distance = 100*last_distance // Set high so we can minimize sequence next_steps = {} for n=0 to nmax do integer n2 = n*n for m=n 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 function make_spiral(integer npoints) // Make a Babylonian spiral of npoints. sequence x = {0,0}, y = {0,1}, deltas atom distance = 1 integer px = 0, py = 1, -- position pdx = 0, pdy = 1 -- previous delta for n=3 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 end for return {x,y} end function sequence {x,y} = make_spiral(10000), first40 = columnize({x[1..40],y[1..40]}) printf(1,"The first 40 Babylonian spiral points are:\n%s\n", {join_by(first40,1,10,fmt:="(%3d,%3d)")}) include pGUI.e include IupGraph.e Ihandle dlg, graph integer p10 = 4 constant mt = {{12,1}, -- 10 {220,20}, -- 100 {2000,400}, -- 1000 {12000,3000}} -- 10000 function get_data(Ihandle /*graph*/) integer n = power(10,p10) IupSetStrAttribute(graph, "GTITLE", "Babylonian spiral (%d)", {n}) sequence xn = x[1..n], yn = y[1..n] integer {m,t} = mt[p10] 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 p10 = max(p10-1,1) end if if c=K_RIGHT then p10 = min(p10+1,round(log10(length(x)))) end if IupUpdate(graph) return IUP_CONTINUE end function IupOpen() graph = IupGraph(get_data,"RASTERSIZE=640x480") 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) 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
<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 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):
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)) 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)
Wren
<lang ecmascript>import "./trait" for Indexed import "./seq" for Lst import "./fmt" for Fmt import "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][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)
// create .csv file for all 10,000 points for display by an external plotter File.create("babylonian_spiral.csv") { |file|
for (p in points10000) {
file.writeBytes("%(p[0]), %(p[1])\n")
}
}</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]