import math
import numpy as np
import matplotlib.pyplot as plt


def f(x):
    """
    Function f(x): Runge function.
    """
    return 1.0 / (1.0 + 25.0 * x * x)


def lagrange_global(xx, xi, yi):
    """
    Global polynomial interpolation in Lagrange form.

    All data points are used to construct and evaluate the
    interpolating polynomial.

    Input:
        xx = abscissa at which the interpolation is evaluated
        xi = array of data abscissas
        yi = array of data ordinates

    Output:
        interpolated value at xx

    Comments:
        The routine works for both equally and unequally spaced points.
        Unlike the local interpolation routine used elsewhere in this
        chapter, all available points are used here so that the effect
        of node placement on a high-order global polynomial is visible.

    Alexander Godunov
    Prepared for the companion website, 2026.
    """

    n = len(xi)
    y = 0.0

    for i in range(n):
        lam = 1.0

        for j in range(n):
            if j != i:
                lam *= (xx - xi[j]) / (xi[i] - xi[j])

        y += yi[i] * lam

    return y


def main():
    """
    Runge phenomenon and Chebyshev nodes

    The Runge function f(x) is sampled using two sets of n base points:
    equally spaced points and Chebyshev points. A global polynomial
    interpolant is then evaluated for both cases at nint points across
    the interval. The program also calculates the average absolute
    interpolation error for each choice of nodes.

    The interpolating polynomial is evaluated in Lagrange form.

    This example illustrates that increasing the polynomial order with
    equally spaced points can produce large oscillations near the ends
    of the interval, while Chebyshev nodes greatly reduce this effect.

    Alexander Godunov
    Prepared for the companion website, 2026.
    """

    n = 13       # base points for interpolation
    nint = 25    # compute interpolation in nint points

    xmin = -1.0
    xmax = 1.0

    xe = np.zeros(n, dtype=float)
    ye = np.zeros(n, dtype=float)
    xc = np.zeros(n, dtype=float)
    yc = np.zeros(n, dtype=float)

    # Step 1: generate equally spaced interpolation points
    step = (xmax - xmin) / (n - 1)

    for i in range(n):
        xe[i] = xmin + step * i
        ye[i] = f(xe[i])

    # Step 2: generate Chebyshev interpolation points
    # The nodes are written in increasing order from left to right.
    for i in range(n):
        xc[i] = math.cos(
            (2.0 * (n - i - 1) + 1.0) * math.pi / (2.0 * n)
        )
        yc[i] = f(xc[i])

    # Step 3: evaluate both global interpolants at nint points
    errav_e = 0.0
    errav_c = 0.0
    step = (xmax - xmin) / (nint - 1)

    print("     Runge Function: Equally Spaced and Chebyshev Nodes")
    print(f"     number of interpolation points = {n:2d}")
    print(
        f"{'x':>13}"
        f"{'equal poly':>13}"
        f"{'equal error':>13}"
        f"{'Cheb poly':>13}"
        f"{'Cheb error':>13}"
    )

    for i in range(nint):
        x = xmin + step * i
        y = f(x)

        yse = lagrange_global(x, xe, ye)
        ysc = lagrange_global(x, xc, yc)

        erre = yse - y
        errc = ysc - y

        print(
            f"{x:13.6f}"
            f"{yse:13.6f}"
            f"{erre:13.6f}"
            f"{ysc:13.6f}"
            f"{errc:13.6f}"
        )

        # Step 4: calculate the average absolute interpolation errors
        errav_e += abs(erre) / nint
        errav_c += abs(errc) / nint

    print()
    print(" Average absolute error:")
    print(f"   equally spaced = {errav_e:12.6f}")
    print(f"   Chebyshev      = {errav_c:12.6f}")

    # Plot the Runge function and both interpolants on a fine grid.
    xfine = np.linspace(xmin, xmax, 600)
    yfine = 1.0 / (1.0 + 25.0 * xfine * xfine)

    ye_fine = np.zeros_like(xfine)
    yc_fine = np.zeros_like(xfine)

    for i in range(len(xfine)):
        ye_fine[i] = lagrange_global(xfine[i], xe, ye)
        yc_fine[i] = lagrange_global(xfine[i], xc, yc)

    plt.figure()
    plt.plot(xfine, yfine, "-", linewidth=1.2, label="Exact function")
    plt.plot(
        xfine,
        ye_fine,
        "--",
        linewidth=1.2,
        label="Equally spaced interpolation",
    )
    plt.plot(
        xfine,
        yc_fine,
        "-.",
        linewidth=1.2,
        label="Chebyshev interpolation",
    )
    plt.plot(xe, ye, "o", markersize=5, label="Equally spaced nodes")
    plt.plot(xc, yc, "s", markersize=5, label="Chebyshev nodes")

    plt.xlabel("x")
    plt.ylabel("f(x)")
    plt.title("Runge function: equally spaced and Chebyshev nodes")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()