Lagrange Interpolation in Haskell

Here is a naive translation of the formula for a Lagrange Polynomial into Haskell:
interp xys x = sum $ zipWith (*) ys $ map f xs where
xs = map fst xys
ys = map snd xys
f xj = product $ map (p xj) xs
p xj xm = if xj == xm then 1 else (x - xm) / (xj - xm)

This is a cool little function- given a set of points, and a new point not in the set (list), it will give the value of that new point as interpolated by a Lagrange Polynomial- the polynomial of the least degree that is the value for y for each value x in the given list.

I've done symbolic regression with Gene Expression Programming, where you find a polynomial of arbitrary degree to fit a set of points, so it is cool to me to see such a simple function to compare against that is guaranteed to be of the lowest possible order.