Interpolation in OCaml

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.