-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Approximation.Chebyshev
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- 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 :: (Double -> Double) -> Double -> Double -> Int -> [Double]
cheby_approx Double -> Double
f Double
a Double
b Int
n = [Double]
f''
    where a' :: Double
a' = Double
0.5 forall a. Num a => a -> a -> a
* (Double
b forall a. Num a => a -> a -> a
- Double
a)
	  b' :: Double
b' = Double
0.5 forall a. Num a => a -> a -> a
* (Double
b forall a. Num a => a -> a -> a
+ Double
a)
	  y :: [Double]
y = [ Double
a' forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos (forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k forall a. Num a => a -> a -> a
+ Double
0.5) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) forall a. Num a => a -> a -> a
+ Double
b' | Int
k <- [Int
0..Int
nforall a. Num a => a -> a -> a
-Int
1] ]
	  f' :: [Double]
f' = forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
f [Double]
y
	  f'' :: [Double]
f'' = [ Double
2 forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) [Double]
f' [ forall a. Floating a => a -> a
cos (forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
j forall a. Num a => a -> a -> a
* (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k forall a. Num a => a -> a -> a
+ Double
0.5) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) | Int
k <- [Int
0..Int
nforall a. Num a => a -> a -> a
-Int
1] ]) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n | Int
j <- [Int
0..Int
nforall a. Num a => a -> a -> a
-Int
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 :: [Double] -> Double -> Double -> Double -> Double
cheby_eval [Double]
f Double
a Double
b Double
x = Double
y forall a. Num a => a -> a -> a
* Array Int Double
dforall i e. Ix i => Array i e -> i -> e
!Int
1 forall a. Num a => a -> a -> a
- Array Int Double
dforall i e. Ix i => Array i e -> i -> e
!Int
2 forall a. Num a => a -> a -> a
+ Double
0.5 forall a. Num a => a -> a -> a
* Array Int Double
cforall i e. Ix i => Array i e -> i -> e
!Int
0
    where y :: Double
y = (Double
2 forall a. Num a => a -> a -> a
* Double
x forall a. Num a => a -> a -> a
- Double
a forall a. Num a => a -> a -> a
- Double
b) forall a. Fractional a => a -> a -> a
/ (Double
b forall a. Num a => a -> a -> a
- Double
a)
          c :: Array Int Double
c = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
n) [Double]
f
	  d :: Array Int Double
d = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (Int
1,Int
nforall a. Num a => a -> a -> a
+Int
2) ((Int
nforall a. Num a => a -> a -> a
+Int
2,Double
0) forall a. a -> [a] -> [a]
: (Int
nforall a. Num a => a -> a -> a
+Int
1,Double
0) forall a. a -> [a] -> [a]
: [ (Int
j,Double
2forall a. Num a => a -> a -> a
*Double
yforall a. Num a => a -> a -> a
*Array Int Double
dforall i e. Ix i => Array i e -> i -> e
!(Int
jforall a. Num a => a -> a -> a
+Int
1) forall a. Num a => a -> a -> a
- Array Int Double
dforall i e. Ix i => Array i e -> i -> e
!(Int
jforall a. Num a => a -> a -> a
+Int
2) forall a. Num a => a -> a -> a
+ Array Int Double
cforall i e. Ix i => Array i e -> i -> e
!Int
j) | Int
j <- [Int
1..Int
n] ])
	  n :: Int
n = forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
f forall a. Num a => a -> a -> a
- Int
1