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


def f(x):
    """
    Function used to generate the interpolation data.
    """
    return math.sin(x)


def spline_coefficients(x, y):
    """
    Calculate the coefficients b[i], c[i], and d[i], i=0,1,...,n-1,
    for cubic spline interpolation

        s(x) = y[i] + b[i]*(x-x[i]) + c[i]*(x-x[i])**2
                    + d[i]*(x-x[i])**3

    for x[i] <= x <= x[i+1].

    Input:
        x = array of data abscissas in strictly increasing order
        y = array of data ordinates

    Output:
        b, c, d = arrays of spline coefficients

    Origin:
        Based on the SPLINE routine presented by
        G. E. Forsythe, M. A. Malcolm, and C. B. Moler,
        Computer Methods for Mathematical Computations,
        Prentice-Hall, 1977.

    This Python version is based on the Fortran implementation by
    Alexander Godunov.
    """

    n = len(x)

    b = np.zeros(n, dtype=float)
    c = np.zeros(n, dtype=float)
    d = np.zeros(n, dtype=float)

    gap = n - 1

    # Check input
    if n < 2:
        return b, c, d

    if n < 3:
        b[0] = (y[1] - y[0]) / (x[1] - x[0])   # linear interpolation
        c[0] = 0.0
        d[0] = 0.0
        b[1] = b[0]
        c[1] = 0.0
        d[1] = 0.0
        return b, c, d

    #
    # Step 1: preparation
    #
    d[0] = x[1] - x[0]
    c[1] = (y[1] - y[0]) / d[0]

    for i in range(1, gap):
        d[i] = x[i + 1] - x[i]
        b[i] = 2.0 * (d[i - 1] + d[i])
        c[i + 1] = (y[i + 1] - y[i]) / d[i]
        c[i] = c[i + 1] - c[i]

    #
    # Step 2: end conditions
    #
    b[0] = -d[0]
    b[n - 1] = -d[n - 2]
    c[0] = 0.0
    c[n - 1] = 0.0

    if n != 3:
        c[0] = c[2] / (x[3] - x[1]) - c[1] / (x[2] - x[0])
        c[n - 1] = (
            c[n - 2] / (x[n - 1] - x[n - 3])
            - c[n - 3] / (x[n - 2] - x[n - 4])
        )

        c[0] = c[0] * d[0] ** 2 / (x[3] - x[0])
        c[n - 1] = -c[n - 1] * d[n - 2] ** 2 / (x[n - 1] - x[n - 4])

    #
    # Step 3: forward elimination
    #
    for i in range(1, n):
        h = d[i - 1] / b[i - 1]
        b[i] = b[i] - h * d[i - 1]
        c[i] = c[i] - h * c[i - 1]

    #
    # Step 4: back substitution
    #
    c[n - 1] = c[n - 1] / b[n - 1]

    for j in range(1, gap + 1):
        i = n - j - 1
        c[i] = (c[i] - d[i] * c[i + 1]) / b[i]

    #
    # Step 5: compute spline coefficients
    #
    b[n - 1] = (
        (y[n - 1] - y[n - 2]) / d[n - 2]
        + d[n - 2] * (c[n - 2] + 2.0 * c[n - 1])
    )

    for i in range(gap):
        b[i] = (
            (y[i + 1] - y[i]) / d[i]
            - d[i] * (c[i + 1] + 2.0 * c[i])
        )

        d[i] = (c[i + 1] - c[i]) / d[i]
        c[i] = 3.0 * c[i]

    c[n - 1] = 3.0 * c[n - 1]
    d[n - 1] = d[n - 2]

    return b, c, d


def spline_eval(u, x, y, b, c, d):
    """
    Evaluate the cubic spline at point u:

        spline_eval = y[i] + b[i]*(u-x[i]) + c[i]*(u-x[i])**2
                           + d[i]*(u-x[i])**3

    where x[i] <= u <= x[i+1].

    Input:
        u       = abscissa at which the spline is evaluated
        x, y    = arrays of given data points
        b, c, d = spline coefficients computed by spline_coefficients

    Output:
        interpolated value at u

    Origin:
        Based on the SEVAL routine presented by
        G. E. Forsythe, M. A. Malcolm, and C. B. Moler,
        Computer Methods for Mathematical Computations,
        Prentice-Hall, 1977.

    This Python version is based on the Fortran implementation by
    Alexander Godunov.
    """

    n = len(x)

    # If u is outside the x interval, return the nearest boundary value.
    if u <= x[0]:
        return y[0]

    if u >= x[n - 1]:
        return y[n - 1]

    #
    # Binary search for i such that x[i] <= u <= x[i+1]
    #
    i = 0
    j = n

    while j > i + 1:
        k = (i + j) // 2

        if u < x[k]:
            j = k
        else:
            i = k

    #
    # Evaluate spline interpolation
    #
    dx = u - x[i]

    return y[i] + dx * (b[i] + dx * (c[i] + dx * d[i]))


def main():
    """
    Cubic spline interpolation

    Function values f(x) are calculated at n base points. The spline
    coefficients are then computed, and the spline is evaluated at
    nint points across the interval. The program also calculates the
    average absolute interpolation error.

    Origin:
        Based on the SPLINE and SEVAL routines presented by
        G. E. Forsythe, M. A. Malcolm, and C. B. Moler,
        Computer Methods for Mathematical Computations,
        Prentice-Hall, 1977.

    This Python version is based on the Fortran implementation by
    Alexander Godunov.
    Revised for the companion website, 2026.
    """

    n = 11       # base points for interpolation
    nint = 21    # compute interpolation in nint points

    xmin = 0.0
    xmax = 2.0

    xi = np.zeros(n, dtype=float)
    yi = np.zeros(n, dtype=float)

    # Step 1: generate xi and yi from f(x), xmin, xmax, and n
    step = (xmax - xmin) / (n - 1)

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

    # Step 2: calculate spline coefficients
    b, c, d = spline_coefficients(xi, yi)

    # Step 3: evaluate the spline at nint points
    errav = 0.0
    step = (xmax - xmin) / (nint - 1)

    xplot = np.zeros(nint, dtype=float)
    ysplot = np.zeros(nint, dtype=float)

    print(f"{'x':>12}{'spline':>12}{'error':>12}")

    for i in range(nint):
        x = xmin + step * i
        y = f(x)
        ys = spline_eval(x, xi, yi, b, c, d)
        error = ys - y

        print(f"{x:12.5f}{ys:12.5f}{error:12.5f}")

        # Step 4: calculate the average absolute interpolation error
        errav += abs(y - ys) / nint

        xplot[i] = x
        ysplot[i] = ys

    print(f"{'Average error':>24}{errav:12.5f}")

    # Plot the exact function, interpolation points, and cubic spline.
    xfine = np.linspace(xmin, xmax, 400)
    yfine = np.sin(xfine)

    plt.figure()
    plt.plot(xfine, yfine, "-", linewidth=1.2, label="Exact function")
    plt.plot(xplot, ysplot, "--", linewidth=1.2, label="Cubic spline")
    plt.plot(xi, yi, "o", markersize=6, label="Interpolation points")

    plt.xlabel("x")
    plt.ylabel("f(x)")
    plt.title("Cubic spline interpolation")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()