PI

/*      pi.c
 *
 *      Computation of pi by a "montecarlo" method:
 *      Every process computes a number of points in the unit square
 *      and returns the number of hits in the quarter circle
 *      pi is  #hits/#total * 4
 *      
 *      usage:
 *              pi [total]
 */

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

#define TAG_WORK         1
#define TAG_RESULT       2
#define DEFAULT_ANZAHL   10000L       /* Defaultwert für Anzahl */

void main(int argc, char *argv[])
{
  int    myid;
  void   master(int argc, char *argv[]);
  void   slave(void);
   
  /* start MPI */
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
  
  if  (myid == 0) {
    master(argc, argv);
  } else {
    slave();
  }

  /* leave MPI */
  MPI_Finalize();
}

void master(int argc, char *argv[])
{
  long        total;                    /* total number of points */
  long        mytotal;                  /* number of points per process */
  long        myhits;                   /* number of hits per process */
  long        totalhits;                /* total number of hits */
  double      pi;                       /* result */
  int         nproc;                    /* number of processes */
  MPI_Status  status;
  int         i;
  long        calc(long total);         /* calculation */
  
  /* compute number of points per process */
  if (argc != 2) {
    total = DEFAULT_ANZAHL;
  } else {
    total = 1000 * atol(argv[1]);
  }
  MPI_Comm_size(MPI_COMM_WORLD, &nproc);
  mytotal = (long) ceil(total/nproc);
  total = nproc * mytotal;              /* correct total */
  
  /* send work to all processes */
  for (i=1; i<nproc; i++) {
    MPI_Send(&mytotal, 1, MPI_LONG, i, TAG_WORK, MPI_COMM_WORLD);
  }
  
  /*  compute partial results */
  myhits = calc(mytotal);
  
  /* combine partial results */
  totalhits = myhits;
  for (i = 1; i < nproc; i++) {
    MPI_Recv(&myhits, 1, MPI_LONG, MPI_ANY_SOURCE, TAG_RESULT,
	     MPI_COMM_WORLD, &status);
    totalhits += myhits;
  }  
  
  /* compute final result */
  pi = totalhits/(double)total*4;
  printf("\nPI = %lf\n", pi);
}

void slave(void)
{
  long   mytotal;                       /* number of points per process */
  long   myhits;                        /* number of hits per process */
  MPI_Status      status;
  long   calc(long total);              /* calculation */
  
  
  /* get work from master */
  MPI_Recv(&mytotal, 1, MPI_LONG, 0, TAG_WORK, MPI_COMM_WORLD, &status);
  
  /*  compute partial results */
  myhits = calc(mytotal);
  
  /* send result to master */
  MPI_Send(&myhits, 1, MPI_LONG, 0, TAG_RESULT, MPI_COMM_WORLD);
}

#include <sys/types.h>

long calc(long total)
{
  /* 
   * compute total random points in the unit square
   * and return the number of hits in the sector (x*x + y*y < 1)
   */
  
  double  x, y;                     /* random coordinates */
  long    hits = 0;                 /* number of hits */
  long     i;
  
  /* initialize random generator (otherwise all return the same result!) */
  srand(getpid());
  
  for(i=0; i<total; i++) {
    x = ((double) rand())/RAND_MAX;
    y = ((double) rand())/RAND_MAX;
      
    if ( x*x + y*y <= 1.0 ) {
      hits++;
    }
  }
  
  return(hits); 
}