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
Ada
pragma Ada_2022;
with Ada.Containers.Vectors;
with Ada.Numerics; use Ada.Numerics;
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;
with Ada.Text_IO; use Ada.Text_IO;
procedure Babylonian_Spiral is
type Fixed_10_4 is delta 10.0 ** (-4) digits 9;
type Point is record
X, Y : Integer;
end record;
type Points_Array is array (1 .. 10_000) of Point;
Points : Points_Array := [1 => (0, 0), 2 => (0, 1), others => (0, 0)];
Ante_Previous_Point, Previous_Point : Point;
package Point_Vectors is new Ada.Containers.Vectors (Positive, Point);
use Point_Vectors;
Point_Vector, PVC : Point_Vectors.Vector;
Min_Radius : Integer := 1;
Previous_Dist : Fixed_10_4; -- The previous distance which must be exceeded
Min_Dist, Tmp_Dist : Fixed_10_4;
Max_Angle, Tmp_Angle : Fixed_10_4;
Tmp_Ix, Cand_Ix : Integer;
Fixed_Pi : constant Fixed_10_4 := 3.1416;
SVG_File : File_Type;
X, Y : Integer;
function Distance (Point1, Point2 : Point) return Fixed_10_4 is
(Fixed_10_4 (Sqrt (Float (Point2.X - Point1.X)**2 +
Float (Point2.Y - Point1.Y)**2)));
function Generate_Square_Of_Points (Centre : Point; Max_Radius : Positive)
return Point_Vectors.Vector is
PV : Point_Vectors.Vector;
begin
for X in Centre.X - Max_Radius .. Centre.X + Max_Radius loop
for Y in Centre.Y - Max_Radius .. Centre.Y + Max_Radius loop
PV.Append ((X, Y), 1);
end loop;
end loop;
return PV;
end Generate_Square_Of_Points;
function Atan2 (Y, X : Float) return Float is
Res : Float;
begin
if X > 0.0 then Res := Arctan (Y / X);
elsif X < 0.0 and then Y >= 0.0 then Res := Arctan (Y / X) + Pi;
elsif X < 0.0 and then Y < 0.0 then Res := Arctan (Y / X) - Pi;
elsif X = 0.0 and then Y > 0.0 then Res := Pi / 2.0;
elsif X = 0.0 and then Y > 0.0 then Res := -Pi / 2.0;
else Res := -Pi / 2.0; -- Technically: Undefined
end if;
return Res;
end Atan2;
function Angle (Centre, P2, P3 : Point) return Fixed_10_4 is
Res : Float;
begin
Res := Atan2 (Float (P3.Y) - Float (Centre.Y), Float (P3.X) - Float (Centre.X)) -
Atan2 (Float (P2.Y) - Float (Centre.Y), Float (P2.X) - Float (Centre.X));
return Fixed_10_4 (Res);
end Angle;
function Collinear (Pt1, Pt2, Pt3 : Point) return Boolean is
begin
if Pt1.X = Pt2.X and then Pt2.X = Pt3.X then
return True;
end if;
if Pt1.Y = Pt2.Y and then Pt2.Y = Pt3.Y then
return True;
end if;
return Pt1.X * (Pt2.Y - Pt3.Y) + Pt2.X * (Pt3.Y - Pt1.Y) + Pt3.X * (Pt1.Y - Pt2.Y) = 0;
end Collinear;
begin
for Point_Ix in 3 .. Points'Last loop
Ante_Previous_Point := Points (Point_Ix - 2);
Previous_Point := Points (Point_Ix - 1);
Previous_Dist := Distance (Points (Point_Ix - 1), Points (Point_Ix - 2));
Min_Radius := Integer (Previous_Dist);
Min_Dist := 99999.999;
Max_Angle := 0.0;
-- examine surrounding points
Point_Vector := Generate_Square_Of_Points (Previous_Point, Min_Radius + 4);
for Pt of Point_Vector loop
if not Collinear (Pt, Previous_Point, Ante_Previous_Point) then
Tmp_Dist := Distance (Pt, Previous_Point);
if Tmp_Dist > Previous_Dist and then Tmp_Dist < Min_Dist then
Min_Dist := Tmp_Dist;
end if;
end if;
end loop;
-- Grab closest Points
PVC.Clear;
for Pt of Point_Vector loop
if Distance (Pt, Previous_Point) = Min_Dist then
PVC.Append (Pt);
end if;
end loop;
Tmp_Ix := 1;
for Pt of PVC loop
Tmp_Angle := Angle (Previous_Point, Ante_Previous_Point, Pt);
if Tmp_Angle < -(Fixed_Pi) then
Tmp_Angle := Tmp_Angle + (2.0 * Fixed_Pi);
end if;
if Tmp_Angle < Fixed_Pi and then Tmp_Angle > Max_Angle then
Cand_Ix := Tmp_Ix;
Max_Angle := Tmp_Angle;
end if;
Tmp_Ix := Tmp_Ix + 1;
end loop;
if Point_Ix < 41 then
Put ("(" & PVC (Cand_Ix).X'Image & "," & PVC (Cand_Ix).Y'Image & ") ");
end if;
Points (Point_Ix) := PVC (Cand_Ix);
end loop;
Create (SVG_File, Out_File, "babylonian_spiral.svg");
Put_Line (SVG_File, "<svg xmlns='http://www.w3.org/2000/svg' width='12000' height='15000'>");
Put_Line (SVG_File, "<rect width='100%' height='100%' fill='white'/>");
X := 500 + Points (1).X;
Y := 10000 - Points (1).Y;
Put (SVG_File, "<path stroke-width='2' stroke='black' fill='none' d='M" &
X'Image & "," & Y'Image);
for Pt_Ix in 2 .. Points'Last loop
X := 500 + Points (Pt_Ix).X;
Y := 10000 - Points (Pt_Ix).Y;
Put (SVG_File, " L" & X'Image & "," & Y'Image);
end loop;
Put_Line (SVG_File, "'/>\n</svg>");
Close (SVG_File);
end Babylonian_Spiral;
- Output:
( 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)
https://static.miraheze.org/rosettacodewiki/c/c7/Babylonian_spiral.svg
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
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
""" 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))
M2000 Interpreter
module Babylonian_spiral (max=40) {
declare math math
const TAU=6.283185307179586
dim sq()
sq()=lambda (m)-> {
for i=0 to m-1:data i*i:next
=array([])
}(max)
xydeltas=((0&, 0&), (0&, 1&))
def isqrt(n)=val(sqrt(n)->long)
def pr(c)="("+c#str$(", ")+")"
long δsquared = 1
for L=0 to max-3
(x,y)=xydeltas#val(len(xydeltas)-1)
θ = @atan2(y, x)
candidates=stack
while len(candidates)=0
δsquared++
for i=0& to max-1
a=sq(i)
if a > δsquared/2& then exit for
for j =isqrt(δsquared) + 1 to 1
b = sq(j)
if a + b < δsquared then exit for
if a + b = δsquared then
stack candidates {
Data (i, j), (-i, j), (i, -j), (-i, -j), (j, i), (-j, i), (j, -i), (-j, -i)
}
end if
next
next
end while
minimum=(,)
minVal=TAU
stack candidates {
while not empty
read candidate()
val = (θ - @atan2(candidate(1), candidate(0))) mod# TAU
if val < minVal then minVal = val : minimum = candidate()
end while
}
append xyDeltas, (minimum,)
next
for i=0 to len(xyDeltas)-2
Return xyDeltas, i+1:=@add(xyDeltas#val(i+1), xyDeltas#val(i))
next
Drawing 12000, 12000 {
Cls, 0
pen 0
(m1, m2, s1)=(6000, 6000, twipsX*10)
move m1, m2
k=each(xyDeltas)
dim b()
while k
b()=array(k)
b=b()
b*=s1
draw to b(0)+m1, m2-b(1): circle fill 1, 30
end while
} as Drw
image drw
clipboard drw as "emf"
Document Doc$
for i=0 to len(xyDeltas)-2
Doc$=pr(xyDeltas#val(i))+", "
next
Doc$=pr(xyDeltas#val(i))
Print Doc$
Save.Doc Doc$, "out.txt"
function add(a(), b())
=a(0)+b(0), a(1)+b(1)
end function
function atan2(y, x)
method math, "atan2", y, x as ret
=ret // radians
end function
}
Babylonian_spiral 40
- Output:
to file
(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)
MATLAB
Directly translate Julia code
% 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)
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
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
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
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
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
Zig
const std = @import("std");
const ArrayList = std.ArrayList;
const print = std.debug.print;
const Vec = struct {
x: i32,
y: i32,
pub fn sum_sq(v: Vec) i32 {
return v.x * v.x + v.y * v.y;
}
pub fn add(a: Vec, b: Vec) Vec {
return .{
.x = a.x + b.x,
.y = a.y + b.y,
};
}
pub fn atan(v: Vec) f32 {
const x: f32 = @floatFromInt(v.x);
const y: f32 = @floatFromInt(v.y);
return std.math.atan2(x, y);
}
pub fn lessThan(ctx: void, a: Vec, b: Vec) std.math.Order {
_ = ctx;
return std.math.order(a.sum_sq(), b.sum_sq());
}
};
pub fn add_candidates(candidates: *ArrayList(Vec), base: Vec) !void {
const i = base.x;
const j = base.y;
const candidate_steps = [_]Vec{
.{ .x = i, .y = j },
.{ .x = -i, .y = j },
.{ .x = i, .y = -j },
.{ .x = -i, .y = -j },
.{ .y = i, .x = j },
.{ .y = -i, .x = j },
.{ .y = i, .x = -j },
.{ .y = -i, .x = -j },
};
try candidates.appendSlice(&candidate_steps);
}
pub fn babylonian_spiral(n_points: usize, allocator: std.mem.Allocator) !ArrayList(Vec) {
const tau: f32 = std.math.tau;
var queue = std.PriorityQueue(Vec, void, Vec.lessThan).init(allocator, {});
defer queue.deinit();
try queue.add(.{ .x = 1, .y = 1 });
var point = Vec{ .x = 0, .y = 1 };
var last = Vec{ .x = 0, .y = 1 };
var last_angle = last.atan();
var points = ArrayList(Vec).init(allocator);
try points.ensureUnusedCapacity(n_points);
points.appendAssumeCapacity(.{ .x = 0, .y = 0 });
points.appendAssumeCapacity(.{ .x = 0, .y = 1 });
var candidates = ArrayList(Vec).init(allocator);
defer candidates.deinit();
var sum_sq: i32 = 1;
var loaded: i32 = 1;
var count: usize = 2;
while (count < n_points) : (count += 1) {
if (loaded * loaded < sum_sq + 1) {
loaded += 1;
var y: i32 = 0;
while (y <= loaded) : (y += 1) {
try queue.add(.{ .x = loaded, .y = y });
}
}
var base = queue.remove();
try add_candidates(&candidates, base);
sum_sq = base.sum_sq();
while (queue.peek().?.sum_sq() == sum_sq) {
base = queue.remove();
try add_candidates(&candidates, base);
}
var min = candidates.items[0];
var min_angle = min.atan() - last_angle;
if (min_angle < 0) min_angle += tau;
for (candidates.items[1..]) |step| {
var angle = step.atan() - last_angle;
if (angle < 0) angle += tau;
if (angle < min_angle) {
min = step;
min_angle = angle;
}
}
point = point.add(min);
points.appendAssumeCapacity(point);
last = min;
last_angle = min.atan();
candidates.clearRetainingCapacity();
}
return points;
}
pub fn main() !void {
var gpa = std.heap.GeneralPurposeAllocator(.{}){};
const allocator = gpa.allocator();
var points = try babylonian_spiral(10_000, allocator);
defer points.deinit();
print("The first 40 Babylonian spiral points are:\n", .{});
for (points.items[0..40], 0..) |p, n| {
print("({},{}) ", .{ p.x, p.y });
if (n % 10 == 9) {
print("\n", .{});
}
}
const file = try std.fs.cwd().createFile("babylonian_spiral.svg", .{});
defer file.close();
var buffered_writer = std.io.bufferedWriter(file.writer());
defer {
buffered_writer.flush() catch {
print("Failed to flush buffer to file", .{});
};
}
var w = buffered_writer.writer();
var min_xi: i32 = 0;
var max_xi: i32 = 0;
var min_yi: i32 = 0;
var max_yi: i32 = 0;
for (points.items) |p| {
if (min_xi > p.x) min_xi = p.x;
if (min_yi > p.y) min_yi = p.y;
if (max_xi < p.x) max_xi = p.x;
if (max_yi < p.y) max_yi = p.y;
}
const min_x: f32 = @floatFromInt(min_xi);
const min_y: f32 = @floatFromInt(min_yi);
const max_x: f32 = @floatFromInt(max_xi);
const max_y: f32 = @floatFromInt(max_yi);
const width = max_x - min_x;
const height = max_y - min_y;
const view_width: u32 = 1000;
const view_height: u32 = 1000;
try w.print(
\\<svg viewBox="0 0 {} {}" xmlns="http://www.w3.org/2000/svg">
\\<polyline fill="none" stroke="black" points="
, .{ view_width, view_height });
for (points.items) |p| {
const x = (@as(f32, @floatFromInt(p.x)) - min_x) * view_width / width;
const y = (@as(f32, @floatFromInt(p.y)) - min_y) * view_height / height;
try w.print("{d:.0},{d:.0} ", .{ x, view_height - y });
}
try w.print("{s}", .{"\"/></svg>"});
}
- 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)
