#include "libscl.h"
using namespace scl;
using namespace std;

class density_base {
public:
  virtual den_val operator() (REAL) = 0;
  virtual ~density_base() { }
};

class uniform : public density_base {
public:
  den_val operator() (REAL x) 
  {
    if (0.0<=x && x<=1.0) return den_val(true,0.0);
    return den_val();
  }
};

class exponential : public density_base {
public:
  den_val operator() (REAL x) 
  {
    if (0.0<=x) return den_val(true,-x);
    return den_val();
  }
};

class normal : public density_base {
public:
  den_val operator() (REAL x) 
  {
    const REAL pi = 4.0*atan(1.0);
    const REAL lr2pi = 0.5*log(2.0*pi);
    return den_val(true,-0.5*x*x - lr2pi);
  }
};

REAL prob(density_base* f, REAL a, REAL b, INTEGER n)
{ //Compute probability with trapezoid rule
  REAL sum = 0.0;
  REAL x = a;
  REAL inc = (b - a)/REAL(n);
  if ((*f)(a).positive) sum += exp((*f)(a).log_den)/2.0;
  for (INTEGER i=1; i<n; ++i) { 
    x += inc;
    if ((*f)(x).positive) sum += exp((*f)(x).log_den);
  }
  if ((*f)(b).positive) sum += exp((*f)(b).log_den)/2.0;
  return sum*inc;
}

int main(int argc, char** argp, char** envp)
{
  char d = 'n';
  REAL a=0.0;
  REAL b=1.0;
  INTEGER g=100;
  switch (argc) {
    case 5: g=atoi(argp[4]);
    case 4: d=argp[1][0], a=atof(argp[2]); b=atof(argp[3]); break;
    default: error(string("Usage: ")+argp[0]+" d a b [g] "); break; 
  }
  density_base* f = 0;
  switch(d) {
    case 'e': f = new exponential; break;
    case 'u': f = new uniform; break; 
    case 'n': f = new normal; break; 
    default: error("Error, d must be e or n or u"); 
  }
  cout << prob(f,a,b,g) << '\n';
  delete f;  // Never forget to delete the pointer
  return 0;
}