/*==============================================================
    Numerical integration of f(x) on [a,b]
    method: Monte-Carlo method
    written by: Alex Godunov 
----------------------------------------------------------------
input:
    f   - a single argument real function (supplied by the user)
    a,b - the two end-points of the interval of integration
    n   - number of random points for xi
output:
    r      - result of integration
    errest - estimated error
Comments:
    be sure that following headers are included
    #include <cstdlib>
    #include <ctime>     
================================================================*/

    double mc_int1(double(*f)(double), double& errest, double a, double b, int n)
{
    double r, x, u, f2s, fs;
/* variables fs and f2s are used to estimate an error of integration */    
    
    srand(time(NULL));               /* initial seed value (use system time) */

    fs  = 0.0;
    f2s = 0.0;

    for (int i = 1; i <= n; i=i+1)
    {
        u = 1.0*rand()/(RAND_MAX+1.0); /* random number between 0.0 and 1.0 */
        x = a + (b-a)*u;
        fs = fs + f(x);
        f2s= f2s+ f(x)*f(x);
    }
    r = fs*(b-a)/n;
/* evaluate integration error */
    fs = fs/n;
    f2s = f2s/n;
    errest = (b-a)*sqrt((f2s - fs*fs)/n);
    return r;
}