content: #!/usr/bin/env runhaskell
  
  import Control.Monad (mapM_)
  import Text.Printf (printf)
  
  type Point = (Double,Double)
  type Interval = (Double,Double)
  
  breakAt :: Int -> [a] -> ([a], a, [a])
  breakAt k as = (take (k-1) as, head $ drop (k-1) as, drop k as)
  
  -- produces the kth Lagrange basis polynomial for the given xs
  basis :: Int -> [Double] -> (Double -> Double)
  basis k xs t = product (map factor xs')
    where (as, xk, bs) = breakAt k xs
          xs' = as ++ bs
          factor xj = (t - xj) / (xk - xj)
  
  -- produces the Lagrange polynomial for the given points
  interpolate :: [Point] -> (Double -> Double)
  interpolate ps x = 
    let (xs,ys) = unzip ps 
     in sum . map ($ x) $ do
       (y,j) <- ys `zip` [1..]
       let b = basis j xs
       return $ \x -> y * b x
  
  -- produces a list of n+1 values from a to b with consecutive values
  -- (b-a)/n apart
  uniform :: Interval -> Int -> [Double]
  uniform (a,b) n = [a,(a+d)..b] where d = (b - a) / (fromIntegral n)
  
  -- produces a list of n+1 values from a to b with the distance between
  -- consecutive values determined by Chebyshev's criterion
  chebyshev :: Interval -> Int -> [Double]
  chebyshev (a,b) n = let c0 = (a+b) / 2
                          c1 = (a-b) / 2
                          n' = fromIntegral n
                       in map (\i -> c0 + c1 * cos(pi/n'*(i-1))) [1..n'+1]
  
  -- produces the Lagrange polynomial approximation of the given
  -- function f with samples drawn from the given xs
  approximate :: (Double -> Double) -> [Double] -> (Double -> Double)
  approximate f xs = \x -> interpolate (xs `zip` map f xs) x
  
  pr (x,y) = printf "%10.3f\t%10.3f\n" x y
  
  main = do
    putStrLn "Interpolate a set of points:"
    putStrLn "\tx\t\ty"
    let points = [(-3,7),(4, -1),(6,2),(8,1)]
        xs = [-10,-9.5..10]
        l = interpolate points
    mapM_ pr (xs `zip` map l xs)
    
    putStrLn "Approximate a function with a polynomial:"
    let f x = abs (cos (2*x))
    let xs' = [-pi,-pi+0.1..pi]
        n = 40
        interval = (-pi,pi)
    putStrLn "using uniform nodes:"
    putStrLn "\tx\t\ty"
    let uniformApproximation = approximate f (uniform interval n)
    mapM_ pr (xs' `zip` map uniformApproximation xs')
    putStrLn "using Chebyshev nodes:"
    putStrLn "\tx\t\ty"
    let chebyshevApproximation = approximate f (chebyshev interval n)
    mapM_ pr (xs' `zip` map chebyshevApproximation xs')
  