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

using namespace std;

double f(double x);
double polynomial_interp(double xx, const double xi[], const double yi[],
                         int ni, int npts);

int main()
{
    //====================================================================
    // General polynomial interpolation
    //
    // Function values f(x) are calculated at n base points. A local
    // polynomial of selected order is then used to interpolate the
    // function at nint points across the interval. The program also
    // calculates the average absolute interpolation error.
    //
    // The interpolating polynomial is evaluated in Lagrange form.
    // For order = 1, 2, 3, ... the routine uses 2, 3, 4, ...
    // neighboring data points, respectively.
    //
    // Original implementation by Alexander Godunov, January 2010.
    // 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];
    double xmin, xmax;
    double x, y, step, ys, error, errav;
    int i, order;

    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: choose the interpolation order
    // order = 1 -> linear
    // order = 2 -> quadratic
    // order = 3 -> cubic
    order = 3;

    cout << "    General Polynomial Interpolation" << endl;
    cout << "     order of interpolation = " << setw(2) << order << endl;

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

    cout << setw(12) << "x"
         << setw(16) << "interpolated"
         << setw(12) << "error" << endl;

    cout << fixed << setprecision(5);

    for (i = 0; i < nint; i++)
    {
        x = xmin + step*i;
        y = f(x);

        ys = polynomial_interp(x, xi, yi, n, order + 1);

        error = ys - y;

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

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

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

    return 0;
}


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


double polynomial_interp(double xx, const double xi[], const double yi[],
                         int ni, int npts)
/*====================================================================
  Local polynomial interpolation of selected order.

  The interpolating polynomial is evaluated in Lagrange form using
  npts neighboring data points. The interpolation order is npts-1.

  Input:
    xx   = abscissa at which the interpolation is evaluated
    xi   = array of data abscissas
    yi   = array of data ordinates
    ni   = number of data points in xi and yi
    npts = number of points used for interpolation
           (interpolation order = npts-1)

  Output:
    polynomial_interp = interpolated value at xx

  Comments:
    The routine works for both equally and unequally spaced xi.
    A binary search locates the interval containing xx, and a local
    group of neighboring data points is selected for interpolation.
    If xx lies outside the data interval, the nearest boundary value
    is returned.

  Original implementation by Alexander Godunov, January 2010.
  Revised for the companion website, 2026.
====================================================================*/
{
    int i, j, k, js, jl, nuse;
    double y;
    vector<double> lambda(ni);

    // Check the number of interpolation points
    nuse = npts;
    if (nuse > ni) nuse = ni;

    // If xx is outside the xi interval, return a boundary value
    if (xx <= xi[0]) return yi[0];
    if (xx >= xi[ni-1]) return yi[ni-1];

    // Binary search to find i such that xi[i] < xx < xi[i+1]
    i = 0;
    j = ni - 1;

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

        if (xx < xi[k])
            j = k;
        else
            i = k;
    }

    // Shift i so that xx lies near the middle of the selected data points
    i = i + 1 - nuse/2;

    // Keep the interpolation points inside the available data range
    if (i < 0) i = 0;
    if (i + nuse - 1 > ni - 1) i = ni - nuse;

    // Evaluate the interpolating polynomial in Lagrange form
    y = 0.0;

    for (js = i; js < i + nuse; js++)
    {
        lambda[js] = 1.0;

        for (jl = i; jl < i + nuse; jl++)
        {
            if (jl != js)
            {
                lambda[js] = lambda[js] *
                    (xx - xi[jl])/(xi[js] - xi[jl]);
            }
        }

        y = y + yi[js]*lambda[js];
    }

    return y;
}