diff --git a/bed-and-breakfast.cabal b/bed-and-breakfast.cabal --- a/bed-and-breakfast.cabal +++ b/bed-and-breakfast.cabal @@ -1,5 +1,5 @@ Name: bed-and-breakfast -Version: 0.2.3 +Version: 0.3 Synopsis: Efficient Matrix operations in 100% Haskell. Description: Efficient Matrix operations in 100% Haskell. . @@ -37,7 +37,10 @@ . [@v0.2.3@] Added 'Read' instance for @Matrix@. Improved on documentation. - + . + [@v0.3@] Added a QuickCheck test suite, fixed a bug in @det@ + (det would crash for singular matrices, where it should + return 0; it does so no). License: MIT License-File: LICENSE @@ -60,4 +63,10 @@ array >= 0.4 Hs-Source-Dirs: src +Test-Suite quickcheck + type: exitcode-stdio-1.0 + main-is: quick-check-tests.hs + build-depends: base >= 4.5 && < 5, + bed-and-breakfast == 0.3, + QuickCheck >= 2.4.2 diff --git a/quick-check-tests.hs b/quick-check-tests.hs new file mode 100644 --- /dev/null +++ b/quick-check-tests.hs @@ -0,0 +1,36 @@ +{-# LANGUAGE Haskell2010, TemplateHaskell #-} + +import Control.Monad + +import Numeric.Matrix + +import Test.QuickCheck +import Test.QuickCheck.All + +import System.Exit + + +dim :: Num a => a +dim = 4 + +instance (MatrixElement e, Arbitrary e) => Arbitrary (Matrix e) where + + arbitrary = sequence (replicate dim (vector dim)) >>= return . fromList + + +prop_det :: Matrix Rational -> Bool +prop_det m1 = let m1' = inv m1 in case m1' of (Just _) -> det m1 /= 0; _ -> det m1 == 0 + +prop_inv :: Matrix Rational -> Bool +prop_inv m1 = let m1' = inv m1 in case m1' of (Just m1') -> m1' * m1 == unit dim; _ -> True + +prop_rank :: Matrix Rational -> Bool +prop_rank m1 = let r = rank m1 in if r < dim then det m1 == 0 else r == dim + + +main = do + success <- $(quickCheckAll) + + (if success then exitSuccess else exitFailure) + + diff --git a/src/Numeric/Matrix.hs b/src/Numeric/Matrix.hs --- a/src/Numeric/Matrix.hs +++ b/src/Numeric/Matrix.hs @@ -290,22 +290,48 @@ -- | The dimensions of a given matrix. dimensions :: Matrix e -> (Int, Int) + + -- | The number of rows in the given matrix. numRows :: Matrix e -> Int + + -- | The number of columns in the given matrix. numCols :: Matrix e -> Int -- | Builds a matrix from a list of lists. -- + -- The innermost lists represent the rows. This function will create a m-n-matrix, + -- where m is the number of rows, which is the minimum length of the row lists + -- and n is the number of columns, i.e. the length of the outer list. + -- -- > fromList [[1,2,3],[2,1,3],[3,2,1]] :: Matrix Rational fromList :: [[e]] -> Matrix e toList :: Matrix e -> [[e]] -- | An identity square matrix of the given size. + -- + -- >> unit 4 + -- 1 0 0 0 + -- 0 1 0 0 + -- 0 0 1 0 + -- 0 0 0 1 unit :: Int -> Matrix e -- | A square matrix of the given size consisting of all zeros. + -- + -- >> zero 2 + -- 0 0 + -- 0 0 zero :: Int -> Matrix e + -- | A square matrix which trace is the given list, all other components + -- set to zero. + -- + -- >> diag [1,4,7,9] + -- 1 0 0 0 + -- 0 4 0 0 + -- 0 0 7 0 + -- 0 0 0 9 diag :: [e] -> Matrix e -- | Check whether the matrix is the empty matrix. @@ -318,9 +344,6 @@ times :: Matrix e -> Matrix e -> Matrix e inv :: Matrix e -> Maybe (Matrix e) --- adjugate :: Matrix e -> Matrix e --- cofactors :: Matrix e -> Matrix e ; cofactors = undefined - -- | Applies Bareiss multistep integer-preserving -- algorithm for finding the determinant of a matrix. -- Returns 0 if the matrix is not a square matrix. @@ -376,7 +399,7 @@ (x:xs) -> (length xs + 1, length x) adjugate = transpose . cofactors - transpose = fromList . L.transpose . toList + transpose m = matrix (dimensions m) (\(i,j) -> m `at` (j,i)) trace = select (uncurry (==)) inv _ = Nothing @@ -703,34 +726,40 @@ pivotR <- newSTRef 1 flip mapM_ [1..size] $ \k -> do + sign <- readSTRef signR + unless (sign == 0) $ do - prev <- readSTRef pivotR - pivot <- read a k k >>= tee (writeSTRef pivotR) + prev <- readSTRef pivotR + pivot <- read a k k >>= tee (writeSTRef pivotR) - when (pivot == 0) $ do - s <- flip mapM [(k+1)..size] $ \r -> do - a_rk <- read a r k - if a_rk == 0 then return 0 else return r - let sf = filter (>0) s + when (pivot == 0) $ do + s <- flip mapM [(k+1)..size] $ \r -> do + a_rk <- read a r k + if a_rk == 0 then return 0 else return r + let sf = filter (>0) s - when (not $ null sf) $ do - let sw = head sf + when (not $ null sf) $ do + let sw = head sf - row <- readArray a sw - readArray a k >>= writeArray a sw - writeArray a k row + row <- readArray a sw + readArray a k >>= writeArray a sw + writeArray a k row - read a k k >>= writeSTRef pivotR - readSTRef signR >>= writeSTRef signR . negate + read a k k >>= writeSTRef pivotR + readSTRef signR >>= writeSTRef signR . negate - pivot' <- readSTRef pivotR - flip mapM [(k+1)..size] $ \i -> do - a_i <- readArray a i - flip mapM [(k+1)..size] $ \j -> do - a_ij <- readArray a_i j - a_ik <- readArray a_i k - a_kj <- read a k j - writeArray a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev) + when (null sf) (writeSTRef signR 0) + + sign' <- readSTRef signR + unless (sign' == 0) $ do + pivot' <- readSTRef pivotR + flip mapM_ [(k+1)..size] $ \i -> do + a_i <- readArray a i + flip mapM [(k+1)..size] $ \j -> do + a_ij <- readArray a_i j + a_ik <- readArray a_i k + a_kj <- read a k j + writeArray a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev) liftM2 (*) (readSTRef pivotR) (readSTRef signR)