/*
 Numerical integration of f(x) on [a,b]
 Three methods:
 1. adaptive Simpson
 2. adaptive Gauss-Kronrod 7-15
 3. adaptive Newton-Cotes Quanc8
 Last update: february 2022
*/
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <ctime>

using namespace std;

//function prototypes
double f(double);
double gauss715a(double(*f)(double),double a,double b, double eps,
                 double& epsout, int& ncall);
double  simpson1(double(*f)(double),double a,double b, double eps,
                 double& epsout, int& nfun);
void      quanc8(double(*fun)(double), double a, double b,
                 double abserr, double relerr,
                 double& result, double& errest, int& nofun,double& flag);

int main()
{
    double a, b, eps;
    double epsout, integral;
    double nc8, errest, flag, abserr, relerr;
    int ncall, nfun, nofun;
    const double pi = 3.1415926535897932;
       
    cout.precision(12);
    cout.setf(ios::fixed | ios::showpoint);

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

/* initial information */
    a = 0.0;             // left endpoint
    b = 1.0*pi;          // right endpoint
    eps = 1.0e-8;        // tolerance
    
/* Results for adaptive Simpson */

    integral =  simpson1(f, a, b, eps, epsout, nfun);
    cout << endl;
    cout << " Adaptive Simpson = " << setw(14) << integral   << endl;
    cout << " Tolerance        = " << setw(14) << eps    << endl;
    cout << " Tol. estimated   = " << setw(14) << epsout << endl;
    cout << " Function calls   = " << setw(14) << nfun  << endl;


/*  Result for adaptive Gauss-Kronrod */
    integral = gauss715a(f, a, b, eps, epsout, ncall);
    
    cout << endl;
    cout << " Adaptive Gauss-K = " << setw(14) << integral   << endl;
    cout << " Tolerance        = " << setw(14) << eps    << endl;
    cout << " Tol. estimated   = " << setw(14) << epsout << endl;
    cout << " Function calls   = " << setw(14) << ncall  << endl;
   
/* Results for 8-panel Newton-Cotes rule */
    abserr=0.0;
    relerr = eps;
    quanc8(f, a, b, abserr, relerr, nc8, errest, nofun, flag);

    cout << endl;
    cout << " Integral Quanc8  = " << setw(14) << nc8   << endl;
    cout << " Tolerance        = " << setw(14) << eps    << endl;
    cout << " Tol. estimated   = " << setw(14) << errest << endl;
    cout << " Function calls   = " << setw(14) << nofun  << endl;
//    cout << " Flag             = " << setw(14) << 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.1415926535897932;
    double y;
   y = sin(x);
//   y = 1.0/(1.0-0.998*pow(x,2));
//   y = x*sin(30.0*x)*cos(50.0*x);
//        y = x/(exp(x)-1.0);
//    y = sin(pi*x*x)/((x-pi)*(x-pi)+1.0);
//    y = sin(20*x)/(x*x + 1.0);
//    y = 1 - cos(32.0*x);
//    y = 4.0/(x*x + 1.0);
//    y = pow(x-1.0,1.0/5.0)/(x*x+1.0);
//    y = (2.0/sqrt(pi))*exp(-x*x);

//    y = exp(0.0-x*x)/(1+x*x);
//    y = x*sin(x)/(1+x*x);

/* from a battery test */
//y = exp(x);   // 1 on [0,1].   1.7182818284590452354
//y = 4.0*pi*pi*x*sin(20.0*pi*x)*cos(2.0*pi*x);
    //on [0,1] -0.63466518254339257343
//y = cos(cos(x)+3.0*sin(x)+2.0*cos(2.0*x)+3.0*sin(2.0*x)+3.0*cos(3.0*x));
    // on [0,pi] exact value=0.83867634269442961454

return y;
}


/*==============================================================
 Integration of f(x) on [a,b]
 Method: Gauss7 and Kronrod15 rules applied to whole interval
 This is a service program for adaptive non-recursive gauss715a
 written by: Alex Godunov 
--------------------------------------------------------------
 IN:
 f   - Function to integrate
 a   - Lower limit of integration
 b   - Upper limit of integration

 OUT:
 g7  - integral Gauss 7 point rule
 k15 - integral Kronrod 15 point rule
============================================================== */
double GK715(double(*f)(double),double a,double b, double& g7, double& k15)

