#include <stdio.h>
#include <math.h>
#include "mpi.h"
#define NP 10000
#define f(x,y) exp(-x*x-0.5*y*y)*(10.0*x*x - x)*(y*y-y)
int main(int argc, char **argv)
{

/* 
   Exercise 5.2: compute integral of f() over 
   star-shaped region: central square is 1 x 1, each of 
   N,S,E,W points are 2 x 1 triangles of the same area 
   as the central square 
                         W. Petersen, 26 May, 2002, ETHZ 
*/
 
   int i,j,size,ip,master,rank;
   int ierr;
   static float seed;
   int seeds[]={331,557,907,1103,1303};  /* seeds for each cpu */
   float x,y,t,tot;
   float buff[1];      /* buff gets partial sums from slaves */
   float sum[4];       /* partial sums from "points" */
   float ggl(float*);  /* random number generator */
   MPI_Status stat;

   MPI_Init(&argc,&argv);               /* initialize MPI */
   MPI_Comm_size(MPI_COMM_WORLD, &size); /* get number of cpus(5) */
   MPI_Comm_rank(MPI_COMM_WORLD, &rank); /* get rank of this one */
   master = 0;                           /* rank of master */

   ip = rank;
   if(ip==master){
/* master computes integral over center square */
     tot = 0.0; 
     seed    = (float) seeds[ip];
     for(i=0;i<NP;i++){
        x = ggl(&seed)-0.5;
        y = ggl(&seed)-0.5;
        tot += f(x,y); 
     } 
     tot *= 1.0/((float) NP);  /* int f(x,y) dx dy over center square */
   } else { 
/* rank != 0 computes E,N,W,S points of star */
     seed  = seeds[ip];
     seed  = 0;
     sum[ip-1] = 0.0;
     for(i=0;i<NP;i++){
        x = ggl(&seed);
        y = ggl(&seed) - 0.5;
        if(y > (0.5-0.25*x)){
           x = 2.0 - x;
           y = -(y-0.5);
        }
        if(y < (-0.5+0.25*x)){
           x = 2.0 - x;
           y = -(y+0.5);
        }
        x += 0.5;
        if(ip==2){
           t = x;
           x = y;
           y = -t;
        } else if(ip==3){
           x = -x; 
           y = -y;
        } else if(ip==4){
           t = x;
           x = -y; 
           y = t;
        }
        sum[ip-1] += f(x,y);
     }
     sum[ip-1] *= 1.0/((float) NP);
     buff[0] = sum[ip-1];
     ierr=MPI_Send(buff,1,MPI_FLOAT,0,0,MPI_COMM_WORLD);
   }
   if(ip==master){   /* get partial sums from other cpus */
      for(i=1;i<5;i++){ 
         ierr = MPI_Recv(buff,1,MPI_FLOAT,
                MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,
                &stat);
         tot += buff[0];
      }
      printf(" Integral = %e\n",tot);
   }
   MPI_Finalize();
}
#include <math.h>
float ggl(float *ds)
{
/* generate u(0,1) distributed random numbers.
   Seed ds must be saved between calls. ggl is
   essentially the same as the IMSL routine RNUM.

   W. Petersen and M. Troyer, 24 Oct. 2002, ETHZ */

   double t,d2=0.2147483647e10;
   t   = (double) *ds;
   t   = fmod(0.16807e5*t,d2);
   *ds = (float) t;
   return((float)((t-1.0e0)/(d2-1.0e0)));
}