//=======================================================================
// ch07_quanc8.cpp
//
// Adaptive numerical integration using QUANC8.
//
// The program evaluates
//
//               pi
//               /
//              |  sin(x) dx = 2
//              /
//               0
//
// with the adaptive 8-panel Newton-Cotes QUANC8 algorithm.
//
// QUANC8 repeatedly subdivides intervals and compares successive
// Newton-Cotes estimates to control the integration error.
//
// Based on the QUANC8 routine presented by
// G. E. Forsythe, M. A. Malcolm, and C. B. Moler,
// Computer Methods for Mathematical Computations,
// Prentice-Hall, 1977.
//
// C++ version adapted by Alexander Godunov, February 2007.
// Revised and reorganized for the companion website, 2026.
//=======================================================================

#include <cmath>
#include <iomanip>
#include <iostream>

using namespace std;


// Function prototypes.
double f(double x);

void quanc8(double (*fun)(double),
            double a, double b,
            double abserr, double relerr,
            double& result, double& errest,
            int& nfun, double& flag);


//=======================================================================
int main()
{
    const double pi = 3.14159265358979323846;

    // ---- Problem setup ------------------------------------------------
    const double a      = 0.0;
    const double b      = pi;
    const double abserr = 0.0;
    const double relerr = 1.0e-8;

    const double exact = 2.0;

    double result, errest, flag;
    int nfun;

    quanc8(f, a, b, abserr, relerr, result, errest, nfun, flag);

    cout << "QUANC8 adaptive integration\n";
    cout << "Integral of sin(x) from 0 to pi\n";

    cout << scientific << setprecision(7);
    cout << "Exact value          = " << setw(15) << exact  << '\n';
    cout << "Calculated value     = " << setw(15) << result << '\n';

    cout << scientific << setprecision(4);
    cout << "Estimated error      = " << setw(12) << errest << '\n';

    cout << defaultfloat;
    cout << "Function calls       = " << setw(12) << nfun << '\n';

    cout << fixed << setprecision(5);
    cout << "Flag                 = " << setw(12) << flag << '\n';

    return 0;
}


//=======================================================================
double f(double x)
{
    //-------------------------------------------------------------------
    // Function to be integrated.
    //
    // To use this program for another problem, replace only the line
    // defining y below and change the integration limits in main()
    // if needed.
    //
    // Input:
    //   x - integration variable
    //
    // Output:
    //   return value - value of the integrand
    //-------------------------------------------------------------------

    const double y = sin(x);

    return y;
}


