//=======================================================================
// ch07_principal_value.cpp
//
// Cauchy principal value integral using analytical subtraction.
//
// The program evaluates
//
//                 1
//                 /
//                |       cos(x)
//      PV        |     --------- dx
//                |      x - 0.2
//                /
//                -1
//
// The pole is removed analytically:
//
//      f(x)/(x-x0)
//        = [f(x)-f(x0)]/(x-x0) + f(x0)/(x-x0).
//
// Therefore
//
//      PV integral
//        = integral [f(x)-f(x0)]/(x-x0) dx
//          + f(x0) log[(b-x0)/(x0-a)].
//
// The remaining regular integral is evaluated with Gauss-Legendre
// quadrature on [a,x0] and [x0,b].
//
// 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 principal_value(double (*fun)(double),
                       double a, double b, double x0, int n);

double gauss_legendre_regular(double (*fun)(double),
                              double a, double b,
                              double x0, double f0, int n);

int main()
{
    const double a  = -1.0;
    const double b  =  1.0;
    const double x0 =  0.2;

    const double reference = -0.5912784964342436;
    const int nvalues[] = {4, 8, 16, 32};

    cout << "Cauchy principal value integral\n";
    cout << "PV integral of cos(x)/(x-0.2) from -1 to 1\n";
    cout << scientific << setprecision(7);
    cout << "Reference value      = " << setw(15) << reference << "\n\n";

    cout << "        n       Principal value\n";
    cout << "--------------------------------\n";

    for (int n : nvalues) {
        const double result = principal_value(f, a, b, x0, n);

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

    return 0;
}


double f(double x)
{
    //-------------------------------------------------------------------
    // Numerator function in PV integral f(x)/(x-x0).
    //
    // To use another problem of the same form, replace only this line
    // and set a, b, and x0 in main().
    //-------------------------------------------------------------------
    return cos(x);
}


double principal_value(double (*fun)(double),
                       double a, double b, double x0, int n)
{
    //-------------------------------------------------------------------
    // Cauchy principal value using analytical subtraction.
    //-------------------------------------------------------------------

    if (x0 <= a || x0 >= b) {
        cerr << "principal_value: x0 must lie inside (a,b).\n";
        return 0.0;
    }

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

    const double f0 = fun(x0);

    const double left =
        gauss_legendre_regular(fun, a, x0, x0, f0, n);

    const double right =
        gauss_legendre_regular(fun, x0, b, x0, f0, n);

    const double singular_part =
        f0*log((b-x0)/(x0-a));

    return left + right + singular_part;
}


double gauss_legendre_regular(double (*fun)(double),
                              double a, double b,
                              double x0, double f0, int n)
{
    //-------------------------------------------------------------------
    // Gauss-Legendre quadrature for the regularized integrand
    //
    //          [f(x)-f(x0)]/(x-x0).
    //
    // In this application x0 is an endpoint of the subinterval, so it
    // is never sampled by Gauss-Legendre quadrature.
    //-------------------------------------------------------------------

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

    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_regular: Newton iteration failed.\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;

        const double gleft  = (fun(xleft)-f0)/(xleft-x0);
        const double gright = (fun(xright)-f0)/(xright-x0);

        if (fabs(z) <= tol) {
            integral += weight*(fun(midpoint)-f0)/(midpoint-x0);
        } else {
            integral += weight*(gleft+gright);
        }
    }

    return halfwidth*integral;
}