#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

double f(double x);
void spline(double x[], double y[], double b[], double c[], double d[], int n);
double spline_eval(double u, double x[], double y[],
                   double b[], double c[], double d[], int n);

int 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 C++ version is based on the Fortran implementation by
    // Alexander Godunov.
    // Revised for the companion website, 2026.
    //====================================================================

    const int n = 11;       // base points for interpolation
    const int nint = 21;    // compute interpolation in nint points

    double xi[n], yi[n], b[n], c[n], d[n];
    double xmin, xmax;
    double x, y, step, ys, error, errav;
    int i;

    xmin = 0.0;
    xmax = 2.0;

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

    for (i = 0; i < n; i++)
    {
        xi[i] = xmin + step*i;
        yi[i] = f(xi[i]);
    }

    // Step 2: calculate spline coefficients
    spline(xi, yi, b, c, d, n);

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

    cout << setw(12) << "x"
         << setw(12) << "spline"
         << setw(12) << "error" << endl;

    cout << fixed << setprecision(5);

    for (i = 0; i < nint; i++)
    {
        x = xmin + step*i;
        y = f(x);
        ys = spline_eval(x, xi, yi, b, c, d, n);
        error = ys - y;

        cout << setw(12) << x
             << setw(12) << ys
             << setw(12) << error << endl;

        // Step 4: calculate the average absolute interpolation error
        errav = errav + fabs(y - ys)/nint;
    }

    cout << setw(24) << "Average error"
         << setw(12) << errav << endl;

    return 0;
}


//
// Function f(x)
//
double f(double x)
{
    return sin(x);
}


void spline(double x[], double y[], double b[], double c[], double d[], int n)
/*======================================================================
  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
    n = number of data points (n >= 2)

  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 C++ version is based on the Fortran implementation by
  Alexander Godunov.
======================================================================*/
{
    int i, j, gap;
    double h;

    gap = n - 1;

    // Check input
    if (n < 2) return;

    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;
    }

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

    for (i = 1; i < gap; i++)
    {
        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]*d[0]/(x[3] - x[0]);
        c[n-1] = -c[n-1]*d[n-2]*d[n-2]/(x[n-1] - x[n-4]);
    }

    //
    // Step 3: forward elimination
    //
    for (i = 1; i < n; i++)
    {
        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 = 1; j <= gap; j++)
    {
        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 = 0; i < gap; i++)
    {
        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];
}


double spline_eval(double u, double x[], double y[],
                   double b[], double c[], double d[], int n)
/*======================================================================
  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
    n       = number of data points

  Output:
    spline_eval = 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 C++ version is based on the Fortran implementation by
  Alexander Godunov.
======================================================================*/
{
    int i, j, k;
    double dx;

    // 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]));
}