#include <iostream.h>
#include <fstream.h>
#include <stdio.h>
#include <stdlib.h>
#include <iomanip.h>
#include <math.h>

// Calculate Monte Carlo integral of f(x) = x^2
// using 4 different sampling strategies
// dump the results in matlab format for plotting

inline double f( double y ) {
  return y*y;
}

int main( int argc, char **argv ) {

  if( argc != 5 ) {
    cerr << "Usage: montecarlo <Nmax> <deltaN> <a> <b>\n";
    exit(1);
  }

  double intvalue = 8.0/3.0;

  double Nmax = atof( argv[1] );
  double deltaN = atof( argv[2] );
  double a = atof( argv[3] );
  double b = atof( argv[4] );
  int i;
  double N = deltaN;

  ofstream un( "uniform.m", ios::out );
  ofstream st( "stratified.m", ios::out );
  ofstream im( "importance.m", ios::out );
  ofstream si( "strimp.m", ios::out );

  if( !un || !st || !im || !si ) {
    cerr << "Error: could not open a file\n";
    exit(1);
  }

  un << "N = [ ";
  st << "N = [ ";
  im << "N = [ ";
  si << "N = [ ";

  // fill in the vector of number of samples
  for( i = (int)N; i <= (int)Nmax; i += (int)deltaN ) {
    un << i << " ...\n";
    st << i << " ...\n";
    im << i << " ...\n";
    si << i << " ...\n";
  }
  un << "];\nS = [ ";
  st << "];\nS = [ ";
  im << "];\nS = [ ";
  si << "];\nS = [ ";

  // fill in the vector of standard deviations
  for( i = (int)N; i <= (int)Nmax; i += (int)deltaN ) {
    int j;
    double yi, yi_1;
    double sum1 = 0.0;
    double sum2 = 0.0;
    double e1, e2, e3, e4;
    un << sqrt( 256.0 / (45.0 * (double)i) ) << " ...\n";
    im << sqrt( 16.0 / (3.0 * (double)i) ) << " ...\n";

    yi = a;
    yi_1 = a + (b-a)/(double)i;
    for( j = 0; j < i; j++ ) {
      // stratified
      e1 = (yi_1 - yi) * ( pow(yi_1,5) - pow(yi,5) ) / 5.0;
      e2 = pow( ( pow(yi_1,3) - pow(yi,3) ), 2 ) / 9.0;
      sum1 += (e1 - e2) / (double)i;

      // stratified imp.
      e3 = pow( yi_1 - yi, 2) * ( pow(yi_1,6) - pow(yi,6) ) /
	( ( yi_1*yi_1 - yi*yi ) * 3.0 );
      e4 = (yi_1 - yi) * ( pow(yi_1,4) - pow(yi,4) ) / ( yi_1*yi_1 - yi*yi );
      e4 = e4 * e4 / 4;
      sum2 += (e3 - e4) / (double)i;
      
      yi = yi_1;
      yi_1 += (b-a)/(double)i;
    }
    st << sqrt( sum1 ) << " ...\n";
    si << sqrt( sum2 ) << " ...\n";
  }
  un << "];\nI = [ ";
  st << "];\nI = [ ";
  im << "];\nI = [ ";
  si << "];\nI = [ ";

  while( N <= Nmax ) {

    // uniform sampling
    double sum = 0.0;
    double p = 1.0 / (b - a);
    double xi, y;
    
    for( i = 0; i < (int)N; i++ ) {
      xi = drand48();
      y = a + (b - a) * xi;
      sum += f(y) / p;
    }
    double integral = sum / N;
    un << setprecision(12) << fabs(integral - intvalue) << " ...\n";

    // stratified sampling
    sum = 0.0;
    double strat_y[ (int)N ];
    double y_i = a;
    double y_i_1 = a + (b-a)/N;
    
    for( i = 0; i < (int)N; i++ ) {
      xi = drand48();

      // fill in the array of sample points between [0,2)
      strat_y[i] = y_i + ( y_i_1 - y_i ) * xi;

      // use the sample points to integrate
      sum += f( strat_y[i] ) / p;

      // set up the next interval/stratum
      y_i += (b-a)/N;
      y_i_1 += (b-a)/N;
    }
    integral = sum / N;
    st << setprecision(12) << fabs(integral - intvalue) << " ...\n";

    // importance sampling
    // pdf = kx  int(pdf) = 1 => k = 1/2
    // P(x) = int( x/2, 0, x ) = x^2/4
    // xi = inverse(P(x))  => x = 2*sqrt(xi)
    sum = 0.0;
    for( i = 0; i < (int)N; i++ ) {
      y = 2.0 * sqrt( drand48() );
      p = y / 2.0;
      sum += f( y ) / p;
    }
    integral = sum / N;
    im << setprecision(12) << fabs(integral - intvalue) << " ...\n";

    // stratified importance
    sum = 0.0;
    for( i = 0; i < (int)N; i++ ) {

      strat_y[i] = 2.0 * sqrt( strat_y[i] / (b-a) );
      p = strat_y[i] / 2.0;

      sum += f( strat_y[i] ) / p;
    }
    integral = sum / N;
    si << setprecision(12) << fabs(integral - intvalue) << " ...\n";

    // ready the next iteration
    N += deltaN;
  }
  un << "];\n";
  st << "];\n";
  im << "];\n";
  si << "];\n";
  
  un.close();
  st.close();
  im.close();
  si.close();

}