pppack a FORTRAN77 code which carries out interpolation and approximation using piecewise polynomials, especially cubic splines, by Carl de Boor.
An original, true, correct version of PPPACK is available at https://www.netlib.org/pppack/index.html Files shown here are modifications made in pursuance of my own interests and needs, and should not be considered in preference to the original versions.
Typically, a set of data ( X(I), Y(I) ) for I=1, L+1 is available which is to be interpolated. A function F(X) is to be found which passes through the given data. F(X) is to be constructed from some order of polynomials, often cubic, and is to be continuously differentiable to all orders except at certain 'break' points, most likely at the same X points given with the data. Thus, to determine the types of functions F we can construct, it is necessary to specify the number of break points L+1, or intervals L, the location of the breakpoints 'BREAK', the values that the function is to assume at the breakpoints 'Y' or some other condition, and the order 'K' of the polynomial pieces. Sometimes auxilliary information, such as the slope or second derivative of the function at the left and right endpoints is part of the prescription.
Given a set of interpolating conditions or other requirements, the piecewise polynomial routines that follow will produce or manipulate a representation of the function F which has the form:
( BREAK, PCOEF, K, L )
These quantities represent the function F in the following way. If a point X is between breakpoints BREAK(I) and BREAK(I+1), then set
H = X - BREAK(I),
and we have
F(X)= PCOEF(1,I)
+ PCOEF(2,I) * H
+ PCOEF(3,I) * H^2 / 2!
+ PCOEF(4,I) * H^3 / 3!
...
+ PCOEF(K,I) * H^(K-1) /(K-1)!
Note that the piecewise polynomial functions make F and its derivatives continuous from the right. Thus at the break point I+1, we use the definition of F appropriate for the interval ( BREAK(I+1), BREAK(I+2) ) and not ( BREAK(I), BREAK(I+1) ).
Note also that the behavior of the function F for X values below the first breakpoint BREAK(1) or above the last breakpoint BREAK(L+1) is not specified. In fact, generally, F is set to zero below BREAK(1), and the definition of F on the last interval ( BREAK(L), BREAK(L+1) ) is extended all the way to the right.
Whenever you have a piecewise polynomial representation of the above form, you can evaluate the function F(X) by calling the function PPVALU. Moreover, other routines like BSPLPP can convert a B-spline representation into a piecewise polynomial representation. Also, KNOTS can use the information in the breakpoint sequence BREAK, with the continuity conditions required, to construct an equivalent knot sequence for a B-spline.
The computer code and data files described and made available on this web page are distributed under the GNU LGPL license.
pppack is available in a FORTRAN77 version and a FORTRAN90 version.
bernstein_polynomial, a FORTRAN77 library which evaluates the Bernstein polynomials, useful for uniform approximation of functions;
DIVDIF, a FORTRAN77 library which computes interpolants by divided differences.
HERMITE, a FORTRAN77 library which computes the Hermite interpolant, a polynomial that matches function values and derivatives.
SLATEC, a FORTRAN90 library which includes PPPACK.
SPLINE, a FORTRAN77 library which includes many routines to construct and evaluate spline interpolants and approximants.
TEST_APPROX, a FORTRAN77 library which defines a number of test problems for approximation and interpolation.
Carl de Boor