/* 
 * 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;
    unsigned long seed = 1UL;
    double mean, sd;
    int i, j;
    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);
        fflush (stdout);
        p *= 2;
    }

    gsl_rng_free (r);

    return 0;
}