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