neville.c

大师写的二代小波经典之作

/*
 *  -*- Mode: ANSI C -*-
 *  (Polynomial Interpolation)
 *  $Id: neville.c,v 1.5 1996/09/17 16:10:41 fernande Exp $
 *  $Source: /sgi.acct/sweldens/cvs/liftpack/Lifting/neville.c,v $
 *  Author: Wim Sweldens, Gabriel Fernandez
 *
 *  Given n points of the form (x,f), this program uses the Neville algorithm
 *  to find the value at xx of the polynomial of degree (n-1) interpolating
 *  the points (x,f).
 *  Ref: Stoer and Bulirsch, Introduction to Numerical Analysis,
 *  Springer-Verlag.
 */
/* do not edit anything above this line */

/* FLWT header files */
#include <flwtdef.h>   /* Flt and Vector definitions */
#include <flwterr.h>   /* Error() function */
#include <mem.h>       /* vector() and free_vector() functions */
#include <neville.h>   /* external declaration of Neville() */
#include <util.h>      /* zero() function */

Flt
Neville ( const Flt *x, const Flt *f, const int n, const Flt xx )
{
    register int i,j;
    Vector vy;
    Flt y;

    vy = vector( 0, (long)(n-1) );   /* allocate a temporary vector */
    for ( i=0; i<n; i++ ) {
        vy[i] = f[i];
        for ( j=i-1; j>=0; j--) {
            Flt den = x[i] - x[j];
            if ( zero(den) ) {
                Error ("Neville", DENOM_NEAR_ZERO, RETURN);
            }
            vy[j] = vy[j+1] + (vy[j+1] - vy[j]) * (xx - x[i]) / den;
        }
    }
    y = vy[0];

    /* Free memory */
    free_vector( vy, 0, (long)(n-1) );

    return y;
}