"""
ch07_principal_value.py

Cauchy principal value integral using analytical subtraction.

The program evaluates

                1
                /
               |       cos(x)
     PV        |     --------- dx
               |      x - 0.2
               /
               -1

The singularity is removed analytically:

    f(x)/(x-x0)
      = [f(x)-f(x0)]/(x-x0) + f(x0)/(x-x0).

Hence

    PV integral
      = integral [f(x)-f(x0)]/(x-x0) dx
        + f(x0) log[(b-x0)/(x0-a)].

The regular part is evaluated with Gauss-Legendre quadrature on
[a,x0] and [x0,b].

Based on the Fortran implementation by Alexander Godunov.
Python version prepared for the companion website, 2026.

Matplotlib is used only for visualization.
"""

import math
import matplotlib.pyplot as plt


def f(x):
    """Numerator function in f(x)/(x-x0)."""
    return math.cos(x)


def legendre_pair(n, z):
    if n == 1:
        return z, 1.0

    p0 = 1.0
    p1 = z

    for j in range(2, n + 1):
        p2 = ((2.0 * j - 1.0) * z * p1 - (j - 1.0) * p0) / j
        p0 = p1
        p1 = p2

    return p1, p0


def gauss_legendre_regular(fun, a, b, x0, f0, n):
    """Gauss-Legendre quadrature for [f(x)-f(x0)]/(x-x0)."""
    tol = 1.0e-14
    max_iter = 100

    midpoint = 0.5 * (a + b)
    halfwidth = 0.5 * (b - a)
    integral = 0.0
    m = (n + 1) // 2

    for i in range(1, m + 1):
        z = math.cos(math.pi * (i - 0.25) / (n + 0.5))
        converged = False

        for _ in range(max_iter):
            pn, pnm1 = legendre_pair(n, z)
            dp = 1.0 if n == 1 else n * (z * pn - pnm1) / (z * z - 1.0)

            zold = z
            z = zold - pn / dp

            if abs(z - zold) <= tol:
                converged = True
                break

        if not converged:
            raise RuntimeError(
                "gauss_legendre_regular: Newton iteration failed."
            )

        pn, pnm1 = legendre_pair(n, z)
        dp = 1.0 if n == 1 else n * (z * pn - pnm1) / (z * z - 1.0)
        weight = 2.0 / ((1.0 - z * z) * dp * dp)

        xleft = midpoint - halfwidth * z
        xright = midpoint + halfwidth * z

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

        if abs(z) <= tol:
            integral += weight * (fun(midpoint) - f0) / (midpoint - x0)
        else:
            integral += weight * (gleft + gright)

    return halfwidth * integral


def principal_value(fun, a, b, x0, n):
    """Cauchy principal value using analytical subtraction."""
    if x0 <= a or x0 >= b:
        raise ValueError("principal_value: x0 must lie inside (a,b).")

    if n < 1:
        raise ValueError("principal_value: n must be positive.")

    f0 = fun(x0)

    left = gauss_legendre_regular(fun, a, x0, x0, f0, n)
    right = gauss_legendre_regular(fun, x0, b, x0, f0, n)

    singular_part = f0 * math.log((b - x0) / (x0 - a))

    return left + right + singular_part


def regularized_value(x, x0):
    """Regularized integrand, including its removable limit at x=x0."""
    if abs(x - x0) < 1.0e-12:
        return -math.sin(x0)
    return (f(x) - f(x0)) / (x - x0)


def plot_integrands(a, b, x0):
    # Original singular integrand, plotted on both sides of the pole.
    gap = 0.01
    nplot = 1000

    xleft = [a + i * (x0 - gap - a) / nplot for i in range(nplot + 1)]
    xright = [
        x0 + gap + i * (b - x0 - gap) / nplot
        for i in range(nplot + 1)
    ]

    yleft = [f(x) / (x - x0) for x in xleft]
    yright = [f(x) / (x - x0) for x in xright]

    plt.figure()
    plt.plot(xleft, yleft, linewidth=1.5)
    plt.plot(xright, yright, linewidth=1.5)
    plt.xlabel("x")
    plt.ylabel("f(x)/(x-x0)")
    plt.title("Original integrand with a pole at x0 = 0.2")
    plt.grid(True)

    # Regularized integrand after analytical subtraction.
    xplot = [a + i * (b - a) / (2 * nplot) for i in range(2 * nplot + 1)]
    gplot = [regularized_value(x, x0) for x in xplot]

    plt.figure()
    plt.plot(xplot, gplot, linewidth=1.5)
    plt.xlabel("x")
    plt.ylabel("[f(x)-f(x0)]/(x-x0)")
    plt.title("Regularized integrand after analytical subtraction")
    plt.grid(True)


def main():
    a = -1.0
    b = 1.0
    x0 = 0.2

    reference = -0.5912784964342436
    nvalues = [4, 8, 16, 32]

    plot_integrands(a, b, x0)

    print("Cauchy principal value integral")
    print("PV integral of cos(x)/(x-0.2) from -1 to 1")
    print(f"Reference value      = {reference:15.7e}\n")

    print("        n       Principal value")
    print("--------------------------------")

    for n in nvalues:
        result = principal_value(f, a, b, x0, n)
        print(f"{n:9d}{result:21.7e}")

    plt.show()


if __name__ == "__main__":
    main()