```-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Approximation.Chebyshev
--
-- Stability   :  experimental
-- Portability :  portable
--
-- Function approximation using Chebyshev polynomials
--
-- @ f(x) = ( sum (k=0..N-1) c_k * T_k(x) ) - 0.5 * c_0 @
--
-- over the interval @ [a,b] @
--
-- Reference: NRiC
--
-----------------------------------------------------------------------------

module Numeric.Approximation.Chebyshev (cheby_approx,
cheby_eval) where

import Data.Array

-- | Calculates the Chebyshev approximation to @f(x)@ over @[a,b]@

cheby_approx :: (Double -> Double) -- ^ f(x)
-> Double             -- ^ a
-> Double             -- ^ b
-> Int                -- ^ N
-> [Double]           -- ^ c_n

cheby_approx f a b n = f''
where a' = 0.5 * (b - a)
b' = 0.5 * (b + a)
y = [ a' * cos (pi * (fromIntegral k + 0.5) / fromIntegral n) + b' | k <- [0..n-1] ]
f' = map f y
f'' = [ 2 * sum (zipWith (*) f' [ cos (pi * fromIntegral j * (fromIntegral k + 0.5) / fromIntegral n) | k <- [0..n-1] ]) / fromIntegral n | j <- [0..n-1] ]

-- | Evaluates the Chebyshev approximation to @f(x)@ over @[a,b]@ at @x@

cheby_eval :: [Double] -- ^ c_n
-> Double   -- ^ a
-> Double   -- ^ b
-> Double   -- ^ x
-> Double   -- ^ f(x)

cheby_eval f a b x = y * d!1 - d!2 + 0.5 * c!0
where y = (2 * x - a - b) / (b - a)
c = listArray (0,n) f
d = array (1,n+2) ((n+2,0) : (n+1,0) : [ (j,2*y*d!(j+1) - d!(j+2) + c!j) | j <- [1..n] ])
n = length f - 1
```