{
double h, c, f0;
int i, j;
double fl[8], fr[8];
double x[8] = {0.000000000000000, 0.207784955007898,
               0.405845151377397, 0.586087235467691,
               0.741531185599394, 0.864864423359769,
               0.949107912342758, 0.991455371120813};
double wg[4]= {0.417959183673469, 0.381830050505119,
               0.279705391489277, 0.129484966168870};
double wk[8]= {0.209482141084728, 0.204432940075298,
               0.190350578064785, 0.169004726639267,
               0.140653259715525, 0.104790010322250,
               0.063092092629979, 0.022935322010529};

/* *** initialization *** */
g7  = 0.0;
k15 = 0.0;

h = (b-a)/2.0;
c = (b+a)/2.0;

/* evaluate function f at points x using the symmetry for -x and +x
   the gauss interval [-1,1] has to be rescaled to [a,b]
   integral {f(x)dx} on [a,b]
   = ((b-a)/2)*integral{f(((b-a)/2)*y + (b+a)/2)}dy on [-1,1] */
f0 = f(c);
for (i = 1; i<=7; i=i+1)
{
  fl[i] = f(h*(-1.0)*x[i]+c);
  fr[i]  = f(h*x[i]+c);
}

/* evaluate Gauss 7 */
g7 = wg[0]*f0;
for (i = 1; i<=3; i=i+1)
{
  g7 = g7 + wg[i]*(fl[2*i]+fr[2*i]);
}
g7 = g7*h;

/* evaluate Kronrod 15 */
k15 = wk[0]*f0;
for (i = 1; i<=7; i=i+1)
{
  k15 = k15 + wk[i]*(fl[i]+fr[i]);
}
k15 = k15*h;
    return 0;
}
/* end of function GK715 */


/*=======================================================================
 Integration of f(x) on [a,b]
 gauss715a(f, a, b, eps, epsout, ncall)
 Method: adaptive non-recursive Gauss-Kronrod 7-15 rule
 written by: Alex Godunov
-------------------------------------------------------------------------
 IN:
 f    - Function to integrate
 a    - Lower limit of integration
 b    - Upper limit of integration
 eps  - Error tolerance (200.0*abs(k15-g7))**1.5 <=eps

 OUT:
 gauss715a - result of integration
 epsout    - estimated tolerance
 ncall     - number of function values used in calculation

 PARAMETERS:
 im - max number of levels (i.e. the smallest interval (b-a)/2^imax)

 Comment:
 the function calls subroutine GK715

 20.07.2012
1. think about relative error estimation (see the condition)
2. program calculates how many subintervals failed to converge (nfail)
   and what fraction of [a,b] failed to converge (xfail)
!========================================================================  */
double gauss715a(double(*f)(double),double a,double b, double eps,
                 double& epsout, int& ncall)

{
const int im=32;
double tol[im], x[im], h[im];
double x0, step, err, sum, xfail;
double g7, k15;
int il[im];
int i, imax, deep;
int nfail;

/* stage 1: initialization for level 1 (top level) */

sum    = 0.0;
epsout = 0.0;
ncall  = 0;
i      = 1;
imax   = im;

nfail  = 0;
xfail  = 0.0;

x[1]   = a;
h[1]   = b-a;
tol[1] = eps;
il[1]  = 1;

/* stage 2: main part: adaptive integration */

while (i > 0) //major loop
{
/* stage 2a: calculate integral on [x(i),(x(i)+h(i))] */
  GK715(f,x[i],x[i]+h[i],g7,k15);
  ncall = ncall + 15;

/* stage 2b: save data at this level */
  x0   = x[i];
  step = h[i];
  err  = tol[i];
  deep = il[i];

/* stage 2c: the current level has been computed */
  i=i-1;

/* stage 2d: local condition for convergence */
  if(pow(200.0*fabs(k15-g7),1.5) <=err)
  {
  sum = sum + k15;
  epsout = epsout + pow(200.0*fabs(k15-g7),1.5);
  }

/* stage 2e: check if the code can continue to subdivide intervals */
  else if( deep >= imax )
    {nfail = nfail + 1;
     xfail = xfail + step;}
/* stage 2f: make smaller intervals (i.e. h=h/2) */
/* data for the right subinterval                */
     else
     {
      i = i+1;
      h[i]  = step/2.0;
      x[i]  = x0 + h[i];
         tol[i]= err/2.0;
      il[i] = deep + 1;
/* data for the left subinterval */
      i = i+1;
      x[i]  = x0;
      h[i]  = h[i-1];
      tol[i]= tol[i-1];
      il[i] = il[i-1];
     };
   
} //end of the major loop

xfail = xfail/fabs(b-a);

return sum;

}

