/*
 Numerical integration of f(x) on [a,b] using 3 following methodss:
 # Trapezoid approximation
 # Simpson's rule
 # 8-panel newton-cotes rule (quanc8.cpp)
 time of integration is calculated and printed
 Update: AG - Februaru 2022
*/
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <ctime>
//#include "trapezoid.cpp"
//#include "simpson.cpp"
//#include "quanc8.cpp"

using namespace std;
//function prototypes
double trapezoid(double(*f)(double), double a, double b, int n);
double simpson(double(*f)(double), 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, trap, simp;
    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
       
    cout.precision(6);
    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
    n = 2;                // initial number of intervals
    ntimes = 18;          // number of interval doublings with nmax=2^(ntimes)

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

/* step 2: integration with trapezoid and simpson */
    for (i=1; i <=ntimes; i=i+1)
    {
        trap = trapezoid(f, a, b, n);
        simp = simpson(f, a, b, n);
        cout << setw(10)  << n << setw(10) << trap << setw(10) << simp <<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   << 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;
  
//    system ("pause");
    return 0;
}
/*
 *  Function for integration
*/
    double f(double x)
{
    const double pi = 3.1415926;
    double y;
    y = sin(x);
//    y = x*sin(30.0*x)*cos(50.0*x);                                   // int (1)
//    y = sin(pi*x*x)/((x-pi)*(x-pi)+1.0);           // int (2)
//    y = sin(20*x)/(x*x + 1.0);
//    y = 1 - cos(32.0*x);
//    y = log(x)*sin(10.0*pi*x);
    return y;
}


/*==============================================================
    Numerical integration of f(x) on [a,b]
    method: trapezoid rule
    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 intervals
output:
    r - result of integration
================================================================*/

    double trapezoid(double(*f)(double), double a, double b, int n)
{
    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;
}

/*==============================================================
    Numerical integration of f(x) on [a,b]
    method: Simpson rule
    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 intervals
output:
    s - result of integration
================================================================*/

    double simpson(double(*f)(double), double a, double b, int n)
{
    double s, dx, x;
// if n is odd - 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=2; i<=n-1; i=i+2)
    {
        x = a+static_cast<float>(i)*dx;
        s = s + 2.0*f(x) + 4.0*f(x+dx);
    }
    s = (s + f(a)+f(b)+4.0*f(a+dx) )*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 (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;
}




/* Test output for f(x) = sin(x) on [0.0,pi]

 Intervals  Trapez.   Simpson
         2  1.570796  2.094395
         4  1.896119  2.004560
         8  1.974232  2.000269
        16  1.993570  2.000017
        32  1.998393  2.000001
        64  1.999598  2.000000
       128  1.999900  2.000000
       256  1.999975  2.000000
       512  1.999994  2.000000
      1024  1.999998  2.000000
      2048  2.000000  2.000000
      4096  2.000000  2.000000

    result from quanc8
    nofun =         33
 integral =   2.000000
est.error =   0.000000
     flag =   0.000000

*/

/* Test output for f(x) = sin(pi*x*x)/((x-pi)*(x-pi)+1.0) on [0.0,2.0*pi]

 Intervals  Trapez.   Simpson
         2 -1.395444 -1.764472
         4  0.026751  0.500817
         8 -0.747037 -1.004966
        16 -0.297934 -0.148234
        32  0.225043  0.399368
        64  0.042146 -0.018819
       128  0.042308  0.042362
       256  0.042339  0.042349
       512  0.042346  0.042348
      1024  0.042348  0.042348
      2048  0.042348  0.042348
      4096  0.042348  0.042348

    result from quanc8
    nofun =       1153
 integral =   0.042348
est.error =   0.000000
     flag =   0.000000

*/