#include <stdio.h>
#include <math.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_sf_gamma.h>

double vol (int dim, gsl_rng *r, long ntot);

double vol (int dim, gsl_rng *r, long ntot)
{
    long i, nin = 0;
    double r2, u;
    int j;

    for (i = 0; i < ntot; i++)
    {
        r2 = 0.;
        for (j = 0; j < dim; j++)
        {
            u = gsl_rng_uniform (r);
            r2 += u*u;
        }
        if (r2 < 1.)
        {
            nin += 1;
        }
    }
    return ( (((double) nin)/ntot)*pow(2., (double) dim));
}

int main (void)
{
    int n, dimmin = 1, dimmax = 12;
    unsigned long seed = 10LU;
    long ntot = 100000000L;

    gsl_rng * r = gsl_rng_alloc (gsl_rng_taus);
    gsl_rng_set (r, seed);

    for (n = dimmin; n <= dimmax; n++)
    {
        double theor = pow(M_PI, n*.5)/gsl_sf_gamma (.5*n + 1.);
        printf ("%5d    %.5f    %5f\n", n, vol(n, r, ntot), theor);
    }

    gsl_rng_free (r);

    return 0;
}