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

double f1(double x){
  double r=1.0;
  double return_val;
  double arg = r*r-x*x;
  if (arg<0.) arg=0.;
  return_val = sqrt(arg);
  return return_val;
}


double f2(double x){
  return x*sin(10*x);
}


double trapez(double a, double b, int n){
  int i;
  double sum=0.;
  double h=(b-a)/(float)n;
  double xa=a;
  double xb=a+h;

  for (i=0; i<n; i++){
    sum += (f2(xa)+f2(xb))/2. * h;
    xa = xb;
    xb += h;
  }
  return sum;
}


void main(int argc, char **argv){

  int myrank,nproz,in;
  double integral,gesamt;
  double a,b;
  double my_a,my_b;
  double result;
  int source;
  int n_ges;
  int my_n;

  MPI_Status status;

  MPI_Init(&argc,&argv);

  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);

  MPI_Comm_size(MPI_COMM_WORLD, &nproz);

  n_ges=400000000;
  a = 0.;
  b = 2*3.141592654;

  my_a = a + (b-a)/(float)nproz*myrank;
  my_b = a + (b-a)/(float)nproz*(myrank+1);

  my_n = n_ges/nproz;

  integral = trapez(my_a,my_b,my_n);
  printf("%d: integral = %f\n",myrank,integral);

  MPI_Send(&integral,1,MPI_DOUBLE,0,0,MPI_COMM_WORLD);
  if ( myrank==0 ){
    for (source=0; source<nproz; source++){
      MPI_Recv(&result,1,MPI_DOUBLE,source,0,MPI_COMM_WORLD,&status);
      gesamt += result;
    }
  }

  if ( myrank==0 )
    printf("Summe aus Teilintegralen bei Prozeß %d: %f\n",myrank,gesamt);

  fflush(stdout);
  MPI_Barrier(MPI_COMM_WORLD);

  MPI_Allreduce(&integral,&gesamt,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);

  if ( myrank==0 )
    printf("Summe aus Teilintegralen bei Prozeß %d: %f\n",myrank,gesamt);

  fflush(stdout);
  MPI_Barrier(MPI_COMM_WORLD);

  MPI_Finalize();
}