/*
 Integration of a function using following approximations:
 # Trapezoid approximation
 # Simpson's rule
 # Monte Carlo method
 # 8-panel newton-cotes rule
 time of integration is also calculated and printed
*/
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <cstdlib>
#include <ctime>
//#include "trapezoid.cpp"
//#include "simpson.cpp"
//#include "quanc8.cpp"
//#include "mc_int1.cpp"
using namespace std;

double trapezoid(double(*f)(double), double a, double b, int n);
double simpson(double(*f)(double), double a, double b, int n);
double mc_int1(double(*f)(double), double& errest, double a, double b, int n);
void quanc8(double(*fun)(double), double a, double b,
            double abserr, double relerr,
            double& result, double& errest, int& nofun,double& flag);


double f(double); 


int main()
{
    double a, b;
    double trap, simp, mc;
    double nc8, errest, flag;
    int i, n, nofun;
    int ntimes;
    const double pi = 3.1415926;
    const double abserr=0.0, relerr=1.0e-12;  //see quanc8.cpp
    double tail, errmc;    
       
    cout.precision(6);
    cout.setf(ios::fixed | ios::showpoint);

// current local time
    time_t rawtime;
    time ( &rawtime );
    cout << " The current local time is:" << endl <<" "<< ctime (&rawtime) << endl;

// current time in seconds (begin calculations)
    time_t seconds_i;
    seconds_i = time (NULL);

/* initial information */    
    a = 0.0;              // left endpoint
    b = pi;               // right endpoint
    n = 2;                // initial number of intervals
    ntimes = 16;          // number of interval doublings with nmax=2^ntimes

/* "tail part for the 0->infinity integral */
//    tail = 4.0/sqrt(b);  // it works ONLY for the last integral in the list
    tail = 0.0;


    cout << " Intervals  "<<"Trapez.   "<<"Simpson   "<< "MC      " 
         << "  MC error " << endl;

/* step 2: integration with trapezoid and simpson */
    for (i=0; i <=ntimes; i=i+1)
    {
        trap = trapezoid(f, a, b, n);
        simp =   simpson(f, a, b, n);
        mc =     mc_int1(f, errmc, a, b, n);
        cout << setw(10)  << n  << setw(10) << trap +tail  << setw(10) << simp+tail 
             << setw(10)  << mc << setw(10) << errmc << endl;
        n = n*2;
    }

/* step 3: 8-panel newton-cotes rule */
    quanc8(f, a, b, abserr, relerr, nc8, errest, nofun, flag);

    cout << endl     << "    result from quanc8 " << endl;
    cout << setw(12) << "nofun = "    << setw(10) << nofun << endl;
    cout << setw(12) << "integral = " << setw(10) << nc8+tail   << endl;
//    cout << setw(12) << "est.error = "<< setw(10) << errest << endl;
//    cout << setw(12) << "flag = "     << setw(10) << flag   << endl;    

// current time in seconds (end of calculations)
    time_t seconds_f;
    seconds_f = time (NULL);
    cout << endl << "total elapsed time = " << seconds_f - seconds_i << " seconds" << endl << endl;
  
    return 0;
}
/*
 *  Function for integration
*/
    double f(double x)
{
    const double pi = 3.1415926;
    double y;
    y = sin(x);
//    y = cos(x);
//    y = 2.0*sqrt(x)/((x+1.0)*(x+1.0));
//    y = 0.2/(pow((x-6.0),2.0) + 0.02);
//    y = x*cos(10.0*pow(x,2.0))/(pow(x,2.0) + 1.0);
//    y = sqrt(x)/(x*x + 1.0);
//    y = sin(pi*x*x)/((x-pi)*(x-pi)+1);
//    y = 1.0/(1.0-0.998*x*x);
//    y = x*sin(30.*x)*cos(50.*x);
//    y = x/(exp(x)-1.0);
    return y;
}

    double trapezoid(double(*f)(double), double a, double b, int n)
/*
    Integration by trapezoid rule of f(x) on [a,b]
    method:     trapezoid rule
    written by: Alex Godunov
-------------------------------------------------
    f - Function to integrate (supplied by a user)
    a - Lower limit of integration
    b - Upper limit of integration
    r - Result of integration (out)
    n - number of intervals
=================================================
*/
{
    double r, dx, x;
    r = 0.0;
    dx = (b-a)/static_cast<float>(n);
    for (int i = 1; i <= n-1; i=i+1)
    {
        x = a + static_cast<float>(i)*dx;
        r = r + f(x);
    }
    r = (r + (f(a)+f(b)/2.0))*dx;
    return r;
}



