PICOLL

/*
 *      picoll.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 DEFAULT_ANZAHL 10000L       /* Defaultwert für Anzahl */


void 
main(int argc, char *argv[])
{
  
  int    myid;
  int    nproc;                         /* number of processes */
  long   total;                         /* total number of points */
  long   mytotal;                       /* number of points per process */
  long   myhits;                        /* number of hits per process */
  long   totalhits = 0;                 /* total number of hits */
  double pi;                            /* result */
  long   calc(long total);              /* calculation */
  
  /* start MPI */
  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &nproc);
  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
  
  if  (myid == 0) {
    /* compute number of points per process */
    if (argc != 2) {
      total = DEFAULT_ANZAHL;
    } else {
      total = 1000 * atol(argv[1]);
    }
    mytotal = (long) ceil(total/nproc);
    total = nproc * mytotal;              /* correct total */
  }       
  
  /* send work to all processes */
  MPI_Bcast(&mytotal, 1, MPI_LONG, 0, MPI_COMM_WORLD);
  
  /*  compute partial results */
  myhits = calc(mytotal);
  
  /* combine partial results */
  MPI_Reduce(&myhits, &totalhits, 1, MPI_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
  
  if (myid == 0) {
    /* compute final result */
    pi = totalhits/(double)total*4;
    printf("\nPI = %lf\n", pi);
  }
  
  /* leave MPI */
  MPI_Finalize();
}

#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); 
}