//================================================== file = integrate.c =====
//=  Monte Carlo simulation to integrate a function f(x) from 0 to x_max    =
//===========================================================================
//=  Notes:                                                                 =
//=    1) This Monte Carlo simulation estimates the area under the curve    =
//=       f(x) from 0 to x_max                                              =
//=    2) The maximum y value (y_max) must be known                         =
//=    3) Must manually set NUM_ITER in the #define                         =
//=    4) Must manually initialize x_max and y_max in main()                =
//=-------------------------------------------------------------------------=
//=  Example execution: (NUM_ITER = 1000000)                                =
//=                                                                         =
//=    *** BEGIN SIMULATION ***                                             =
//=    -------------------------------------------------------------        =
//=    --      *** Results from Monte Carlo integration ***       --        =
//=    -------------------------------------------------------------        =
//=    - Number of iterations = 1000000                                     =
//=    -------------------------------------------------------------        =
//=    - Estimated area (from 0 to 3.000000) = 9.000072                     =
//=    -------------------------------------------------------------        =
//=    *** END SIMULATION ***                                               =
//=-------------------------------------------------------------------------=
//=  Build: bcc32 integrate.c                                               =
//=-------------------------------------------------------------------------=
//=  Execute: integrate.c                                                   =
//=-------------------------------------------------------------------------=
//=  Author: Ken Christensen                                                =
//=          University of South Florida                                    =
//=          WWW: http://www.csee.usf.edu/~christen                         =
//=          Email: christen@csee.usf.edu                                   =
//=-------------------------------------------------------------------------=
//=  History: KJC (06/08/09) - Genesis (from findPi.c)                      =
//===========================================================================
//----- Include files -------------------------------------------------------
#include <stdio.h>               // Needed for printf()
#include <math.h>                // Needed for pow()

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

//----- Function prototypes -------------------------------------------------
double f(double x);              // Function to integrate
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   x_max;                // Integration bound for x for f(x)
  double   y_max;                // Integration bound for y for f(x)
  double   hitCount;             // Count of hits within the quarter circle
  double   area_est;             // Estimated area
  int      i;                    // Loop counter

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

  // Initialize x_max and y_max for f(x)
  x_max = 3;
  y_max = f(x_max);

  // 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 within square X-by-Y)
    x = rand_val(0) * x_max;
    y = rand_val(0) * y_max;

    // Determine if the point is under the f(x) curve
    if (y < f(x)) hitCount++;
  }

  // Estimate the area
  area_est = (hitCount / NUM_ITER) * (x_max * y_max);

  // Output results
  printf("------------------------------------------------------------- \n");
  printf("--      *** Results from Monte Carlo integration ***       -- \n");
  printf("------------------------------------------------------------- \n");
  printf("- Number of iterations = %d \n", NUM_ITER);
  printf("------------------------------------------------------------- \n");
  printf("- Estimated area (from 0 to %f) = %f \n", x_max, area_est);
  printf("------------------------------------------------------------- \n");
  printf("*** END SIMULATION *** \n");

  return(0);
}

//=========================================================================
//= Function to integrate                                                 =
//=  - Input is x, output is f(x)                                         =
//=========================================================================
double f(double x)
{
  double   y;                    // y = f(x) value

  // Function is f(x) = x^2
  y = pow(x, 2);

  // Return y = f(x)
  return(y);
}

//=========================================================================
//= 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);
}