Data.Vec.LinAlg

Documentation

dot :: (Num a, Num v, Fold v a) => v -> v -> a

dot / inner / scalar product

normSq :: (Num a, Num v, Fold v a) => v -> a

vector norm, squared

norm :: (Num v, Floating a, Fold v a) => v -> a

vector / L2 / Euclidean norm

normalize :: (Floating a, Num v, Fold v a, Map a a v v) => v -> v

normalize v is a unit vector in the direction of v. v is assumed non-null.

cross :: Num a => Vec3 a -> Vec3 a -> Vec3 a

3d cross product.

homPoint :: (Snoc v a v', Num a) => v -> v'

lift a point into homogenous coordinates

homVec :: (Snoc v a v', Num a) => v -> v'

point-at-infinity in homogenous coordinates

project :: (Reverse' () t1 v', Fractional t1, Vec a t t1, Reverse' () v (t :. t1)) => v -> v'

project a vector from homogenous coordinates. Last vector element is assumed non-zero.

multvm :: (Transpose m mt, Map v a mt v', Fold v a, Num a, Num v) => v -> m -> v'

row vector * matrix

multmv :: (Map v a m v', Num v, Fold v a, Num a) => m -> v -> v'

matrix * column vector

multmm :: (Map v v' m1 m3, Map v a b v', Transpose m2 b, Fold v a, Num v, Num a) => m1 -> m2 -> m3

matrix * 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

apply a translation to a projective transformation matrix

column :: (Transpose m mt, Access n v mt) => n -> m -> v

get the n-th column as a vector. n is a type-level natural.

row :: Access n a v => n -> v -> a

get the n-th row as a vector. n is a type-level natural.

class Transpose a b | a -> b, b -> a where

matrix transposition

Methods

transpose :: a -> b

Instances

Transpose () ()

(Vec (Succ n) s (:. s ra), Vec (Succ m) (:. s ra) (:. (:. s ra) a), Vec (Succ m) s (:. s rb), Vec (Succ n) (:. s rb) (:. (:. s rb) b), Transpose' (:. (:. s ra) a) (:. (:. s rb) b)) => Transpose (:. (:. s ra) a) (:. (:. s rb) b)

scale :: (GetDiagonal' N0 () m r, Num r, Vec n a r, Vec n r m, SetDiagonal' N0 r m) => r -> m -> m

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.

diagonal :: (Vec n a v, Vec n v m, SetDiagonal v m, Num m) => v -> m

diagonal v is a square matrix with the vector v as the diagonal, and 0 elsewhere.

identity :: (Vec n a v, Vec n v m, Num v, Num m, SetDiagonal v m) => m

identity matrix (square)

det :: forall n a r m. (Vec n a r, Vec n r m, Det' m a) => m -> a

Determinant by minor expansion, i.e. Laplace's formula. Unfolds into a closed form expression. This should be the fastest way for 4x4 and smaller, but snd . gaussElim works too.

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' m a1, Det' a a1) => m -> v -> v

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.

class GaussElim a m | m -> a where

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.

Methods

gaussElim :: m -> (m, a)

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.

Instances

(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))

(Num a, Pivot a (:. r ())) => GaussElim a (:. r ())

invert :: forall n a r m r' m'. (Num r, Num m, Vec n a r, Vec n r m, Append r r r', ZipWith r r r' m m m', Drop n r' r, Map r' r m' m, SetDiagonal r m, GaussElim a m', BackSubstitute m') => m -> Maybe m

invert m returns Just the inverse of m or Nothing if m is singular.

invertAndDet :: forall n a r m r' m'. (Num a, Num r, Num m, Vec n a r, Vec n r m, Append r r r', ZipWith r r r' m m m', Drop n r' r, Map r' r m' m, SetDiagonal r m, GaussElim a m', BackSubstitute m') => m -> (m, a)

inverse and determinant. If det = 0, inverted matrix is garbage.

solve :: forall n a v r m r' m'. (Num r, Num m, Vec n a r, Vec n r m, Snoc r a r', ZipWith r a r' m r m', Drop n r' (a :. ()), Map r' a m' r, GaussElim a m', BackSubstitute m') => m -> r -> Maybe r

Solution of linear system by Gaussian elimination. Returns Nothing if no solution.