Babylonian spiral: Difference between revisions
m
→{{header|Wren}}: Made displayed image a bit smaller.
m (→{{header|Wren}}: Wren-trait -> Wren-iterate) |
m (→{{header|Wren}}: Made displayed image a bit smaller.) |
||
| (19 intermediate revisions by 4 users not shown) | |||
Line 1:
{{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}}==
Line 297 ⟶ 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 743 ⟶ 1,102:
{{libheader|Wren-plot}}
Generates an image similar to the OEIS one.
<syntaxhighlight lang="
import "graphics" for Canvas, Color
import "./iterate" for Indexed, Stepped
Line 847 ⟶ 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]]
| |||