On my last
graphics assignment, I spent far more effort on some non-requirements than on requirements. Among those was a
Vector class, for dealing with vectors and matrices of any size. Much of it went unused, but rather than let it just go to waste, I figure I’ll share it instead. It’s hypothetically possible that it could be helpful to
somebody, so here it is, cleaned up a bit and with more documentation.
(This post is in
Literate Haskell format. This means that, in addition to looking pretty, you can copy’n’paste the whole thing right into
Vector.lhs and any Haskell compiler will be able to read it.)
I use fancy powers like multi-parameter classes here, so let’s inform the compiler that we’re using extensions beyond the
Haskell 98 specification. (Think
#pragma.)
> {-# OPTIONS -fallow-incoherent-instances -fallow-undecidable-instances #-}> {-# OPTIONS -fglasgow-exts #-}Legalese blah.
> -- Copyright (c) 2007, 林子傑> -- All rights reserved.> --> -- Redistribution and use in source and binary forms, with or without> -- modification, are permitted provided that the following conditions are met:> --> -- 1. Redistributions of source code must retain the above copyright notice,> -- this list of conditions and the following disclaimer.> -- 2. 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.Alright, let’s get to it!
Vector.lhs should contain a module named
Vector. It’s always good style to have an explicit export list. I’m leaving one out here just because I’m lazy… please don’t imitate me.
> module Vector where> > import Control.Arrow> import Control.Monad> import Data.Ix> import qualified Data.List as List> import qualified Prelude> import Prelude hiding> ( foldl, foldr, head, last, length, map, negate, scanl, scanr> , reverse, unzip, zip, zip3, zipWith, zipWith3 )> First up, here are some tools we will use.
> -- types that may be rearranged> class Transpose a b where transpose :: a (b t) -> b (a t)> > -- groups, like in number theory> class Group a where> infixl 6 .+, .-; (.+), (.-) :: a -> a -> a; (.-) = (negate .) . (.+)> negate :: a -> a; negate = (.-) identity; identity :: a> instance Num a => Group a where> (.+) = (+); (.-) = (-); negate = Prelude.negate; identity = 0> > -- types with a "multiplication" defined between them> class Multiply a b c | a b -> c where infixl 7 .*; (.*) :: a -> b -> c> Oh my… no automatic transitivity. Why does this overflow the typechecker’s context reduction stack? Because by giving the compiler a way to chain
Casts, it can generate an infinite number of
casts between any two
Castable types
a and
b (
a→b,
a→a→b,
a→b→b,
a→b→a→b, …).
> -- safely coerceable types> class Cast u v where cast :: u -> v> instance Cast self self where cast = id> -- oops, this overflows the typechecker's context reduction stack> --instance (Cast a b, Cast b c) => Cast a c where cast = cast . cast> Does
0/0 = 0? Not always, but later on, we’ll use this to let the
0 Vector normalize to
0 instead of
NaN.
> -- division lazy in the second argument> infixl 7 /?; (/?) :: Fractional a => a -> a -> a; 0 /? _ = 0; x /? y = x/y> The ternary operator is an evil, evil thing. Never forget that. Haskell doesn’t have one. …I made one.
> -- define a ternary operator> infixl 1 ?; (?) :: Bool -> a -> a -> a; (?) True = const; (?) False = flip const> It almost seems like
\n -> unfoldr (Just . splitAt n) is all we need (don’t the types work out surprisingly well?), but unfortunately that results in a never-ending tail of
[]. We need to teach
unfoldr how to stop…
> -- take n elements at a time> groupInto :: Int -> [a] -> [[a]]> groupInto n = List.unfoldr maybeSplitAt> where maybeSplitAt l = case splitAt n l of ([], _) -> Nothing; l2 -> Just l2> Here is our definition of what a
Vector is. Hmm, looks like a
[] list…
> -- a special Functor with nice abilities> class Functor v => Vector v where> zipWith :: (a -> b -> c) -> v a -> v b -> v c> zipWith = (. zip) . (.) . map . uncurry> foldl :: (a -> b -> a) -> a -> v b -> a> foldl f k = Prelude.foldl f k . elems> foldr :: (a -> b -> b) -> b -> v a -> b> foldr f k = Prelude.foldr f k . elems> head :: v a -> a; head = Prelude.head . elems> last :: v a -> a; last = Prelude.last . elems> reverse :: v a -> v a; reverse = new . Prelude.reverse . elems> length :: v a -> Int; length = Prelude.length . elems> elems :: v a -> [a]; elems = foldr (:) []> new :: [a] -> v a> In fact, a normal
[] list indeed
is a
Vector.
> -- generic, arbitrary-length vector> instance Vector v => Transpose v [] where> transpose = Prelude.map new . List.transpose . elems> instance Vector v => Transpose [] v where> transpose = new . List.transpose . Prelude.map elems> instance (Vector u, Vector v) => Transpose u v where> transpose = map new . transpose . elems> instance Transpose [] [] where transpose = List.transpose> instance Vector [] where zipWith = Prelude.zipWith; elems = id; new = id> Let’s say that a
Matrix is a
Vector of numeric
Vectors,
Transposeable.
> newtype (Num a, Num (v1 a), Num (v2 a), Transpose v1 v2, Transpose v2 v1) =>> Matrix v1 v2 a = Matrix (v1 (v2 a)) deriving (Eq, Ord, Show)> By defining all the standard numeric operations on
Vectors, they can be treated just like all other
Numbers. Code-wise, this isn’t very interesting, though. We’re just doing everything component-wise.
> -- allow a Vector to pretend to be a number, acting component-wise> instance (Num a, Vector v, Eq (v a), Show (v a)) => Num (v a) where> (+) = zipWith (+); (*) = zipWith (*); (-) = zipWith (-)> negate = map Prelude.negate; abs = map abs; signum = map signum> fromInteger = new . repeat . fromInteger> instance (Fractional a, Vector v, Num (v a)) => Fractional (v a) where> (/) = zipWith (/); recip = map recip> fromRational = new . repeat . fromRational> instance (Floating a, Vector v, Fractional (v a)) => Floating (v a) where> pi = new $ repeat pi; exp = map exp; sqrt = map sqrt; log = map log> (**) = zipWith (**); logBase = zipWith logBase> sin = map sin; tan = map tan; cos = map cos> asin = map asin; atan = map atan; acos = map acos> sinh = map sinh; tanh = map tanh; cosh = map cosh> asinh = map asinh; atanh = map atanh; acosh = map acosh> instance (Real a, Vector v, Num (v a), Ord (v a)) => Real (v a) where> toRational = toRational . Prelude.head . elems> instance (RealFrac a, Vector v, Real (v a), Fractional (v a)) => RealFrac (v a)> where> properFraction x = (y, x - fromIntegral y) where y = truncate x> truncate = truncate . head; round = round . head> ceiling = ceiling . head; floor = floor . head> instance (Ix a, Vector v, Ord (v a), Num (v Int)) => Ix (v a) where> range (a, b) = map new . foldr accum [[]] $ zipWith (curry range) a b> where accum l = (=<<) $ flip Prelude.map l . flip (:)> index (a, b) c => sum . Prelude.zipWith (*) (elems $ zipWith3 (curry index) a b c) .> scanl (*) 1 $ zipWith (curry rangeSize) a b> inRange (a, b) = and . elems . zipWith inRange (zip a b)> rangeSize = product . elems . uncurry (zipWith (curry rangeSize))> Haskell’s
Num class is like a Ring from number theory, and non-square matrices do not form a ring. Matrices do form a
Group, though. We’ll define that later, but as long as we eventually do, Haskell lets us results
now by saying that
((+), (-), negate) are just like the versions for non-square matrices.
> -- similarly for matrices, act like numbers> -- unfortunately, the numeric classes set up in Haskell's Prelude place addition> -- and multiplication in the same class, so we only handle squares here> instance (Num a, Vector v, Num (v a), Eq (v (v a)), Show (v (v a))) =>> Num (Matrix v v a) where> (+) = (.+); (*) = (.*); (-) = (.-)> negate = negate; abs = cast . det; signum = cast . signum . det> fromInteger = (cast :: a -> Matrix v v a) . fromInteger> instance (Fractional a, Ord a, Vector v, Num (v a), Num (Matrix v v a)) =>> Fractional (Matrix v v a) whereDivision by
Gaussian elimination. It may look ugly, but it’s still short and sweet compared to implementations I’ve seen in other languages.
> Matrix a / Matrix b = Matrix . new . sort' . backward . forward $ zip b a> where> forward = snd . mapAccumL accumF (id, [])> backward = snd . mapAccumR accumB id> sort' = map snd . List.sortBy ((. fst) . compare . fst)> accumF (f, d) p = ((elim n p' . f, n:d), (n, p'))> where> (x, y) = f p; p' = (map (/e) x, map (/e) y)> (n, e) = max' . filter (not . (`elem` d) . fst) $ mkIndex x> max' = List.maximumBy ((. abs . snd) . compare . abs . snd)> accumB f (n, p) = (f . elim n p', (n, y)) where p'@(_, y) = f p> elim n (x, y) (x', y') = (x' - map (*e) x, y' - map (*e) y)> where e = elems x' !! n> fromRational = (cast :: a -> Matrix v v a) . fromRational> instance (Real a, Vector v, Num (v a), Num (Matrix v v a),> Ord (Matrix v v a)) => Real (Matrix v v a) where> toRational = toRational . det> Addition and subtraction of matrices are of course component-wise.
Matrix⋅
Matrix product can be defined in terms of
Vector⋅
Matrix product which can be defined in terms of
Vector⋅
Vector product.
> -- generic matrix addition/subtraction/multiplication (non-squares)> instance (Num a, Vector v1, Vector v2, Num (v1 a), Num (v2 a)) =>> Group (Matrix v1 v2 a) where> Matrix a .+ Matrix b = Matrix $ zipWith (+) a b> Matrix a .- Matrix b = Matrix $ zipWith (-) a b> negate (Matrix a) = Matrix $ map Prelude.negate a> identity = Matrix . new . repeat . new $ repeat 0> instance (Num a, Vector v1, Vector v2, Num (v1 a), Num (v2 a)) =>> Multiply (v1 a) (Matrix v1 v2 a) (v2 a) where> v .* Matrix a = map (v `dot`) $ transpose a> instance (Num a, Vector v1, Vector v2, Num (v1 a), Num (v2 a)) =>> Multiply (Matrix v1 v2 a) (v2 a) (v1 a) where> Matrix a .* v = map (`dot` v) a> instance (Num a, Vector v1, Vector v2, Vector v3, Num (v1 a), Num (v2 a),> Num (v3 a)) => Multiply (Matrix v1 v2 a) (Matrix v2 v3 a) (Matrix v1 v3 a)> where Matrix a .* b = Matrix $ map (.* b) a> This is far from the most efficient way of calculating the determinant, but it was so easy to write! Just expand every single product, determine the sign, and add ’em all up.
> -- matrix determinant> det :: (Num a, Vector v, Num (v a)) => Matrix v v a -> a> det (Matrix a) = sum . map snd . foldr accum [([], 1)] $ map mkIndex a> where> accum l1 l2 = [(n:d, s $ e*m) | (n, e) <- l1, (d, m) <- l2, s <- sign n d]> sign n d = deranges False n d >>= return . ($ id) . (? Prelude.negate)> deranges k x (y:l) = case compare x y of> LT -> deranges k x l; EQ -> mzero; GT -> deranges (not k) x l> deranges k _ _ = return k> Boring, aside from that a plain
Number
casts to a nice diagonal
Matrix.
> -- coersions> instance (Vector u, Vector v) => Cast (u a) (v a) where cast = new . elems> instance (Num a, Vector t, Vector u, Vector v, Vector w, Num (t a), Num (u a),> Num (v a), Num (w a)) => Cast (Matrix t u a) (Matrix v w a) where> cast (Matrix a) = Matrix . new . map (new . elems) $ elems a> instance (Num a, Vector v1, Vector v2, Num (v1 a), Num (v2 a)) =>> Cast a (Matrix v1 v2 a) where> cast x = Matrix . new . map new . iterate (0:) $ x:repeat 0> instance (Num a, Vector v1, Vector v2, Num (v1 a), Num (v2 a)) =>> Cast (Matrix v1 v2 a) [a] where> cast (Matrix a) = (>>= elems) $ elems a> instance (Num a, Vector v1, Vector v2, Num (v1 a), Num (v2 a)) =>> Cast [a] (Matrix v1 v2 a) where> cast l = Matrix $ new vs> where vs@(~(v:_)) = map new $ groupInto (length v) l> This was used in both the division and the determinant algorithm. It’s the equivalent of
zip [0..] on
Vectors.
> -- pair ascending index numbers to each element of a vector> mkIndex :: Vector v => v a -> [(Int, a)]> mkIndex = snd . mapAccumL (\n e -> (succ n, (n, e))) 0> These should be pretty self-explanatory. My
norm/
magnitude terminology is probably wrong though.
> -- dot product> infixl 7 `dot`> dot :: (Num a, Vector v, Num (v a)) => v a -> v a -> a> dot a b = foldr (+) 0 $ a*b> > -- norm> norm :: (Num a, Vector v, Num (v a)) => v a -> a; norm a = a `dot` a> > -- magnitude> mag :: (Floating a, Vector v, Num (v a)) => v a -> a; mag = sqrt . norm> > -- normalize> normalize :: (Floating a, Vector v, Num (v a)) => v a -> v a> normalize a = map (/? mag a) a> Some fixed-length vectors. This is pretty boring, but we want to hard-wire a few of these both for performance reasons and so that we can ensure that we are always working on
Vectors of a specific, known small size when we want to.
> data Vector0 a = Vector0 deriving (Eq, Ord, Show)> data Vector1 a = Vector1 !a deriving (Eq, Ord, Show)> data Vector2 a = Vector2 !a !a deriving (Eq, Ord, Show)> data Vector3 a = Vector3 !a !a !a deriving (Eq, Ord, Show)> data Vector4 a = Vector4 !a !a !a !a deriving (Eq, Ord, Show)> type Matrix0 = Matrix Vector0 Vector0> type Matrix1 = Matrix Vector1 Vector1> type Matrix2 = Matrix Vector2 Vector2> type Matrix3 = Matrix Vector3 Vector3> type Matrix4 = Matrix Vector4 Vector4> > -- vector specialization: 0-D> instance Functor Vector0 where fmap _ _ = Vector0> instance Transpose [] Vector0 where transpose _ = Vector0> instance Transpose Vector0 [] where transpose _ = []> instance Vector Vector0 where> zipWith _ _ _ = Vector0; foldl _ k _ = k; foldr _ k _ = k> reverse _ = Vector0; length _ = 0; elems _ = []; new _ = Vector0> > -- vector specialization: 1-D> instance Functor Vector1 where fmap f (Vector1 x) = Vector1 $ f x> instance Transpose [] Vector1 where> transpose = Vector1 . map (\(Vector1 x) -> x)> instance Transpose Vector1 [] where transpose (Vector1 xs) = map Vector1 xs> instance Vector Vector1 where> zipWith f (Vector1 x0) (Vector1 x1) = Vector1 $ f x0 x1> foldl f k (Vector1 x) = f k x; foldr f k (Vector1 x) = f x k> head (Vector1 x) = x; last (Vector1 x) = x; reverse = id; length _ = 1> elems (Vector1 x) = [x]> new (x:_) = Vector1 x; new _ = error "Vector.new: short list"> > -- vector specialization: 2-D> instance Functor Vector2 where fmap f (Vector2 x y) = Vector2 (f x) (f y)> instance Transpose [] Vector2 where> transpose = uncurry Vector2 . Prelude.unzip . map (\(Vector2 x y) -> (x, y))> instance Transpose Vector2 [] where> transpose (Vector2 xs ys) = Prelude.zipWith Vector2 xs ys> instance Vector Vector2 where> zipWith f (Vector2 x0 y0) (Vector2 x1 y1) = Vector2 (f x0 x1) (f y0 y1)> foldl f k (Vector2 x y) = k `f` x `f` y> foldr f k (Vector2 x y) = f x $ f y k> head (Vector2 x _) = x; last (Vector2 _ y) = y> reverse (Vector2 x y) = Vector2 y x; length _ = 2> elems (Vector2 x y) = [x, y]> new (x:y:_) = Vector2 x y; new _ = error "Vector.new: short list"> > -- vector specialization: 3-D> instance Functor Vector3 where> fmap f (Vector3 x y z) = Vector3 (f x) (f y) (f z)> instance Transpose [] Vector3 where> transpose as = Vector3 xs ys zs> where (xs, ys, zs) = unzip3 $ map (\(Vector3 x y z) -> (x, y, z)) as> instance Transpose Vector3 [] where> transpose (Vector3 xs ys zs) = Prelude.zipWith3 Vector3 xs ys zs> instance Vector Vector3 where> zipWith f (Vector3 x0 y0 z0) (Vector3 x1 y1 z1) => Vector3 (f x0 x1) (f y0 y1) (f z0 z1)> foldl f k (Vector3 x y z) = k `f` x `f` y `f` z> foldr f k (Vector3 x y z) = f x . f y $ f z k> head (Vector3 x _ _) = x; last (Vector3 _ _ z) = z> reverse (Vector3 x y z) = Vector3 z y x; length _ = 3> elems (Vector3 x y z) = [x, y, z]> new (x:y:z:_) = Vector3 x y z; new _ = error "Vector.new: short list"> > -- vector specialization: 4-D> instance Functor Vector4 where> fmap f (Vector4 x y z t) = Vector4 (f x) (f y) (f z) (f t)> instance Transpose [] Vector4 where> transpose as = Vector4 xs ys zs ts> where (xs, ys, zs, ts) = List.unzip4 $> map (\(Vector4 x y z t) -> (x, y, z, t)) as> instance Transpose Vector4 [] where> transpose (Vector4 xs ys zs ts) = List.zipWith4 Vector4 xs ys zs ts> instance Vector Vector4 where> zipWith f (Vector4 x0 y0 z0 t0) (Vector4 x1 y1 z1 t1) => Vector4 (f x0 x1) (f y0 y1) (f z0 z1) (f t0 t1)> foldl f k (Vector4 x y z t) = k `f` x `f` y `f` z `f` t> foldr f k (Vector4 x y z t) = f x . f y . f z $ f t k> head (Vector4 x _ _ _) = x; last (Vector4 _ _ _ t) = t> reverse (Vector4 x y z t) = Vector4 t z y x; length _ = 4> elems (Vector4 x y z t) = [x, y, z, t]> new (x:y:z:t:_) = Vector4 x y z t; new _ = error "Vector.new: short list"> The only definition of a
Vector cross-product that I’m aware of operates only in 3-D, so this is the only function in this entire module which works only on one specific type of
Vector.
> -- cross product> infixl 7 `cross`> cross :: Num a => Vector3 a -> Vector3 a -> Vector3 a> cross (Vector3 x0 y0 z0) (Vector3 x1 y1 z1) => Vector3 (y0*z1-y1*z0) (z0*x1-z1*x0) (x0*y1-x1*y0)> If you’re me, you’ve got a few bits unhinged when it comes to writing APIs. I re-implement about every useful
[] list function, on
Vectors.
> -- we're used to typing map rather than fmap> map :: Functor v => (a -> b) -> v a -> v b; map = fmap> > -- map+foldl> mapAccumL :: Vector v => (acc -> x -> (acc, y)) -> acc -> v x -> (acc, [y])> mapAccumL f acc = second ($[]) . foldl accum (acc, id)> where accum (a, c) = second ((.) c . (:)) . f a> > -- map+foldr> mapAccumR :: Vector v => (acc -> x -> (acc, y)) -> acc -> v x -> (acc, [y])> mapAccumR f acc = foldr (flip accum) (acc, [])> where accum (b, c) = second (:c) . f b> > -- accumulations of foldl> scanl :: Vector v => (a -> b -> a) -> a -> v b -> [a]> scanl f k = flip snd [] . foldl accum (k, (k:))> where accum (a, c) b = let d = f a b in (d, c . (d:))> > -- accumulations of foldr> scanr :: Vector v => (a -> b -> b) -> b -> v a -> [b]> scanr f = foldr accum . (:[]) where accum a b = f a (Prelude.head b):b> > -- synonym> unzip :: Transpose a b => a (b t) -> b (a t); unzip = transpose> > -- zips> zip :: Vector v => v a -> v b -> v (a, b); zip = zipWith (,)> zip3 :: Vector v => v a -> v b -> v c -> v (a, b, c); zip3 = zipWith3 (,,)> zip4 :: Vector v => v a -> v b -> v c -> v d -> v (a, b, c, d)> zip4 = zipWith4 (,,,)> zip5 :: Vector v => v a -> v b -> v c -> v d -> v e -> v (a, b, c, d, e)> zip5 = zipWith5 (,,,,)> zip6 :: Vector v =>> v a -> v b -> v c -> v d -> v e -> v f -> v (a, b, c, d, e, f)> zip6 = zipWith6 (,,,,,)> zip7 :: Vector v =>> v a -> v b -> v c -> v d -> v e -> v f -> v g -> v (a, b, c, d, e, f, g)> zip7 = zipWith7 (,,,,,,)> > -- extended zipWith> zipWith3 :: Vector v => (a -> b -> c -> d) -> v a -> v b -> v c -> v d> zipWith3 f a = zipWith ($) . zipWith f a> zipWith4 :: Vector v => (a -> b -> c -> d -> e) ->> v a -> v b -> v c -> v d -> v e> zipWith4 f a b = zipWith ($) . zipWith3 f a b> zipWith5 :: Vector v => (a -> b -> c -> d -> e -> f) ->> v a -> v b -> v c -> v d -> v e -> v f> zipWith5 f a b c = zipWith ($) . zipWith4 f a b c> zipWith6 :: Vector v => (a -> b -> c -> d -> e -> f -> g) ->> v a -> v b -> v c -> v d -> v e -> v f -> v g> zipWith6 f a b c d = zipWith ($) . zipWith5 f a b c d> zipWith7 :: Vector v => (a -> b -> c -> d -> e -> f -> g -> h) ->> v a -> v b -> v c -> v d -> v e -> v f -> v g -> v h> zipWith7 f a b c d e = zipWith ($) . zipWith6 f a b c d e
Post a Comment