Here's my favorite way to fit a sequence to a polynomial, in working Python. (As I mentioned before, this is O(2^len(ys)), but can be cured of that.) def fit(ys, x): """ The x at the first of the ys is 0. """ d = len(ys) - 1 # (max x) - (min x) if d == 0: return ys[0] else: left = fit(ys[:-1], x) # Up to the 2nd-last y. right = fit(ys[1:], x - 1) # Start with the 2nd y. return ((d - x) * left + x * right) / d If x is a number, the function just processes numbers, as in any sensible language: ys = [0, 4, 16] fit(ys, -1 ) => 4 fit(ys, 0 ) => 0 fit(ys, 0.5 ) => 1.0 fit(ys, 1 ) => 4 fit(ys, 1.5 ) => 9.0 fit(ys, 2 ) => 16 fit(ys, 2.5 ) => 25.0 But if the language happens to let you define objects that act like polynomials, then the same code processes polynomials: X = UniPoly([0, 1]) P = fit(ys, X) => UniPoly([0, 0, 4]) P( -1 ) => 4 P( 0 ) => 0 P( 0.5 ) => 1 P( 1 ) => 4 P( 1.5 ) => 9 P( 2 ) => 16 P( 2.5 ) => 25 --Steve