{-# LANGUAGE FlexibleContexts #-}
-- |
-- Module    : Numeric.Polynomial.Chebyshev
-- Copyright : (c) 2009, 2011 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Chebyshev polynomials.
module Numeric.Polynomial.Chebyshev (
    -- * Chebyshev polinomials
    -- $chebyshev
    chebyshev
  , chebyshevBroucke
    -- * References
    -- $references
  ) where

import qualified Data.Vector.Generic as G



-- $chebyshev
--
-- A Chebyshev polynomial of the first kind is defined by the
-- following recurrence:
--
-- \[\begin{aligned}
-- T_0(x)     &= 1 \\
-- T_1(x)     &= x \\
-- T_{n+1}(x) &= 2xT_n(x) - T_{n-1}(x) \\
-- \end{aligned}
-- \]

data C = C {-# UNPACK #-} !Double {-# UNPACK #-} !Double

-- | Evaluate a Chebyshev polynomial of the first kind. Uses
-- Clenshaw's algorithm.
chebyshev :: (G.Vector v Double) =>
             Double      -- ^ Parameter of each function.
          -> v Double    -- ^ Coefficients of each polynomial term, in increasing order.
          -> Double
chebyshev :: Double -> v Double -> Double
chebyshev Double
x v Double
a = C -> Double
fini (C -> Double) -> (v Double -> C) -> v Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> C -> C) -> C -> v Double -> C
forall (v :: * -> *) a b.
Vector v a =>
(a -> b -> b) -> b -> v a -> b
G.foldr' Double -> C -> C
step (Double -> Double -> C
C Double
0 Double
0) (v Double -> C) -> (v Double -> v Double) -> v Double -> C
forall b c a. (b -> c) -> (a -> b) -> a -> c
. v Double -> v Double
forall (v :: * -> *) a. Vector v a => v a -> v a
G.tail (v Double -> Double) -> v Double -> Double
forall a b. (a -> b) -> a -> b
$ v Double
a
    where step :: Double -> C -> C
step Double
k (C Double
b0 Double
b1) = Double -> Double -> C
C (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b1) Double
b0
          fini :: C -> Double
fini   (C Double
b0 Double
b1) = v Double -> Double
forall (v :: * -> *) a. Vector v a => v a -> a
G.head v Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b1
          x2 :: Double
x2               = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
2
{-# INLINE chebyshev #-}

data B = B {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double

-- | Evaluate a Chebyshev polynomial of the first kind. Uses Broucke's
-- ECHEB algorithm, and his convention for coefficient handling. It
-- treat 0th coefficient different so
--
-- > chebyshev x [a0,a1,a2...] == chebyshevBroucke [2*a0,a1,a2...]
chebyshevBroucke :: (G.Vector v Double) =>
             Double      -- ^ Parameter of each function.
          -> v Double    -- ^ Coefficients of each polynomial term, in increasing order.
          -> Double
chebyshevBroucke :: Double -> v Double -> Double
chebyshevBroucke Double
x = B -> Double
fini (B -> Double) -> (v Double -> B) -> v Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> B -> B) -> B -> v Double -> B
forall (v :: * -> *) a b.
Vector v a =>
(a -> b -> b) -> b -> v a -> b
G.foldr' Double -> B -> B
step (Double -> Double -> Double -> B
B Double
0 Double
0 Double
0)
    where step :: Double -> B -> B
step Double
k (B Double
b0 Double
b1 Double
_) = Double -> Double -> Double -> B
B (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b1) Double
b0 Double
b1
          fini :: B -> Double
fini   (B Double
b0 Double
_ Double
b2) = (Double
b0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
0.5
          x2 :: Double
x2                 = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
2
{-# INLINE chebyshevBroucke #-}



-- $references
--
-- * Broucke, R. (1973) Algorithm 446: Ten subroutines for the
--   manipulation of Chebyshev series. /Communications of the ACM/
--   16(4):254&#8211;256.  <http://doi.acm.org/10.1145/362003.362037>
--
-- * Clenshaw, C.W. (1962) Chebyshev series for mathematical
--   functions. /National Physical Laboratory Mathematical Tables 5/,
--   Her Majesty's Stationery Office, London.
--