/* 
 * Monte Carlo calculation Pi
 */
#include <stdio.h>

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

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

int main(void) {
    double pi[NEXP];      /* calculated values of pi */
    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 */
            pi[j] = mc_pi(p, r);
        }
        mean = gsl_stats_mean(pi, 1, NEXP);
        sd = gsl_stats_sd_m(pi, 1, NEXP, mean);
        printf("%15ld   %10.8f   %10.8f\n", p, mean, sd);
        p *= 2;
    }

    gsl_rng_free(r);

    return(0);
}