#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();
  }
};

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)
{
  REAL a=0.0;
  REAL b=1.0;
  INTEGER g=100;
  uniform u;
  exponential e;
  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; 
  }
  cout << prob(u,a,b,g) << '\n';
  cout << prob(e,a,b,g) << '\n';
  return 0;
}