numeric-quest 0.1.1.3 → 0.2.1
raw patch · 12 files changed
Files
- Eigensystem.hs +2/−0
- EigensystemNum.hs +10/−8
- Fraction.hs +2/−0
- Makefile +11/−0
- Orthogonals.lhs +15/−29
- QuantumVector.lhs +3/−1
- Roots.hs +2/−10
- RowEchelon.hs +412/−0
- Tensor.lhs +3/−0
- numeric-quest.cabal +38/−15
- test/Test/Main.hs +10/−0
- test/Test/RowEchelon.hs +212/−0
Eigensystem.hs view
@@ -54,6 +54,8 @@ import QuantumVector import LinearAlgorithms (triangular, tridiagonal, triangular2) import Data.List (findIndex)+import Prelude2010+import Prelude () ---------------------------------------------------------------------------- -- Category: Eigensystem for QuantumVector
EigensystemNum.hs view
@@ -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"
Fraction.hs view
@@ -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)
Makefile view
@@ -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
Orthogonals.lhs view
@@ -304,7 +304,7 @@ by some algorithms. So far we have been able to avoid this. <pre> -> 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 </pre> <hr>@@ -783,7 +779,7 @@ <pre> -> 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 @@ <pre> -> 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 @@ <pre> -> 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 @@ <pre> -> 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 =
QuantumVector.lhs view
@@ -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 ") "
Roots.hs view
@@ -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) ----------------------------------------------------------------------------- --
+ RowEchelon.hs view
@@ -0,0 +1,412 @@+{- |+<https://en.wikipedia.org/wiki/Row_echelon_form>++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
Tensor.lhs view
@@ -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
numeric-quest.cabal view
@@ -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
+ test/Test/Main.hs view
@@ -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
+ test/Test/RowEchelon.hs view
@@ -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]]))"]]