//=======================================================================
// ch07_gauss_legendre_2d.cpp
//
// Two-dimensional integration using a Gauss-Legendre product rule.
//
// The program evaluates
//
//          1   1
//          /   /
//         |   |        1
//     I = |   |   -------------  dy dx
//         |   |   1 + x^2 + y^2
//         /   /
//          0   0
//
// An n-point Gauss-Legendre rule is used in each direction. The
// resulting tensor-product rule requires n^2 function evaluations.
//
// Based on the Fortran implementation by Alexander Godunov.
// C++ version prepared for the companion website, 2026.
//=======================================================================

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

using namespace std;

double f(double x, double y);

double gauss_legendre_2d(double (*fun)(double, double),
                         double ax, double bx,
                         double ay, double by,
                         int n, int& nfun);

bool gauss_legendre_nodes_weights(int n, double a, double b,
                                  vector<double>& node,
                                  vector<double>& weight);

int main()
{
    const double ax = 0.0;
    const double bx = 1.0;
    const double ay = 0.0;
    const double by = 1.0;

    const double reference = 6.3951035187031102e-01;
    const int nvalues[] = {2, 4, 8, 16, 32};

    cout << "Two-dimensional Gauss-Legendre integration\n";
    cout << "Integral of 1/(1+x^2+y^2) over [0,1] x [0,1]\n";
    cout << scientific << setprecision(7);
    cout << "Reference value      = " << setw(15) << reference << "\n\n";

    cout << "        n     Function calls          Integral       Absolute error\n";
    cout << "-------------------------------------------------------------------\n";

    for (int n : nvalues) {
        int nfun = 0;
        const double integral =
            gauss_legendre_2d(f, ax, bx, ay, by, n, nfun);
        const double error = fabs(integral-reference);

        cout << setw(9) << defaultfloat << n
             << setw(19) << nfun
             << scientific << setprecision(7)
             << setw(20) << integral
             << setw(20) << error << '\n';
    }

    return 0;
}


double f(double x, double y)
{
    //-------------------------------------------------------------------
    // Function to be integrated over the rectangular domain.
    //
    // To use another problem, replace only the line below and change
    // the limits in main() if needed.
    //-------------------------------------------------------------------
    return 1.0/(1.0 + x*x + y*y);
}


double gauss_legendre_2d(double (*fun)(double, double),
                         double ax, double bx,
                         double ay, double by,
                         int n, int& nfun)
{
    //-------------------------------------------------------------------
    // Two-dimensional Gauss-Legendre tensor-product rule.
    //-------------------------------------------------------------------

    if (n < 1) {
        cerr << "gauss_legendre_2d: n must be positive.\n";
        nfun = 0;
        return 0.0;
    }

    vector<double> xnode, xweight;
    vector<double> ynode, yweight;

    if (!gauss_legendre_nodes_weights(n, ax, bx, xnode, xweight) ||
        !gauss_legendre_nodes_weights(n, ay, by, ynode, yweight)) {
        nfun = 0;
        return 0.0;
    }

    double integral = 0.0;
    nfun = 0;

    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            integral += xweight[i]*yweight[j]*
                        fun(xnode[i], ynode[j]);
            ++nfun;
        }
    }

    return integral;
}


bool gauss_legendre_nodes_weights(int n, double a, double b,
                                  vector<double>& node,
                                  vector<double>& weight)
{
    //-------------------------------------------------------------------
    // Generate n Gauss-Legendre nodes and weights on [a,b].
    //-------------------------------------------------------------------

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

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

    node.assign(n, 0.0);
    weight.assign(n, 0.0);

    const double midpoint = (a+b)/2.0;
    const double halfwidth = (b-a)/2.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 root iteration did not converge.\n";
            return false;
        }

        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 w = 2.0/((1.0-z*z)*dp*dp);

        const int left  = i-1;
        const int right = n-i;

        node[left]  = midpoint - halfwidth*z;
        node[right] = midpoint + halfwidth*z;

        weight[left]  = halfwidth*w;
        weight[right] = halfwidth*w;
    }

    return true;
}