{-# LANGUAGE DataKinds #-} {-# LANGUAGE TemplateHaskell #-} module Numeric.QuaterFloatTest (runTests) where import Numeric.DataFrame.Arbitraries () import Numeric.Quaternion import Numeric.Scalar import Numeric.Vector import Test.QuickCheck type T = Float eps :: T eps = 1e-3 -- | Some non-linear function are very unstable; -- it would be a downting task to determine the uncertainty precisely. -- Instead, I just make sure the numbers are not so big to find at least -- the most obvious bugs. smaller :: Quater T -> Quater T smaller q@(Quater x y z w) = if square q >= 1 then Quater (f x) (f y) (f z) (f w) else q where f :: T -> T f a | abs a > 2 = signum a * log (abs a) | otherwise = a (=~=) :: Quater T -> Quater T -> Property (=~=) = approxEq eps infix 4 =~= -- Even lower precision - for naturally fragile operations approxEq :: T -> Quater T -> Quater T -> Property approxEq e a b = property \$ err <= c * scalar e where v1 = toVec4 a v2 = toVec4 b c = 1 `max` normL1 v1 `max` normL1 v2 err = normL1 \$ v1 - v2 prop_Eq :: Quater T -> Bool prop_Eq q = and [ q == q , Quater (takei q) (takej q) (takek q) (taker q) == q , fromVecNum (imVec q) (taker q) == q , im q + re q == q , fromVec4 (toVec4 q) == q ] prop_DoubleConjugate :: Quater T -> Property prop_DoubleConjugate q = property \$ conjugate (conjugate q) == q prop_Square :: Quater T -> Property prop_Square q = q * conjugate q =~= realToFrac (square q) prop_RotScale :: Quater T -> Vector T 3 -> Property prop_RotScale q v = fromVecNum (rotScale q v) 0 =~= q * fromVecNum v 0 * conjugate q prop_GetRotScale :: Vector T 3 -> Vector T 3 -> Property prop_GetRotScale a b = a /= 0 ==> approxEq eps' (fromVecNum (rotScale q a) 0) (fromVecNum b 0) where q = getRotScale a b -- when a and b are almost opposite, precision of getRotScale suffers a lot -- compensate it by increasing eps: eps' = if cross a b == 0 then eps else eps / (max eps . min 1 . unScalar \$ 1 + normalized a `dot` normalized b) prop_InverseRotScale :: Quater T -> Vector T 3 -> Property prop_InverseRotScale q v = square q > eps ==> fromVecNum (rotScale (1/q) (rotScale q v)) 0 =~= fromVecNum v 0 prop_NegateToMatrix33 :: Quater T -> Bool prop_NegateToMatrix33 q = toMatrix33 q == toMatrix33 (negate q) prop_NegateToMatrix44 :: Quater T -> Bool prop_NegateToMatrix44 q = toMatrix44 q == toMatrix44 (negate q) prop_FromToMatrix33 :: Quater T -> Property prop_FromToMatrix33 q = q /= 0 ==> fromMatrix33 (toMatrix33 q) =~= q .||. fromMatrix33 (toMatrix33 q) =~= negate q prop_FromToMatrix44 :: Quater T -> Property prop_FromToMatrix44 q = q /= 0 ==> fromMatrix44 (toMatrix44 q) =~= q .||. fromMatrix44 (toMatrix44 q) =~= negate q prop_RotationArg :: Quater T -> Property prop_RotationArg q | q == 0 = axisRotation (imVec q) (qArg q) =~= 1 | otherwise = axisRotation (imVec q) (qArg q) =~= signum q prop_UnitQ :: Quater T -> Property prop_UnitQ q = square q > eps ==> realToFrac (square (q / q)) =~= 1 prop_ExpLog :: Quater T -> Property prop_ExpLog q | square q < eps = log (exp q) =~= q | otherwise = exp (log q) =~= q prop_SinAsin :: Quater T -> Property prop_SinAsin q = sin (asin q') =~= q' where q' = signum q prop_CosAcos :: Quater T -> Property prop_CosAcos q = cos (acos q') =~= q' where q' = signum q prop_TanAtan :: Quater T -> Property prop_TanAtan q = tan (atan q') =~= q' where q' = sqrt (smaller q) -- because this inverse fun is fragile prop_SinhAsinh :: Quater T -> Property prop_SinhAsinh q = sinh (asinh q') =~= q' where q' = sqrt (smaller q) -- because this inverse fun is fragile prop_CoshAcosh :: Quater T -> Property prop_CoshAcosh q = cosh (acosh q) =~= q prop_TanhAtanh :: Quater T -> Property prop_TanhAtanh q = tanh (atanh q') =~= q' where q' = sqrt (smaller q) -- because this inverse fun is fragile prop_SqrtSqr :: Quater T -> Property prop_SqrtSqr q = sqrt q * sqrt q =~= q prop_SinCos :: Quater T -> Property prop_SinCos q = sin q' * sin q' + cos q' * cos q' =~= 1 where q' = sqrt (smaller q) -- because it involves exponents prop_SinhCosh :: Quater T -> Property prop_SinhCosh q = cosh q' * cosh q' - sinh q' * sinh q' =~= 1 where q' = sqrt (smaller q) -- because it involves exponents prop_ReadShow :: Quater T -> Bool prop_ReadShow q = q == read (show q) return [] runTests :: Int -> IO Bool runTests n = \$forAllProperties \$ quickCheckWithResult stdArgs { maxSuccess = n }