Add Chebyshev approximations

bench/Main.hs

@@ -107,18 +107,6 @@ vdot' :: forall f k p a.
          => a -> f a -> f a -> a
 vdot' r xs ys = r + sum (liftA2' (unapply (uncurry (+) . unprod)) (xs, ys))
                 \\ (proveCartesian @k :: (Obj k a, Obj k a) :- Obj k (p a a))
--- vdot' r xs ys = r + sum (liftA2' add (xs, ys))
---                    where {-# INLINE add #-}
---                          add :: p a a `k` a
---                          add = unapply (uncurry (+) . unprod)
---                                \\ (proveCartesian @k ::
---                                        (Obj k a, Obj k a) :- Obj k (p a a))
---                          -- add :: (a, a) -> a
---                          -- add (x, y) = x + y
---                          -- add :: (a *#* a) -#> a
---                          -- add = UFun $ \(UProd (x, y)) -> x + y
---                          --       \\ (proveCartesian @k ::
---                          --               (Obj k a, Obj k a) :- Obj k (p a a))
 
 -- bench_vdot1 :: Double -> Double
 -- bench_vdot1 r = let xs :: V.Vector Double

lib/Chebyshev.hs

@@ -0,0 +1,75 @@
+-- | Adapted from [math-functions-0.2.1.0] Numeric.Polynomial.Chebyshev
+module Chebyshev ( chebyshev
+                 , chebyshevApprox
+                 , chebyshevApprox'
+                 ) where
+
+import Prelude hiding ( id, (.), curry, uncurry
+                      , Functor(..)
+                      , Foldable(..)
+                      )
+
+import Data.Maybe
+
+import Category
+import Foldable
+import Functor
+import Unfoldable
+
+
+
+-- $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 a = C !a !a
+
+
+
+-- | Evaluate a Chebyshev polynomial of the first kind, using the
+-- Clenshaw algorithm; see
+-- <https://en.wikipedia.org/wiki/Clenshaw_algorithm>.
+chebyshev :: (Foldable v, Obj (Dom v) a, Num a)
+             => v a   -- ^ Polynomial coefficients in increasing order
+             -> a     -- ^ Position
+             -> a
+chebyshev cs' x =
+    if null cs then 0 else fini . foldr step (C 0 0) . tail $ cs
+    where cs = toList cs'
+          step k (C b0 b1) = C (k + 2 * x * b0 - b1) b0
+          fini   (C b0 b1) = head cs + x * b0 - b1
+{-# INLINE chebyshev #-}
+
+
+
+-- | Approximate a given function via Chebyshev polynomials.
+-- See <http://mathworld.wolfram.com/ChebyshevApproximationFormula.html>.
+chebyshevApprox :: ( Unfoldable v, Obj (Dom v) b
+                   , RealFloat a, Fractional b
+                   ) => Int -> (a -> b) -> v b
+chebyshevApprox n = chebyshevApprox' (2 * n) n
+
+chebyshevApprox' :: forall f a b.
+                    ( Unfoldable f, Obj (Dom f) b
+                    , RealFloat a, Fractional b
+                    ) => Int -> Int -> (a -> b) -> f b
+chebyshevApprox' np nc f =
+    (fromJust . fromList) [coeff i | i <- [0..nc-1]]
+    where coeff j = rf ((if j == 0 then 1 else 2) / fi np) *
+                    sum [f (x k) * rf (cheb j k) | k <- [0..np-1]]
+          x :: Int -> a
+          x k = cos (pi * (fi k + 0.5) / fi np)
+          cheb :: Int -> Int -> a
+          cheb j k = cos (pi * fi j * (fi k + 0.5) / fi np)
+          fi = fromIntegral
+          rf = realToFrac

stack.yaml

@@ -15,7 +15,7 @@
 # resolver:
 #  name: custom-snapshot
 #  location: "./custom-snapshot.yaml"
-resolver: lts-11.6
+resolver: lts-11.7
 
 # User packages to be built.
 # Various formats can be used as shown in the example below.

test/ChebyshevSpec.hs

@@ -0,0 +1,50 @@
+module ChebyshevSpec where
+
+import Test.QuickCheck
+import Test.QuickCheck.Instances()
+
+import Chebyshev
+
+
+
+-- | Approximate floating-point equality
+approx :: (RealFrac a, Show a) => a -> a -> a -> Property
+approx eps x y = counterexample (show x ++ " /= " ++ show y) (abs (x - y) < eps)
+
+-- | n-th Chebyshev polynomial
+cheb :: Num a => Int -> a -> a
+cheb n = chebyshev (replicate n 0 ++ [1])
+
+-- | Recurrence relation for Chebyshev polynomials
+recur :: Num a => Int -> a -> a
+recur n x = if | n == 0 -> 1
+               | n == 1 -> x
+               | otherwise -> 2 * x * cheb (n - 1) x - cheb (n - 2) x
+
+
+
+prop_Chebyshev_recurrence ::
+    NonNegative (Small Int) -> Integer -> Property
+prop_Chebyshev_recurrence (NonNegative (Small n)) x =
+    cheb n x === recur n x
+
+
+
+-- prop_Chebyshev_approx :: NonNegative (Small Int) -> Property
+-- prop_Chebyshev_approx (NonNegative (Small n)) =
+--     conjoin (zipWith (approx 1.0e-13) coeffs' coeffs)
+--     where phi :: Double -> Double
+--           phi = cheb n
+--           coeffs = replicate n 0 ++ [1] ++ replicate n 0
+--           coeffs' = chebyshevApprox (2 * n + 1) phi
+
+
+
+prop_Chebyshev_approx :: [Double] -> Property
+prop_Chebyshev_approx cs' =
+    conjoin (zipWith (approx (1.0e-13 * maxc)) cs fcs)
+    where cs = cs' ++ [1]
+          maxc = maximum (fmap abs cs)
+          n = length cs
+          f = chebyshev cs
+          fcs = chebyshevApprox (2 * n) f