//=======================================================================
void quanc8(double (*fun)(double),
            double a, double b,
            double abserr, double relerr,
            double& result, double& errest,
            int& nfun, double& flag)
{
    //-------------------------------------------------------------------
    // Adaptive integration using the 8-panel Newton-Cotes QUANC8
    // algorithm.
    //
    // Input:
    //   fun    - function to integrate
    //   a      - lower integration limit
    //   b      - upper integration limit; b may be less than a
    //   abserr - requested absolute error tolerance, abserr >= 0
    //   relerr - requested relative error tolerance, relerr >= 0
    //
    // Output:
    //   result - numerical approximation to the integral
    //   errest - estimated magnitude of the numerical error
    //   nfun   - number of function evaluations
    //   flag   - reliability indicator
    //
    // Reliability flag:
    //   flag = 0 indicates that the requested tolerance was probably met.
    //   A nonzero value indicates that one or more intervals did not
    //   converge normally or that the function-evaluation limit was
    //   approached.
    //
    // Method:
    //   QUANC8 is an automatic adaptive quadrature routine based on an
    //   8-panel Newton-Cotes formula. The interval is subdivided where
    //   the difference between successive estimates is too large.
    //
    //   The local error estimate is based on
    //
    //       |Q_new - Q_old| / 1023.
    //
    //   Previously computed function values are reused whenever an
    //   interval is subdivided.
    //
    // Based on the QUANC8 routine presented by
    // G. E. Forsythe, M. A. Malcolm, and C. B. Moler,
    // Computer Methods for Mathematical Computations,
    // Prentice-Hall, 1977.
    //
    // C++ version adapted by Alexander Godunov, February 2007.
    //-------------------------------------------------------------------

    // Initialize outputs so the zero-length interval case is defined.
    result = 0.0;
    errest = 0.0;
    flag   = 0.0;
    nfun   = 0;

    if (abserr < 0.0 || relerr < 0.0) {
        cerr << "quanc8: error tolerances must be non-negative.\n";
        return;
    }

    if (a == b) {
        return;
    }

    // ---- Storage ------------------------------------------------------
    double qright[32] = {};
    double fval[17]   = {};
    double x[17]      = {};
    double fsave[9][31] = {};
    double xsave[9][31] = {};

    // ---- General initialization --------------------------------------
    int levmin = 1;
    int levmax = 30;
    int levout = 6;
    int nomax  = 5000;

    // Trouble section is entered when nfun approaches this limit.
    int nofin = nomax - 8*(levmax - levout + 128);

    // Newton-Cotes coefficients.
    const double w0 =   3956.0 / 14175.0;
    const double w1 =  23552.0 / 14175.0;
    const double w2 =  -3712.0 / 14175.0;
    const double w3 =  41984.0 / 14175.0;
    const double w4 = -18160.0 / 14175.0;

    double cor11 = 0.0;
    double area  = 0.0;

    // ---- Initialize first interval -----------------------------------
    int lev = 0;
    int nim = 1;

    double x0    = a;
    x[16]        = b;
    double qprev = 0.0;
    double f0    = fun(x0);

    const double 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 (int j = 2; j <= 16; j += 2) {
        fval[j] = fun(x[j]);
    }

    nfun = 9;

    // ---- Main adaptive loop ------------------------------------------
    while (nfun <= nomax) {

        // Complete the 17-point grid on the current interval.
        x[1] = (x0 + x[2])/2.0;
        fval[1] = fun(x[1]);

        for (int j = 3; j <= 15; j += 2) {
            x[j] = (x[j-1] + x[j+1])/2.0;
            fval[j] = fun(x[j]);
        }

        nfun += 8;

        const double step = (x[16]-x0)/16.0;

        // Newton-Cotes estimates on the left and right halves.
        const double qleft =
            ( w0*(f0       + fval[8])
            + w1*(fval[1]  + fval[7])
            + w2*(fval[2]  + fval[6])
            + w3*(fval[3]  + fval[5])
            + w4*fval[4] )*step;

        qright[lev+1] =
            ( w0*(fval[8]  + fval[16])
            + w1*(fval[9]  + fval[15])
            + w2*(fval[10] + fval[14])
            + w3*(fval[11] + fval[13])
            + w4*fval[12] )*step;

        const double qnow  = qleft + qright[lev+1];
        const double qdiff = qnow - qprev;

        area += qdiff;

        // ---- Local convergence test ----------------------------------
        const double esterr = fabs(qdiff)/1023.0;

        double tolerr = (abserr >= relerr*fabs(area))
                      ? abserr
                      : relerr*fabs(area);

        tolerr *= step/stone;

        int key;

        if (lev < levmin) {
            key = 1;
        }
        else if (lev >= levmax) {
            key = 2;
        }
        else if (nfun > nofin) {
            key = 3;
        }
        else if (esterr <= tolerr) {
            key = 4;
        }
        else {
            key = 1;
        }

        switch (key) {

            case 1:
                // No convergence: subdivide the current interval.
                nim = 2*nim;
                lev = lev + 1;

                // Save the right half for later.
                for (int i = 1; i <= 8; ++i) {
                    fsave[i][lev] = fval[i+8];
                    xsave[i][lev] = x[i+8];
                }

                // Continue immediately with the left half.
                qprev = qleft;

                for (int i = 1; i <= 8; ++i) {
                    const int j = -i;
                    fval[2*j+18] = fval[j+9];
                    x[2*j+18]    = x[j+9];
                }

                continue;

            case 2:
                // Maximum subdivision level reached.
                flag += 1.0;
                break;

            case 3:
                // Function-evaluation limit is being approached.
                nofin  = 2*nofin;
                levmax = levout;
                flag  += (b-x0)/(b-a);
                break;

            case 4:
                // Current interval satisfies the convergence test.
                break;
        }

        // ---- Accept current interval ---------------------------------
        result += qnow;
        errest += esterr;
        cor11  += qdiff/1023.0;

        // Locate the next interval that still has to be processed.
        while (nim != 2*(nim/2)) {
            nim = nim/2;
            lev = lev - 1;
        }

        nim = nim + 1;

        if (lev <= 0) {
            break;
        }

        // Restore saved data for the next interval.
        qprev = qright[lev];
        x0    = x[16];
        f0    = fval[16];

        for (int i = 1; i <= 8; ++i) {
            fval[2*i] = fsave[i][lev];
            x[2*i]    = xsave[i][lev];
        }
    }

    // ---- Final correction and error estimate -------------------------
    result += cor11;

    if (errest == 0.0) {
        return;
    }

    // Make sure errest is not below the representable roundoff level.
    // This is the intended QUANC8 roundoff guard.
    double temp = fabs(result) + errest;

    while (temp == fabs(result)) {
        errest *= 2.0;
        temp = fabs(result) + errest;
    }
}