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

/*
  You need an INCLUDE file here.
*/

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

/******************************************************************************/

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

/******************************************************************************/
/*
  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 D_NUM (try 1, 2, 4, 8)

    The function to be integrated is

      F(X) = 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 = ( double * ) malloc ( dim_num * sizeof ( double ) );

  printf ( "\n" );
  printf ( "MC\n" );
  printf ( "  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 );

  printf ( "\n" );
  printf ( "   Samples Dimension      Estimate          Error\n" );
  printf ( "\n" );
  printf ( "  %8d  %8d  %16f  %16e\n", sample_num, dim_num, q, q_error );

  free ( x );

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

double f ( int dim_num, double x[] )

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

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

/******************************************************************************/
/*
  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 )
  {
    printf ( "\n" );
    printf ( "R8VEC_UNIFORM_01 - Fatal error!\n" );
    printf ( "  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;
}