#include <chrono>
#include <ratio>
#include <thread>
#include <iostream>
#include <cmath>

#define EWALD_F   1.12837917
#define EWALD_P   0.3275911
#define A1        0.254829592
#define A2       -0.284496736
#define A3        1.421413741
#define A4       -1.453152027
#define A5        1.061405429

double erfca(const double grij) {
   double expm2 = exp(-grij*grij);
   double t 	= 1.0 / (1.0 + EWALD_P*grij);
   return t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
}


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

  if (argc != 3) {
    std::cerr << "Wrong number of parameters. Synopsis: erfc <<# of evaluations>> <<mode (0 or 1)>>" << std::endl;
    return -1;
  }


  int N 	= atoi(argv[1]);   	// number of evaluations
  int approx    = atoi(argv[2]);	// use approximation (0/1)
  int L 	= 1.;			// interval


  std::cout << "Settings:  " << N  << ((approx) ? " approximate" : " exact") << " evaluations." << std::endl;
  

  double * y = new double[N];		// y_i = erfc(x_i)
  double * x = new double[N];		// x_i

  // create random input values x_1, ..., x_N
  srand(100);
  for (int i=0; i<N; i++) {
    x[i] = (rand()%N + 1) *L/N;
  }

  // measure the time
  auto t1 = std::chrono::high_resolution_clock::now();

  // evaluate function N times
  if (approx) {
    for (int i=0; i<N; i++) 
      y[i] = erfca(x[i]); 
  }
  else {
    for (int i=0; i<N; i++) 
      y[i] = erfc(x[i]); 
  }

  // measure the time again
  auto t2 = std::chrono::high_resolution_clock::now();
 
  // integral duration: requires duration_cast
  auto int_ms = std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1);
 
  // fractional duration: no duration_cast needed
  std::chrono::duration<double, std::milli> fp_ms = t2 - t1;
  std::cout << "It took " << fp_ms.count() << " ms." << std::endl;


  delete [] y;
  delete [] x;
  return 0;
}