#include "interpolate.h"
/* returns value of f(x) at X */
double lagrange
(
const double *x, /* array of interpolation points */
const double *fx, /* array of known values of f(x) */
int n, /* number of values in x and fx */
double X /* value of x at which to compute f(x) */
)
{
double l, p = 0.0;
int i, k;
for (k = 0; k < n; k++) {
l = 1.0;
/*
** Obtain lk(x) for the next term to add to the
** interpolating polynomial.
*/
for (i = 0; i < n; i++) {
if (i != k)
l = l * ((X - x[i]) / (x[k] - x[i]));
}
/*
** Add the next term computed from lk(x) to the
** interpolating polynomial.
*/
p = p + (fx[k] * l);
}
return p;
}