double simpson(double(*f)(double), double a, double b, int n)
/*
    Integration by integration by Simpson rule of f(x) on [a,b]
    method:     Simpson rule
    written by: Alex Godunov
-------------------------------------------------
    f - Function to integrate (supplied by a user)
    a - Lower limit of integration
    b - Upper limit of integration
    s - Result of integration (out)
    n - number of intervals
=================================================
*/
{
    double s, dx, x;
// if n is odd we add +1 interval to make it even
    if((n/2)*2 != n) {n=n+1;}
    s = 0.0;
    dx = (b-a)/static_cast<float>(n);
    for ( int i=1; i<=n-1; i=i+2)
    {
        x = a+static_cast<float>(i)*dx;
        s = s + 4.0*f(x) + 2.0*f(x+dx);
    }
    s = (s + f(a)+f(b))*dx/3.0;
    return s;
}



void quanc8(double(*fun)(double), double a, double b,
            double abserr, double relerr,
            double& result, double& errest, int& nofun,double& flag)
/*
   estimate the integral of fun(x) from a to b to a user provided tolerance.
   an automatic adaptive routine based on the 8-panel newton-cotes rule.

   input ...
   fun     the name of the integrand function subprogram fun(x).
   a       the lower limit of integration.
   b       the upper limit of integration.(b may be less than a.)
   relerr  a relative error tolerance. (should be non-negative)
   abserr  an absolute error tolerance. (should be non-negative)

   output ...
   result  an approximation to the integral hopefully satisfying the
           least stringent of the two error tolerances.
   errest  an estimate of the magnitude of the actual error.
   nofun   the number of function values used in calculation of result.
   flag    a reliability indicator.  if flag is zero, then result
           probably satisfies the error tolerance.  if flag is
           xxx.yyy , then  xxx = the number of intervals which have
           not converged and  0.yyy = the fraction of the interval
           left to do when the limit on  nofun  was approached.

   comments ...
   Alex Godunov, Last update 18 February 2007
   the program is based on a fortran version of program quanc8.f
*/
{
    double w0,w1,w2,w3,w4,area,x0,f0,stone,step,cor11,temp;
    double qprev,qnow,qdiff,qleft,esterr,tolerr;
    double qright[32], f[17], x[17], fsave[9][31], xsave[9][31];
    double dabs,dmax1;
    int    levmin,levmax,levout,nomax,nofin,lev,nim,i,j;
    int    key;

//  ***   stage 1 ***   general initialization

    levmin = 1;
    levmax = 30;
    levout = 6;
    nomax = 5000;
    nofin = nomax - 8*(levmax - levout + 128);
//  trouble when nofun reaches nofin

    w0 =   3956.0 / 14175.0;
    w1 =  23552.0 / 14175.0;
    w2 =  -3712.0 / 14175.0;
    w3 =  41984.0 / 14175.0;
    w4 = -18160.0 / 14175.0;

//  initialize running sums to zero.

    flag   = 0.0;
    result = 0.0;
    cor11  = 0.0;
    errest = 0.0;
    area   = 0.0;
    nofun  = 0;
    if (a == b) return;

//  ***   stage 2 ***   initialization for first interval

    lev = 0;
    nim = 1;
    x0 = a;
    x[16] = b;
    qprev  = 0.0;
    f0 = fun(x0);
    stone = (b - a) / 16.0;
    x[8]  =  (x0    + x[16])   / 2.0;
    x[4]  =  (x0    + x[8])    / 2.0;
    x[12] =  (x[8]  + x[16])   / 2.0;
    x[2]  =  (x0    + x[4])    / 2.0;
    x[6]  =  (x[4]  + x[8])    / 2.0;
    x[10] =  (x[8]  + x[12])   / 2.0;
    x[14] =  (x[12] + x[16])   / 2.0;
    for (j=2; j<=16; j = j+2)
    {
      f[j] = fun(x[j]);
    }
    nofun = 9;

//  ***   stage 3 ***   central calculation

    while(nofun <= nomax)
    {
      x[1] = (x0 + x[2]) / 2.0;
      f[1] = fun(x[1]);
      for(j = 3; j<=15; j = j+2)
        {
          x[j] = (x[j-1] + x[j+1]) / 2.0;
          f[j] = fun(x[j]);
        }
      nofun = nofun + 8;
      step  = (x[16] - x0) / 16.0;
      qleft  = (w0*(f0 + f[8])  + w1*(f[1]+f[7])  + w2*(f[2]+f[6])
             + w3*(f[3]+f[5])  +  w4*f[4]) * step;
      qright[lev+1] = (w0*(f[8]+f[16])+w1*(f[9]+f[15])+w2*(f[10]+f[14])
                    + w3*(f[11]+f[13]) + w4*f[12]) * step;
      qnow  = qleft + qright[lev+1];
      qdiff = qnow  - qprev;
      area  = area  + qdiff;

//  ***   stage 4 *** interval convergence test

      esterr = fabs(qdiff) / 1023.0;
      if(abserr >= relerr*fabs(area))
        tolerr = abserr;
        else
        tolerr = relerr*fabs(area);
      tolerr = tolerr*(step/stone);

// multiple logic conditions for the convergence test
      key = 1;
      if (lev < levmin) key = 1;
        else if (lev >= levmax)
        key = 2;
        else if (nofun > nofin)
        key = 3;
        else if (esterr <= tolerr)
        key = 4;
        else
        key = 1;

      switch (key) {
// case 1 ********************************* (mark 50)
      case 1:
//      ***   stage 5   ***   no convergence
//      locate next interval.
        nim = 2*nim;
        lev = lev+1;

//      store right hand elements for future use.
        for(i=1; i<=8; i=i+1)
        {
          fsave[i][lev] = f[i+8];
          xsave[i][lev] = x[i+8];
        }

//      assemble left hand elements for immediate use.
        qprev = qleft;
        for(i=1; i<=8; i=i+1)
        {
          j = -i;
          f[2*j+18] = f[j+9];
          x[2*j+18] = x[j+9];
        }
        continue;  // go to start of stage 3 "central calculation"
      break;

// case 2 ********************************* (mark 62)
      case 2:
        flag = flag + 1.0;
      break;
// case 3 ********************************* (mark 60)
      case 3:
//    ***   stage 6   ***   trouble section
//    number of function values is about to exceed limit.
        nofin = 2*nofin;
        levmax = levout;
        flag = flag + (b - x0) / (b - a);
      break;
// case 4 ********************************* (continue mark 70)
      case 4:
      break;
// default ******************************** (continue mark 70)
      default:
      break;
// end case section ***********************
}

//   ***   stage 7   ***   interval converged
//   add contributions into running sums.

    result = result + qnow;
    errest = errest + esterr;
    cor11  = cor11  + qdiff / 1023.0;

//  locate next interval

    while (nim != 2*(nim/2))
    {
      nim = nim/2;
      lev = lev-1;
    }
    nim = nim + 1;
    if (lev <= 0) break;  // may exit futher calculation

//  assemble elements required for the next interval.

    qprev = qright[lev];
    x0 = x[16];
    f0 = f[16];
    for (i =1; i<=8; i=i+1)
    {
      f[2*i] = fsave[i][lev];
      x[2*i] = xsave[i][lev];
    }
}
//  *** end stage 3 ***   central calculation

//  ***   stage 8   ***   finalize and return

    result = result + cor11;

//  make sure errest not less than roundoff level.

    if (errest == 0.0) return;
    do
    {
     temp = fabs(result) + errest;
     errest = 2.0*errest;
    }
    while (temp == fabs(result));

    return;
}

/*==============================================================
    Numerical integration of f(x) on [a,b]
    method: Monte-Carlo method
    written by: Alex Godunov 
----------------------------------------------------------------
input:
    f   - a single argument real function (supplied by the user)
    a,b - the two end-points of the interval of integration
    n   - number of random points for xi
output:
    r      - result of integration
    errest - estimated error
Comments:
    be sure that following headers are included
    #include <cstdlib>
    #include <ctime>     
================================================================*/

    double mc_int1(double(*f)(double), double& errest, double a, double b, int n)
{
    double r, x, u, f2s, fs;
/* variables fs and f2s are used to estimate an error of integration */    
    
    srand(time(NULL));               /* initial seed value (use system time) */

    fs  = 0.0;
    f2s = 0.0;

    for (int i = 1; i <= n; i=i+1)
    {
        u = 1.0*rand()/(RAND_MAX+1.0); /* random number between 0.0 and 1.0 */
        x = a + (b-a)*u;
        fs = fs + f(x);
        f2s= f2s+ f(x)*f(x);
    }
    r = fs*(b-a)/n;
/* evaluate integration error */
    fs = fs/n;
    f2s = f2s/n;
    errest = (b-a)*sqrt((f2s - fs*fs)/n);
    return r;
}