#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

double f(double x)
{
    return 1.0 / (1.0 + x + pow(x,2));
}

double getRandomValue(double min, double max)
{
    double f = (double)rand() / RAND_MAX;
    return min + f * (max - min);
}


double monteCarlo(double a, double b, int N)
{
    std::vector<double> xi(N), yi(N);

    double step = (b - a) / N;

    double MMax = 0;
    for (double i = a; i <= b; i += step)
    {
        if (f(i) > MMax)
            MMax = f(i);
    }

    int Ncounter = 0;

    double S = (b - a) * MMax;

    for (size_t i = 0; i < xi.size(); ++i)
    {
        xi[i] = a + getRandomValue(0, 1) * (b - a);
        yi[i] = f(xi[i]);

        double Y = getRandomValue(0, 1) * MMax;
        if ( Y <= yi[i] )
            Ncounter++;
    }

    return S * Ncounter / static_cast<double>(N);
}
int main()
{
    cout << monteCarlo(0.0, 1.0, 5000) << endl;
    return 0;
}