Total circles area: Difference between revisions
Content added Content deleted
(+ C entry) |
(+ D entry) |
||
Line 130: | Line 130: | ||
{{out}} |
{{out}} |
||
<pre>Approximated area: 21.56262288</pre> |
<pre>Approximated area: 21.56262288</pre> |
||
=={{header|D}}== |
|||
This version converges much faster than both the ordered grid and Montecarlo sampling solutions. |
|||
<lang d>import std.stdio, std.math, std.algorithm, std.typecons, std.range; |
|||
alias real Fp; |
|||
struct Circle { Fp x, y, r; } |
|||
void removeInternalDisks(ref Circle[] circles) { |
|||
static bool isFullyInternal(in Circle c1, in Circle c2) |
|||
pure nothrow { |
|||
if (c1.r > c2.r) // quick exit |
|||
return false; |
|||
return (c1.x - c2.x) ^^ 2 + (c1.y - c2.y) ^^ 2 < |
|||
(c2.r - c1.r) ^^ 2; |
|||
} |
|||
// Heuristics for performance: large radii first. |
|||
circles.sort!q{a.r > b.r}(); |
|||
// What circles will be kept. |
|||
auto keep = new bool[circles.length]; |
|||
keep[] = true; |
|||
foreach (i1, c1; circles) |
|||
if (keep[i1]) |
|||
foreach (i2, c2; circles) |
|||
if (keep[i2]) |
|||
if (i1 != i2 && isFullyInternal(c2, c1)) |
|||
keep[i2] = false; |
|||
// Pack circles array, removing fully internal circles. |
|||
size_t pos = 0; |
|||
foreach (i, k; keep) |
|||
if (k) |
|||
circles[pos++] = circles[i]; |
|||
circles.length = pos; |
|||
// Alternative implementation of the packing: |
|||
// circles = zip(circles, keep) |
|||
// .filter!(ck => ck[1])() |
|||
// .map!(ck => ck[0])() |
|||
// .array(); |
|||
} |
|||
void main() { |
|||
Circle[] circles = [ |
|||
{ 1.6417233788, 1.6121789534, 0.0848270516}, |
|||
{-1.4944608174, 1.2077959613, 1.1039549836}, |
|||
{ 0.6110294452, -0.6907087527, 0.9089162485}, |
|||
{ 0.3844862411, 0.2923344616, 0.2375743054}, |
|||
{-0.2495892950, -0.3832854473, 1.0845181219}, |
|||
{ 1.7813504266, 1.6178237031, 0.8162655711}, |
|||
{-0.1985249206, -0.8343333301, 0.0538864941}, |
|||
{-1.7011985145, -0.1263820964, 0.4776976918}, |
|||
{-0.4319462812, 1.4104420482, 0.7886291537}, |
|||
{ 0.2178372997, -0.9499557344, 0.0357871187}, |
|||
{-0.6294854565, -1.3078893852, 0.7653357688}, |
|||
{ 1.7952608455, 0.6281269104, 0.2727652452}, |
|||
{ 1.4168575317, 1.0683357171, 1.1016025378}, |
|||
{ 1.4637371396, 0.9463877418, 1.1846214562}, |
|||
{-0.5263668798, 1.7315156631, 1.4428514068}, |
|||
{-1.2197352481, 0.9144146579, 1.0727263474}, |
|||
{-0.1389358881, 0.1092805780, 0.7350208828}, |
|||
{ 1.5293954595, 0.0030278255, 1.2472867347}, |
|||
{-0.5258728625, 1.3782633069, 1.3495508831}, |
|||
{-0.1403562064, 0.2437382535, 1.3804956588}, |
|||
{ 0.8055826339, -0.0482092025, 0.3327165165}, |
|||
{-0.6311979224, 0.7184578971, 0.2491045282}, |
|||
{ 1.4685857879, -0.8347049536, 1.3670667538}, |
|||
{-0.6855727502, 1.6465021616, 1.0593087096}, |
|||
{ 0.0152957411, 0.0638919221, 0.9771215985}]; |
|||
writeln("Input Circles: ", circles.length); |
|||
removeInternalDisks(circles); |
|||
writeln("Circles left: ", circles.length); |
|||
immutable Fp xMin = reduce!((acc, c) => min(acc, c.x - c.r)) |
|||
(Fp.max, circles[]); |
|||
immutable Fp xMax = reduce!((acc, c) => max(acc, c.x + c.r)) |
|||
(cast(Fp)0, circles[]); |
|||
alias Tuple!(Fp,"y0", Fp,"y1") yRange; |
|||
auto yRanges = new yRange[circles.length]; |
|||
Fp computeTotalArea(in Fp nSlicesX) { |
|||
Fp total = 0; |
|||
// Adapted from an idea by Cosmologicon. |
|||
foreach (p; cast(int)(xMin * nSlicesX) .. |
|||
cast(int)(xMax * nSlicesX) + 1) { |
|||
immutable Fp x = p / nSlicesX; |
|||
size_t nPairs = 0; |
|||
// Look for the circles intersecting the current |
|||
// vertical secant: |
|||
foreach (ref const Circle c; circles) { |
|||
immutable Fp d = c.r ^^ 2 - (c.x - x) ^^ 2; |
|||
immutable Fp sd = sqrt(d); |
|||
if (d > 0) |
|||
// And keep only the intersection chords. |
|||
yRanges[nPairs++] = yRange(c.y - sd, c.y + sd); |
|||
} |
|||
// Merge the ranges, counting the overlaps only once. |
|||
yRanges[0 .. nPairs].sort(); |
|||
Fp y = -Fp.max; |
|||
foreach (r; yRanges[0 .. nPairs]) |
|||
if (y < r.y1) { |
|||
total += r.y1 - max(y, r.y0); |
|||
y = r.y1; |
|||
} |
|||
} |
|||
return total / nSlicesX; |
|||
} |
|||
// Iterate to reach some precision. |
|||
enum Fp epsilon = 1e-9; |
|||
Fp nSlicesX = 1_000; |
|||
Fp oldArea = -1; |
|||
while (true) { |
|||
immutable Fp newArea = computeTotalArea(nSlicesX); |
|||
if (abs(oldArea - newArea) < epsilon) { |
|||
writeln("N. vertical slices: ", nSlicesX); |
|||
writefln("Approximate area: %.17f", newArea); |
|||
return; |
|||
} |
|||
oldArea = newArea; |
|||
nSlicesX *= 2; |
|||
} |
|||
}</lang> |
|||
{{out}} |
|||
<pre>Input Circles: 25 |
|||
Circles left: 14 |
|||
N. vertical slices: 256000 |
|||
Approximate area: 21.56503660593628004</pre> |
|||
Runtime is about 7.6 seconds. DMD 2.061, -O -release -inline -noboundscheck. |
|||
=={{header|Python}}== |
=={{header|Python}}== |