# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <cmath>

//
//  An INCLUDE file is needed here.
//
using namespace std;

int main ( int argc, char *argv[] );
double f ( int dim_num, double x[] );
void r8vec_uniform_01 ( int n, int *seed, double r[] );

//****************************************************************************80

int main ( int argc, char *argv[] )

//****************************************************************************80
//
//  Purpose:
//
//    MAIN is the main program for MC.
//
//  Discussion:
//
//    Things you can vary include:
//
//    * the function to be integrated (see the list below)
//    * the number of samples used, SAMPLE_NUM.
//    * the spatial dimension DIM_NUM (try 1, 2, 4, 8)
//
//    The function to be integrated is.
//
//      F = Product ( 1 <= DIM <= DIM_NUM ) ( 2 * abs ( 2 * X(DIM) - 1 ) ) 
//      Exact integral is 1.
//
{
  int dim_num = 4;
  double q;
  double q_error;
  double q_exact = 1.0;
  int sample;
  int sample_num = 1000;
  int seed;
  double *x;

  x = new double[dim_num];

  cout << "\n";
  cout << "MC\n";
  cout << "  C++ version\n";

  seed = 123456789;

  q = 0.0;

  for ( sample = 1; sample <= sample_num; sample++ )
  {
    r8vec_uniform_01 ( dim_num, &seed, x );

    q = q + f ( dim_num, x );
  }

  q = q / ( double ) ( sample_num );
  q_error = fabs ( q - q_exact );

  cout << "\n";
  cout << "   Samples Dimension      Estimate          Error\n";
  cout << "\n";

  cout << "  " << setw(8) << sample_num
       << "  " << setw(8)  << dim_num
       << "  " << setw(16) << q
       << "  " << setw(16) << q_error << "\n";

  delete [] x;

  return 0;
}
//****************************************************************************80

double f ( int dim_num, double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    F evaluates the function F(X) which we are integrating.
//
//  Modified:
//
//    04 September 2008
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the number of spatial dimensions.
//
//    Input, double X[DIM_NUM], the point at which to evaluate F.
//
//    Output, double F, the value of F(X).
//
{
  int dim;
  double value;

  value = 1.0;
  for ( dim = 0; dim < dim_num; dim++ )
  {
    value = value * 2.0 * fabs ( 2.0 * x[dim] - 1.0 );
  }

  return value;
}
//****************************************************************************80

void r8vec_uniform_01 ( int n, int *seed, double r[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_UNIFORM_01 returns a unit pseudorandom vector.
//
//  Parameters:
//
//    Input, int N, the number of entries in the vector.
//
//    Input/output, int *SEED, a seed for the random number generator.
//
//    Output, double R[N], the vector of pseudorandom values.
//
{
  int i;
  int i4_huge = 2147483647;
  int k;

  if ( *seed == 0 )
  {
    cerr << "\n";
    cerr << "R8VEC_UNIFORM_01 - Fatal error!\n";
    cerr << "  Input value of SEED = 0.\n";
    exit ( 1 );
  }

  for ( i = 0; i < n; i++ )
  {
    k = *seed / 127773;

    *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;

    if ( *seed < 0 )
    {
      *seed = *seed + i4_huge;
    }

    r[i] = ( double ) ( *seed ) * 4.656612875E-10;
  }

  return;
}