Chebyshev approximation of arbitrary functions

Chebyshev.hs

import Math.Polynomial
import Math.Polynomial.Chebyshev

coeff f n j = (2 / nf) * (sum [f(xk) * yk | (xk, yk) <- zip x y])
    where
        x = map ((\x -> cos (pi * (x + 0.5) / nf)) . fromIntegral) [0..n-1]
        y = map ((\x -> cos (jf * pi * (x + 0.5) / nf)) . fromIntegral) [0..n-1]
        nf = fromIntegral n
        jf = fromIntegral j

approx f n = foldl1 addPoly $ (constPoly (-0.5 * c0)) : [scalePoly ck (t k) | (ck, k) <- zip c [0..n-1]]
    where
        c  = map ((coeff f n) . fromIntegral) [0..n-1]
        c0 = head c