#include <iostream>
#include <chrono>
#include <cmath>
 
 
 
#define PI 3.14159265358979323846  /* pi */
 
double integral(double boundaryLeft, double boundaryRight, long divisions);
double f(double x);
 
int main(int argc, char **argv)
{
    int number_of_processes = 1;
    int my_rank = 0;
    const int repetitions = 1;
 
    double result = 0.0;
    double local_integral = 0.0;
    const auto start_chrono = std::chrono::high_resolution_clock::now();
 
    const long number_of_intervals = 5040 * 100000;
 
    const double global_a = 0.0;
    const double global_b = 1.0;
 
    for (int round = 0; round < repetitions; round++)
    {
 
        const long local_number_of_intervals = number_of_intervals / number_of_processes;
        const double local_interval_width = (global_b - global_a) / (double)number_of_processes;
        const double local_a = global_a + my_rank * local_interval_width;
        const double local_b = local_a + local_interval_width;
 
        local_integral = integral(local_a, local_b, local_number_of_intervals);
 
        // send / receive
        
         result += local_integral;
 
    }
 
    result = result / repetitions;
 
    const auto finish_chrono = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff = finish_chrono - start_chrono;
    if (my_rank == 0)
    {
        std::cout << "Runtime: " << diff.count() << " s. " << std::endl;
        printf("Result: %0.20f \n", result * 4.);
        printf("Error: %0.20f \n", fabs(result * 4 - PI));
    }
 
    return 0;
}
 
double integral(double boundaryLeft, double boundaryRight, long divisions)
{
    double tmp = 0.0;
    double locIntervalWidth = (boundaryRight - boundaryLeft) / (double)divisions;
    for (int i = 0; i < divisions; i++)
    {
        tmp += f(boundaryLeft + locIntervalWidth * ((double)i + 0.5)) * locIntervalWidth;
    }
    return tmp;
}
 
double f(double x)
{
    return (std::sqrt(1. - x * x));
}