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