#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>

// Conversion from degrees to radians
double degtorad(double degang) {
  return ((M_PI * degang) / 180.0);
}

// Function to be integrated (replace with your actual function)
double func(double x) {
  return tan(x); // Example function
}

// Trapezoidal rule function with OpenMP parallelization
double trap(int n, double h, double *fx) {
  double area = 0.0;

  #pragma omp parallel for reduction(+:area) private(i)
  for (int i = 0; i < n; i++) {
    area += h * func(a + i * h); // Use actual function here
  }

  return area;
}

int main() {
  int rank, size, n_local, n_global = 2; // Initial guess for global n
  double a, b, h, area, local_area;
  double *x, *fx;
  double diff, tol = 1.0e-6;
  MPI_Status status;

  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Comm_size(MPI_COMM_WORLD, &size);

  // Enter a and b in process 0
  if (rank == 0) {
    printf("Enter a and b\n");
    scanf("%lf %lf", &a, &b);
  }

  // Broadcast a and b to all processes
  MPI_Bcast(&a, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&b, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  // Find optimal n in parallel using a doubling strategy
  while (1) {
    // Calculate local sub-interval size on each process
    h = (b - a) / (n_global - 1);
    n_local = n_global / size + (rank < n_global % size); // Account for uneven distribution

    // Allocate memory for x and fx on each process
    x = (double*)malloc(n_local * sizeof(double));
    fx = (double*)malloc(n_local * sizeof(double));

    // Calculate x and fx values (can be parallelized within each process)
    for (int i = 0; i < n_local; i++) {
      x[i] = a + rank * (b - a) / size + i * h;
      fx[i] = func(degtorad(x[i])); // Convert to radians before function call
    }

    // Call trap function (parallelized) to find local area
    local_area = trap(n_local, h, fx);

    // Free memory for x and fx
    free(x);
    free(fx);

    // Gather local areas from all processes and sum on process 0
    MPI_Reduce(&local_area, &area, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);

    if (rank == 0) {
      diff = fabs(area - log(2.0));
      printf("n = %d, Area = %lf, Difference = %lf\n", n_global, area, diff);
      if (diff < tol) {
        break;
      }
    }

    // Double the global number of intervals for the next iteration
    n_global *= 2;

    // Check for potential overflow of n_global
    int ret = MPI_Bcast(&n_global, 1, MPI_INT, 0, MPI_COMM_WORLD);
    if (ret != MPI_SUCCESS) {
      printf("Error: MPI_Bcast failed (%d)\n", ret);
      MPI_Finalize();
      return 1;
    }
  }

  MPI_Finalize();
  return 0;
}