I was reading a bit about random numbers and remembered that I wrote a simple Monte Carlo integrator in C++ a few years ago. I took a few minutes to cleanup and comment the code, which is presented below.

Monte Carlo integration is simple, but surprisingly powerful:

$I={\int}_{a}^{b}f\left(x\right)dx$

$I\approx \sum _{i=1}^{n}\frac{f\left({x}_{i}\right)}{p\left({x}_{i}\right)}$

**x**_{i} is a random value within the range [a,b]

**p(x)** represents the distribution of random values, for a uniform distribution:

$p\left(x\right)=\frac{1}{b\u2013a}$
This presentation by Fabio Pellacini provides a lot more details.

The test code in the main() method computes the integral of **sin**^{2}(x) in the interval **[3,5]**.

`#include <iostream>`

#include <cmath>

#include <ctime>

using namespace std;

// Functor base class for encapsulating 1-dimensional function to be integrated

class Function1d

{

public:

virtual double operator()(double x)

{

return x;

}

};

// Functor for sine squared function

class SineSquared : public Function1d

{

public:

double operator()(double x)

{

return pow(sin(x), 2);

}

};

// Monte Carlo integrator class declaration

class MonteCarloIntegrator

{

public:

double Run(int numSamples, Function1d& func, double intervalMin, double intervalMax);

};

// Monte Carlo integrator implementation

// ::Run() method implementation

double MonteCarloIntegrator::Run(int numSamples, Function1d& func, double intervalMin, double intervalMax)

{

double sum = 0;

double div = 1.0 / (double)RAND_MAX;

double intervalScale = 1.0 / (intervalMax-intervalMin);

for(int i=0; i<numSamples; i++)

{

double rnd1 = intervalMin + ( ((double)rand() * div) * (intervalMax-intervalMin) );

sum += (func)(rnd1) / intervalScale;

}

return (1.0/(double)numSamples) * sum;

}

// main() function with test code

void main()

{

MonteCarloIntegrator integrator;

double output = integrator.Run(5000000, SineSquared(), 3, 5);

srand(time(NULL));

std::cout << output << endl;

getchar();

}