#include <iostream>
#include <random>
#include <chrono>

float lowerX = 0;
float upperX = 20;

double randomFloat()
{
    return (double)(rand()) / (double)(RAND_MAX);
}

class Point {
private:
    double x{};
    double y{};

public:
    Point(double x, double y): x(x), y(y) {}

    [[nodiscard]] double get_x() const {
        return this->x;
    }

    [[nodiscard]] double get_y() const {
        return this->y;
    }

    Point() {};
};

class LinearBounds {
private:
    float lower;
    float higher;

public:
    LinearBounds(float lower, float higher): lower(lower), higher(higher) {}
    [[nodiscard]] double length() const {
        return std::abs(higher - lower);
    }

    [[nodiscard]] double getRandomValue() const {
        return this->lower + (randomFloat() * this->length());
    }
    [[nodiscard]] float get_lower() const {
        return lower;
    }
    [[nodiscard]] float get_upper() const {
        return higher;
    }
};


class Bounds {
protected:
    LinearBounds x;
    LinearBounds y;
public:
    Bounds(LinearBounds x, LinearBounds y): x(x), y(y) {}
    double area() {
        return x.length() * y.length();
    }
    Point getRandomPoint() {
        return {x.getRandomValue(), y.getRandomValue()};
    }
    LinearBounds get_x() {
        return x;
    }

};




double function(double x) {
    return pow(x,4)-4*pow(x,3)+18*pow(x,2)-12*x-69;
}
double indefiniteIntegral(double x) {
    return (1.0/5) * pow(x, 5) - pow(x, 4) + 6*pow(x, 3) - 6 * pow(x, 2) -69*x;
}

bool getIsInside(Point p) {
    double y_0 = function(p.get_x());
    return std::abs(y_0) > std::abs(p.get_y()) && y_0* p.get_y() >= 0;
}

double getRealIntegral(LinearBounds bounds) {
    return indefiniteIntegral(bounds.get_upper())- indefiniteIntegral(bounds.get_lower());
}
float getUpperLimit() {
    float upperLimit = 0;
    for(float i = lowerX; i< upperX; i+= 0.025) {
        if(function((double)i) > upperLimit) {
            upperLimit = (float) (function(i) * 1.025);
        }
    }
    return upperLimit;
}
float getLowerLimit() {
    float lowerLimit = 0;
    for(float i = lowerX; i< upperX; i+= 0.025F) {
        if(function(i) < lowerLimit) {
            lowerLimit = (float) (function(i) * 1.025);
        }
    }
    return lowerLimit;
}

int main() {


    auto startTime = std::chrono::high_resolution_clock::now();
    double pointsInside = 0;
    double samples = 10000000;
    Bounds bounds = Bounds(LinearBounds(lowerX, upperX), LinearBounds(getLowerLimit(), getUpperLimit()));
    Point toTestPoint;
    for (int i = 0; i< samples; i++) {
        toTestPoint = bounds.getRandomPoint();
        if (getIsInside(toTestPoint)) {
            pointsInside++;
        }
    }

    std::cout << std::fixed;
    std::cout.precision(5);
    double integral = (pointsInside / samples) * bounds.area();
    std::cout << "The approximated Integral of the function is: " << integral << "\n";

    double real_value = getRealIntegral(bounds.get_x());
    std::cout << "The real Integral of the function is: " << real_value << "\n";

    double error = std::abs(real_value-integral);
    std::cout << "That's an error of " << error <<" or " << (error/real_value)*100 << "%\n";
    std::chrono::duration<double> duration = std::chrono::duration_cast<std::chrono::duration<double>>(std::chrono::high_resolution_clock::now()-startTime);
    std::cout << "And the whole thing took " << duration.count() << " Seconds for ";
    std::cout.precision(0);
    std::cout << samples << " samples";
    return 0;
}