#ifndef _utl_PolynomialInterpolation_h_
#define _utl_PolynomialInterpolation_h_

#include <cmath>
#include <utl/ErrorLogger.h>


namespace utl {

  /** \brief Perform polynomial interpolation or extrapolation

    Neville's algorithm.

    With x = { -1, 0, 1, 2, 3 } and y = { 0, 1, 0, 2, -1 } the following results
    are obtained for interpolation and extrapolation with different n (see image):

    \image html utl_PolynomialInterpolation.png

    \author Darko Veberic
    \date 28 Apr 2007
    \ingroup math

    \param n: number of points to take from the x and y array
    \param px: array of x coordinates of the interpolating points
    \param py: array of y coordinates of the interpolating points
    \param x: polynomial evaluation point
    \param dy: returns the accuracy of the approximation (used in utl::Integrator)
    \return interpolated value for y in point x
  */

  inline
  double
  PolynomialInterpolation(const unsigned int n, const double px[], const double py[],
                          const double x, double& dy)
  {
    if (!n)
      return 0;

    double minDist = std::fabs(x - px[0]);
    unsigned int k = 0;
    double c[n];
    double d[n];
    c[0] = py[0];
    d[0] = py[0];
    double oldx = px[0];
    for (unsigned int i = 1; i < n; ++i) {
      const double newx = px[i];
      if (oldx == newx) {
        ERROR("x values should be distinct.");
        return 0;
      }
      oldx = newx;
      const double dist = std::fabs(x - newx);
      if (dist < minDist) {
        minDist = dist;
        k = i;
      }
      c[i] = d[i] = py[i];
    }

    double y = py[k];

    for (unsigned int m = 1; m < n; ++m) {
      const unsigned int nm = n - m;
      for (unsigned int i = 0; i < nm; ++i) {
        const double xa = px[i];
        const double xb = px[i+m];
        const double dx = xb - xa;
        const double cd = (c[i+1] - d[i]) / dx;
        c[i] = (x-xa)*cd;
        d[i] = (x-xb)*cd;
      }
      dy = (2*k < nm) ? c[k] : d[--k];
      y += dy;
    }

    return y;
  }


  /** \brief Perform polynomial interpolation or extrapolation with approximation control

    \author Darko Veberic
    \date 28 Apr 2007
    \ingroup math

    \param n: number of points to take from the px and py array
    \param px: array of x coordinates of the interpolating points
    \param py: array of y coordinates of the interpolating points
    \param x: polynomial evaluation point
  */

  inline
  double
  PolynomialInterpolation(const unsigned int n, const double px[], const double py[],
                          const double x)
  {
    double dummy;
    return PolynomialInterpolation(n, px, py, x, dummy);
  }

}

#endif