#include <cmath>
#include <iostream>

#define LOGGING
#ifdef LOGGING
#define Logg(message) std::clog << message << std::endl;
#else // LOGGING
#define Logg(message)
#endif /* NO_ARG_CHECK */


template<class T, class functor>
T integrate_trapez_simple(functor& F, double a, double b, int N)
{// Trapezoid rule for fixed number of function evaluations N at points a, a+(b-a)/(N-1),...
  double dh = (b-a)/(N-1);
  T sum = 0.5*F(a);
  double x = a+dh;
  for (int i=1; i<N-1; i++, x += dh) sum += F(x);
  sum += 0.5*F(x);
  return sum*dh;
}

template <class T, class functor>
T integrate4(functor& F, double a, double b, int Np)
{ // Forth order integration routine with equidistant mesh
  static double coeff[4] = {17./48., 59./48.,43./48.,49./48.};
  double dh = (b-a)/(Np-1), x = a;
  T sum = 0;
  for (int i=0; i<4; i++, x+=dh) sum += coeff[i]*F(x);
  for (int i=4; i<Np-4; i++, x+=dh) sum += F(x);
  for (int i=Np-4; i<Np; i++, x+=dh) sum += coeff[Np-1-i]*F(x);
  return sum*dh;
}

template<class T, class functor>
T integrate_trapez(functor& F, double a, double b, double precision=1e-5, int max_recursion=22)
{
  T olds = 0, s = 0.5*(F(a)+F(b))*(b-a); // lowest order approximation
  int m=1;                          // number of subdivisions
  for (int i=0; i<max_recursion; i++, m*=2) { // m==2^(i-2)
    double dh = (b-a)/m;         // This is the distance between points to be added at this stage
    double x = a + 0.5*dh;       // first point added at
    T sum = 0;
    for (int j=0; j<m; j++, x+=dh) sum += F(x); // Evaluates sum S(n)
    s = 0.5*s + 0.5*dh*sum;     // The best approximation so far
    if (i > 4  &&  fabs(s-olds) <= precision*fabs(olds)){ //Avoid spurious early convergence.
      Logg("Trapezoid: The integral step used =" << 0.5*dh << " with " << m << " divisions");
      return s;
    }
    olds=s;
  }
  std::cerr<<"Too many steps in routine integrate_trapez"<<std::endl;
  return s; // The best approximation
}

template<class T, class functor>
T integrate_simps(functor& F, double a, double b, double precision=1e-5, int max_recursion=22)
{
  T oldsm = 0, sm = 0; // simpson approximation
  T olds = 0, s = 0.5*(F(a)+F(b))*(b-a); // trapezoid approximation, lowest order approximation
  int m=1;                          // number of subdivisions
  for (int i=0; i<max_recursion; i++, m*=2) { // m==2^(i-2)
    double dh = (b-a)/m;         // This is the distance between points to be added at this stage
    double x = a + 0.5*dh;       // first point added at
    T sum = 0;
    for (int j=0; j<m; j++, x+=dh) sum += F(x); 
    s = 0.5*s + 0.5*dh*sum;     // The best approximation so far
    sm = (4*s-olds)/3.;         // This is the first difference with trapezoid
    if (i > 4  &&  fabs(sm-oldsm) <= precision*fabs(oldsm)){ //Avoid spurious early convergence.
      Logg("Simpson:   The integral step used ="<<0.5*dh<<" with "<<m<<" divisions");
      return sm; 
    }
    olds=s;
    oldsm=sm;
  }
  std::cerr<<"Too many steps in routine integrate_simps"<<std::endl;
  return sm; // The best approximation
}


// Integral is 0.5*(sqrt(pi)-2*Gamma(0.75))=-0.3391897770124197
double F0(double x)
{ return (1-2*sqrt(fabs(x)))*exp(-x*x);}
double F1(double x)
{ return cos(x);}
double F2(double x)
{return 1;}

using namespace std;

int main()
{
  double IF0_exact = -0.3391897770124197*2;
  cout.precision(16);
  cout<<" Founction F0=(1-2*sqrt(fabs(x)))*exp(-x*x)"<<endl;
  double v = integrate_trapez_simple<double>(F0, -5, 5, 200);
  cout<<"Precision of simple trapezoid rule with 200 points "<<v-IF0_exact<<endl;
  double v1 = integrate_trapez<double>(F0, -5, 5,1e-4);
  cout<<"Precision of adaptive trapezoid rule "<<v1-IF0_exact<<endl;
  double v2 = integrate_simps<double>(F0, -5, 5,1e-4);
  cout<<"Precision of adaptive simpon's rule "<<v2-IF0_exact<<endl;
  
  cout<<" Founction F1=cos(x)"<<endl;
  double p = integrate_trapez_simple<double>(F1, 0.0, 51*M_PI/2., 200);
  cout<<"Precision of simple trapezoid rule with 200 points "<<p+1<<endl;
  double p1 = integrate_trapez<double>(F1, 0.0, 51*M_PI/2., 1e-4);
  cout<<"Precision of adaptive trapezoid rule "<<p1+1<<endl;
  double p2 = integrate_simps<double>(F1, 0.0, 51*M_PI/2., 1e-4);
  cout<<"Precision of adaptive simpon's rule "<<p2+1<<endl;
  double p3 = integrate4<double>(F1, 0.0, 51*M_PI/2., 200);
  cout<<"Precision of forth order routine with 200 points "<<p3+1<<endl;
  
  return 0;
}