/*==========================================================
 Integration of f(x) on [a,b]
 Method: adaptive non-recursive Simpson rule
 written by: Alex Godunov 
----------------------------------------------------------
 IN:
 f    - Function to integrate
 a    - Lower limit of integration
 b    - Upper limit of integration
 eps  - Error tolerance abs(s1+s2-s0)<=15*eps (using 10*eps recommended)

 OUT:
 sum  - the result of integration
 nfun - the number of function values used in calculation

 PARAMETERS:
 im - max number of levels (i.e. the smallest interval (b-a)/2**imax)
 nmax - max number of function calls

 Comment:
 reference points for two Simpson's intervals
 numbers      0   1   2   3   4
 in the code  a---1---m---3---b

 14.07.2012
 so far nmax has not been used
 think about relative error estimation (see the condition)
========================================================== */
double  simpson1(double(*f)(double),double a,double b, double eps,
                 double& epsout, int& nfun)
{
const int im=32;
double tol[im], x[im], h[im], fa[im], fm[im], fb[im], s[im];
double sum, step, err;
double x0, f0, f1, f2, f3, f4, s0, s1, s2;
int il[im];
int i, imax, deep;
int nfail;

/* stage 1: initialization for level 1 (top level) */

sum = 0.0;
epsout = 0.0;
i = 1;
imax = im;
nfail = 0;

x[1]   = a;
h[1]   = (b-a)/2.0;
fa[1]  = f(a);
fm[1]  = f(a+h[1]);
fb[1]  = f(b);
tol[1] = 15.0*eps;
il[1]  = 1;

/* Simpson's method for [a,b] */
s[1] = h[1]*(fa[1]+4*fm[1]+fb[1])/3.0;
nfun = 3;

/* stage 2: main part: adaptive integration */

while (i > 0)
{
/* stage 2a: calculate function values at h/2 and 3h/2 */
  f1 = f(x[i] +     h[i]/2.0);
  f3 = f(x[i] + 3.0*h[i]/2.0);
  nfun = nfun + 2;

/* Simpson's integrals for the left and right intervals */
  s1 = h[i]*(fa[i]+4.0*f1+fm[i])/6.0;
  s2 = h[i]*(fm[i]+4.0*f3+fb[i])/6.0;

/* stage 2b: save data at this level */
  x0 = x[i];
  f0 = fa[i];
  f2 = fm[i];
  f4 = fb[i];
  step = h[i];
  err  = tol[i];
  s0   = s[i];
  deep = il[i];

/* stage 2c: the current level has been computed */
  i=i-1;
  
/* stage 2d: local condition for convergence */
  if(abs(s1+s2-s0) <= err )
    {sum = sum + (s1+s2);
     epsout = epsout + abs(s1+s2-s0)/15.0;}

/* stage 2e: check if the code can continue to subdivide intervals */
      else if( deep >= imax )
      {
       nfail = nfail + 1;
      }

/* stage 2f: make smaller intervals (i.e. h=h/2) */
//    data for the right subinterval
      else
      {
        i = i+1;
        x[i]  = x0 + step;
        fa[i] = f2;
        fm[i] = f3;
        fb[i] = f4;
        h[i]  = step/2.0;
        tol[i]= err/2.0;
        s[i]  = s2;
        il[i] = deep + 1;
//     data for the left subinterval
        i = i+1;
        x[i]  = x0;
        fa[i] = f0;
        fm[i] = f1;
        fb[i] = f2;
        h[i]  = h[i-1];
        tol[i]= tol[i-1];
        s[i]  = s1;
        il[i] = il[i-1];
        }
  
} // end major loop while

return sum;
        
}

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