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

Peter Junglas 11.5.2000