"""
ch07_gauss_legendre.py

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.
Python version prepared for the companion website, 2026.

The numerical routine uses only the Python standard library.
Matplotlib is used only for visualization.
"""

import math
import matplotlib.pyplot as plt


def f(x):
    """Test function."""
    return math.sin(x)


def legendre_pair(n, z):
    """Return P_n(z) and P_(n-1)(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(fun, a, b, n):
    """n-point Gauss-Legendre quadrature on [a,b]."""
    tol = 1.0e-14
    max_iter = 100

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

    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: Newton iteration did not converge."
            )

        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

        if abs(z) <= tol:
            integral += weight * fun(midpoint)
        else:
            integral += weight * (fun(xleft) + fun(xright))

    return halfwidth * integral


def plot_integrand(a, b):
    nplot = 1000
    xplot = [a + i * (b - a) / nplot for i in range(nplot + 1)]
    yplot = [f(x) for x in xplot]

    plt.figure()
    plt.plot(xplot, yplot, linewidth=1.5)
    plt.xlabel("x")
    plt.ylabel("f(x)")
    plt.title("Integrand for Gauss-Legendre quadrature")
    plt.grid(True)


def main():
    a = 0.0
    b = math.pi
    exact = 2.0
    nvalues = [2, 4, 8, 16]

    plot_integrand(a, b)

    print("Gauss-Legendre quadrature")
    print("Integral of sin(x) from 0 to pi")
    print(f"Exact value          = {exact:15.7e}\n")

    print("        n       Gauss-Legendre")
    print("--------------------------------")

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

    plt.show()


if __name__ == "__main__":
    main()