//=======================================================================
// ch07_gauss_legendre.cpp
//
// Gauss-Legendre quadrature for numerical integration.
//
// The program evaluates
//
//       integral sin(x) dx,   0 <= x <= pi,
//
// using several n-point Gauss-Legendre rules. The exact value is 2.
//
// The Gauss-Legendre nodes and weights are generated numerically.
// Roots of the Legendre polynomial P_n(x) are found by Newton iteration,
// and symmetry is used so that only half of the roots are computed.
//
// Based on the original 8-point and 16-point Gauss quadrature programs
// written by Alexander Godunov, October 2009.
// C++ version prepared for the companion website, 2026.
//=======================================================================

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

using namespace std;

double f(double x);
double gauss_legendre(double (*fun)(double),
                      double a, double b, int n);

int main()
{
    const double pi = 3.14159265358979323846;
    const double a = 0.0;
    const double b = pi;
    const double exact = 2.0;
    const int nvalues[] = {2, 4, 8, 16};

    cout << "Gauss-Legendre quadrature\n";
    cout << "Integral of sin(x) from 0 to pi\n";
    cout << scientific << setprecision(7);
    cout << "Exact value          = " << setw(15) << exact << "\n\n";

    cout << "        n       Gauss-Legendre\n";
    cout << "--------------------------------\n";

    for (int n : nvalues) {
        const double integral = gauss_legendre(f, a, b, n);

        cout << setw(9) << defaultfloat << n
             << scientific << setprecision(7)
             << setw(21) << integral << '\n';
    }

    return 0;
}


double f(double x)
{
    //-------------------------------------------------------------------
    // Test function for the integration example.
    //-------------------------------------------------------------------
    return sin(x);
}


double gauss_legendre(double (*fun)(double),
                      double a, double b, int n)
{
    //-------------------------------------------------------------------
    // n-point Gauss-Legendre quadrature on [a,b].
    //
    // The nodes are the roots of P_n(x), found by Newton iteration.
    // The corresponding weights are calculated from P'_n(x).
    // Symmetry is used so that only half of the roots are computed.
    //-------------------------------------------------------------------

    const double pi = 3.14159265358979323846;
    const double tol = 1.0e-14;
    const int max_iter = 100;

    if (n < 1) {
        cerr << "gauss_legendre: n must be positive.\n";
        return 0.0;
    }

    const double midpoint = (a+b)/2.0;
    const double halfwidth = (b-a)/2.0;

    double integral = 0.0;
    const int m = (n+1)/2;

    for (int i = 1; i <= m; ++i) {

        double z = cos(pi*(static_cast<double>(i)-0.25) /
                       (static_cast<double>(n)+0.5));

        double p0 = 0.0, p1 = 0.0, p2 = 0.0, dp = 0.0;
        bool converged = false;

        for (int iter = 0; iter < max_iter; ++iter) {

            p0 = 1.0;
            p1 = z;

            if (n == 1) {
                p2 = p1;
            } else {
                for (int j = 2; j <= n; ++j) {
                    p2 = ((2.0*j-1.0)*z*p1 - (j-1.0)*p0)/j;
                    p0 = p1;
                    p1 = p2;
                }
            }

            dp = (n == 1) ? 1.0
                 : n*(z*p1-p0)/(z*z-1.0);

            const double zold = z;
            z = zold - p1/dp;

            if (fabs(z-zold) <= tol) {
                converged = true;
                break;
            }
        }

        if (!converged) {
            cerr << "gauss_legendre: Newton iteration did not converge.\n";
            return 0.0;
        }

        // Re-evaluate P_n and P_(n-1) at the converged root.
        p0 = 1.0;
        p1 = z;

        if (n == 1) {
            dp = 1.0;
        } else {
            for (int j = 2; j <= n; ++j) {
                p2 = ((2.0*j-1.0)*z*p1 - (j-1.0)*p0)/j;
                p0 = p1;
                p1 = p2;
            }

            dp = n*(z*p1-p0)/(z*z-1.0);
        }

        const double weight = 2.0/((1.0-z*z)*dp*dp);

        const double xleft  = midpoint - halfwidth*z;
        const double xright = midpoint + halfwidth*z;

        if (fabs(z) <= tol) {
            integral += weight*fun(midpoint);
        } else {
            integral += weight*(fun(xleft) + fun(xright));
        }
    }

    return halfwidth*integral;
}