Commits

Aleksey Khudyakov committed 8611699

Rewrite tests for Chebyshev polynomials

Reference implementation now uses Rational for intermediate
calculations. It's more robust since they do not suffer from
rounding errors.

At the moment tests fails occasionly. One typically need to run
test 1e4 times in order t get a failure. It's not a bug but an
effect of rounding errors.

Distribution of errors would be much more enlightening but
QC do not allow to build it.

  • Participants
  • Parent commits b395cc6

Comments (0)

Files changed (1)

tests/Tests/Chebyshev.hs

 import Data.Vector.Unboxed                  (fromList)
 import Test.Framework
 import Test.Framework.Providers.QuickCheck2
-import Test.QuickCheck                      (Arbitrary(..))
+import Test.QuickCheck                      (Arbitrary(..),printTestCase,Property)
 
 import Tests.Helpers
 import Numeric.Polynomial.Chebyshev
 tests :: Test
 tests = testGroup "Chebyshev polynomials"
   [ testProperty "Chebyshev 0" $ \a0 (Ch x) ->
-      (ch0 x * a0) ≈ (chebyshev x $ fromList [a0])
+      testCheb [a0] x
   , testProperty "Chebyshev 1" $ \a0 a1 (Ch x) ->
-      (a0*ch0 x + a1*ch1 x) ≈  (chebyshev x $ fromList [a0,a1])
+      testCheb [a0,a1] x
   , testProperty "Chebyshev 2" $ \a0 a1 a2 (Ch x) ->
-       (a0*ch0 x + a1*ch1 x + a2*ch2 x) ≈ (chebyshev x $ fromList [a0,a1,a2])
+      testCheb [a0,a1,a2] x
   , testProperty "Chebyshev 3" $ \a0 a1 a2 a3 (Ch x) ->
-       (a0*ch0 x + a1*ch1 x + a2*ch2 x + a3*ch3 x) ≈ (chebyshev x $ fromList [a0,a1,a2,a3])
+      testCheb [a0,a1,a2,a3] x
   , testProperty "Chebyshev 4" $ \a0 a1 a2 a3 a4 (Ch x) ->
-       (a0*ch0 x + a1*ch1 x + a2*ch2 x + a3*ch3 x + a4*ch4 x) ≈ (chebyshev x $ fromList [a0,a1,a2,a3,a4])
+       testCheb [a0,a1,a2,a3,a4] x
   ]
-  where (≈) = eq 1e-12
+  where
 
+testCheb :: [Double] -> Double -> Property
+testCheb as x
+  = printTestCase (">>> Exact   = " ++ show exact)
+  $ printTestCase (">>> Numeric = " ++ show num  )
+  $ printTestCase (">>> rel.err.= " ++ show err  )
+  $ eq 1e-12 num exact
+  where
+    exact = evalCheb as x
+    num   = chebyshev x (fromList as)
+    err   = abs (num - exact) / abs exact
+
+evalCheb :: [Double] -> Double -> Double
+evalCheb as x
+  = realToFrac
+  $ sum
+  $ zipWith (*) (map realToFrac as)
+  $ map ($ realToFrac x) cheb
 
 -- Chebyshev polynomials of low order
-ch0,ch1,ch2,ch3,ch4 :: Double -> Double
-ch0 _ = 1
-ch1 x = x
-ch2 x = 2*x^2 - 1
-ch3 x = 4*x^3 - 3*x
-ch4 x = 8*x^4 - 8*x^2 + 1
-
+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