Rosetta Code

Total circles area: Difference between revisions

Page Discussion

Total circles area (view source)

Revision as of 23:48, 25 September 2012

1,403 bytes removed ,  11 years ago
→‎C: Scanline Method: removed separate sorting pass to have fewer things to look at; some comments
m (→‎C: Scanline Method: make it not compile, so users trying it won't be puzzled by incorrect-looking output)
(→‎C: Scanline Method: removed separate sorting pass to have fewer things to look at; some comments)
Line 167:
This version performs about 5 million scanlines in about a second, result should be accurate to maybe 10 decimal points.
<lang c>#include <stdio.h>
#include <mathstring.h>
#include <stdlib.h>
#include <stdboolmath.h>
 
typedef double Fpflt;
typedef struct { Fp x, y, r, r2, y0, y1; } Circle;
flt x, y, r, r2;
flt y0, y1; // extent of circle y+r and y-r
flt x0, x1; // where scanline intersects circle
} circle_t;
#define SZ sizeof(circle_t)
 
Circlecircle_t circles[] = {
{ 1.6417233788, 1.6121789534, 0.0848270516},
#error data snipped, copyfor dataspace; copy from previous C example.
};
{ 0.0152957411, 0.0638919221, 0.9771215985}};
 
inlineflt Fp minmax(Fpflt ax, Fpflt by) { return ax < by ? ay : bx; }
inlineflt Fp maxmin(Fpflt ax, Fpflt by) { return ax > by ? ay : bx; }
inlineflt Fp squaresq(Fpflt x) { return x * x; }
flt cdist(circle_t *c1, circle_t *c2) {
 
return sqrt(sq(c1->x - c2->x) + sq(c1->y - c2->y));
inline Fp dist(const Fp x0, const Fp y0, const Fp x1, const Fp y1) {
return sqrt(square(x0 - x1) + square(y0 - y1));
}
 
inline void swap_c(circle_t *c)
typedef struct { Fp x0, x1; } Sect;
{
circle_t tmp = c[0];
c[0] = c[1], c[1] = tmp;
}
 
flt area(circle_t *circs, int n_circ, flt ymin, flt ymax, flt step)
Fp spans(Circle *c, size_t n_circles,
{
const Fp ymin, const Fp step, const Fp ymax) {
int i, n = n_circ;
// Safe if n_circles is not very large.
Sect sects[n_circles];
 
circle_t *c = malloc(SZ * n);
Fp total = 0;
memcpy(c, circs, SZ * n);
 
while (n--)
// GNU C99 extension.
for (i = 0; i < n; i++)
inline void swap_s(const size_t i, const size_t j) {
if (c[i].y1 < c[i+1].y1) swap_c(c + i);
const Sect ss = sects[i];
sects[i] = sects[j];
sects[j] = ss;
}
 
flt total = 0;
inline void swap_c(const size_t i, const size_t j) {
const Circle cc = c[i];
c[i] = c[j];
c[j] = cc;
}
 
int row = 1 + ceil((ymax - ymin) / step);
for (size_t i = 0; i < n_circles; i++)
while (row--) {
for (size_t j = i + 1; j < n_circles; j++)
flt y = ymin + step * row;
if (c[i].y1 < c[j].y1)
for (n = 0; n < n_circ; n++)
swap_c(i, j);
if (y >= c[n].y1) // rest of circles below scanline, ignore
break;
else if (y > c[n].y0) {
flt dx = sqrt(c[n].r2 - sq(y - c[n].y));
c[n].x0 = c[n].x - dx;
c[n].x1 = c[n].x + dx;
 
// keep circles sorted by left intersection
for (int row = 1 + ceil((ymax - ymin) / step); row; row--) {
for (i = n; i-- && c[i].x0 > c[i+1].x0; swap_c(c + i));
// Could probably do "y -= step", not sure.
const Fp y = ymin + step * row;
 
} else {// remove a circle when scanline has passed it
size_t n = n_circles;
memmove(c + n, c + n + 1, forSZ * (size_t--n_circ i = 0; i <- n)); ) {
n--;
if (y >= c[i].y1)
}
n = i;
else if (y <= c[i].y0) {
const Circle cc = c[i];
for (size_t j = i + 1; j < n_circles; j++)
c[j - 1] = c[j];
n_circles--;
c[n_circles] = cc;
n--;
} else {
const Fp dx = sqrt(c[i].r2 - square(y - c[i].y));
sects[i].x0 = c[i].x - dx;
sects[i].x1 = c[i].x + dx;
i++;
}
}
 
if (!n) continue;
continue;
 
flt right = c->x1;
for (size_t i = n; i--; ) { // Bubble sort. Important.
total += c->x1 - c->x0;
bool swapped = false;
for (size_t j = 0; j < i; j++)
if (sects[j].x0 > sects[j + 1].x0) {
swap_s(j, j + 1);
swap_c(j, j + 1);
swapped = true;
}
if (!swapped)
break;
}
 
for (i = 1; i < n; i++) {
Sect *cur = sects;
if (c[i].x1 <= right) continue;
total += cur->x1 - cur->x0;
total += c[i].x1 - max(c[i].x0, right);
for (size_t i = 1; i < n; i++) {
right = c[i].x1;
if (sects[i].x1 <= cur->x1)
}
continue;
}
total += sects[i].x1 - max(sects[i].x0, cur->x1);
cur = sects + i;
}
}
 
free(c);
return total;
return total * step;
}
 
int main(void) {
{
size_t n_circles = sizeof(circles) / sizeof(Circle);
int n_circ = sizeof(circles) / SZ;
flt ymin = INFINITY, ymax = -INFINITY;
 
circle_t *c1, *c2;
// Remove fully covered circles.
for (Circle *c1 = circles + n_circlesn_circ; c1-- > circles; ) {
for (Circle *c2 = circles + n_circlesn_circ; c2-- > circles; )
// Removethrow c1out if it'scircles inside c2. Helps, notanother important.circle
if (c1 != c2 && cdist(c1, c2) + if (c1->r !<= c2->r) &&{
*c1 = circles[--n_circ];
dist(c1->x, c1->y, c2->x, c2->y) + c1->r <= c2->r) {
break;
n_circles--;
}
*c1 = circles[n_circles];
ymin = min(ymin, c1->y0 = c1->y - c1->r);
break;
ymax = max(ymax, c1->y1 = c1->y + c1->r);
}
c1->r2 = sq(c1->r);
}
 
flt s = 1. / (1 << 20);
Fp miny = INFINITY, maxy = -INFINITY;
flt y0 = floor(ymin / s) * s;
for (size_t i = 0; i < n_circles; i++) {
flt a = area(circles, n_circ, y0, ymax, s);
circles[i].r2 = square(circles[i].r);
int nlines = (ymax - y0) / s;
circles[i].y0 = circles[i].y - circles[i].r;
circles[i].y1 = circles[i].y + circles[i].r;
miny = min(miny, circles[i].y0);
maxy = max(maxy, circles[i].y1);
}
 
// roughly, it cease to make sense if sqrt(nlines) * 1e-14 >> s * s
const size_t s = 1U << 20; // Definitely overkill. Definitely.
printf("area = %.10f\tat %d scanlines\n", a, nlines);
const Fp sp = spans(circles, n_circles,
floor(miny * s) / s, 1.0 / s, maxy);
printf("Area = %.10f at\t%d scanlines.\n",
sp / s, (int)((maxy - floor(miny * s) / s) * s));
 
return 0;
}</lang>
{{out}}
Areaarea = 21.5650366037 at 5637290 scanlines.
 
=={{header|D}}==
Anonymous user
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/219432"