1{-# OPTIONS_GHC -fno-warn-type-defaults #-}
2
3module Tests.Chebyshev ( tests ) where
4
5import Data.Vector.Unboxed                  (fromList)
6import Test.Tasty
7import Test.Tasty.QuickCheck                (testProperty)
8import Test.QuickCheck                      (Arbitrary(..),counterexample,Property)
9
10import Tests.Helpers
11import Numeric.Polynomial.Chebyshev
12import Numeric.MathFunctions.Comparison
13
14
15tests :: TestTree
16tests = testGroup "Chebyshev polynomials"
17  [ testProperty "Chebyshev 0" $ \a0 (Ch x) ->
18      testCheb [a0] x
19  , testProperty "Chebyshev 1" $ \a0 a1 (Ch x) ->
20      testCheb [a0,a1] x
21  , testProperty "Chebyshev 2" $ \a0 a1 a2 (Ch x) ->
22      testCheb [a0,a1,a2] x
23  , testProperty "Chebyshev 3" $ \a0 a1 a2 a3 (Ch x) ->
24      testCheb [a0,a1,a2,a3] x
25  , testProperty "Chebyshev 4" $ \a0 a1 a2 a3 a4 (Ch x) ->
26      testCheb [a0,a1,a2,a3,a4] x
27  , testProperty "Broucke" $ testBroucke
28  ]
29  where
30
31testBroucke :: Ch -> [Double] -> Property
32testBroucke _      []     = counterexample "" True
33testBroucke (Ch x) (c:cs)
34  = counterexample (">>> Chebyshev  = " ++ show c1)
35  $ counterexample (">>> Brouke     = " ++ show cb)
36  $ counterexample (">>> rel.err.   = " ++ show (relativeError c1 cb))
37  $ counterexample (">>> diff. ulps = " ++ show (ulpDistance   c1 cb))
38  $ within 64 c1 cb
39  where
40    c1 = chebyshev        x (fromList $ c : cs)
41    cb = chebyshevBroucke x (fromList $ c*2 : cs)
42
43
44testCheb :: [Double] -> Double -> Property
45testCheb as x
46  = counterexample (">>> Exact      = " ++ show exact)
47  $ counterexample (">>> Numeric    = " ++ show num  )
48  $ counterexample (">>> rel.err.   = " ++ show err  )
49  $ counterexample (">>> diff. ulps = " ++ show (ulpDistance num exact))
50  $ eq 1e-12 num exact
51  where
52    exact = evalCheb as x
53    num   = chebyshev x (fromList as)
54    err   = relativeError num exact
55
56evalCheb :: [Double] -> Double -> Double
57evalCheb as x
58  = realToFrac
59  $ sum
60  $ zipWith (*) (map realToFrac as)
61  $ map ($ realToFrac x) cheb
62
63-- Chebyshev polynomials of low order
64cheb :: [Rational -> Rational]
65cheb =
66  [ \_ -> 1
67  , \x -> x
68  , \x -> 2*x^2 - 1
69  , \x -> 4*x^3 - 3*x
70  , \x -> 8*x^4 - 8*x^2 + 1
71  ]
72
73-- Double in the [-1 .. 1] range
74newtype Ch = Ch Double
75             deriving Show
76instance Arbitrary Ch  where
77  arbitrary = do x <- arbitrary
78                 return $ Ch $ 2 * (abs . snd . properFraction) x - 1
79