//=======================================================================
// ch07_improper_infinite.cpp
//
// Improper integral with an infinite upper limit.
//
// The program evaluates
//
//              infinity
//                 /
//                |  exp(-x) sin(x) dx = 1/2
//                /
//                0
//
// by transforming [0,infinity) to [0,1):
//
//        x = t/(1-t),        dx = dt/(1-t)^2 .
//
// Gauss-Legendre quadrature is then applied on [0,1]. Because
// Gauss-Legendre nodes do not include the endpoints, t = 1 is never
// sampled.
//
// Based on the Fortran implementation by Alexander Godunov.
// C++ version prepared for the companion website, 2026.
//=======================================================================

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

using namespace std;

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

int main()
{
    const double a = 0.0;
    const double b = 1.0;
    const double exact = 0.5;
    const int nvalues[] = {8, 16, 32, 64};

    cout << "Improper integral: infinite interval\n";
    cout << "Integral of exp(-x) sin(x) from 0 to infinity\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(g, a, b, n);

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

    return 0;
}


double f(double x)
{
    //-------------------------------------------------------------------
    // Original integrand on the semi-infinite interval.
    //
    // To use another problem, replace only the line below.
    //-------------------------------------------------------------------
    return exp(-x)*sin(x);
}


double g(double t)
{
    //-------------------------------------------------------------------
    // Transformed integrand after x = t/(1-t).
    //-------------------------------------------------------------------

    const double one_minus_t = 1.0 - t;
    const double x = t/one_minus_t;

    return f(x)/(one_minus_t*one_minus_t);
}


double gauss_legendre(double (*fun)(double),
                      double a, double b, int n)
{
    //-------------------------------------------------------------------
    // n-point Gauss-Legendre quadrature 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: 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;
        }

        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;
}