Вычисление интеграла геометрическим методом Монте-Карло

#include "stdafx.h"
#include <iostream>
#include <random>
 
int main()
{
    double x, y;
    double Integr;
    double Integr_0 = 0.0;  // начальное значение для сравнения
    double summ = 0.0;      // для подсчета суммы
    int n = 0;              // для подсчета количества точек, попавших в заданную область
    int N = 500;            // первоначальное количество точек для вычисления
 
    std::mt19937 gen{ std::random_device()() };
    // область интегрирования будет квадрат 1x1, 
    // поэтому задаем диапазон генерации случайных чисел от 0 до 1
    std::uniform_real_distribution<double> dist(0.0, 1.0);
 
    do {
        for (int i = 0; i < N; i++) {
            x = dist(gen);
            y = dist(gen);
            // если наша точка попадает в необходимую область, ниже прямой y = 1 - x
            // то вычисляем и накапливаем сумму, а также количество точек, попавших в данную область
            if (y < 1. - x) {
                summ += 1. - x - y;
                n++;
            }
        }
 
        // вычисляем наш интеграл. 0.5 - это площадь интегрирования. 
        Integr = 0.5 * summ / n;
        //меняем местами для сохранения в дальнейших вычислениях
        std::swap(Integr, Integr_0);
        // увеличиваем количество точек в 2 раза
        N *= 2;
        //можно обнулить значения и вычислять по новой, а можно и оставить для накопления результата
//      n = 0;
//      summ = 0.;
    } while (abs(Integr - Integr_0) > 0.0001);
 
    // так как был обмен, то последний результат хранится в Integr_0
    std::cout << Integr_0 << std::endl;
 
    system ("pause");
}