I just wrote some code to do polynomial interpolation in O'Caml. The full darcs repository is at http://web.mit.edu/farr/www/onumerics/; here's the guts of the code. Eventually (depending on my time and level of continued interest), this may turn into a more general O'Caml numerical library. As usual, it's released under the GPL.

(* interp.ml: Interpolation.
Copyright (C) 2008 Will M. Farr
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see .
*)
let interp xs ys x =
let n = Array.length xs in
assert(Array.length ys = n);
let ytab = Array.make n ys in
for i = 1 to n - 1 do
let m = n - i in
let last_col = ytab.(i-1) and
this_col = Array.make m 0.0 in
for j = 0 to m - 1 do
let xlow = xs.(j) and
xhigh = xs.(j+i) and
ylow = last_col.(j) and
yhigh = last_col.(j+1) in
this_col.(j) <- (x-.xlow)*.yhigh/.(xhigh-.xlow) +. (x-.xhigh)*.ylow/.(xlow-.xhigh)
done;
ytab.(i) <- this_col
done;
ytab.(n-1).(0)

Here's how it works. Initially, we have two arrays of points, `xs` and `ys`, and a point at which we want to evaluate the interpolated polynomial, `x`. This means we have an array of 0-order interpolating polynomials (the constants `ys`). We can use a linear interpolation between neighboring points to construct the values of first-order interpolating polynomials at `x`. Then we can linearly interpolate between the neighboring first-order polynomials to construct second-order interpolating polynomials. We iterate this process until we have a single value, which is the value of the n-point interpolating polynomial at `x`.
In the case of four points, here is a graphical representation of the procedure:

x0: y0
y01
x1: y1 y012
y12 y0123
x2: y2 y123
y23
x3: y3

The

`y`'s subscripts indicate the points involved in the interpolation. Each

`y` is a linear interpolation of the

`y` values above and below it to the left at the corresponding

`x` points.