//=======================================================================
// ch07_adaptive_gauss_kronrod.cpp
//
// Adaptive numerical integration using embedded Gauss-Kronrod rules.
//
// Two quadrature pairs are included:
//   1. Gauss 7-point / Kronrod 15-point (G7-K15)
//   2. Gauss 10-point / Kronrod 21-point (G10-K21)
//
// The same non-recursive adaptive subdivision algorithm is used with
// either pair.
//
// Example:
//   Integral of 1/[1+400(x-0.2)^2] from 0 to 1.
//
// Exact value:
//   [atan(16) + atan(4)] / 20
//
// The integrand has a narrow peak near x = 0.2, making it a useful
// example for adaptive integration.
//
// Based on integration codes written by Alexander Godunov, July 2012.
// Revised and reorganized for the companion website, 2026.
//=======================================================================

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

using namespace std;


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

void G7K15(double (*fun)(double), double a, double b,
           double& g7, double& k15);

void G10K21(double (*fun)(double), double a, double b,
            double& g10, double& k21);

void GKinteg(double (*fun)(double), double a, double b, double eps,
             double& result, double& errest, int& nfun, int key);


//=======================================================================
int main()
{
    const double a   = 0.0;
    const double b   = 1.0;
    const double eps = 1.0e-8;

    const double exact = (atan(16.0) + atan(4.0))/20.0;

    double result, errest;
    int nfun;

    cout << "Adaptive Gauss-Kronrod integration\n";
    cout << "Integral of 1/[1+400(x-0.2)^2] from 0 to 1\n";

    cout << scientific << setprecision(10);
    cout << "Exact value          = " << setw(18) << exact << "\n\n";

    GKinteg(f, a, b, eps, result, errest, nfun, 1);

    cout << "G7-K15 result        = " << setw(18) << result << '\n';
    cout << scientific << setprecision(4);
    cout << "Estimated error      = " << setw(12) << errest << '\n';
    cout << defaultfloat;
    cout << "Function calls       = " << setw(10) << nfun << "\n\n";

    GKinteg(f, a, b, eps, result, errest, nfun, 2);

    cout << scientific << setprecision(10);
    cout << "G10-K21 result       = " << setw(18) << result << '\n';
    cout << scientific << setprecision(4);
    cout << "Estimated error      = " << setw(12) << errest << '\n';
    cout << defaultfloat;
    cout << "Function calls       = " << setw(10) << nfun << '\n';

    return 0;
}


//=======================================================================
double f(double x)
{
    //-------------------------------------------------------------------
    // Test function for the integration example.
    //
    // Input:
    //   x - point at which the function is evaluated
    //
    // Output:
    //   return value - 1/[1+400(x-0.2)^2]
    //
    // To use another integrand, replace only the line defining y below
    // and change the integration limits in main() if needed.
    //-------------------------------------------------------------------

    const double y = 1.0/(1.0 + 400.0*(x-0.2)*(x-0.2));

    return y;
}


//=======================================================================
void G7K15(double (*fun)(double), double a, double b,
           double& g7, double& k15)
{
    //-------------------------------------------------------------------
    // Gauss-Kronrod quadrature pair G7-K15 on one interval [a,b].
    //
    // Input:
    //   fun - function to integrate
    //   a   - lower integration limit
    //   b   - upper integration limit
    //
    // Output:
    //   g7  - Gauss 7-point approximation
    //   k15 - Kronrod 15-point approximation
    //
    // The 7 Gauss points are contained within the 15 Kronrod points.
    // Both estimates therefore reuse the same function evaluations.
    // Their difference is used by GKinteg as a local error indicator.
    //
    // The tabulated nodes and weights refer to [-1,1] and are mapped
    // linearly to [a,b].
    //-------------------------------------------------------------------

    const double x[8] = {
        0.000000000000000,
        0.207784955007898,
        0.405845151377397,
        0.586087235467691,
        0.741531185599394,
        0.864864423359769,
        0.949107912342758,
        0.991455371120813
    };

    const double wg[4] = {
        0.417959183673469,
        0.381830050505119,
        0.279705391489277,
        0.129484966168870
    };

    const double wk[8] = {
        0.209482141084728,
        0.204432940075298,
        0.190350578064785,
        0.169004726639267,
        0.140653259715525,
        0.104790010322250,
        0.063092092629979,
        0.022935322010529
    };

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

    const double f0 = fun(c);

    double fl[7], fr[7];

    for (int i = 1; i <= 7; ++i) {
        fl[i-1] = fun(c - h*x[i]);
        fr[i-1] = fun(c + h*x[i]);
    }

    // Gauss 7-point estimate.
    g7 = wg[0]*f0;

    for (int i = 1; i <= 3; ++i) {
        g7 += wg[i]*(fl[2*i-1] + fr[2*i-1]);
    }

    g7 *= h;

    // Kronrod 15-point estimate.
    k15 = wk[0]*f0;

    for (int i = 1; i <= 7; ++i) {
        k15 += wk[i]*(fl[i-1] + fr[i-1]);
    }

    k15 *= h;
}


