Rosetta Code

Babylonian spiral: Difference between revisions

Page Discussion

Babylonian spiral (view source)

Revision as of 10:14, 19 February 2024

15,464 bytes added ,  2 months ago
m
→‎{{header|Wren}}: Made displayed image a bit smaller.
m (→‎{{header|Julia}}: iterator version)
m (→‎{{header|Wren}}: Made displayed image a bit smaller.)
 
(25 intermediate revisions by 6 users not shown)
Line 1:
{{draft task}}
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 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 78 ⟶ 165:
 
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}}==
Line 149 ⟶ 435:
 
=== iterator version ===
See Python iteratorgenerator version.
<syntaxhighlight lang="julia">using PlotsDataStructures
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)[1][12]
pushenqueue!(sums, (i^2, i, 0) => i^2)
i += 1
end
nextsum, xy = peek(sums)[2], Tuple{Int, Int}[]
sort!(sums)
while !isempty(sums) && peek(sums)[2] == nextsum # pop all vectors with same length
nextsum, xy = sums[1][1], Tuple{Int, Int}[]
while !isempty(sums) && sums[1][1] == nextsum, #a, popb all= vectors with same lengthdequeue!(sums)
nextsum, a, b = popfirst!(sums)
push!(xy, (a, b))
if a > b
push!(sums,hsq = (a^2 + (b + 1)^2, a, b + 1))
enqueue!(sums, (hsq, a, b + 1) => hsq)
end
end
Line 173 ⟶ 460:
 
function babylonian_spiral(N)
sqsumsdx, dy, points = SquareSums0, 1, [(0, 0)]
dx,for dy,xys xydeltasin = 0, 1, [SquareSums(0, 0)]
for xys in sqsums
for i in 1:length(xys)
a, b = xys[i]
Line 186 ⟶ 472:
_, idx = findmax(p -> p[1] * dx + p[2] * dy, xys)
dx, dy = xys[idx]
push!(xydeltaspoints, (points[end][1] + dx, points[end][2] + dy))
length(xydeltaspoints) >= N && break
end
return accumulate((a, b) -> (a[1] + b[1], a[2] + b[2]), @view xydeltaspoints[begin:N])
end
 
Line 199 ⟶ 485:
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}}==
Line 640 ⟶ 1,097:
{{trans|Python}}
{{libheader|DOME}}
{{libheader|Wren-traititerate}}
{{libheader|Wren-seq}}
{{libheader|Wren-fmt}}
{{libheader|Wren-plot}}
Generates an image similar to the OEIS one.
<syntaxhighlight lang="ecmascriptwren">import "dome" for Window
import "graphics" for Canvas, Color
import "./traititerate" for Indexed, Stepped
import "./seq" for Lst
import "./fmt" for Fmt
Line 749 ⟶ 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]]
PureFox
9,476

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/330011...355574"