#include <float.h>
#include <math.h>
#include <stdlib.h>
#include "interpolate.h"
/* returns value of f(x) at X, or DBL_MAX on errors */
double newton
(
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 tolerance, /* delta for successive pk(x) at which */
/* to stop */
int *degree /* returned -- degree of interpolating */
/* polynomial */
)
{
double *table, delta, m = 1.0, p;
int i, k;
/*
** Allocate enough storage for the longest diagonal required
** in the divided-diffrence table.
*/
if ((table = (double *)malloc(sizeof(double) * n)) == NULL)
return DBL_MAX;
/*
** Set the initial interpolating polynomial and first element
** of the divided-difference table to f(x0).
*/
p = table[0] = fx[0];
*degree = 0;
for (k = 0; k < n - 1; k++) {
(*degree)++;
table[k + 1] = fx[k + 1];
/*
** Compute the next diagonal placed in the
** divided-difference table to obtain the next coefficient.
*/
for (i = k; i >= 0; i--)
table[i] = (table[i+1] - table[i]) / (x[k+1] - x[i]);
/*
** Compute the next pk(x), one degree larger than the
** last, and determine whether the specified tolerance will
** be met.
*/
m = m * (X - x[k]);
p = p + (delta = table[0] * m);
if (fabs(delta) < tolerance) {
free(table);
return p;
}
}
/*
** The specified tolerance was never met for differences
** between successive interpolating polynomials.
*/
free(table);
return DBL_MAX;
}