Babylonian spiral: Difference between revisions
m
→{{header|Wren}}: Made displayed image a bit smaller.
(→{{header|Wren}}: Converted to a DOME script so we can use Wren-plot to plot the points.) |
m (→{{header|Wren}}: Made displayed image a bit smaller.) |
||
| (40 intermediate revisions by 10 users not shown) | |||
Line 1:
{{
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,
Line 26:
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: [[File:From_oies_A@97346_A297347_plot2a.png]]
; See also
Line 36 ⟶ 35:
=={{header|C++}}==
<syntaxhighlight lang="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;
}
}
}
</syntaxhighlight>
{{ out }}
<pre>
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)
</pre>
=={{header|J}}==
Line 43 ⟶ 129:
Implementation:
<syntaxhighlight lang="j">
require'stats'
bspir=: {{
Line 57 ⟶ 143:
end.
}}
</syntaxhighlight>
Task example:
<
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
</syntaxhighlight>
Also:
<syntaxhighlight lang="j">
require'plot'
plot bspir 40
plot bspir 1e4
</syntaxhighlight>
There's an online example of
This could be made significantly faster with a better estimator (or with a better implementation of J compiled to javascript).
=={{header|Java}}==
<syntaxhighlight lang="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();
}
}
</syntaxhighlight>
{{ out }}
[[Media:BabylonianSpiralJava.svg]]
<pre>
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)
</pre>
=={{header|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.
<syntaxhighlight lang=jq>
# 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);
</syntaxhighlight>
'''Print the first 40 in groups of 10 per line'''
<syntaxhighlight lang=jq>
"The first 40 Babylonian spiral points are:",
([Points(40)] | _nwise(10) | map(tostring) | join(" ") )
</syntaxhighlight>
'''For plotting 10000 points'''
<syntaxhighlight lang=jq>
Points(10000) | join(" ")
</syntaxhighlight>
'''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:
<pre>
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
)
</pre>
=={{header|Julia}}==
{{trans|Python}}
<
using GLMakie
Line 139 ⟶ 424:
lines(babylonianspiral(10_000), linewidth = 1)
</
<pre>
The first 40 Babylonian spiral points are:
Line 146 ⟶ 431:
(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)
</pre>
[[File:babylonuanspiral.png]]
=== iterator version ===
See Python generator version.
<syntaxhighlight lang="julia">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))
</syntaxhighlight>
=={{header|MATLAB}}==
=== Directly translate Julia code ===
{{trans|Julia}}
<syntaxhighlight lang="MATLAB}}">
% 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
</syntaxhighlight>
{{out}}
<pre>
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)
</pre>
[[File:BabylonianSpiral.png ]]
=={{header|Nim}}==
{{trans|Python}}
{{libheader|gnuplotlib.nim}}
<syntaxhighlight lang="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")
</syntaxhighlight>
{{out}}
<pre>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)
</pre>
=={{header|Perl}}==
{{trans|Raku}}
<
use warnings;
use feature <say state>;
Line 196 ⟶ 705:
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;</
{{out}}
<pre>The first 40 Babylonian spiral points are:
Line 208 ⟶ 717:
{{libheader|Phix/online}}
You can run this online [http://phix.x10.mx/p2js/Babylonian_spiral.htm here]. Use left/right arrow keys to show less/more edges.
<!--<
<span style="color: #000080;font-style:italic;">--
-- demo/rosetta/Babylonian_spiral.exw
Line 353 ⟶ 862:
<span style="color: #7060A8;">IupClose</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<!--</
{{out}}
<pre>
Line 364 ⟶ 873:
=={{header|Python}}==
<
from itertools import accumulate
Line 418 ⟶ 927:
# stretch portion of task
plot(*zip(*points10000), color="navy", linewidth=0.2)
axis('scaled')
show()
</
<pre>
The first 40 Babylonian spiral points are:
Line 429 ⟶ 938:
(5, 9) (7, 17) (13, 23) (21, 26) (28, 21) (32, 13) (32, 4) (31, -5) (29, -14) (24, -22)
</pre>
[[File:babspiral.png]]
=== 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.
<syntaxhighlight lang="python">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()</syntaxhighlight>
=={{header|Raku}}==
Line 434 ⟶ 993:
===Translation===
<syntaxhighlight lang="raku"
my @squareCache = (0..nsteps).hyper.map: *²;
my @dxys =
my $dsq = 1;
Line 445 ⟶ 1,004:
until @candidates.elems {
$dsq++;
last if
next if
}
}
Line 478 ⟶ 1,036:
],
],
);</
{{out}}
<pre>
Line 493 ⟶ 1,051:
Exact same output; about one tenth the execution time.
<syntaxhighlight lang="raku"
sub next-interval (Int $int) {
Line 533 ⟶ 1,091:
],
],
);</
{{out|Same output}}
Line 539 ⟶ 1,097:
{{trans|Python}}
{{libheader|DOME}}
{{libheader|Wren-
{{libheader|Wren-seq}}
{{libheader|Wren-fmt}}
{{libheader|Wren-plot}}
Generates an image similar to the OEIS one.
<
import "graphics" for Canvas, Color
import "./
import "./seq" for Lst
import "./fmt" for Fmt
Line 638 ⟶ 1,196:
}
var Game = Main.new()</
{{out}}
Line 648 ⟶ 1,206:
[5, 9] [7, 17] [13, 23] [21, 26] [28, 21] [32, 13] [32, 4] [31, -5] [29, -14] [24, -22]
</pre>
[[File:Wren-Babylonian_spiral.png|600px]]
| |||