// Monte Carlo simulation (method 2)

#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <getopt.h>
#include <vector>
#include <algorithm>
#include <ctime>

using namespace std;

double X       = -1.0; // if negative, use (Start,End,Step) range
double Start   = -1.0; // 0.01;
double End     = 3.0;
double Step    = 0.01;
int    N       = 10000;
int    J       = 100;
int    C       = 10000;
int    Prec    = 400;
bool   PrnFuncValOnly = false;
bool   Quiet   = false;
bool   Verbose = false;
bool   Debug   = false;
bool   help    = false;

int toint (const char* const s) {
  size_t len = strlen(s);
  char  *ptr;
  long   val;
  val = strtol(s, &ptr, 0);
  if (ptr != (s+len)) {
    fprintf(stderr,"Not an integer argument.");
    exit(-1);
  }
  return val;
}
double todouble(const char*const s) {
  size_t len = strlen(s);
  char  *ptr;
  double val;
  val = strtod(s, &ptr);
  if (ptr != (s+len)) {
    fprintf(stderr,"Not an floating point argument.");
    exit(-1);
  }
  return val;
}

// Calculate Vn = \sqrt(N) (\max(I_n/N - U_n) - \min ((I_n -1) - U_n)
double calcVn(int N) {
  double sqrtN = sqrtl((double)N);
  //  Generate a vector of random uniform variables
  vector<double> uvec;
  for (int n = 0; n < N; ++n) {
    double u = (double)rand()/RAND_MAX;
    uvec.push_back(u);
  }
  //  Sort in ascending order
  sort(uvec.begin(),uvec.end());
  // Find the max and min
  double max = 0.0;
  double min = 1.0;
  for (int n = 0; n < N; ++n) {
    double t1  = (double)(n+1) / N - uvec[n];
    double t2  = (double)(n)   / N - uvec[n];
    if(t1 > max) max = t1;
    if(t2 < min) min = t2;
  }
  // M_3
  double vn = sqrtN * (max - min);
  return vn;
}

// Calculate approximation to F_3(x)
double fc(double x, int C, int J, int N) {
  vector<double> msvec;
  for(int c = 0; c < C; ++c) {
    double produ = 1.0;
    double ms = 0.0;
    for (int j = 0; j < J; ++j) {
      double vn = calcVn(N);
      double u = (double)rand()/RAND_MAX;
      double l = produ * u;
      produ *= (1.0-u);
      double m = sqrtl(l) * vn;
      if(m > ms) ms = m;
    }
    msvec.push_back(ms);
  }
  sort(msvec.begin(),msvec.end());

  int c = 0;
  if (Start > 0) {
    double t = Start;
    for(c = 0; (c < C) && (t <= End); ) {
      if(msvec[c] > t && t <= End) {
        printf("%.3f\t%e\n",t,(double)c/C);
        t += Step;
      } else {
        ++c;
      }
    }
    if(t <= End) {
      for(double i = t; i <= End; i += Step) {
        printf("%.3f\t%lf\n",i,1.0);
      }
    }
  } else {
    for(c = 0; (c < C) && (msvec[c] <= x); ++c) {
      ;
    }
  }
  return ((double)c) / C;
}

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

  srand(time(NULL));

  while (1) {
    static struct option long_options[] = {
      {"N",       required_argument,  0,  'n' },
      {"X",       required_argument,  0,  'x' },
      {"Start",   required_argument,  0,  '0' },
      {"End",     required_argument,  0,  'f' },
      {"Step",    required_argument,  0,  'h' },
      {"C",       required_argument,  0,  'c' },
      {"J",       required_argument,  0,  'j' },
      {"Prec",    required_argument,  0,  'p' },
      {"prnFuncValOnly",   optional_argument,  0,  'y' },
      {"!prnFuncValOnly",   optional_argument,  0,  'z' },
      {"quiet",   optional_argument,  0,  'q' },
      {"verbose", optional_argument,  0,  'v' },
      {"debug",   optional_argument,  0,  'd' },
      {0,0,0,0} /* This is a filler for -1 */
    };
    int option_index = 0;
    int c = getopt_long (argc, argv, "0:f:h:n:x:c:j:p:vdqyz", long_options, &option_index);
    if (c == -1) break;
    switch (c) {
    case 'n':
      N   = toint(optarg);
      break;
    case 'x':
      X   = todouble(optarg);
      break;
    case '0':
      Start = todouble(optarg);
      break;
    case 'f':
      End = todouble(optarg);
      break;
    case 'h':
      Step = todouble(optarg);
      break;
    case 'c':
      C   = toint(optarg);
      break;
    case 'j':
      J = toint(optarg);
      break;
    case 'p':
      Prec = toint(optarg);
      break;
    case 'y':
      PrnFuncValOnly = true;
    case 'z':
      PrnFuncValOnly = false;
    case 'q' :
      Quiet = true;
      break;
    case 'v' :
      Verbose = true;
      break;
    case 'd' :
      Debug   = true;
      break;
    default:
      fprintf(stderr,"Unknown option '%c'\n",c);
      return 1;
      break;
    }
  }

  if(X >= 0) {
    if(!PrnFuncValOnly) {
      printf("%.2f\t",X);
    }
    printf("%lf",fc(X,C,J,N));
    printf("\n");
  } else if(Start > 0 && (End > Start && Step > 0.0) ) {
    fc(1.0,C,J,N);
  } else {
    fprintf(stderr,"Need to supply -x param (X positive) or -0 param (Start positivet)\n");
  }

  return 0;
}