Total circles area: Difference between revisions

Content added Content deleted
(Used INFINITY as in the second C entry, for uniformity + tags)
(Small changes in second C entry: C99, VLAs, less commas, 4 spaces indented, variables scope narrowing)
Line 169: Line 169:
#include <math.h>
#include <math.h>
#include <stdlib.h>
#include <stdlib.h>
#include <stdbool.h>


typedef double F;
typedef double Fp;
typedef struct { F x, y, r, r2, y0, y1; } circle;
typedef struct { Fp x, y, r, r2, y0, y1; } Circle;
circle circles[] = {
{ 1.6417233788, 1.6121789534, 0.0848270516},...
/* data snipped for space, copy from previous example */
};


Circle circles[] = {
inline F min(F a, F b) { return a < b ? a : b; }
{ 1.6417233788, 1.6121789534, 0.0848270516},
inline F max(F a, F b) { return a > b ? a : b; }
// ... data snipped for space, copy from previous example.
inline F sq(F x) { return x * x; }
{ 0.0152957411, 0.0638919221, 0.9771215985}};
inline F dist(F x, F y, F x1, F y1) { return sqrt(sq(x - x1) + sq(y - y1)); }


inline Fp min(Fp a, Fp b) { return a < b ? a : b; }
typedef struct { F x0, x1; } sect_t;
inline Fp max(Fp a, Fp b) { return a > b ? a : b; }
F spans(circle *c, int n_circs, F ymin, F step, F ymax)
inline Fp square(Fp x) { return x * x; }
{
F y, total = 0;
int i, j, n, row = 1 + ceil((ymax - ymin) / step);


inline Fp dist(const Fp x0, const Fp y0, const Fp x1, const Fp y1) {
sect_t *sect = malloc(n_circs * sizeof(sect_t));
return sqrt(square(x0 - x1) + square(y0 - y1));
inline void swap_s(int i, int j) {
}
sect_t ss;
ss = sect[i], sect[i] = sect[j], sect[j] = ss;
}


typedef struct { Fp x0, x1; } Sect;
inline void swap_c(int i, int j) {
circle cc;
cc = c[i], c[i] = c[j], c[j] = cc;
}


Fp spans(Circle *c, size_t n_circles,
for (i = 0; i < n_circs; i++)
const Fp ymin, const Fp step, const Fp ymax) {
for (j = i + 1; j < n_circs; j++)
// Safe if n_circles is not very large.
if (c[i].y1 < c[j].y1) swap_c(i, j);
Sect sects[n_circles];


Fp total = 0;
while (row--) {
y = ymin + step * row; // could probably do "y -= step"; not sure


// GNU C99 extension.
for (n = n_circs, i = 0; i < n; ) {
inline void swap_s(const int i, const int j) {
if (y >= c[i].y1) n = i;
const Sect ss = sects[i];
else if (y <= c[i].y0) {
sects[i] = sects[j];
circle cc = c[i];
for (j = i + 1; j < n_circs; j++) c[j-1] = c[j];
sects[j] = ss;
}
c[--n_circs] = cc;
--n;
} else {
F dx = sqrt(c[i].r2 - sq(y - c[i].y));
sect[i].x0 = c[i].x - dx;
sect[i].x1 = c[i].x + dx;
i++;
}
}


inline void swap_c(const int i, const int j) {
if (!n) continue;
const Circle cc = c[i];
c[i] = c[j];
c[j] = cc;
}


for (i = n; i--; ) { // Bubble sort. Important.
for (size_t i = 0; i < n_circles; i++)
for (size_t j = i + 1; j < n_circles; j++)
int swapped = 0;
for (j = 0; j < i; j++)
if (c[i].y1 < c[j].y1)
swap_c(i, j);
if (sect[j].x0 > sect[j+1].x0) {
swap_s(j, j + 1);
swap_c(j, j + 1);
swapped = 1;
}
if (!swapped) break;
}


for (int row = 1 + ceil((ymax - ymin) / step); row; row--) {
sect_t *cur = sect;
// Could probably do "y -= step", not sure.
total += cur->x1 - cur->x0;
const Fp y = ymin + step * row;
for (i = 1; i < n; i++) {
if (sect[i].x1 <= cur->x1) continue;
total += sect[i].x1 - max(sect[i].x0, cur->x1);
cur = sect + i;
}
}
free(sect);
return total;
}


size_t n = n_circles;
int main(void)
for (size_t i = 0; i < n; ) {
{
if (y >= c[i].y1)
int i, n_circs = sizeof(circles) / sizeof(circle);
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;

for (size_t i = n; i--; ) { // Bubble sort. Important.
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;
}

Sect *cur = sects;
total += cur->x1 - cur->x0;
for (size_t i = 1; i < n; i++) {
if (sects[i].x1 <= cur->x1)
continue;
total += sects[i].x1 - max(sects[i].x0, cur->x1);
cur = sects + i;
}
}

return total;
}


int main() {
F miny = INFINITY, maxy = -INFINITY;
size_t n_circles = sizeof(circles) / sizeof(Circle);
for (i = 0; i < n_circs; i++) {
circles[i].r2 = sq(circles[i].r);
circles[i].y0 = circles[i].y - circles[i].r;
circles[i].y1 = circles[i].y + circles[i].r;
miny = min(circles[i].y0, miny);
maxy = max(circles[i].y1, maxy);
}


// Remove fully covered circles.
circle *c1, *c2;
for (c1 = circles + n_circs; c1-- > circles; )
for (Circle *c1 = circles + n_circles; c1-- > circles; )
for (c2 = circles + n_circs; c2-- > circles; )
for (Circle *c2 = circles + n_circles; c2-- > circles; )
// Remove c1 if it's inside c2. Helps, not important.
// Remove c1 if it's inside c2. Helps, not important.
if (c1 != c2 && dist(c1->x, c1->y, c2->x, c2->y) + c1->r <= c2->r) {
if (c1 != c2 &&
dist(c1->x, c1->y, c2->x, c2->y) + c1->r <= c2->r) {
*c1 = circles[--n_circs];
n_circles--;
break;
*c1 = circles[n_circles];
}
break;
}


Fp miny = INFINITY, maxy = -INFINITY;
int s = 1 << 20; // Definitely overkill. Definitely.
for (size_t i = 0; i < n_circles; i++) {
circles[i].r2 = square(circles[i].r);
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);
}


const size_t s = 1U << 20; // Definitely overkill. Definitely.
printf("area = %.10f at\t%d scanlines\n",
const Fp sp = spans(circles, n_circles,
spans(circles, n_circs, floor(miny * s) / s, 1./s, maxy) / s,
(int)((maxy - floor(miny * s) / s) * s));
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;
return 0;
}</lang>
}</lang>
{{out}}
{{out}}
Retrieved from "https://rosettacode.org/wiki/Total_circles_area"