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

using namespace std;

double f(double x);
double lagrange_global(double xx, const double xi[], const double yi[], int n);

int 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.
    //====================================================================

    const int n = 13;       // base points for interpolation
    const int nint = 25;    // compute interpolation in nint points
    const double pi = 3.14159265358979323846;

    double xe[n], ye[n], xc[n], yc[n];
    double xmin, xmax;
    double x, y, step;
    double yse, ysc, erre, errc;
    double errav_e, errav_c;
    int i;

    xmin = -1.0;
    xmax =  1.0;

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

    for (i = 0; i < n; i++)
    {
        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 = 0; i < n; i++)
    {
        xc[i] = cos((2.0*(n-i-1) + 1.0)*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);

    cout << "     Runge Function: Equally Spaced and Chebyshev Nodes" << endl;
    cout << "     number of interpolation points = "
         << setw(2) << n << endl;

    cout << setw(13) << "x"
         << setw(13) << "equal poly"
         << setw(13) << "equal error"
         << setw(13) << "Cheb poly"
         << setw(13) << "Cheb error" << endl;

    cout << fixed << setprecision(6);

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

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

        erre = yse - y;
        errc = ysc - y;

        cout << setw(13) << x
             << setw(13) << yse
             << setw(13) << erre
             << setw(13) << ysc
             << setw(13) << errc << endl;

        // Step 4: calculate the average absolute interpolation errors
        errav_e += fabs(erre)/nint;
        errav_c += fabs(errc)/nint;
    }

    cout << endl;
    cout << " Average absolute error:" << endl;
    cout << "   equally spaced = " << setw(12) << errav_e << endl;
    cout << "   Chebyshev      = " << setw(12) << errav_c << endl;

    return 0;
}


//
// Function f(x): Runge function
//
double f(double x)
{
    return 1.0/(1.0 + 25.0*x*x);
}


double lagrange_global(double xx, const double xi[],
                       const double yi[], int n)
/*====================================================================
  Global polynomial interpolation in Lagrange form.

  All n 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
    n  = number of interpolation points

  Output:
    lagrange_global = 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.
====================================================================*/
{
    double lambda, y;
    int i, j;

    y = 0.0;

    for (i = 0; i < n; i++)
    {
        lambda = 1.0;

        for (j = 0; j < n; j++)
        {
            if (j != i)
            {
                lambda = lambda*(xx - xi[j])/(xi[i] - xi[j]);
            }
        }

        y += yi[i]*lambda;
    }

    return y;
}