Pi Solution

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

// function prototypes

double myRand(void);

//  main program start

int main(int argc, char *argv[]) {
//
int   	    i, j;
int   	    ir;
int		    nbatch = 10;
int		    batch_size = 100;
int		    nseed = 1;
unsigned long int   batch_success;
double x, y, r;
double area_batch;
double area_accum, area_accum2, area;
double var, std;
//
printf("\n ***** Program to Calulate PI by Monte Carlo Methods *****\n\n");
//
if (argc < 4) {
printf(" ***** ERROR: Not enough command line arguments! \n");
printf("\n      Please enter in the following order:\n");
printf("      Random Number Seed\n");
printf("      Number of Batches\n");
printf("      Histories per Batch\n");
exit(0);
}
//
printf("Program Name is    : %s\n", argv[0]);
printf("Random Number seed : %s\n", argv[1]);
printf("Number of Batches  : %s\n", argv[2]);
printf("Histories per Batch: %s\n", argv[3]);
//
//  Re-assign command line input
//
nseed      = atoi(argv[1]);
nbatch     = atoi(argv[2]);
batch_size = atoi(argv[3]);
//
//  Initialize the random number generator
//
srand(nseed);
//
//  Initialize the variance accumulators for x and x-squared
//
area_accum  = 0.0;
area_accum2 = 0.0;
//
for (j=0; j<nbatch; j++)  {
batch_success = 0;
for (i=0; i<batch_size; i++)  {
x = myRand();
y = myRand();
r = x * x + y * y;
if (r < 1.0)  {
batch_success++;
}
}
area_batch = ((double) batch_success) / ((double) batch_size);
//
//  Increment the variance accumulators
//
area_accum  += area_batch;
area_accum2 += area_batch * area_batch;
}
//
//  Normalize results
//
area = area_accum / nbatch;
//
//  Calculate variance and standard deviation
//
var  = (area_accum2 / nbatch) - (area * area);
var  = var / (nbatch - 1.0);
std  = sqrt(var);
//
//  Print results
//
printf("\n *** RESULTS ***\n");
printf("Area: %f +/- %f (%f)\n", area, std, std/area);
printf("Pi  : %f +/- %f (%f)\n", area*4.0, std*4.0, std/area);
//
}			// end main function
//
double myRand()  {
//
//  Return a random number between 0 and 1
//
double randx;
randx = ((double) rand()) / (RAND_MAX + 1.0);
return randx;
}

```