/* 
 * Monte Carlo calculation of area with error estimation
 */
#include <stdio.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_statistics.h>
     
double area (long m, gsl_rng *r);

#define POINTS 1024       /* 2^10, initial number of random points to generate */
#define NEXP 64           /* number of experiments for each value of # points */
#define M 12              /* number of different # of points values */

int main(void) {
    double ar[NEXP];      /* calculated areas */
    gsl_rng *r;
    double mean, sd;
    int i, j;
    unsigned long seed = 1UL;
    long p = POINTS;
     
    /* allocate random number generator */
    r = gsl_rng_alloc(gsl_rng_taus2);

    /* seed the random number generator */
    gsl_rng_set(r, seed);

    for (i = 0; i < M; i++) {
        for (j = 0; j < NEXP; j++) {
            /* calculate pi using Monte Carlo algorithm */
            ar[j] = area(p, r);
        }
        mean = gsl_stats_mean(ar, 1, NEXP);
        sd = gsl_stats_sd_m(ar, 1, NEXP, mean);
        printf("%15ld   %10.8f   %10.8f\n", p, mean, sd);
        p *= 2;
    }

    gsl_rng_free(r);

    return(0);
}