Vec (empty) → 0.9.0
raw patch · 9 files changed
+1685/−0 lines, 9 filesdep +basesetup-changed
Dependencies added: base
Files
- Data/Vec.hs +137/−0
- Data/Vec/Base.hs +416/−0
- Data/Vec/Instances.hs +117/−0
- Data/Vec/LinAlg.hs +726/−0
- Data/Vec/Nat.hs +69/−0
- Data/Vec/Packed.hs +152/−0
- LICENSE +22/−0
- Setup.hs +4/−0
- Vec.cabal +42/−0
+ Data/Vec.hs view
@@ -0,0 +1,137 @@+{-+Copyright (c) 2008, Scott E. Dillard+All rights reserved.++Redistribution and use in source and binary forms, with or without modification,+are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright notice,+ this list of conditions and the following disclaimer.+ * Redistributions in binary form must reproduce the above copyright notice,+ this list of conditions and the following disclaimer in the documentation+ and/or other materials provided with the distribution.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON+ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.+-}++{- |++Vec : a library for fixed-length lists and low-dimensional linear algebra++Scott E. Dillard <sedillard@gmail.com>++darcs : <http://graphics.cs.ucdavis.edu/~sdillard/Vec>++/Synopsis/++Vectors are represented by lists with type-encoded lengths. The constructor is+@:.@, which acts like a cons both at the value and type levels, with @()@+taking the place of nil. So @x:.y:.z:.()@ is a 3d vector. The library provides+a set of common list-like functions (map, fold, etc) for working with vectors.+Built up from these functions are a small but useful set of linear algebra+operations: matrix multiplication, determinants, solving linear systems,+inverting matrices.++/Design/+++* Simplicity : +Beyond the initial complexities of type-level lists and+numbers, I've tried to keep the API simple. There is no vector-space+class, nor a complicated hierarchy of linear\/affine\/projective+transformations. These can be added on top of the library easily.++* Purity :+The library is written in the functional style. For most+functions this does not hinder performance at all, but some I am still+working on (Gaussian elimination) so if this library is a bottleneck you+can easily drop down to C. ++* Low Dimension :+Although the dimensionality is limited only by what GHC+will handle, the library is meant for 2,3 and 4 dimensions. For general+linear algebra, check out the excellent hmatrix library and blas bindings.++To the point of simplicity, vectors and matrices are instances of Num and+Fractional. All arithmetic is done component-wise and literals construct+uniform vectors and matrices. There are many interesting projects aiming to+overhaul Haskell's number classes, but for now the type of @(*)@ is @a -> a ->+a@ so that's what we're working with. It is easy to incorporate this library+into a more mathematically consistent class hierarchy (provided you can design+one.) ++The rule is simple : + If the method is unary, it's a map. + If it's binary, it's a zipWith.++/Performance/++@(:.)@ is strict in both arguments, but it is also polymorphic, so at runtime+vectors will be realized as linked lists, albeit with less pattern matching.+However the library provides packed representations for 2,3 and 4d vectors of+Ints, Floats and Doubles. @'Vec3F' x y z@ constructs a packed vector of+unboxed Floats. Functions @'pack'@ and @'unpack'@ convert between packed and+unpacked types. When vector operations are bracketed by 'pack' and 'unpack',+GHC can unfold them into very efficient code. The 'Storable' instances for+vectors also generate fast code. Without optimizations, the code falls back+into linked-list mode. The optimizations depend on inlining, so you may need+to increase your unfolding threshold in certain situations.++/GHC Extensions/++This library makes heavy use of functional dependencies. I have tried to+tweak things so that they \"just work.\" However, every now and then you will+get incomprehensible error messages, usually about how this isn't an+instance of that. These are how type errors typically manifest, so first+double check to make sure you aren't trying to mix vectors of different+dimension or component types. If you still get these errors, manual type+annotations usually make them go away.+++/Related Work/++See previous work by David Menendez,+ <http://haskell.org/pipermail/haskell/2005-May/015815.html>++and of course Oleg Kiselyov,+ <http://okmij.org/ftp/papers/number-parameterized-types.pdf>++Other vector and linear algebra packages :++vector-space, by Conal Elliott : + <http://hackage.haskell.org/cgi-bin/hackage-scripts/package/vector-space>++hmatrix, by Alberto Ruiz :+ <http://hackage.haskell.org/cgi-bin/hackage-scripts/package/hmatrix>++blas bindings, by Patrick Perry :+ <http://hackage.haskell.org/cgi-bin/hackage-scripts/package/blas>++templatized geometry library (C++), by Oliver Kreylos :+ <http://graphics.cs.ucdavis.edu/~okreylos/ResDev/Geometry/index.html>+-}++module Data.Vec + (module Data.Vec.Base+ ,module Data.Vec.LinAlg+ ,module Data.Vec.Packed+ ,module Data.Vec.Nat+ )+where++import Data.Vec.Base+import Data.Vec.LinAlg+import Data.Vec.Packed+import Data.Vec.Nat+import Data.Vec.Instances++
+ Data/Vec/Base.hs view
@@ -0,0 +1,416 @@+{- Copyright (c) 2008, Scott E. Dillard. All rights reserved. -}++{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE EmptyDataDecls #-}+{-# LANGUAGE ExistentialQuantification #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE UndecidableInstances #-}++{-# HADDOCK_OPTIONS prune #-}++module Data.Vec.Base where++import Data.Vec.Nat++import Prelude hiding (map,zipWith,foldl,foldr,reverse,+ take,drop,head,tail,sum,last,product,+ minimum,maximum)+import qualified Prelude as P++++-- | The vector constructor. @(:.)@ for vectors is like @(:)@ for lists, and+-- @()@ takes the place of @[]@. (The list of instances here is not meant to be+-- readable.)++data a :. b = !a :. !b+ deriving (Eq,Ord,Read)++infixr :.++--derived show outputs in prefix notation+instance (Show a, ShowVec v) => Show (a:.v) where+ show (a:.v) = "(" ++ show a ++ ":." ++ showVec v ++ ")"+++-- | Helper to keep parentheses at bay. Just use @show@ as usual.+class ShowVec v where+ showVec :: v -> String++instance ShowVec () where+ showVec = show+ {-# INLINE showVec #-}++instance (Show a, ShowVec v) => ShowVec (a:.v) where+ showVec (a:.v) = show a ++ ":." ++ showVec v+ {-# INLINE showVec #-}+++-- * Vector Types+type Vec2 a = a :. a :. ()+type Vec3 a = a :. (Vec2 a)+type Vec4 a = a :. (Vec3 a)+type Vec5 a = a :. (Vec4 a)+type Vec6 a = a :. (Vec5 a)+type Vec7 a = a :. (Vec6 a)+type Vec8 a = a :. (Vec7 a)+type Vec9 a = a :. (Vec8 a)+type Vec10 a = a :. (Vec9 a)+type Vec11 a = a :. (Vec10 a)+type Vec12 a = a :. (Vec11 a)+type Vec13 a = a :. (Vec12 a)+type Vec14 a = a :. (Vec13 a)+type Vec15 a = a :. (Vec14 a)+type Vec16 a = a :. (Vec15 a)+type Vec17 a = a :. (Vec16 a)+type Vec18 a = a :. (Vec17 a)+type Vec19 a = a :. (Vec18 a)+++++-- | The type constraint @Vec n a v@ infers the vector type @v@ from the+-- length @n@, a type-level natural, and underlying component type @a@. +-- So @x :: Vec N4 a v => v@ declares @x@ to be a 4-vector of @a@s.++class Vec n a v | n a -> v, v -> n a where+ -- | Make a uniform vector of a given length. @n@ is a type-level natural.+ -- Use `vec` when the length can be inferred.+ mkVec :: n -> a -> v++ -- | turn a list into a vector of inferred length+ fromList :: [a] -> v++ -- | get a vector element, which one is determined at runtime+ getElem :: Int -> v -> a++ -- | set a vector element, which one is determined at runtime+ setElem :: Int -> a -> v -> v++instance Vec N1 a ( a :. () ) where+ mkVec _ a = a :. ()+ fromList (a:_) = a :. ()+ fromList [] = error "fromList: list too short"+ getElem !i (a :. _) + | i == 0 = a+ | otherwise = error "getElem: index out of bounds"+ setElem !i a _ + | i == 0 = a :. ()+ | otherwise = error "setElem: index out of bounds"+ {-# INLINE setElem #-}+ {-# INLINE getElem #-}+ {-# INLINE mkVec #-}+ {-# INLINE fromList #-}++instance Vec (Succ n) a (a':.v) => Vec (Succ (Succ n)) a (a:.a':.v) where+ mkVec _ a = a :. (mkVec undefined a)+ fromList (a:as) = a :. (fromList as)+ fromList [] = error "fromList: list too short"+ getElem !i (a :. v)+ | i == 0 = a+ | otherwise = getElem (i-1) v+ setElem !i a (x :. v)+ | i == 0 = a :. v+ | otherwise = x :. (setElem (i-1) a v)+ {-# INLINE setElem #-}+ {-# INLINE getElem #-}+ {-# INLINE mkVec #-}+ {-# INLINE fromList #-}+++-- | Make a uniform vector. The length is inferred.+vec :: (Vec n a v) => a -> v+vec = mkVec undefined+{-# INLINE vec #-}+++-- | get or set a vector element, known at compile+--time. Use the Nat types to access vector components. For instance, @get n0@+--gets the x component, @set n2 44@ sets the z component to 44. +++class Access n a v | v -> a where+ get :: n -> v -> a+ set :: n -> a -> v -> v++instance Access N0 a (a :. v) where+ get _ (a :. _) = a+ set _ a (_ :. v) = a :. v+ {-# INLINE set #-}+ {-# INLINE get #-}++instance Access n a v => Access (Succ n) a (a :. v) where+ get _ (_ :. v) = get (undefined::n) v+ set _ a' (a :. v) = a :. (set (undefined::n) a' v)+ {-# INLINE set #-}+ {-# INLINE get #-}+++-- * List-like functions++-- | The first element.++class Head v a | v -> a where + head :: v -> a++instance Head (a :. as) a where + head (a :. _) = a+ {-# INLINE head #-}+++-- | All but the first element. ++class Tail v v_ | v -> v_ where + tail :: v -> v_++instance Tail (a :. as) as where + tail (_ :. as) = as+ {-# INLINE tail #-}+++++-- | Apply a function over each element in a vector. Constraint @Map a b u v@+-- states that @u@ is a vector of @a@s, @v@ is a vector of @b@s with the same+-- length as @u@, and the function is of type @a -> b@.++class Map a b u v | u -> a, v -> b, b u -> v, a v -> u where+ map :: (a -> b) -> u -> v++instance Map a b (a :. ()) (b :. ()) where+ map f (x :. ()) = (f x) :. ()+ {-# INLINE map #-}++instance Map a b (a':.u) (b':.v) => Map a b (a:.a':.u) (b:.b':.v) where+ map f (x:.v) = (f x):.(map f v)+ {-# INLINE map #-}+++++-- | Combine two vectors using a binary function. The length of the result is+-- the min of the lengths of the arguments. The constraint @ZipWith a b c u v+-- w@ states that @u@ is a vector of @a@s, @v@ is a vector of @b@s, @w@ is a+-- vector of @c@s, and the binary function is of type @a -> b -> c@.++class ZipWith a b c u v w | u->a, v->b, w->c, u v c -> w where+ zipWith :: (a -> b -> c) -> u -> v -> w++instance ZipWith a b c (a:.()) (b:.()) (c:.()) where+ zipWith f (x:._) (y:._) = f x y :.()+ {-# INLINE zipWith #-}++instance ZipWith a b c (a:.()) (b:.b:.bs) (c:.()) where+ zipWith f (x:._) (y:._) = f x y :.()+ {-# INLINE zipWith #-}++instance ZipWith a b c (a:.a:.as) (b:.()) (c:.()) where+ zipWith f (x:._) (y:._) = f x y :.()+ {-# INLINE zipWith #-}++instance + ZipWith a b c (a':.u) (b':.v) (c':.w) + => ZipWith a b c (a:.a':.u) (b:.b':.v) (c:.c':.w) + where+ zipWith f (x:.u) (y:.v) = f x y :. zipWith f u v+ {-# INLINE zipWith #-}+++-- | Fold a function over a vector. ++class Fold a v | v -> a where+ fold :: (a -> a -> a) -> v -> a+ foldl :: (b -> a -> b) -> b -> v -> b+ foldr :: (a -> b -> b) -> b -> v -> b++instance Fold a (a:.()) where+ fold f (a:._) = a + foldl f z (a:._) = (f $! z) $! a+ foldr f z (a:._) = (f $! a) $! z+ {-# INLINE fold #-}+ {-# INLINE foldl #-}+ {-# INLINE foldr #-}++instance Fold a (a':.u) => Fold a (a:.a':.u) where+ fold f (a:.v) = (f $! a) $! (fold f v)+ foldl f z (a:.v) = (f $! (foldl f z v)) $! a+ foldr f z (a:.v) = (f $! a) $! (foldr f z v)+ {-# INLINE fold #-}+ {-# INLINE foldl #-}+ {-# INLINE foldr #-}++-- | Reverse a vector +reverse v = reverse' () v+{-# INLINE reverse #-}++-- Reverse helper function : builds the reversed list as its first argument+class Reverse' p v v' | p v -> v' where+ reverse' :: p -> v -> v'+ +instance Reverse' p () p where+ reverse' p () = p+ {-# INLINE reverse' #-}++instance Reverse' (a:.p) v v' => Reverse' p (a:.v) v' where+ reverse' p (a:.v) = reverse' (a:.p) v + {-# INLINE reverse' #-}+++-- | Append two vectors ++class Append v1 v2 v3 | v1 v2 -> v3, v1 v3 -> v2 where + append :: v1 -> v2 -> v3++instance Append () v v where+ append _ = id+ {-# INLINE append #-}++instance Append (a:.()) v (a:.v) where+ append (a:.()) v = a:.v+ {-# INLINE append #-}++instance (Append (a':.v1) v2 v3) => Append (a:.a':.v1) v2 (a:.v3) where+ append (a:.u) v = a:.(append u v)+ {-# INLINE append #-}++++-- | @take n v@ constructs a vector from the first @n@ elements of @v@. @n@ is a+-- type-level natural. For example @take n3 v@ makes a 3-vector of the first+-- three elements of @v@.++class Take n v v' | n v -> v', n v' -> v where+ take :: n -> v -> v'++instance Take N0 v () where+ take _ _ = ()+ {-# INLINE take #-}++instance Take n v v' => Take (Succ n) (a:.v) (a:.v') where+ take _ (a:.v) = a:.(take (undefined::n) v)+ {-# INLINE take #-}+++-- | @drop n v@ strips the first @n@ elements from @v@. @n@ is a type-level+-- natural. For example @drop n2 v@ drops the first two elements.++class Drop n v v' | n v -> v', n v' -> v where+ drop :: n -> v -> v'+ +instance Drop N0 v v where+ drop _ = id+ {-# INLINE drop #-}++instance (Tail v' v'', Drop n v v') => Drop (Succ n) v v'' where+ drop _ = tail . drop (undefined::n)+ {-# INLINE drop #-}+++-- | Get the last element, usually significant for some reason (quaternions,+-- homogenous coordinates, whatever)+class Last v a | v -> a where+ last :: v -> a++instance Last (a:.()) a where + last (a:._) = a+ {-# INLINE last #-}++instance Last (a':.v) a => Last (a:.a':.v) a where+ last (a:.v) = last v+ {-# INLINE last #-}++-- | @snoc v a@ appends the element a to the end of v. ++class Snoc v a v' | v a -> v', v' -> v a where + snoc :: v -> a -> v'++instance Snoc () a (a:.()) where+ snoc _ a = (a:.())+ {-# INLINE snoc #-}++instance Snoc v a (a:.v) => Snoc (a:.v) a (a:.a:.v) where+ snoc (b:.v) a = b:.(snoc v a)+ {-# INLINE snoc #-}++++-- | sum of vector elements+sum :: (Fold a v, Num a) => v -> a+sum x = fold (+) x+{-# INLINE sum #-}++-- | product of vector elements+product :: (Fold a v, Num a) => v -> a+product x = fold (*) x+{-# INLINE product #-}++-- | maximum vector element+maximum :: (Fold a v, Ord a) => v -> a+maximum x = fold max x+{-# INLINE maximum #-}++-- | minimum vector element+minimum :: (Fold a v, Ord a) => v -> a+minimum x = fold min x+{-# INLINE minimum #-}++toList :: (Fold a v) => v -> [a]+toList = foldr (:) [] +{-# INLINE toList #-}++++++++-- * Matrix Types++type Mat22 a = Vec2 (Vec2 a)+type Mat23 a = Vec2 (Vec3 a)+type Mat24 a = Vec2 (Vec4 a)++type Mat32 a = Vec3 (Vec2 a)+type Mat33 a = Vec3 (Vec3 a)+type Mat34 a = Vec3 (Vec4 a)+type Mat35 a = Vec3 (Vec5 a)+type Mat36 a = Vec3 (Vec6 a)++type Mat42 a = Vec4 (Vec2 a)+type Mat43 a = Vec4 (Vec3 a)+type Mat44 a = Vec4 (Vec4 a)+type Mat45 a = Vec4 (Vec5 a)+type Mat46 a = Vec4 (Vec6 a)+type Mat47 a = Vec4 (Vec7 a)+type Mat48 a = Vec4 (Vec8 a)++-- | convert a matrix to a list-of-lists+matToLists :: (Fold a v, Fold v m) => m -> [[a]]+matToLists = (P.map toList) . toList+{-# INLINE matToLists #-}++-- | convert a matrix to a list in row-major order+matToList :: (Fold a v, Fold v m) => m -> [a]+matToList = concat . matToLists+{-# INLINE matToList #-}++-- | convert a list-of-lists into a matrix+matFromLists :: (Vec j a v, Vec i v m) => [[a]] -> m+matFromLists = fromList . (P.map fromList)+{-# INLINE matFromLists #-}++-- | convert a list into a matrix. (row-major order)+matFromList :: forall i j v m a. (Vec i v m, Vec j a v, Nat i) => [a] -> m+matFromList = matFromLists . groupsOf (nat(undefined::i))+ where groupsOf n xs = let (a,b) = splitAt n xs in a:(groupsOf n b)+{-# INLINE matFromList #-}+++
+ Data/Vec/Instances.hs view
@@ -0,0 +1,117 @@+{- Copyright (c) 2008, Scott E. Dillard. All rights reserved. -}++{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE EmptyDataDecls #-}+{-# LANGUAGE ExistentialQuantification #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE UndecidableInstances #-}++module Data.Vec.Instances where++import Data.Vec.Base as V+import Data.Vec.Nat+import Foreign.Storable+import Foreign.Ptr++-- Storable instances. ++instance Storable a => Storable (a:.()) where+ sizeOf _ = sizeOf (undefined::a)+ alignment _ = alignment (undefined::a)+ peek !p = peek (castPtr p) >>= \a -> return (a:.())+ peekByteOff !p !o = peek (p`plusPtr`o)+ peekElemOff !p !i = peek (p`plusPtr`(i*sizeOf(undefined::a)))+ poke !p (a:._) = poke (castPtr p) a+ pokeByteOff !p !o !x = poke (p`plusPtr`o) x+ pokeElemOff !p !i !x = poke (p`plusPtr`(i*sizeOf(undefined::a))) x+ {-# INLINE sizeOf #-}+ {-# INLINE alignment #-}+ {-# INLINE peek #-}+ {-# INLINE peekByteOff #-}+ {-# INLINE peekElemOff #-}+ {-# INLINE poke #-}+ {-# INLINE pokeByteOff #-}+ {-# INLINE pokeElemOff #-}++instance (Vec (Succ (Succ n)) a (a:.a:.v), Storable a, Storable (a:.v)) + => Storable (a:.a:.v) + where+ sizeOf _ = sizeOf (undefined::a) + sizeOf (undefined::(a:.v))+ alignment _ = alignment (undefined::a)+ peek !p = + peek (castPtr p) >>= \a -> + peek (castPtr (p`plusPtr`sizeOf(undefined::a))) >>= \v -> + return (a:.v)+ peekByteOff !p !o = peek (p`plusPtr`o)+ peekElemOff !p !i = peek (p`plusPtr`(i*sizeOf(undefined::(a:.a:.v))))+ poke !p (a:.v) = + poke (castPtr p) a >> + poke (castPtr (p`plusPtr`sizeOf(undefined::a))) v+ pokeByteOff !p !o !x = poke (p`plusPtr`o) x+ pokeElemOff !p !i !x = poke (p`plusPtr`(i*sizeOf(undefined::(a:.a:.v)))) x+ {-# INLINE sizeOf #-}+ {-# INLINE alignment #-}+ {-# INLINE peek #-}+ {-# INLINE peekByteOff #-}+ {-# INLINE peekElemOff #-}+ {-# INLINE poke #-}+ {-# INLINE pokeByteOff #-}+ {-# INLINE pokeElemOff #-}+++-- Num and Fractional instances : All arithmetic is done component-wise and+-- literals construct uniform vectors and matrices. +--+-- The rule is simple : +-- If the method is unary, it's a map. +-- If it's binary, it's a zipWith.+--+-- You are free to ignore these instances if the definition of (*) offends you.++instance+ (Eq (a:.u)+ ,Show (a:.u)+ ,Num a+ ,Map a a (a:.u) (a:.u) + ,ZipWith a a a (a:.u) (a:.u) (a:.u)+ ,Vec (Succ l) a (a:.u)+ )+ => Num (a:.u) + where+ (+) u v = V.zipWith (+) u v + (-) u v = V.zipWith (-) u v+ (*) u v = V.zipWith (*) u v+ abs u = V.map abs u+ signum u = V.map signum u+ fromInteger i = vec (fromInteger i)+ {-# INLINE (+) #-}+ {-# INLINE (-) #-}+ {-# INLINE (*) #-}+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ {-# INLINE fromInteger #-}+++instance + (Fractional a+ ,Ord (a:.u)+ ,ZipWith a a a (a:.u) (a:.u) (a:.u)+ ,Map a a (a:.u) (a:.u)+ ,Vec (Succ l) a (a:.u)+ ,Show (a:.u)+ ) + => Fractional (a:.u) + where+ (/) u v = V.zipWith (/) u v+ recip u = V.map recip u+ fromRational r = vec (fromRational r)+ {-# INLINE (/) #-}+ {-# INLINE recip #-}+ {-# INLINE fromRational #-}
+ Data/Vec/LinAlg.hs view
@@ -0,0 +1,726 @@+{- Copyright (c) 2008, Scott E. Dillard. All rights reserved. -}++{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE EmptyDataDecls #-}+{-# LANGUAGE ExistentialQuantification #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE PatternSignatures #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE UndecidableInstances #-}++{-# OPTIONS_HADDOCK ignore-exports,prune #-}++module Data.Vec.LinAlg + (dot+ ,normSq+ ,norm+ ,normalize+ ,cross+ ,homPoint+ ,homVec+ ,project+ ,multvm+ ,multmv+ ,multmm+ ,translate+ ,column+ ,row+ ,Transpose(transpose)+ ,SetDiagonal(setDiagonal)+ ,GetDiagonal(getDiagonal)+ ,scale+ ,diagonal+ ,identity+ ,Det(det)+ ,cramer'sRule+ ,NearZero(nearZero)+ ,GaussElim(gaussElim)+ ,BackSubstitute(backSubstitute)+ ,BackSubstitute'(backSubstitute')+ ,invert+ ,invertAndDet+ ,solve+ ) where++import Prelude hiding (map,zipWith,foldl,foldr,reverse,take,drop,+ head,tail,sum,length,last)+import qualified Prelude as P+import Data.Vec.Base+import Data.Vec.Nat+import Data.Vec.Instances++import Control.Monad+import Data.Maybe+++-- | dot / inner / scalar product+dot :: (Num a, Num v, Fold a v) => v -> v -> a+dot u v = sum (u*v)+{-# INLINE dot #-}++-- | vector norm, squared+normSq :: (Num a, Num v, Fold a v) => v -> a+normSq v = dot v v+{-# INLINE normSq #-}++-- | vector / L2 / Euclidean norm+norm :: (Num v, Floating a, Fold a v) => v -> a+norm v = sqrt (dot v v)+{-# INLINE norm #-}++-- | @normalize v@ is a unit vector in the direction of @v@. @v@ is assumed+-- non-null.+normalize :: (Floating a, Num v, Fold a v, Map a a v v) => v -> v+normalize v = map (/(norm v)) v+{-# INLINE normalize #-}++-- | 3d cross product.+cross :: Num a => Vec3 a -> Vec3 a -> Vec3 a+cross (ux:.uy:.uz:.()) (vx:.vy:.vz:.()) =+ (uy*vz-uz*vy):.(uz*vx-ux*vz):.(ux*vy-uy*vx):.()+{-# INLINE cross #-}++-- | lift a point into homogenous coordinates+homPoint :: (Snoc v a v', Num a) => v -> v'+homPoint v = snoc v 1+{-# INLINE homPoint #-}++-- | point-at-infinity in homogenous coordinates+homVec :: (Snoc v a v', Num a) => v -> v'+homVec v = snoc v 0+{-# INLINE homVec #-}++-- | project a vector from homogenous coordinates. Last vector element is+-- assumed non-zero.+project :: + ( Reverse' () t1 v'+ , Fractional t1+ , Vec a t t1+ , Reverse' () v (t :. t1)+ ) => v -> v'+project v = case reverse v of (w:.u) -> reverse (u/vec w)+{-# INLINE project #-}+++-- | row vector * matrix+multvm :: + ( Transpose m mt+ , Map v a mt v'+ , Fold a v+ , Num a+ , Num v+ ) => v -> m -> v'+multvm v m = map (dot v) (transpose m)+{-# INLINE multvm #-}++-- | matrix * column vector+multmv :: + ( Map v a m v'+ , Num v+ , Fold a v+ , Num a+ ) => m -> v -> v'+multmv m v = map (dot v) m+{-# INLINE multmv #-}++-- | matrix * matrix +multmm :: + (Map v v' m1 m3+ ,Map v a b v'+ ,Transpose m2 b+ ,Fold a v+ ,Num v+ ,Num a+ ) => m1 -> m2 -> m3+multmm a b = map (\v -> map (dot v) (transpose b)) a+{-# INLINE multmm #-}++-- | apply a translation to a projective transformation matrix+translate :: + (Transpose m mt+ ,Reverse' () mt (v' :. t)+ ,Reverse' (v' :. ()) t v'1+ ,Transpose v'1 m+ ,Num v'+ ,Num a+ ,Snoc v a v'+ ) => v -> m -> m+translate v m = + case reverse (transpose m) of+ (h:.t) -> transpose (reverse (((homVec v) + h) :. t))+{-# INLINE translate #-}++-- | get the @n@-th column as a vector. @n@ is a type-level natural.+column :: (Transpose m mt, Access n v mt) => n -> m -> v+column n = get n . transpose +{-# INLINE row #-}++-- | get the @n@-th row as a vector. @n@ is a type-level natural.+row :: (Access n a v) => n -> v -> a+row n = get n+{-# INLINE column #-}+++-- Matrix transpose wrapper class: infers type of one argument from the other,+-- because Transpose` can't do it, the fundeps there can't be bijective++-- | matrix transposition+class Transpose a b | a -> b, b -> a where + transpose :: a -> b++instance Transpose () () where+ transpose = id++instance + (Vec (Succ n) s (s:.ra) --(s:ra) is an n-vector of s'es (row of a)+ ,Vec (Succ m) (s:.ra) ((s:.ra):.a) --a is an m-vector of ra's+ ,Vec (Succ m) s (s:.rb) --rb is an m-vector of s'es (row of b)+ ,Vec (Succ n) (s:.rb) ((s:.rb):.b) --b is an n-vector of rb's+ ,Transpose' ((s:.ra):.a) ((s:.rb):.b)+ )+ => Transpose ((s:.ra):.a) ((s:.rb):.b)+ where+ transpose = transpose'+ {-# INLINE transpose #-}++++class Transpose' a b | a->b+ where transpose' :: a -> b++instance Transpose' () () where + transpose' = id+ {-# INLINE transpose' #-}++instance + (Transpose' vs vs') => Transpose' ( () :. vs ) vs'+ where+ transpose' (():.vs) = transpose' vs+ {-# INLINE transpose' #-}++instance Transpose' ((x:.()):.()) ((x:.()):.()) where+ transpose' = id++instance + (Head xss_h xss_hh+ ,Map xss_h xss_hh (xss_h:.xss_t) xs'+ ,Tail xss_h xss_ht+ ,Map xss_h xss_ht (xss_h:.xss_t) xss_+ ,Transpose' (xs :. xss_) xss'+ )+ => Transpose' ((x:.xs):.(xss_h:.xss_t)) ((x:.xs'):.xss') + where+ transpose' ((x:.xs):.xss) =+ (x :. (map head xss)) :. (transpose' (xs :. (map tail xss) :: (xs:.xss_)))+ {-# INLINE transpose' #-}++++++class SetDiagonal v m | m -> v, v -> m where+ -- |set the diagonal of an n-by-n matrix to a given n-vector+ setDiagonal :: v -> m -> m++instance (Vec n a v, Vec n r m, SetDiagonal' N0 v m) => SetDiagonal v m where+ setDiagonal v m = setDiagonal' (undefined::N0) v m+ {-# INLINE setDiagonal #-}++class SetDiagonal' n v m where+ setDiagonal' :: n -> v -> m -> m++instance SetDiagonal' n () m where+ setDiagonal' _ _ m = m+ {-# INLINE setDiagonal' #-}++instance + ( SetDiagonal' (Succ n) v m+ , Access n a r+ ) => SetDiagonal' n (a:.v) (r:.m) + where+ setDiagonal' _ (a:.v) (r:.m) = + (set (undefined::n) a r) :. (setDiagonal' (undefined::Succ n) v m)+ {-# INLINE setDiagonal' #-}++++class GetDiagonal m v | m -> v, v -> m where+ -- |get the diagonal of an n-by-n matrix as a vector+ getDiagonal :: m -> v++instance (Vec n a v, Vec n v m, GetDiagonal' N0 () m v) => GetDiagonal m v where+ getDiagonal m = getDiagonal' (undefined::N0) () m+ {-# INLINE getDiagonal #-}++class GetDiagonal' n p m v where+ getDiagonal' :: n -> p -> m -> v++instance + (Access n a r+ ,Append p (a:.()) (a:.p)+ ) => GetDiagonal' n p (r:.()) (a:.p) + where+ getDiagonal' _ p (r:.()) = append p ((get (undefined::n) r) :. ())+ {-# INLINE getDiagonal' #-}++instance + (Access n a r+ ,Append p (a:.()) p'+ ,GetDiagonal' (Succ n) p' (r:.m) v+ ) + => GetDiagonal' n p (r:.r:.m) v+ where+ getDiagonal' _ p (r:.m) = + getDiagonal' (undefined::Succ n) (append p ((get (undefined::n) r):.())) m+ {-# INLINE getDiagonal' #-}+++-- | @scale v m@ multiplies the diagonal of matrix @m@ by the vector @s@, component-wise. So+-- @scale 5 m@ multiplies the diagonal by 5, whereas @scale 2:.1 m@+-- only scales the x component.+scale :: + ( GetDiagonal' N0 () m r+ , Num r+ , Vec n a r+ , Vec n r m+ , SetDiagonal' N0 r m+ ) => r -> m -> m+scale s m = setDiagonal (s * (getDiagonal m)) m+{-# INLINE scale #-}+++-- | @diagonal v@ is a square matrix with the vector v as the diagonal, and 0+-- elsewhere.+diagonal :: (Vec n a v, Vec n v m, SetDiagonal v m, Num m) => v -> m+diagonal v = setDiagonal v 0+{-# INLINE diagonal #-}+++-- | identity matrix (square)+identity :: (Vec n a v, Vec n v m, Num v, Num m, SetDiagonal v m) => m+identity = diagonal 1 +{-# INLINE identity #-}+++-- DropConsec: this is a helper function for computing determinants. Given an+-- n-vector v, drop each element from v and collect the remaning (n-1)-vectors+-- into an n-vector (ie an n-by-(n-1) matrix)+class DropConsec v vv | v -> vv where+ dropConsec :: v -> vv++instance + (Vec n a v+ ,Pred n n_+ ,Vec n_ a v_+ ,Vec n v_ vv+ ,DropConsec' () v vv+ ) => DropConsec v vv+ where+ dropConsec v = dropConsec' () v + {-# INLINE dropConsec #-}++class DropConsec' p v vv where+ dropConsec' :: p -> v -> vv+ +instance DropConsec' p (a:.()) (p:.()) where+ dropConsec' p (a:.()) = (p:.())+ {-# INLINE dropConsec' #-}++instance + (Append p (a:.v) x+ ,Append p (a:.()) y+ ,DropConsec' y (a:.v) z+ ) + => DropConsec' p (a:.a:.v) (x:.z)+ where+ dropConsec' p (a:.v) = + (append p v) :. (dropConsec' (append p (a:.())) v)+ {-# INLINE dropConsec' #-}++++--Alternating: vector of alternating positive/negative values. This is also a+--helper for computing determinants+class Alternating n a v | v -> n a where+ alternating :: n -> a -> v++instance Alternating N1 a (a:.()) where+ alternating _ !a = a:.()+ {-# INLINE alternating #-}++instance (Num a, Alternating n a (a:.v)) => Alternating (Succ n) a (a:.a:.v) where+ alternating _ !a = a:.(alternating (undefined::n) (negate $! a))+ {-# INLINE alternating #-}+++-- The Determinant of a square matrix, by minor expansion. +class Det' a m | m -> a where+ det' :: m -> a++instance Num a => Det' a ((a:.a:.()):.(a:.a:.()):.()) where+ det' ( (a:.b:.()) :. (c:.d:.()) :. () ) = a*d-b*c+ {-# INLINE det' #-}++--this instance is particularly ugly in order to avoid overlapping with the one above+instance+ (Num a+ ,Num (a:.a:.a:.v)+ ,Fold a (a:.a:.a:.v)+ ,Alternating (Succ (Succ (Succ n))) a (a:.a:.a:.v)+ ,DropConsec (a:.a:.a:.v) vv+ ,Map (a:.a:.a:.v) vv ((a:.a:.a:.v):.(a:.a:.a:.v):.m) vmt+ ,Transpose vmt vm+ ,Map ((a:.a:.v):.(a:.a:.v):.m_) a vm (a:.a:.a:.v)+ ,Det' a ((a:.a:.v):.(a:.a:.v):.m_)+ ,Vec (Succ (Succ (Succ n))) a (a:.a:.a:.v)+ ,Vec (Succ (Succ (Succ n))) (a:.a:.a:.v) ((a:.a:.a:.v):.(a:.a:.a:.v):.(a:.a:.a:.v):.m)+ )+ => + Det' a ((a:.a:.a:.v):.(a:.a:.a:.v):.(a:.a:.a:.v):.m)+ where+ det' (mh:.mt) =+ sum ((alternating undefined 1) * mh *+ (map det' (transpose (map dropConsec mt :: vmt))))+ {-# INLINE det' #-}+++-- For now, use wrapper class to allow type inference. I think maybe the+-- squareness of the matrix is keeping Det' from inferring properly, so we'll+-- enforce that here. But really I have no clue.+++class Det n a m | m -> a where+ -- | Determinant by minor expansion. Unfolds into a closed form expression.+ -- This should be the fastest way for 4x4 and smaller, but @snd . gaussElim@+ -- works too.+ det :: m -> a++instance (Vec n a r, Vec n r m, Det' a m) => Det n a m where+ det = det'+ {-# INLINE det #-}++++--ReplConsec : this is a helper for implementing Cramer's rule. Given an+--n-vector v and a value r, replace each consecutive element from v with r,+--and collect the resulting n-vectors into an n-vector (ie an n-by-n matrix)++class ReplConsec a v vv | v->a, v->vv, vv->v, vv->a where+ replConsec :: a -> v -> vv++instance + (Vec n a v+ ,Vec n v vv+ ,ReplConsec' a () v vv+ ) => ReplConsec a v vv+ where+ replConsec a v = replConsec' a () v :: vv+ {-# INLINE replConsec #-}++class ReplConsec' a p v vv where+ replConsec' :: a -> p -> v -> vv++instance ReplConsec' a p () () where+ replConsec' _ _ () = ()+ {-# INLINE replConsec' #-}++instance + (Append p (a:.v) x+ ,Append p (a:.()) y+ ,ReplConsec' a y v z+ ) + => ReplConsec' a p (a:.v) (x:.z)+ where+ replConsec' r p (a:.v) = + (append p (r:.v)) :. (replConsec' r (append p (a :. ())) v)+ {-# INLINE replConsec' #-}+++++-- | @cramer'sRule m v@ computes the solution to @m\`multmv\`x=v@ using the+-- eponymous method. For larger than 3x3 you will want to use 'solve', which+-- uses 'gaussElim'. Cramer's rule, however, unfolds into a closed-form+-- expression, with no branches or allocations (other than the result). You may+-- need to increase the unfolding threshold to see this.++cramer'sRule :: + (Map a a1 b1 v+ ,Transpose w b1+ ,ZipWith a2 b vv v m w+ ,ReplConsec' a2 () b vv+ ,Vec n b vv+ ,Vec n a2 b+ ,Fractional a1+ ,Det' a1 m+ ,Det' a1 a+ ) => m -> v -> v+cramer'sRule m b =+ case map (\m' -> (det' m')/(det' m)) + (transpose (zipWith replConsec b m)) + of b' -> b' `asTypeOf` b +{-# INLINE cramer'sRule #-}+++++++mapFst f (a,b) = (f a,b)+{-# INLINE mapFst #-}+++class Num a => NearZero a where+ -- | @nearZero x@ should be true when x is close enough to 0 to cause+ -- significant error in division. + nearZero :: a -> Bool+ nearZero 0 = True+ nearZero _ = False+ {-# INLINE nearZero #-}++instance NearZero Float where+ nearZero x = abs x < 1e-6+ {-# INLINE nearZero #-}++instance NearZero Double where+ nearZero x = abs x < 1e-14+ {-# INLINE nearZero #-}++instance NearZero Rational+++++-- Pivot1 : find a non-zero pivot column and put a 1 there. Second return+-- argument tracks value of determinant. Returns nothing if no pivot in the+-- first row. Does not try to find the 'best' pivot, only an acceptable one:+-- matrices are assumed small, roundoff error should be negligible. ++class Pivot1 a m | m -> a where+ pivot1 :: m -> Maybe (m,a)++instance Pivot1 a () where+ pivot1 _ = Nothing++instance + ( Fractional a, NearZero a+ ) => Pivot1 a ((a:.()):.()) + where+ pivot1 ((p:._):._) + | nearZero p = Nothing+ | otherwise = Just (1,p)+ {-# INLINE pivot1 #-}++instance + ( Fractional a, NearZero a + , Map a a (a:.r) (a:.r)+ ) => Pivot1 a ((a:.(a:.r)):.()) + where+ pivot1 ((p:.r):._) + | nearZero p = Nothing+ | otherwise = Just ((1 :. (map (/p) r)):.(), p)+ {-# INLINE pivot1 #-}++instance + ( Fractional a, NearZero a+ , Map a a (a:.r) (a:.r)+ , ZipWith a a a (a:.r) (a:.r) (a:.r) + , Map (a:.r) (a:.r) ((a:.r):.rs) ((a:.r):.rs)+ , Pivot1 a ((a:.r):.rs) + ) => Pivot1 a ((a:.r):.(a:.r):.rs) + where+ pivot1 (row@(p:._):.rows) + | nearZero p = pivot1 rows >>= \(r:.rs,p)-> Just(r:.row:.rs,p)+ | otherwise = Just ( first:.(map add rows) , p)+ where first = map (/p) row+ add r@(x:._) = zipWith (-) r . map (*x) $ first + {-# INLINE pivot1 #-}+++-- Pivot : find a pivot. Second return argument tracks determinant.+-- Returns Nothing if no pivot anywhere.++class Pivot a m | m -> a where+ pivot :: m -> Maybe (m,a)++instance Pivot a (():.v) where+ pivot _ = Nothing+ {-# INLINE pivot #-}++instance + ( Fractional a+ , NearZero a+ , Pivot1 a rs + , Tail (a:.r) r+ , Map (a:.r) r ((a:.r):.rs) (r:.rs') + , Map r (a:.r) (r:.rs') ((a:.r):.rs)+ , Pivot1 a ((a:.r):.rs)+ , Pivot a (r:.rs')+ ) => Pivot a ((a:.r):.rs) + where+ pivot m = + mplus (pivot1 m) + (pivot (map tail m) >>= return . mapFst (map (0:.)) )+ {-# INLINE pivot #-}++++-- | Gaussian elimination, adapted from Mirko Rahn:+-- <http://www.haskell.org/pipermail/glasgow-haskell-users/2007-May/012648.html>+--+-- This is more of a proof of concept. Using a foreign C function will run+-- slightly faster, and compile much faster. But where is the fun in that?+-- Set your unfolding threshold as high as possible.++class GaussElim a m | m -> a where+ -- | @gaussElim m@ returns a pair @(m',d)@ where @m'@ is @m@ in row echelon+ -- form and @d@ is the determinant of @m@. The determinant of @m'@ is 1 or 0,+ -- i.e., the leading coefficient of each non-zero row is 1. + + gaussElim :: m -> (m,a)++instance (Num a, Pivot a (r:.())) => GaussElim a (r:.())+ where+ gaussElim m = fromMaybe (m,1) (pivot m) + {-# INLINE gaussElim #-}++instance + ( Fractional a+ , Map (a:.r) r ((a:.r):.rs) rs_+ , Map r (a:.r) rs_ ((a:.r):.rs) + , Pivot a ((a:.r):.(a:.r):.rs)+ , GaussElim a rs_+ ) => GaussElim a ((a:.r):.(a:.r):.rs)+ where+ gaussElim m =+ flip (maybe (m,1)) (pivot m) $ \(row:.rows,p) ->+ case gaussElim (map tail rows)+ of (rows',p') -> ( row:.(map (0:.) rows') , p*p')+ {-# INLINE gaussElim #-}++++class BackSubstitute m where+ -- | backSubstitute takes a full rank matrix from row echelon form to reduced+ -- row echelon form. Returns @Nothing@ if the matrix is rank deficient. + backSubstitute :: m -> Maybe m ++instance BackSubstitute ((a:.r):.()) where+ backSubstitute = Just . id+ {-# INLINE backSubstitute #-}++instance + ( Map (a:.r) r ((a:.r):.rs) rs_ --map tail+ , Map r (a:.r) rs_ ((a:.r):.rs) --map cons+ , Fold (a,a:.r) aas+ , ZipWith a a a (a:.r) (a:.r) (a:.r)+ , Map a a (a:.r) (a:.r)+ , ZipWith a (a:.r) (a,a:.r) r ((a:.r):.rs) aas+ , Num a, NearZero a+ , BackSubstitute rs_+ ) => BackSubstitute ((a:.r):.(a:.r):.rs)+ where+ backSubstitute (r@(rh:.rt):.rs) + | nearZero (1-rh) = + liftM (map (0:.)) (backSubstitute . map tail $ rs) >>= \rs' -> + return . (:.rs') . foldl (\v (a,w) -> sub v a w) r $ + zipWith (,) rt rs'+ | otherwise = Nothing -- rank deficient+ where sub v a = zipWith (-) v . map (*a)+ {-# INLINE backSubstitute #-}+++++class BackSubstitute' m where+ -- | backSubstitute' takes a full rank matrix from row echelon form to reduced+ -- row echelon form. Returns garbage is matrix is rank deficient.+ backSubstitute' :: m -> m ++instance BackSubstitute' ((a:.r):.()) where+ backSubstitute' = id+ {-# INLINE backSubstitute' #-}++instance + ( Map (a:.r) r ((a:.r):.rs) rs_ --map tail+ , Map r (a:.r) rs_ ((a:.r):.rs) --map cons+ , Fold (a,a:.r) aas+ , ZipWith a a a (a:.r) (a:.r) (a:.r)+ , Map a a (a:.r) (a:.r)+ , ZipWith a (a:.r) (a,a:.r) r ((a:.r):.rs) aas+ , Num a+ , BackSubstitute' rs_+ ) => BackSubstitute' ((a:.r):.(a:.r):.rs)+ where+ backSubstitute' (r@(_:.rt):.rs) = + case map (0:.) (backSubstitute' . map tail $ rs) + of rs' -> (:.rs') $ foldl (\ v (a,w) -> sub v a w) r + (zipWith (,) rt rs')+ where sub v a = zipWith (-) v . map (*a)+ {-# INLINE backSubstitute' #-}+++-- | @invert m@ returns @Just@ the inverse of @m@ or @Nothing@ if @m@ is singular.+invert :: forall n a r m r' m'. + ( Num r, Num m+ , Vec n a r -- r is row type+ , Vec n r m -- m is matrix type+ , Append r r r' -- r' is a row of augmented matrix+ , ZipWith r r r' m m m' -- m' is the augmented matrix+ , Drop n r' r -- get the right half of an augmented matrix row+ , Map r' r m' m -- get the right half of the augmented matrix+ , SetDiagonal r m -- needed to make identity matrix+ , GaussElim a m'+ , BackSubstitute m'+ ) => m -> Maybe m+invert m = + return i >>= backSubstitute . fst . gaussElim . zipWith append m + >>= return . map dropn+ where dropn = drop (undefined::n)+ i = identity :: m+{-# INLINE invert #-}++-- | inverse and determinant. If det = 0, inverted matrix is garbage.+invertAndDet :: forall n a r m r' m'. + ( Num r, Num m+ , Vec n a r -- r is row type+ , Vec n r m -- m is matrix type+ , Append r r r' -- r' is a row of augmented matrix+ , ZipWith r r r' m m m' -- m' is the augmented matrix+ , Drop n r' r -- get the right half of an augmented matrix row+ , Map r' r m' m -- get the right half of the augmented matrix+ , SetDiagonal r m -- needed to make identity matrix+ , GaussElim a m'+ , BackSubstitute' m'+ ) => m -> (m,a)+invertAndDet m = + mapFst ( (map dropn) . backSubstitute') . gaussElim . zipWith append m $ i+ where dropn = drop (undefined::n)+ i = identity :: m+{-# INLINE invertAndDet #-}+++-- | Solution of linear system by Gaussian elimination. Returns @Nothing@+-- if no solution. +solve :: forall n a v r m r' m'. + ( Num r, Num m+ , Vec n a r -- r is row type+ , Vec n r m -- m is matrix type+ , Snoc r a r' -- a row of the extended matrix is one longer+ , ZipWith r a r' m r m' -- m' is the augmented matrix+ , Drop n r' (a:.()) -- get the right part of an augmented matrix row+ , Map r' a m' r -- get the right part of the augmented matrix+ , GaussElim a m'+ , BackSubstitute m'+ ) => m -> r -> Maybe r+solve m v = + return v >>= backSubstitute . fst . gaussElim . zipWith snoc m + >>= return . map (head . drop (undefined::n)) +{-# INLINE solve #-}+
+ Data/Vec/Nat.hs view
@@ -0,0 +1,69 @@+{- Copyright (c) 2008, Scott E. Dillard. All rights reserved. -}++{-# LANGUAGE EmptyDataDecls #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UndecidableInstances #-}++-- | Type level naturals. @Ni@ is a type, @ni@ an undefined value of that type,+-- for @i <- [0..19]@+module Data.Vec.Nat where+++data N0+data Succ a++type N1 = Succ N0+type N2 = Succ N1+type N3 = Succ N2+type N4 = Succ N3+type N5 = Succ N4+type N6 = Succ N5+type N7 = Succ N6+type N8 = Succ N7+type N9 = Succ N8+type N10 = Succ N9+type N11 = Succ N10+type N12 = Succ N11+type N13 = Succ N12+type N14 = Succ N13+type N15 = Succ N14+type N16 = Succ N15+type N17 = Succ N16+type N18 = Succ N17+type N19 = Succ N18++n0 :: N0 ; n0 = undefined+n1 :: N1 ; n1 = undefined+n2 :: N2 ; n2 = undefined+n3 :: N3 ; n3 = undefined+n4 :: N4 ; n4 = undefined+n5 :: N5 ; n5 = undefined+n6 :: N6 ; n6 = undefined+n7 :: N7 ; n7 = undefined+n8 :: N8 ; n8 = undefined+n9 :: N9 ; n9 = undefined+n10 :: N10 ; n10 = undefined+n11 :: N11 ; n11 = undefined+n12 :: N12 ; n12 = undefined+n13 :: N13 ; n13 = undefined+n14 :: N14 ; n14 = undefined+n15 :: N15 ; n15 = undefined+n16 :: N16 ; n16 = undefined+n17 :: N17 ; n17 = undefined+n18 :: N18 ; n18 = undefined+n19 :: N19 ; n19 = undefined++-- | @nat n@ yields the @Int@ value of the type-level natural @n@.+class Nat n where nat :: n -> Int+instance Nat N0 where nat _ = 0+instance Nat a => Nat (Succ a) where nat _ = 1+(nat (undefined::a))++class Pred x y | x -> y, y -> x+instance Pred (Succ N0) N0+instance Pred (Succ n) p => Pred (Succ (Succ n)) (Succ p)+
+ Data/Vec/Packed.hs view
@@ -0,0 +1,152 @@+{- Copyright (c) 2008, Scott E. Dillard. All rights reserved. -}++{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE UndecidableInstances #-}++-- | Packed vectors : use these whenever possible. The regular vector type is+-- just a gussied up linked list, but when vector functions are applied to+-- these types, bracketed by @'pack'@ and @'unpack'@, then things unfold into+-- perfectly optimized code.++module Data.Vec.Packed where++import Data.Vec.Base as V++-- * Packed Vector Types+data Vec2I = Vec2I {-#UNPACK#-} !Int + {-#UNPACK#-} !Int ++data Vec3I = Vec3I {-#UNPACK#-} !Int + {-#UNPACK#-} !Int + {-#UNPACK#-} !Int++data Vec4I = Vec4I {-#UNPACK#-} !Int + {-#UNPACK#-} !Int + {-#UNPACK#-} !Int+ {-#UNPACK#-} !Int++data Vec2F = Vec2F {-#UNPACK#-} !Float + {-#UNPACK#-} !Float ++data Vec3F = Vec3F {-#UNPACK#-} !Float + {-#UNPACK#-} !Float + {-#UNPACK#-} !Float++data Vec4F = Vec4F {-#UNPACK#-} !Float + {-#UNPACK#-} !Float + {-#UNPACK#-} !Float+ {-#UNPACK#-} !Float++data Vec2D = Vec2D {-#UNPACK#-} !Double + {-#UNPACK#-} !Double ++data Vec3D = Vec3D {-#UNPACK#-} !Double + {-#UNPACK#-} !Double + {-#UNPACK#-} !Double++data Vec4D = Vec4D {-#UNPACK#-} !Double + {-#UNPACK#-} !Double + {-#UNPACK#-} !Double+ {-#UNPACK#-} !Double++-- * Packed Matrix Types. +type Mat22I = Vec2 Vec2I +type Mat23I = Vec2 Vec3I +type Mat33I = Vec3 Vec3I +type Mat34I = Vec3 Vec4I +type Mat44I = Vec4 Vec3I ++type Mat22F = Vec2 Vec2F +type Mat23F = Vec2 Vec3F +type Mat33F = Vec3 Vec3F +type Mat34F = Vec3 Vec4F +type Mat44F = Vec4 Vec3F ++type Mat22D = Vec2 Vec2D +type Mat23D = Vec2 Vec3D +type Mat33D = Vec3 Vec3D +type Mat34D = Vec3 Vec4D +type Mat44D = Vec4 Vec4D +++-- | pack a matrix+packMat :: (Map v pv m pm, PackedVec pv v) => m -> pm+packMat = V.map pack +{-# INLINE packMat #-}++-- | unpack a matrix+unpackMat :: (Map pv v pm m, PackedVec pv v) => pm -> m+unpackMat = V.map unpack+{-# INLINE unpackMat #-}++-- | PackedVec class : relates a packed vector type to its unpacked type For+-- now, the fundep is not bijective -- It may be advantageous to have multiple+-- packed representations for a canonical vector type. This may change. In the+-- meantime, you may have to annotate return types.+class PackedVec pv v | pv -> v where+ pack :: v -> pv+ unpack :: pv -> v++instance PackedVec Vec2I (Vec2 Int) where+ pack (x:.y:.()) = Vec2I x y + unpack (Vec2I x y) = x:.y:.()+ {-# INLINE pack #-}+ {-# INLINE unpack #-}++instance PackedVec Vec3I (Vec3 Int) where+ pack (x:.y:.z:.()) = Vec3I x y z+ unpack (Vec3I x y z) = x:.y:.z:.()+ {-# INLINE pack #-}+ {-# INLINE unpack #-}++instance PackedVec Vec4I (Vec4 Int) where+ pack (x:.y:.z:.w:.()) = Vec4I x y z w+ unpack (Vec4I x y z w) = x:.y:.z:.w:.()+ {-# INLINE pack #-}+ {-# INLINE unpack #-}+++instance PackedVec Vec2F (Vec2 Float) where+ pack (x:.y:.()) = Vec2F x y + unpack (Vec2F x y) = x:.y:.()+ {-# INLINE pack #-}+ {-# INLINE unpack #-}++instance PackedVec Vec3F (Vec3 Float) where+ pack (x:.y:.z:.()) = Vec3F x y z+ unpack (Vec3F x y z) = x:.y:.z:.()+ {-# INLINE pack #-}+ {-# INLINE unpack #-}++instance PackedVec Vec4F (Vec4 Float) where+ pack (x:.y:.z:.w:.()) = Vec4F x y z w+ unpack (Vec4F x y z w) = x:.y:.z:.w:.()+ {-# INLINE pack #-}+ {-# INLINE unpack #-}+++instance PackedVec Vec2D (Vec2 Double) where+ pack (x:.y:.()) = Vec2D x y + unpack (Vec2D x y) = x:.y:.()+ {-# INLINE pack #-}+ {-# INLINE unpack #-}++instance PackedVec Vec3D (Vec3 Double) where+ pack (x:.y:.z:.()) = Vec3D x y z+ unpack (Vec3D x y z) = x:.y:.z:.()+ {-# INLINE pack #-}+ {-# INLINE unpack #-}++instance PackedVec Vec4D (Vec4 Double) where+ pack (x:.y:.z:.w:.()) = Vec4D x y z w+ unpack (Vec4D x y z w) = x:.y:.z:.w:.()+ {-# INLINE pack #-}+ {-# INLINE unpack #-}+++
+ LICENSE view
@@ -0,0 +1,22 @@+Copyright (c) 2008, Scott E. Dillard+All rights reserved.++Redistribution and use in source and binary forms, with or without modification,+are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright notice,+ this list of conditions and the following disclaimer.+ * Redistributions in binary form must reproduce the above copyright notice,+ this list of conditions and the following disclaimer in the documentation+ and/or other materials provided with the distribution.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON+ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,4 @@+#! /usr/bin/env runhaskell++import Distribution.Simple+main = defaultMain
+ Vec.cabal view
@@ -0,0 +1,42 @@+Name: Vec+Version: 0.9.0+License: BSD3+License-file: LICENSE+Author: Scott E. Dillard+Maintainer: Scott E. Dillard <sedillard@gmail.com>+Stability: Experimental+Synopsis: Fixed-length lists and low-dimensional linear algebra.+Description: + Vectors are represented by lists with type-encoded lengths. The constructor+ is @:.@, which acts like a cons both at the value and type levels, with @()@+ taking the place of nil. So @x:.y:.z:.()@ is a 3d vector. The library+ provides a set of common list-like functions (map, fold, etc) for working+ with vectors. Built up from these functions are a small but useful set of+ linear algebra operations: matrix multiplication, determinants, solving+ linear systems, inverting matrices.+Cabal-version: >=1.2+Build-type: Simple+Category: Math,Graphics++library+ Build-Depends: base++ Exposed-modules: Data.Vec + Data.Vec.Base,+ Data.Vec.LinAlg,+ Data.Vec.Nat,+ Data.Vec.Instances+ Data.Vec.Packed+ Extensions: + BangPatterns,+ EmptyDataDecls,+ ExistentialQuantification,+ FlexibleInstances, + FlexibleContexts,+ FunctionalDependencies,+ MultiParamTypeClasses, + NoMonomorphismRestriction,+ ScopedTypeVariables,+ TypeOperators, + TypeSynonymInstances,+ UndecidableInstances