/*
 * Calculate f(x) = \int_0^x (sin(1/t))^2 dt using MC integration
 */
#include <stdio.h>
#include <math.h>

#include <gsl/gsl_rng.h>

#include "integral-mc.h"
     
int main() {
    int i, nsteps = 1000, np = 500000;
    unsigned long int seed = 1L;
    double xmin = 0., xmax = 2., f = 0., xleft, xright, step;
    gsl_rng *r;

    step = (xmax - xmin)/nsteps;
    xleft = 0.;
    xright = step;
    printf("%10.6f    %10.6f\n", xleft, f);

    r = gsl_rng_alloc(gsl_rng_taus2);
    gsl_rng_set(r, seed);

    for (i = 0; i < nsteps; i++) {
        f = f + mc_integral(r, xleft, step, np);

        printf("%10.6f    %10.6f\n", xright, f);
        xleft = xright;
        xright += step;
    }

    gsl_rng_free (r);

    return 0;
}

static double mc_integral(gsl_rng *r, double xleft, double step, int np) {
    double f, x, y;
    int i, count = 0;

    for (i = 0; i < np; i++) {
        x = xleft + step*gsl_rng_uniform_pos(r);
        y = gsl_rng_uniform(r);
        if (pow(sin(1./x), 2) >= y) {
            count += 1;
        } 
    }
   
    f = ((double) count)/((double) np)*step;

    return(f);
}