    {-# 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 :: forall (v :: * -> *).
    Vector v Double =>
    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 :: forall (v :: * -> *).
    Vector v Double =>
    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–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.
    --