//===================================================== file = findPi.c =====
//=  Classic Monte Carlo simulation to estimate value of pi                 =
//===========================================================================
//=  Notes:                                                                 =
//=    1) This Monte Carlo simulation estimate pi by "throwing darts" at    =
//=       a quarter circle of radius = 1 contained within a square of       =
//=       side lenfth = 1. The area of the quarter circle is pi/4 and the   =
//=       are of the square is 1.  The percentage of darts that fall        =
//=       inside the circle is used to estimate pi.                         =
//=    2) The value of pi is 3.1415926535 to 10 places                      =
//=    3) Must manually set NUM_ITER in the #define                         =
//=-------------------------------------------------------------------------=
//=  Example execution: (NUM_ITER = 100000)                                 =
//=                                                                         =
//=   *** BEGIN SIMULATION ***                                              =
//=   -------------------------------------------------------------         =
//=   --     *** Results from Monte Carlo pi estimator ***       --         =
//=   -------------------------------------------------------------         =
//=   - Number of iterations = 100000                                       =
//=   -------------------------------------------------------------         =
//=   - Estimated pi = 3.141120                                             =
//=   -------------------------------------------------------------         =
//=   *** END SIMULATION ***                                                =
//=-------------------------------------------------------------------------=
//=  Build: bcc32 findPi.c                                                  =
//=-------------------------------------------------------------------------=
//=  Execute: findPi                                                        =
//=-------------------------------------------------------------------------=
//=  Author: Ken Christensen                                                =
//=          University of South Florida                                    =
//=          WWW: http://www.csee.usf.edu/~christen                         =
//=          Email: christen@csee.usf.edu                                   =
//=-------------------------------------------------------------------------=
//=  History: KJC (03/29/09) - Genesis                                      =
//=           KJC (06/08/09) - Some minor clean-up                          =
//===========================================================================
//----- Include files -------------------------------------------------------
#include <stdio.h>               // Needed for printf()
#include <math.h>                // Needed for pow() and sqrt()

//----- Defines -------------------------------------------------------------
#define  NUM_ITER       100000   // Number of iterations to run

//----- Function prototypes -------------------------------------------------
double rand_val(int seed);       // RNG for uniform(0.0, 1.0) from Jain

//===========================================================================
//=  Main program                                                           =
//===========================================================================
int main()
{
  double   x, y;                 // X and Y values
  double   len;                  // Length from (0,0) of X,Y point
  double   hitCount;             // Count of hits within the quarter circle
  double   pi_est;               // Estimated pi
  int      i;                    // Loop counter

  // Output banner
  printf("*** BEGIN SIMULATION *** \n");

  // Seed the RNG and initialize hitCount to zero
  rand_val(1);
  hitCount = 0.0;

  // Do for NUM_ITER iterations
  for (i=0; i<NUM_ITER; i++)
  {
    // Throw the dart (dart lands at an X,Y coordinate)
    x = rand_val(0);
    y = rand_val(0);

    // Compute the length to X,Y using Pythagorian
    len = sqrt(pow(x, 2.0) + pow(y, 2.0));

    // Test if dart is inside quarter circle and increment hitCount
    if (len <= 1.0) hitCount++;
  }

  // Estimate pi
  pi_est = 4.0 * (hitCount / NUM_ITER);

  // Output results
  printf("------------------------------------------------------------- \n");
  printf("--     *** Results from Monte Carlo pi estimator ***       -- \n");
  printf("------------------------------------------------------------- \n");
  printf("- Number of iterations = %d \n", NUM_ITER);
  printf("------------------------------------------------------------- \n");
  printf("- Estimated pi = %f \n", pi_est);
  printf("------------------------------------------------------------- \n");
  printf("*** END SIMULATION *** \n");

  return(0);
}

//=========================================================================
//= Multiplicative LCG for generating uniform(0.0, 1.0) random numbers    =
//=   - x_n = 7^5*x_(n-1)mod(2^31 - 1)                                    =
//=   - With x seeded to 1 the 10000th x value should be 1043618065       =
//=   - From R. Jain, "The Art of Computer Systems Performance Analysis," =
//=     John Wiley & Sons, 1991. (Page 443, Figure 26.2)                  =
//=-----------------------------------------------------------------------=
//= Input: Seed value of positive integer (0 for no seed)                 =
//= Output: Returns unif RV if input is 0                                 =
//=========================================================================
double rand_val(int seed)
{
  const long  a =      16807;    // Multiplier
  const long  m = 2147483647;    // Modulus
  const long  q =     127773;    // m div a
  const long  r =       2836;    // m mod a
  static long x;                 // Random int value that is seeded
  long        x_div_q;           // x divided by q
  long        x_mod_q;           // x modulo q
  long        x_new;             // New x value

  // Seed the RNG if seed is not equal to zero
  if (seed != 0)
    x = seed;

  // RNG using integer arithmetic
  x_div_q = x / q;
  x_mod_q = x % q;
  x_new = (a * x_mod_q) - (r * x_div_q);
  if (x_new > 0)
    x = x_new;
  else
    x = x_new + m;

  // Return a random value between 0.0 and 1.0
  return((double) x / m);
}