diff --git a/Eigensystem.hs b/Eigensystem.hs
--- a/Eigensystem.hs
+++ b/Eigensystem.hs
@@ -54,6 +54,8 @@
import QuantumVector
import LinearAlgorithms (triangular, tridiagonal, triangular2)
import Data.List (findIndex)
+import Prelude2010
+import Prelude ()
----------------------------------------------------------------------------
-- Category: Eigensystem for QuantumVector
diff --git a/EigensystemNum.hs b/EigensystemNum.hs
--- a/EigensystemNum.hs
+++ b/EigensystemNum.hs
@@ -10,19 +10,20 @@
matSqr x = mult x x
powerIter :: (Fractional a, Ord a) => [[a]] -> [([[a]],[[a]])]
-powerIter x = tail (iterate
+powerIter x =
+ tail $ iterate
(\(_,z)->let s=normalize (matSqr z) in (s,(mult x s)))
([],x)
- )
normalize :: (Fractional a, Ord a) => [[a]] -> [[a]]
-normalize x = map (map (/(matnorm1 x))) x
+normalize x = map (map (/ matnorm1 x)) x
getGrowth :: (Fractional a, Ord a) => ([[a]],[[a]]) -> a
-getGrowth (x,y) = uncurry (/) (maximumBy
+getGrowth (x,y) =
+ uncurry (/) $
+ maximumBy
(\(_,xc) (_,xa) -> compare (abs xc) (abs xa))
(concat (zipWith zip y x))
- )
specRadApprox :: (Fractional a, Ord a) => [[a]] -> [a]
specRadApprox = map getGrowth . powerIter
@@ -31,7 +32,8 @@
eigenValuesApprox = map diagonals . iterate similar_to
limit :: (Num a, Ord a) => a -> [a] -> a
-limit tol (x0:x1:xs) = if abs (x1-x0) < tol * abs x0
- then x0
- else limit tol (x1:xs)
+limit tol (x0:x1:xs) =
+ if abs (x1-x0) < tol * abs x0
+ then x0
+ else limit tol (x1:xs)
limit _ _ = error "Only infinite sequences are allowed"
diff --git a/Fraction.hs b/Fraction.hs
--- a/Fraction.hs
+++ b/Fraction.hs
@@ -266,6 +266,8 @@
--
-- Require: No Taylor coefficient is zero
--
+ -- It is an application of Euler's continued fraction formula.
+ --
zero:one:[higher m | m <- [2..]]
where
zero = (s!!0, s!!1 * x)
diff --git a/Makefile b/Makefile
--- a/Makefile
+++ b/Makefile
@@ -3,3 +3,14 @@
%.html: %.lhs
ln -s $< $@
+
+
+
+run-test: update-test
+ runhaskell Setup.lhs configure --user --enable-tests
+ runhaskell Setup.lhs build
+ runhaskell Setup.lhs haddock
+ ./dist/build/numeric-quest-test/numeric-quest-test
+
+update-test:
+ doctest-extract-0.1 -o test/ --executable-main=Test/Main.hs RowEchelon
diff --git a/Orthogonals.lhs b/Orthogonals.lhs
--- a/Orthogonals.lhs
+++ b/Orthogonals.lhs
@@ -304,7 +304,7 @@
by some algorithms. So far we have been able to avoid this.
-> class Scalar a where
+> class Eq a => Scalar a where
> coupled :: a->a
> norm :: [a] -> a
> almostZero :: a -> Bool
@@ -317,10 +317,10 @@
> scaled = scaled'
> instance Scalar Float where
-> coupled x = x
-> norm u = sqrt (bra_ket u u)
-> almostZero x = (abs x) < 1.0e-8
-> scaled = scaled'
+> coupled x = x
+> norm u = sqrt (bra_ket u u)
+> almostZero x = (abs x) < 1.0e-8
+> scaled = scaled'
> instance (Integral a) => Scalar (Ratio a) where
> coupled x = x
@@ -380,15 +380,11 @@
> normalized :: (Scalar a, Fractional a) => [a] -> [a]
> normalized u =
-> [uk/n | uk <- u]
-> where
-> n = norm u
+> map (/norm u) u
> scaled' :: (Fractional t, Ord t) => [t] -> [t]
> scaled' u =
-> [uk/um | uk <- u]
-> where
-> um = maximum [abs uk| uk <- u]
+> map (/norminf u) u
@@ -783,7 +779,7 @@
-> one_ket_solution :: (Fractional a, Scalar a) => [[a]] -> [a] -> [a]
+> one_ket_solution :: (Scalar a, Fractional a) => [[a]] -> [a] -> [a]
> one_ket_solution a b =
> --
> -- List representing vector |x>, which is
@@ -951,7 +947,7 @@
-> factors_QR :: (Fractional a, Scalar a) => [[a]] -> ([[a]],[[a]])
+> factors_QR :: (Scalar a, Fractional a) => [[a]] -> ([[a]],[[a]])
> factors_QR a =
> --
> -- A pair of matrices (Q, R), such that
@@ -979,7 +975,7 @@
-> determinant :: (Fractional a, Scalar a) => [[a]] -> a
+> determinant :: (Scalar a, Fractional a) => [[a]] -> a
> determinant a =
> let (q,r) = factors_QR a
> -- matrix Q is not normed so we have to respect the norms of its rows
@@ -1188,7 +1184,7 @@
-> similar_to :: (Fractional a, Scalar a) => [[a]] -> [[a]]
+> similar_to :: (Scalar a, Fractional a) => [[a]] -> [[a]]
> similar_to a =
> --
> -- List of columns of matrix A1 similar to A
@@ -1201,7 +1197,7 @@
> where
> (q,r) = factors_QR a
-> iterated_eigenvalues :: (Scalar a1, Fractional a1, Num a) => [[a1]] -> a -> [[a1]]
+> iterated_eigenvalues :: (Scalar a1, Fractional a1, Eq a, Num a) => [[a1]] -> a -> [[a1]]
> iterated_eigenvalues a n
> --
> -- List of vectors representing
@@ -1215,7 +1211,7 @@
> | otherwise = (diagonals a)
> : iterated_eigenvalues (similar_to a) (n-1)
-> eigenvalues :: (Scalar a1, Fractional a1, Num a) => [[a1]] -> a -> [a1]
+> eigenvalues :: (Scalar a1, Fractional a1, Eq a, Num a) => [[a1]] -> a -> [a1]
> eigenvalues a n
> --
> -- Eigenvalues of matrix A
@@ -1789,24 +1785,14 @@
> --
> -- Unit square matrix of with dimensions m x m
> --
-> [g 0 k | k <- [0..(m-1)]]
-> where
-> g i k
-> | i == m = []
-> | i == k = 1:(g (i+1) k)
-> | otherwise = 0:(g (i+1) k)
->
+> [ [ if j==k then 1 else 0 | j <- [0 .. m-1] ] | k <- [0 .. m-1]]
> unit_vector :: Num a => Int -> Int -> [a]
> unit_vector i m =
> --
> -- Unit vector of length m
> -- with 1 at position i, zero otherwise
-> [g i k| k <- [0..(m-1)]]
-> where
-> g j k
-> | j == k = 1
-> | otherwise = 0
+> map (\k -> if k==i then 1 else 0) [0 .. m-1]
> diagonals :: [[a]] -> [a]
> diagonals a =
diff --git a/QuantumVector.lhs b/QuantumVector.lhs
--- a/QuantumVector.lhs
+++ b/QuantumVector.lhs
@@ -157,6 +157,8 @@
> module QuantumVector where
> import Data.Complex -- our Scalar is Complex Double
> import Data.List (nub)
+> import Prelude2010
+> import Prelude ()
> infixl 7 *> -- tensor product of two kets
> infixl 7 <* -- tensor product of two bras
@@ -1129,7 +1131,7 @@
> showsPrec n (j :<+ k) = showsPrec n j . showString " + " . showsPrec n k
-> showsScalar :: (RealFloat t) => Int -> Complex t -> String -> String
+> showsScalar :: (Show t, RealFloat t) => Int -> Complex t -> String -> String
> showsScalar n x@(a :+ b)
> | b == 0 = showsPrec n a . showString " "
> | otherwise = showString "(" .showsPrec n x . showString ") "
diff --git a/Roots.hs b/Roots.hs
--- a/Roots.hs
+++ b/Roots.hs
@@ -77,20 +77,12 @@
u y a b = a + b*y
polynomial_derivative :: Num a => [a] -> [a]
-polynomial_derivative as
+polynomial_derivative as =
--
-- List of coefficients for derivative of polynomial
-- a0 + a1 x + a2 x^2 ...
--
- | as == [] = []
- | otherwise = deriv 1 (drop 1 as) []
- where
- deriv n bs cs
- | bs == [] = reverse2 cs
- | otherwise = deriv (n+1) (tail bs) ((n*(head bs)):cs)
- reverse2 cs
- | cs == [] = []
- | otherwise = reverse cs
+ zipWith (*) (iterate (1+) 1) (drop 1 as)
-----------------------------------------------------------------------------
--
diff --git a/RowEchelon.hs b/RowEchelon.hs
new file mode 100644
--- /dev/null
+++ b/RowEchelon.hs
@@ -0,0 +1,412 @@
+{- |
+
+
+Duplicate of @htam@ package.
+-}
+module RowEchelon (
+ Matrix(Matrix), matrixRows, matrixWidth, matrixHeight,
+ shortBesidesTall, (|||),
+ matrixValid, matrixFromRows,
+ identity, matrixProduct,
+ nullspace, reducedRowEchelon, rowReduction, layoutEchelonBlocks,
+ scatter,
+ Zipper0(Zipper0), Zipper1(Zipper1), altOuter,
+ ) where
+
+import Orthogonals (matrix_matrix)
+
+import qualified Data.Foldable as Fold
+import qualified Data.NonEmpty as NonEmpty
+import qualified Data.List.HT as ListHT
+import qualified Data.List as List
+import Data.NonEmpty ((!:))
+import Data.Maybe (fromMaybe)
+
+
+{- $setup
+>>> import RowEchelon as Matrix
+>>> import qualified Test.QuickCheck as QC
+>>> import Test.QuickCheck ((===))
+>>> import Control.Monad (replicateM)
+>>> import qualified Data.Foldable as Fold
+>>> import qualified Data.NonEmpty.Class as NonEmptyC
+>>> import qualified Data.NonEmpty as NonEmpty
+>>> import qualified Data.List.Match as Match
+>>> import qualified Data.List.HT as ListHT
+>>> import Data.NonEmpty ((!:))
+>>> import Data.Ratio ((%))
+>>>
+>>> genElementUniform, genElementNearZero :: QC.Gen Integer
+>>> genElementUniform = QC.choose (-10,10)
+>>> genElementNearZero = fmap (flip rem 11) $ QC.arbitrary
+>>>
+>>> genMatrixForSize :: Int -> Int -> QC.Gen (Matrix Integer)
+>>> genMatrixForSize m n = do
+>>> fmap (Matrix n) $ replicateM m $ replicateM n genElementNearZero
+>>>
+>>> genMatrix :: QC.Gen (Matrix Integer)
+>>> genMatrix = do
+>>> m <- QC.choose (0,10)
+>>> n <- QC.choose (0,10)
+>>> genMatrixForSize m n
+>>>
+>>> shrinkMatrix :: (Eq a) => Matrix a -> [Matrix a]
+>>> shrinkMatrix matrix@(Matrix width rows) =
+>>> filter (matrix/=) $
+>>> map (Matrix width . snd) (ListHT.removeEach rows)
+>>> ++
+>>> [Matrix (width-1) $ map (drop 1) rows]
+>>>
+>>> forMatrix :: (QC.Testable test) => (Matrix Rational -> test) -> QC.Property
+>>> forMatrix prop =
+>>> QC.forAllShrink genMatrix shrinkMatrix (prop . rationalMatrix)
+>>>
+>>> rationalMatrix :: Matrix Integer -> Matrix Rational
+>>> rationalMatrix = fmap (%1)
+-}
+
+
+data Zipper0 a = Zipper0 [a] [a]
+ deriving (Eq, Show)
+
+data Zipper1 a = Zipper1 [a] a [a]
+ deriving (Eq, Show)
+
+instance Functor Zipper0 where
+ fmap f (Zipper0 xs ys) = Zipper0 (fmap f xs) (fmap f ys)
+
+instance Functor Zipper1 where
+ fmap f (Zipper1 xs y zs) = Zipper1 (fmap f xs) (f y) (fmap f zs)
+
+zipper0FromList :: [a] -> Zipper0 a
+zipper0FromList = Zipper0 []
+
+
+type List0 = []
+type List1 = NonEmpty.T List0
+
+-- data Matrix a = Matrix {matrixWidth :: Int, matrixRows :: [[a]]}
+data Matrix a = Matrix Int [[a]]
+ deriving (Eq, Show)
+
+instance Functor Matrix where
+ fmap f (Matrix width rows) = Matrix width $ map (map f) rows
+
+instance Foldable Matrix where
+ foldMap f (Matrix _width rows) = foldMap (foldMap f) rows
+
+matrixValid :: Matrix a -> Bool
+matrixValid (Matrix width rows) = all ((width==) . length) rows
+
+matrixWidth :: Matrix a -> Int
+matrixWidth (Matrix width _) = width
+
+matrixHeight :: Matrix a -> Int
+matrixHeight = length . matrixRows
+
+matrixRows :: Matrix a -> [[a]]
+matrixRows (Matrix _ rows) = rows
+
+matrixFromRows :: [[a]] -> Matrix a
+matrixFromRows rows =
+ Matrix (ListHT.switchL 0 (const . length) rows) rows
+
+infixr 3 |||
+
+-- | requires that both matrices have the same height
+(|||) :: Matrix a -> Matrix a -> Matrix a
+Matrix widthA a ||| Matrix widthB b =
+ Matrix (widthA+widthB) (zipWith (++) a b)
+
+{- |
+The expression @shortBesidesTall a b@ means:
+
+> /A B_upper\
+> | |
+> \0 B_lower/
+
+Matrix @a@ must be at most as tall as @b@.
+-}
+shortBesidesTall :: (Num a) => Matrix a -> Matrix a -> Matrix a
+shortBesidesTall a b =
+ Matrix (matrixWidth a + matrixWidth b) $
+ zipWith (++)
+ (matrixRows a ++ repeat (replicate (matrixWidth a) 0))
+ (matrixRows b)
+
+matrixProduct :: (Num a) => Matrix a -> Matrix a -> Matrix a
+matrixProduct (Matrix _ x) (Matrix width y) =
+ Matrix width $ matrix_matrix x $ List.transpose y
+
+{- |
+>>> nullspace (Matrix 2 []) :: Matrix Rational
+Matrix 2 [[1 % 1,0 % 1],[0 % 1,1 % 1]]
+
+
+prop> forMatrix $ matrixValid . nullspace
+
+
+prop> :{
+ forMatrix $ \matrix ->
+ matrixWidth matrix == matrixHeight (nullspace matrix)
+:}
+
+prop> :{
+ forMatrix $ \matrix ->
+ matrixWidth (nullspace matrix) <= matrixWidth matrix
+:}
+
+max 0 (width matrix - height matrix) <= width nullspace
+
+prop> :{
+ forMatrix $ \matrix ->
+ matrixWidth matrix <= matrixWidth (nullspace matrix) + matrixHeight matrix
+:}
+
+prop> :{
+ forMatrix $ \matrix ->
+ Fold.all (0==) $ matrixProduct matrix (nullspace matrix)
+:}
+-}
+nullspace :: (RealFrac a) => Matrix a -> Matrix a
+nullspace matrix =
+ let echelon = reducedRowEchelon matrix in
+ let flat = layoutEchelonBlocks echelon in
+ let nullDim = matrixWidth flat in
+ Matrix nullDim $
+ scatter
+ (fmap matrixWidth $ altOuter echelon)
+ (matrixRows $ fmap negate flat) $
+ matrixRows $ identity nullDim
+
+identity :: (Num a) => Int -> Matrix a
+identity n = Matrix n $ take n $ map (take n) $ iterate (0:) $ 1 : repeat 0
+
+
+-- cf. event-list
+data AlternatingList a b = AlternatingList a [(b,a)]
+ deriving (Eq)
+
+instance (Show a, Show b) => Show (AlternatingList a b) where
+ showsPrec p xs =
+ showParen (p>=5) $
+ flip
+ (altFoldr
+ (\a -> showsPrec 5 a . showString " ./ ")
+ (\b -> showsPrec 5 b . showString " /. "))
+ xs
+ .
+ showString "[]"
+
+altSingleton :: a -> AlternatingList a b
+altSingleton a = AlternatingList a []
+
+infixr 5 /. , ./
+
+(./) = altConsOuter
+(/.) = altConsInner
+
+altConsOuter, (./) :: a -> List0 (b,a) -> AlternatingList a b
+altConsOuter = AlternatingList
+
+altConsInner, (/.) :: b -> AlternatingList a b -> List0 (b,a)
+altConsInner b ~(AlternatingList a bas) = (b,a) : bas
+
+altOuter :: AlternatingList a b -> List1 a
+altOuter (AlternatingList a bas) = a !: map snd bas
+
+altFoldr :: (a -> c -> d) -> (b -> d -> c) -> c -> AlternatingList a b -> d
+altFoldr f g c (AlternatingList a0 bas) =
+ f a0 $ foldr (\(b,a) -> g b . f a) c bas
+
+
+{- |
+>>> scatter (3 !: 1 : 4 : 1 : []) "012" ['a'..'k']
+"abc0d1efgh2ijk"
+-}
+scatter :: List1 Int -> [a] -> [a] -> [a]
+scatter (NonEmpty.Cons k_ ks_) =
+ let go _k [] [] dst = dst
+ go k0 (k1:ks) (a:as) dst =
+ case splitAt k0 dst of
+ (prefix,suffix) -> prefix ++ a : go k1 ks as suffix
+ go _ _ _ _ = error "scatter: inconsistent lengths"
+ in go k_ ks_
+
+
+{- |
+> mapM_ print $ layoutEchelonBlocks $ reducedRowEchelon $ matrixFromRows [[1,3,0,5],[0,0,1,7::Rational]]
+
+>>> layoutEchelonBlocks $ reducedRowEchelon $ matrixFromRows [[0::Rational]]
+Matrix 1 []
+
+>>> layoutEchelonBlocks $ reducedRowEchelon $ matrixFromRows [[1,3,0,5::Rational]]
+Matrix 3 [[3 % 1,0 % 1,5 % 1]]
+
+>>> layoutEchelonBlocks $ reducedRowEchelon $ matrixFromRows [[0,1,3,5::Rational]]
+Matrix 3 [[0 % 1,3 % 1,5 % 1]]
+
+>>> layoutEchelonBlocks $ reducedRowEchelon $ matrixFromRows [[1,3,0,5],[0,0,1,7::Rational]]
+Matrix 2 [[3 % 1,5 % 1],[0 % 1,7 % 1]]
+
+prop> :{
+ forMatrix $ \matrix ->
+ (layoutEchelonBlocks $ reducedRowEchelon $
+ identity (matrixHeight matrix) ||| matrix)
+ ===
+ matrix
+:}
+
+Construct a matrix that is already in Echelon form
+and check that it is preserved by Echelon decomposition.
+
+prop> :{
+ QC.forAll (fmap (flip mod 10) . NonEmpty.mapTail (take 9) <$> QC.arbitrary) $
+ \blockWidths ->
+ QC.forAll (traverse (uncurry genMatrixForSize) $
+ NonEmptyC.zip (0!:[1..]) blockWidths) $
+ \blocksInt ->
+
+ let blocks :: NonEmpty.T [] (Matrix Rational)
+ blocks = fmap (fmap (%1)) blocksInt
+ unitBesidesBlock block =
+ Matrix (matrixWidth block + 1) $
+ zipWith (:)
+ (drop 1 $ Match.replicate (matrixRows block) 0 ++ [1])
+ (matrixRows block)
+ echelonWithInterleavedEye =
+ Fold.foldr1 shortBesidesTall $
+ NonEmpty.mapTail (map unitBesidesBlock) blocks
+ in
+
+ layoutEchelonBlocks (reducedRowEchelon echelonWithInterleavedEye)
+ ===
+ Fold.foldr1 shortBesidesTall blocks
+:}
+-}
+layoutEchelonBlocks :: (Num a) => AlternatingList (Matrix a) a -> Matrix a
+layoutEchelonBlocks = Fold.foldr1 shortBesidesTall . altOuter
+
+
+{-
+/I B\ /0\
+| | * P * x = | |
+\0 0/ \0/
+
+B is a rectangular upper block triangular matrix with non-uniform block size.
+
+The solution is
+
+ /-B\
+P*X = | |
+ \ I/
+
+such that
+
+ /-B\
+(I B) * | | = I*(-B) + B*I = 0
+ \ I/
+
+This is the nullspace:
+
+ /-B\
+P^T * | |
+ \ I/
+
+The null-space vectors have the form
+/* * *\
+|1 0 0|
+|* * *|
+|* * *|
+|* * *|
+|0 1 0|
+|* * *|
+|0 0 1|
+|* * *|
+\* * */
+-}
+{- |
+>>> reducedRowEchelon $ matrixFromRows [[2,1,3::Rational]]
+Matrix 0 [] ./ 2 % 1 /. Matrix 2 [[1 % 2,3 % 2]] ./ []
+
+>>> reducedRowEchelon $ matrixFromRows [[2],[1],[3::Rational]]
+Matrix 0 [] ./ 3 % 1 /. Matrix 0 [[]] ./ []
+
+>>> reducedRowEchelon $ matrixFromRows [[1,0,2],[0,1,3::Rational]]
+Matrix 0 [] ./ 1 % 1 /. Matrix 0 [[]] ./ 1 % 1 /. Matrix 1 [[2 % 1],[3 % 1]] ./ []
+
+>>> reducedRowEchelon $ matrixFromRows [[1,3,0,5],[0,0,1,7::Rational]]
+Matrix 0 [] ./ 1 % 1 /. Matrix 1 [[3 % 1]] ./ 1 % 1 /. Matrix 1 [[5 % 1],[7 % 1]] ./ []
+
+prop> :{
+ forMatrix $ \matrix ->
+ Fold.all matrixValid $ altOuter $ reducedRowEchelon matrix
+:}
+-}
+reducedRowEchelon :: (RealFrac a) => Matrix a -> AlternatingList (Matrix a) a
+reducedRowEchelon matrix@(Matrix _ []) = altSingleton matrix
+reducedRowEchelon (Matrix _ rows) =
+ let go matrix@(Zipper0 upper _lower) =
+ case rowReduction matrix of
+ (n, Nothing) -> altSingleton (Matrix n $ reverse upper)
+ (n, Just (hd, block, remaining)) ->
+ Matrix n block ./ hd /. go remaining
+ in go $ zipper0FromList rows
+
+{- |
+>>> rowReduction $ Zipper0 [[1,2,3,4,5],[6,7,8,9,10::Rational]] [[0,2,0,0,0],[3,0,0,0,0]]
+(0,Just (3 % 1,[[],[]],Zipper0 [[0 % 1,0 % 1,0 % 1,0 % 1],[2 % 1,3 % 1,4 % 1,5 % 1],[7 % 1,8 % 1,9 % 1,10 % 1]] [[2 % 1,0 % 1,0 % 1,0 % 1]]))
+-}
+rowReduction ::
+ (RealFrac a) => Zipper0 [a] -> (Int, Maybe (a, [[a]], Zipper0 [a]))
+rowReduction (Zipper0 upper []) =
+ (ListHT.switchL 0 (\row _rows -> length row) upper, Nothing)
+rowReduction (Zipper0 upper (focus:lower)) =
+ case shiftUntilNonZeroColumn 0 $ Zipper1 upper focus lower of
+ (n, Nothing) -> (n, Nothing)
+ (n, Just (Zipper1 upperNE (NonEmpty.Cons maxRowHead maxRow) lowerNE)) ->
+ (n, Just $
+ let maxRowNormalized = map (/maxRowHead) maxRow in
+ (maxRowHead,
+ map (take n) $ reverse upper,
+ Zipper0
+ (maxRowNormalized :
+ map (cancelRow maxRowNormalized) upperNE)
+ (map (cancelRow maxRowNormalized) lowerNE)))
+
+cancelRow :: (Num a) => List0 a -> List1 a -> List0 a
+cancelRow xs (NonEmpty.Cons y0 ys) = zipWith (-) ys $ map (y0*) xs
+
+
+shiftUntilNonZeroColumn ::
+ (Real a) => Int -> Zipper1 [a] -> (Int, Maybe (Zipper1 (List1 a)))
+shiftUntilNonZeroColumn n matrix =
+ case checkHeads matrix of
+ Nothing -> (n, Nothing)
+ Just matrixNE@(Zipper1 upperNE foucsNE lowerNE) ->
+ case maxHead (foucsNE!:lowerNE) of
+ (maxRow, lowerWithoutMax) ->
+ if NonEmpty.head maxRow == 0
+ then (shiftUntilNonZeroColumn $! (n+1)) $
+ fmap NonEmpty.tail matrixNE
+ else (n, Just (Zipper1 upperNE maxRow lowerWithoutMax))
+
+maxHead :: (Real a) => List1 (List1 a) -> (List1 a, [List1 a])
+maxHead = NonEmpty.maximumKey (abs . NonEmpty.head . fst) . NonEmpty.removeEach
+
+checkHeads :: Zipper1 [a] -> Maybe (Zipper1 (List1 a))
+checkHeads (Zipper1 upper focus lower) =
+ case NonEmpty.fetch focus of
+ Nothing ->
+ if all null (upper++lower) then Nothing else error "row lengths differ"
+ Just xs -> Just $
+ Zipper1
+ (fromMaybe (error "upper row lengths differ") $ allNonEmpty upper)
+ xs
+ (fromMaybe (error "lower row lengths differ") $ allNonEmpty lower)
+
+allNonEmpty :: [List0 a] -> Maybe [List1 a]
+allNonEmpty xs =
+ case ListHT.partitionMaybe NonEmpty.fetch xs of
+ (ne,[]) -> Just ne
+ _ -> Nothing
diff --git a/Tensor.lhs b/Tensor.lhs
--- a/Tensor.lhs
+++ b/Tensor.lhs
@@ -122,6 +122,9 @@
> module Tensor where
> import Data.Array(inRange)
+> import Prelude2010
+> import Prelude ()
+>
> infixl 9 # -- used for tensor indexing
> infixl 9 ## -- used for indices expressed as lists
> infixl 7 <*> -- inner product with one bound
diff --git a/numeric-quest.cabal b/numeric-quest.cabal
--- a/numeric-quest.cabal
+++ b/numeric-quest.cabal
@@ -1,5 +1,5 @@
Name: numeric-quest
-Version: 0.1.1.3
+Version: 0.2.1
License: GPL
License-File: LICENSE
Author: Jan Skibinski
@@ -7,9 +7,17 @@
Homepage: http://www.haskell.org/haskellwiki/Numeric_Quest
Category: Math
Synopsis: Math and quantum mechanics
-Description: List based linear algebra, similtaneous linear equations, eigenvalues and eigenvectors, roots of polynomials, transcendent functions with arbitrary precision implemented by continued fractions, quantum operations, tensors
-Tested-With: GHC==6.4.1, GHC==6.6.1, GHC==6.8.2, GHC==6.10.4
-Cabal-Version: >=1.6
+Description:
+ List based linear algebra,
+ similtaneous linear equations,
+ eigenvalues and eigenvectors,
+ roots of polynomials,
+ transcendent functions with arbitrary precision
+ implemented by continued fractions,
+ quantum operations,
+ tensors
+Tested-With: GHC==6.*, GHC==7.*, GHC==8.*, GHC==9.*
+Cabal-Version: >=1.10
Build-Type: Simple
Data-Files:
@@ -17,7 +25,7 @@
README
Source-Repository this
- Tag: 0.1.1.3
+ Tag: 0.2.1
Type: darcs
Location: http://code.haskell.org/~thielema/numeric-quest/
@@ -25,18 +33,15 @@
Type: darcs
Location: http://code.haskell.org/~thielema/numeric-quest/
-Flag splitBase
- description: Choose the new smaller, split-up base package.
-
Library
- If flag(splitBase)
- Build-Depends:
- array >=0.1 && <0.5,
- base >= 2 && <5
- Else
- Build-Depends:
- base >= 1.0 && < 2
+ Build-Depends:
+ non-empty >=0.3.2 && <0.4,
+ utility-ht >=0.0.2 && <0.1,
+ prelude-compat >=0.0.0.1 && <0.1,
+ array >=0.1 && <0.6,
+ base >=3 && <5
+ Default-Language: Haskell98
GHC-Options: -Wall
Hs-source-dirs: .
Exposed-modules:
@@ -47,4 +52,22 @@
Orthogonals
QuantumVector
Roots
+ RowEchelon
Tensor
+
+Test-Suite numeric-quest-test
+ Type: exitcode-stdio-1.0
+ Default-Language: Haskell98
+ GHC-Options: -Wall
+ Hs-Source-Dirs: test
+ Other-Modules: Test.RowEchelon
+ Main-Is: Test/Main.hs
+
+ Build-Depends:
+ numeric-quest,
+ doctest-exitcode-stdio >=0.0 && <0.1,
+ doctest-lib >=0.1 && <0.2,
+ QuickCheck >=2 && <3,
+ non-empty,
+ utility-ht,
+ base
diff --git a/test/Test/Main.hs b/test/Test/Main.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Main.hs
@@ -0,0 +1,10 @@
+-- Do not edit! Automatically created with doctest-extract.
+module Main where
+
+import qualified Test.RowEchelon
+
+import qualified Test.DocTest.Driver as DocTest
+
+main :: IO ()
+main = DocTest.run $ do
+ Test.RowEchelon.test
diff --git a/test/Test/RowEchelon.hs b/test/Test/RowEchelon.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/RowEchelon.hs
@@ -0,0 +1,212 @@
+-- Do not edit! Automatically created with doctest-extract from RowEchelon.hs
+{-# LINE 26 "RowEchelon.hs" #-}
+
+module Test.RowEchelon where
+
+import Test.DocTest.Base
+import qualified Test.DocTest.Driver as DocTest
+
+{-# LINE 27 "RowEchelon.hs" #-}
+import RowEchelon as Matrix
+import qualified Test.QuickCheck as QC
+import Test.QuickCheck ((===))
+import Control.Monad (replicateM)
+import qualified Data.Foldable as Fold
+import qualified Data.NonEmpty.Class as NonEmptyC
+import qualified Data.NonEmpty as NonEmpty
+import qualified Data.List.Match as Match
+import qualified Data.List.HT as ListHT
+import Data.NonEmpty ((!:))
+import Data.Ratio ((%))
+
+genElementUniform, genElementNearZero :: QC.Gen Integer
+genElementUniform = QC.choose (-10,10)
+genElementNearZero = fmap (flip rem 11) $ QC.arbitrary
+
+genMatrixForSize :: Int -> Int -> QC.Gen (Matrix Integer)
+genMatrixForSize m n = do
+ fmap (Matrix n) $ replicateM m $ replicateM n genElementNearZero
+
+genMatrix :: QC.Gen (Matrix Integer)
+genMatrix = do
+ m <- QC.choose (0,10)
+ n <- QC.choose (0,10)
+ genMatrixForSize m n
+
+shrinkMatrix :: (Eq a) => Matrix a -> [Matrix a]
+shrinkMatrix matrix@(Matrix width rows) =
+ filter (matrix/=) $
+ map (Matrix width . snd) (ListHT.removeEach rows)
+ ++
+ [Matrix (width-1) $ map (drop 1) rows]
+
+forMatrix :: (QC.Testable test) => (Matrix Rational -> test) -> QC.Property
+forMatrix prop =
+ QC.forAllShrink genMatrix shrinkMatrix (prop . rationalMatrix)
+
+rationalMatrix :: Matrix Integer -> Matrix Rational
+rationalMatrix = fmap (%1)
+
+test :: DocTest.T ()
+test = do
+ DocTest.printPrefix "RowEchelon:142: "
+{-# LINE 142 "RowEchelon.hs" #-}
+ DocTest.example(
+{-# LINE 142 "RowEchelon.hs" #-}
+ nullspace (Matrix 2 []) :: Matrix Rational
+ )
+ [ExpectedLine [LineChunk "Matrix 2 [[1 % 1,0 % 1],[0 % 1,1 % 1]]"]]
+ DocTest.printPrefix "RowEchelon:146: "
+{-# LINE 146 "RowEchelon.hs" #-}
+ DocTest.property(
+{-# LINE 146 "RowEchelon.hs" #-}
+ forMatrix $ matrixValid . nullspace
+ )
+ DocTest.printPrefix "RowEchelon:149: "
+{-# LINE 149 "RowEchelon.hs" #-}
+ DocTest.property(
+{-# LINE 149 "RowEchelon.hs" #-}
+
+ forMatrix $ \matrix ->
+ matrixWidth matrix == matrixHeight (nullspace matrix)
+ )
+ DocTest.printPrefix "RowEchelon:154: "
+{-# LINE 154 "RowEchelon.hs" #-}
+ DocTest.property(
+{-# LINE 154 "RowEchelon.hs" #-}
+
+ forMatrix $ \matrix ->
+ matrixWidth (nullspace matrix) <= matrixWidth matrix
+ )
+ DocTest.printPrefix "RowEchelon:161: "
+{-# LINE 161 "RowEchelon.hs" #-}
+ DocTest.property(
+{-# LINE 161 "RowEchelon.hs" #-}
+
+ forMatrix $ \matrix ->
+ matrixWidth matrix <= matrixWidth (nullspace matrix) + matrixHeight matrix
+ )
+ DocTest.printPrefix "RowEchelon:166: "
+{-# LINE 166 "RowEchelon.hs" #-}
+ DocTest.property(
+{-# LINE 166 "RowEchelon.hs" #-}
+
+ forMatrix $ \matrix ->
+ Fold.all (0==) $ matrixProduct matrix (nullspace matrix)
+ )
+ DocTest.printPrefix "RowEchelon:224: "
+{-# LINE 224 "RowEchelon.hs" #-}
+ DocTest.example(
+{-# LINE 224 "RowEchelon.hs" #-}
+ scatter (3 !: 1 : 4 : 1 : []) "012" ['a'..'k']
+ )
+ [ExpectedLine [LineChunk "\"abc0d1efgh2ijk\""]]
+ DocTest.printPrefix "RowEchelon:240: "
+{-# LINE 240 "RowEchelon.hs" #-}
+ DocTest.example(
+{-# LINE 240 "RowEchelon.hs" #-}
+ layoutEchelonBlocks $ reducedRowEchelon $ matrixFromRows [[0::Rational]]
+ )
+ [ExpectedLine [LineChunk "Matrix 1 []"]]
+ DocTest.printPrefix "RowEchelon:243: "
+{-# LINE 243 "RowEchelon.hs" #-}
+ DocTest.example(
+{-# LINE 243 "RowEchelon.hs" #-}
+ layoutEchelonBlocks $ reducedRowEchelon $ matrixFromRows [[1,3,0,5::Rational]]
+ )
+ [ExpectedLine [LineChunk "Matrix 3 [[3 % 1,0 % 1,5 % 1]]"]]
+ DocTest.printPrefix "RowEchelon:246: "
+{-# LINE 246 "RowEchelon.hs" #-}
+ DocTest.example(
+{-# LINE 246 "RowEchelon.hs" #-}
+ layoutEchelonBlocks $ reducedRowEchelon $ matrixFromRows [[0,1,3,5::Rational]]
+ )
+ [ExpectedLine [LineChunk "Matrix 3 [[0 % 1,3 % 1,5 % 1]]"]]
+ DocTest.printPrefix "RowEchelon:249: "
+{-# LINE 249 "RowEchelon.hs" #-}
+ DocTest.example(
+{-# LINE 249 "RowEchelon.hs" #-}
+ layoutEchelonBlocks $ reducedRowEchelon $ matrixFromRows [[1,3,0,5],[0,0,1,7::Rational]]
+ )
+ [ExpectedLine [LineChunk "Matrix 2 [[3 % 1,5 % 1],[0 % 1,7 % 1]]"]]
+ DocTest.printPrefix "RowEchelon:252: "
+{-# LINE 252 "RowEchelon.hs" #-}
+ DocTest.property(
+{-# LINE 252 "RowEchelon.hs" #-}
+
+ forMatrix $ \matrix ->
+ (layoutEchelonBlocks $ reducedRowEchelon $
+ identity (matrixHeight matrix) ||| matrix)
+ ===
+ matrix
+ )
+ DocTest.printPrefix "RowEchelon:263: "
+{-# LINE 263 "RowEchelon.hs" #-}
+ DocTest.property(
+{-# LINE 263 "RowEchelon.hs" #-}
+
+ QC.forAll (fmap (flip mod 10) . NonEmpty.mapTail (take 9) <$> QC.arbitrary) $
+ \blockWidths ->
+ QC.forAll (traverse (uncurry genMatrixForSize) $
+ NonEmptyC.zip (0!:[1..]) blockWidths) $
+ \blocksInt ->
+
+ let blocks :: NonEmpty.T [] (Matrix Rational)
+ blocks = fmap (fmap (%1)) blocksInt
+ unitBesidesBlock block =
+ Matrix (matrixWidth block + 1) $
+ zipWith (:)
+ (drop 1 $ Match.replicate (matrixRows block) 0 ++ [1])
+ (matrixRows block)
+ echelonWithInterleavedEye =
+ Fold.foldr1 shortBesidesTall $
+ NonEmpty.mapTail (map unitBesidesBlock) blocks
+ in
+
+ layoutEchelonBlocks (reducedRowEchelon echelonWithInterleavedEye)
+ ===
+ Fold.foldr1 shortBesidesTall blocks
+ )
+ DocTest.printPrefix "RowEchelon:329: "
+{-# LINE 329 "RowEchelon.hs" #-}
+ DocTest.example(
+{-# LINE 329 "RowEchelon.hs" #-}
+ reducedRowEchelon $ matrixFromRows [[2,1,3::Rational]]
+ )
+ [ExpectedLine [LineChunk "Matrix 0 [] ./ 2 % 1 /. Matrix 2 [[1 % 2,3 % 2]] ./ []"]]
+ DocTest.printPrefix "RowEchelon:332: "
+{-# LINE 332 "RowEchelon.hs" #-}
+ DocTest.example(
+{-# LINE 332 "RowEchelon.hs" #-}
+ reducedRowEchelon $ matrixFromRows [[2],[1],[3::Rational]]
+ )
+ [ExpectedLine [LineChunk "Matrix 0 [] ./ 3 % 1 /. Matrix 0 [[]] ./ []"]]
+ DocTest.printPrefix "RowEchelon:335: "
+{-# LINE 335 "RowEchelon.hs" #-}
+ DocTest.example(
+{-# LINE 335 "RowEchelon.hs" #-}
+ reducedRowEchelon $ matrixFromRows [[1,0,2],[0,1,3::Rational]]
+ )
+ [ExpectedLine [LineChunk "Matrix 0 [] ./ 1 % 1 /. Matrix 0 [[]] ./ 1 % 1 /. Matrix 1 [[2 % 1],[3 % 1]] ./ []"]]
+ DocTest.printPrefix "RowEchelon:338: "
+{-# LINE 338 "RowEchelon.hs" #-}
+ DocTest.example(
+{-# LINE 338 "RowEchelon.hs" #-}
+ reducedRowEchelon $ matrixFromRows [[1,3,0,5],[0,0,1,7::Rational]]
+ )
+ [ExpectedLine [LineChunk "Matrix 0 [] ./ 1 % 1 /. Matrix 1 [[3 % 1]] ./ 1 % 1 /. Matrix 1 [[5 % 1],[7 % 1]] ./ []"]]
+ DocTest.printPrefix "RowEchelon:341: "
+{-# LINE 341 "RowEchelon.hs" #-}
+ DocTest.property(
+{-# LINE 341 "RowEchelon.hs" #-}
+
+ forMatrix $ \matrix ->
+ Fold.all matrixValid $ altOuter $ reducedRowEchelon matrix
+ )
+ DocTest.printPrefix "RowEchelon:357: "
+{-# LINE 357 "RowEchelon.hs" #-}
+ DocTest.example(
+{-# LINE 357 "RowEchelon.hs" #-}
+ rowReduction $ Zipper0 [[1,2,3,4,5],[6,7,8,9,10::Rational]] [[0,2,0,0,0],[3,0,0,0,0]]
+ )
+ [ExpectedLine [LineChunk "(0,Just (3 % 1,[[],[]],Zipper0 [[0 % 1,0 % 1,0 % 1,0 % 1],[2 % 1,3 % 1,4 % 1,5 % 1],[7 % 1,8 % 1,9 % 1,10 % 1]] [[2 % 1,0 % 1,0 % 1,0 % 1]]))"]]