//  Monte Carlo Integration
//
//  User input:
//      f - function to be integrated ( f(x) >= 0 )
//      [a,b] - integration interval  ( a < b )
//      max_f - any upper bound on the values of f
//      N - number of random points to be generated
//
#include <iostream>
#include <cstdlib>   // For srand(), my_rand(), EXIT_SUCCESS
#include <ctime>   // For time()
#include <cmath>
using namespace std;

double f( double x )
{
    return( sin( x ) );
}

double my_rand( double a = 0.0, double b = 1.0 )
{
    double x = double( rand() ) / RAND_MAX;
    return( a + ( b - a ) * x );
}

int main()
{
    // User input /////////////////////////////
    double a = 0.0, b = M_PI, max_f = 1.0;
    int N = 1000;
    ///////////////////////////////////////////
    double x, y, approx;
    int P = 0;
    
    srand( time(NULL) );   // Set the seed for the RNG
    
    for ( int i = 1; i <=N; i++ ) {
        x = my_rand( a, b );
        y = my_rand( 0.0, max_f );
        if ( y <= f(x) ) P++;
    }
    approx = P * max_f * ( b - a ) / N;
    cout << "The approximate value is " << approx << endl;
    
    cout << "\n\n"; system( "PAUSE" );
    return( EXIT_SUCCESS );
}