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