MonteCarlo

#include <iostream>
#include <cmath>
#include <vector>
#include <random>
#include <chrono>
#define pb push_bac
#define MONTECARLOITERATIONS 1000000
using namespace std;

namespace OboznyiMaksymArea {

const double PI = acos(-1);

struct Point {
    int x, y;
};

int R, n;
vector<Point> polygon;

double s(Point a, Point b, Point c) {
    return (a.x - b.x) * (a.y + b.y) + (b.x - c.x) * (b.y + c.y) + (c.x - a.x) * (c.y + a.y);
}

bool is_inside(Point p) {
    for (int i = 0; i + 1 < polygon.size(); i++)
        if (s(polygon[i], p, polygon[i + 1]) < 0)
            return 0;
    return 1;
}

int f0 = 1, f1 = 1;

double MonteCarloOurArea() {

    if (MONTECARLOITERATIONS * n > 100'000'000) {
        cout << "The polygon is too big for Monte-Carlo method\n";
        return 0;
    }
    int good = 0;
    for (int k = 0; k < MONTECARLOITERATIONS; k++) {
        Point cur;
        cur.x = f0;
        f1 = f1 + f0;
        f0 = f1 - f0;
        f1 %= (2 * R + 1);
        f0 %= (2 * R + 1);

        cur.y = f0;
        f1 = f1 + f0;
        f0 = f1 - f0;
        f1 %= (2 * R + 1);
        f0 %= (2 * R + 1);


        if (is_inside(cur))
            good++;
    }
    return good * 1.0 / MONTECARLOITERATIONS * (4 * R * R);

}

double MonteCarloStdArea() {
    if (MONTECARLOITERATIONS * n > 100'000'000) {
        cout << "The polygon is too big for Monte-Carlo method\n";
        return 0;
    }
    mt19937 gen;
    gen.seed(time(0));
    int good = 0;
    for (int k = 0; k < MONTECARLOITERATIONS; k++) {
        Point cur;
        cur.x = gen() % (2 * R + 1) - R;
        cur.y = gen() % (2 * R + 1) - R;
        if (is_inside(cur))
            good++;
    }
    return good * 1.0 / MONTECARLOITERATIONS * (4 * R * R);

}

double GeometryArea() {
    return n * 0.5 * R * R * sin(2 * PI / n);
}

double CellMethodArea() {
    if (R * R * n > 20'000'000) {
        cout << "The polygon is too big for cell method\n";
        return 0;
    }
    int full = 0, part = 0;
    for (int i = -R; i < R; i++)
        for (int j = -R; j < R; j++) {
            if (is_inside({i, j}) && is_inside({i, j + 1}) && is_inside({i + 1, j}) && is_inside({i + 1, j + 1}))
                full++;
            else if (is_inside({i, j}) || is_inside({i, j + 1}) || is_inside({i + 1, j}) || is_inside({i + 1, j + 1}))
                part++;
        }
    return full + part * 0.5;
}
void run() {
   cout << "Enter the number of points(from 3 to 100): ";
    while (1) {
        cin >> n;
        if (n < 3 || n > 100)
            cout << "Incorrect value\n";
        else
            break;
    }

    cout << "Enter the radius(from 1 to 1000): ";
    while (1) {
        cin >> R;
        if (R < 1 || R > 1000)
            cout << "Incorrect value\n";
        else
            break;
    }

    double alpha = 2 * PI / n;
    for (int i = 0; i < n; i++)
        polygon.push_back({R * sin(i * alpha), R * cos(i * alpha)});
    polygon.push_back(polygon[0]);

    double ideal = GeometryArea();
    cout << "Geometry: " << ideal << endl;
    double area = MonteCarloOurArea();
    cout << "Our generator: " << area <<  ' ' << (area - ideal) / ideal << endl;
    area = MonteCarloStdArea();
    cout << "Mt19937 generator: " << area <<  ' ' << (area - ideal) / ideal << endl;
    area = CellMethodArea();
    cout << "Cell: " << area <<  ' ' << (area - ideal) / ideal << endl;

}
}

int main()
{

    OboznyiMaksymArea::run();
    return 0;
}