Talk:Total circles area

Monte Carlo

A semi-important point and semi-nitpicking for doing Monte Carlo: you need to estimate the error in a credible way given the number of samples, and should not output more significant digits than the confidence allows. The current C code sample size of 100M can't possibly allow 8 decimal point precision. --Ledrug 04:29, 16 September 2012 (UTC)

Right, in the C/Python outputs only the first few digits are correct. Feel free to add a bit of error estimation code.

This is basically what I used on reddit, but I don't feel like changing the example on the task page because MC really isn't the right tool for this task.

<lang c>#include <stdio.h>

1. include <stdlib.h>
2. include <tgmath.h>
3. include <time.h>

typedef float flt; flt data[][3] = { { 1.6417233788, 1.6121789534, 0.0848270516}, ... };

1. define N sizeof(data)/sizeof(data[0])

int main(void) { flt x0 = 1/0., x1 = -1/0., y0 = 1/0., y1 = -1/0.;

inline flt min(flt x, flt y) { return x < y ? x : y; } inline flt max(flt x, flt y) { return x > y ? x : y; } inline flt sq(flt x) { return x * x; } inline flt rndf(flt x, flt d) { return x + d * rand() / RAND_MAX; }

inline int covered(flt x, flt y) { for (int i = 0; i < N; i++) if (sq(x - data[i][0]) + sq(y - data[i][1]) < data[i][2]) return 1; return 0; }

for (int i = 0; i < N; i++) { flt *p = data[i]; x0 = min(x0, p[0] - p[2]); x1 = max(x1, p[0] + p[2]); y0 = min(y0, p[1] - p[2]); y1 = max(y1, p[1] + p[2]);

p[2] *= p[2]; }

x1 -= x0, y1 -= y0;

flt total = x1 * y1; unsigned long to_try = 0x10000, tries = 0, hit = 0; srand(time(0));

while (1) { hit += covered(rndf(x0, x1), rndf(y0, y1));

if (++tries == to_try) { flt area = total * hit / tries; flt r = (flt) hit / tries; flt s = area * sqrt(r * (1 - r) / tries); printf("%.4f +/- %.4f (%lu samples)\n", area, s, tries); // stop at 3 sigmas if (s * 3 <= 1e-3) break; to_try *= 2; } } return 0; }</lang> --Ledrug 22:09, 16 September 2012 (UTC)

A bad solution algorithm is quite useful for didactic purposes (but I don't know how much Rosettacode cares for didactic).