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

# include "mpi.h"

double f ( int m, double x[] );

int main ( int argc, char *argv[] )
{
  double error;
  double f1;
  double f1_part;
  double f2;
  double f2_part;
  int i;
  int id;
  int j;
  int m = 2;
  int n;
  int np;
  int p;
  unsigned int seed;
  double stdev;
  double sterr;
  double value;
  double x[2];
//
//  Initialize MPI.
//  Get this processor's ID.
//  Get the number of processes.
//
  MPI_Init ( &argc, &argv );
  MPI_Comm_rank ( MPI_COMM_WORLD, &id );
  MPI_Comm_size ( MPI_COMM_WORLD, &p );
//
//  There are P processes available, and P-1 of them are workers.
//  There are N function evaluations to carry out.
//  The master process computes the value NP = N / ( P - 1 ), and then
//  adjusts N so it equals NP * ( P - 1 ).
//
  if ( id == 0 )
  {
    n = 1000000;
    np = n / ( p - 1 );
    n = ( p - 1 ) * np;
  }
//
//  The Broadcast command says that the master (process 0) is sending 
//  an integer "NP" to all the processes, and that each process should
//  receive this value.
//
  MPI_Bcast ( &np, 1, MPI_INT, 0, MPI_COMM_WORLD );
//
//  At this point, each process knows the value NP.
//  The workers should work, and the master can be idle.
//
  f1_part = 0.0;
  f2_part = 0.0;

  if ( 0 < id )
  {
    seed = 12345 + id;
    srand ( seed );

    for ( i = 0; i < np; i++)
    {
      for ( j = 0; j < m; j++ )
      {
        x[j] = ( double ) rand ( ) / ( double ) RAND_MAX;
      }
      value  = f ( m, x );
      f1_part = f1_part + value;
      f2_part = f2_part + value * value;
    }
  }
//
//  The master zeros out F1 and F2 so they can be used to gather the results.
//
  if ( id == 0 )
  {
    f1 = 0.0;
    f2 = 0.0;
  }
//
//  The partial results are gathered into F1 and F2.
//
  MPI_Reduce ( &f1_part, &f1, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
  MPI_Reduce ( &f2_part, &f2, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
//
//  Print the results.
//
  if ( id == 0 )
  {
    f1 = f1 / ( double ) n;
    f2 = f2 / ( double ) ( n - 1 );
    stdev = sqrt ( f2 - f1 * f1 );
    sterr = stdev / sqrt ( ( double ) n );
    error = fabs ( f1 - 1.0 );
    printf ( "\n" );
    printf ( "      N          F           stdev   sterr      err\n" );
    printf ( "\n" );
    printf ( "%7d  %.15g  %6.4f  %6.4f  %8g\n", n, f1, stdev, sterr, error );
  }
//
//  Terminate MPI.
//
  MPI_Finalize ( );

  return 0;
}
double f ( int m, double x[] )
{
  int i;
  double value;

  value = 1.0;
  for ( i = 0; i < m; i++ )
  {
    value = value * fabs ( 4.0 * x[i] - 2.0 );
  }

  return value;
}