#include <cstdlib>
#include <iostream>
#include <ctime>


using namespace std;

inline double frand() { return (double)rand() / (double)RAND_MAX; }

class Function {
public:
  virtual double eval(double x) = 0;
  virtual ~Function() {}
  
  // Finds minimum.
  virtual double min(double a, double b) {
    double delta = (b-a)/0.01;
    double min = eval(a);
    for(double i = a+delta; i <= b; i += delta) {
       double y = eval(i);
       if(y < min) min = y;        
    }     
    return min;
  }
  
  // Finds maximum.
  virtual double max(double a, double b) {
    double delta = (b-a)*0.01;
    double max = eval(a);
    for(double i = a+delta; i <= b; i += delta) {
       double y = eval(i);
       if(y > max) max = y;        
    }     
    return max;          
  }
  
  // Generic, using MC.
  virtual double integrate(unsigned int n, double a, double b) {
    double _min = min(a,b), _max = max(a,b);
    
    // We estimate error 10% on both min and max.
    _min = _min > 0.0 ? _min*0.9 : _min*1.1;
    _max = _max > 0.0 ? _max*1.1 : _max*0.9;
    
    unsigned int P = 0;
    double bma = b-a, mmm = _max - _min;
    for(unsigned int i = 0; i < n; i++) {
       double x = frand()*bma + a;  
       double y = frand()*mmm + _min;
       
       if(eval(x) > y) P++;
    }
    
    return bma*mmm*(double)P/(double)n + _min*bma;
  } 
};


class Polinom : public Function {
public:
  unsigned int n;
  double* data;
  
  // Constructors.
  Polinom(unsigned int _n) : n(_n) { data = new double[n+1]; }
  Polinom(unsigned int _n, double* d) : n(_n), data(d) {}
  virtual ~Polinom() { delete [] data; }
  
  // Evaluates at x.
  virtual double eval(double x) {
    double res = data[0];
    double _x = x;
    for(unsigned int i = 0; i < n; i++) {
      res += data[i+1]*_x;
      _x*=x; 
    }
    return res;
  }    
};

int main(int argc, char *argv[])
{
  srand(time(NULL));
  Polinom* p = new Polinom(5);
  p->data[0] = 10;
  p->data[1] = 7;
  p->data[2] = 5;
  p->data[3] = 2;
  p->data[4] = 3;
  p->data[5] = 4;
  
  cout << "Integral (0, 2) je " << p->integrate(1000000, 0, 1) << "\n";
  system("PAUSE");
  return 0;
}