//=======================================================================
void G10K21(double (*fun)(double), double a, double b,
            double& g10, double& k21)
{
    //-------------------------------------------------------------------
    // Gauss-Kronrod quadrature pair G10-K21 on one interval [a,b].
    //
    // Input:
    //   fun - function to integrate
    //   a   - lower integration limit
    //   b   - upper integration limit
    //
    // Output:
    //   g10 - Gauss 10-point approximation
    //   k21 - Kronrod 21-point approximation
    //
    // The 10 Gauss points are contained within the 21 Kronrod points.
    // Both estimates therefore reuse the same function evaluations.
    // Their difference is used by GKinteg as a local error indicator.
    //
    // The tabulated nodes and weights refer to [-1,1] and are mapped
    // linearly to [a,b].
    //-------------------------------------------------------------------

    const double x[11] = {
        0.0000000000000000000000000,
        0.1488743389816312108848260,
        0.2943928627014601981311266,
        0.4333953941292471907992659,
        0.5627571346686046833390001,
        0.6794095682990244062343274,
        0.7808177265864168970637176,
        0.8650633666889845107320967,
        0.9301574913557082260012072,
        0.9739065285171717200779640,
        0.9956571630258080807355273
    };

    const double wk[11] = {
        0.1494455540029169056649365,
        0.1477391049013384913748415,
        0.1427759385770600807970943,
        0.1347092173114733259280540,
        0.1234919762620658510779581,
        0.1093871588022976418992106,
        0.0931254545836976055350655,
        0.0750396748109199527670431,
        0.0547558965743519960313813,
        0.0325581623079647274788190,
        0.0116946388673718742780644
    };

    const double wg[5] = {
        0.2955242247147528701738930,
        0.2692667193099963550912269,
        0.2190863625159820439955349,
        0.1494513491505805931457763,
        0.0666713443086881375935688
    };

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

    const double f0 = fun(c);

    double fl[10], fr[10];

    for (int i = 1; i <= 10; ++i) {
        fl[i-1] = fun(c - h*x[i]);
        fr[i-1] = fun(c + h*x[i]);
    }

    // Gauss 10-point estimate.
    g10 = 0.0;

    for (int i = 0; i < 5; ++i) {
        const int j = 2*i + 1;  // x indices 1,3,5,7,9
        g10 += wg[i]*(fl[j-1] + fr[j-1]);
    }

    g10 *= h;

    // Kronrod 21-point estimate.
    k21 = wk[0]*f0;

    for (int i = 1; i <= 10; ++i) {
        k21 += wk[i]*(fl[i-1] + fr[i-1]);
    }

    k21 *= h;
}


//=======================================================================
void GKinteg(double (*fun)(double), double a, double b, double eps,
             double& result, double& errest, int& nfun, int key)
{
    //-------------------------------------------------------------------
    // Adaptive non-recursive Gauss-Kronrod integration on [a,b].
    //
    // Input:
    //   fun - function to integrate
    //   a   - lower integration limit
    //   b   - upper integration limit
    //   eps - requested tolerance used in the local convergence test
    //   key - quadrature pair:
    //           1 : Gauss 7-point  / Kronrod 15-point
    //           2 : Gauss 10-point / Kronrod 21-point
    //
    // Output:
    //   result - numerical approximation to the integral
    //   errest - accumulated absolute Gauss-Kronrod differences
    //   nfun   - total number of function evaluations
    //
    // Method:
    //   The current interval is evaluated with an embedded
    //   Gauss-Kronrod pair. If the Gauss-Kronrod difference satisfies
    //   the local error test, the Kronrod estimate is accepted.
    //   Otherwise the interval is divided into two equal subintervals
    //   and both are placed on an explicit stack for later evaluation.
    //
    //   This is the non-recursive adaptive strategy used in the
    //   original integration codes.
    //-------------------------------------------------------------------

    const int im = 32;

    double tol[im], x[im], h[im];
    int level[im];

    result = 0.0;
    errest = 0.0;
    nfun   = 0;

    int i = 0;

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

    while (i >= 0) {

        double gauss, kronrod;

        if (key == 1) {
            G7K15(fun, x[i], x[i] + h[i], gauss, kronrod);
            nfun += 15;
        }
        else if (key == 2) {
            G10K21(fun, x[i], x[i] + h[i], gauss, kronrod);
            nfun += 21;
        }
        else {
            cerr << "GKinteg error: key must be 1 or 2.\n";
            result = 0.0;
            errest = 0.0;
            nfun = 0;
            return;
        }

        const double x0   = x[i];
        const double step = h[i];
        const double err  = tol[i];
        const int deep    = level[i];

        --i;

        const double difference = fabs(kronrod-gauss);

        if (pow(200.0*difference, 1.5) <= err) {
            result += kronrod;
            errest += difference;
        }
        else {
            if (deep >= im) {
                cerr << "GKinteg error: maximum subdivision depth reached.\n";
                return;
            }

            // Right subinterval.
            ++i;
            h[i]     = step/2.0;
            tol[i]   = err/2.0;
            x[i]     = x0 + h[i];
            level[i] = deep + 1;

            // Left subinterval.
            ++i;
            x[i]     = x0;
            h[i]     = h[i-1];
            tol[i]   = tol[i-1];
            level[i] = level[i-1];
        }
    }
}