#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)
{
  uniform u;
  exponential e;
  normal n;
  REAL a=0.0;
  REAL b=1.0;
  INTEGER g=100;
  switch (argc) {
    case 4: g=atoi(argp[3]);
    case 3: a=atof(argp[1]); b=atof(argp[2]); break;
    default: error(string("Usage: ")+argp[0]+" a b [g] "); break; 
  }
  density_base* f;
  f = &u;
  cout << prob(f,a,b,g) << '\n';
  f = &e;
  cout << prob(f,a,b,g) << '\n';
  f = &n;
  cout << prob(f,a,b,g) << '\n';
  return 0;
}