I like vectors! A long time ago my maths teacher introduced me to them and I just like the way three numbers (typically three) can express the notion of a position in space and combining these under different operations can produce other interesting properties.
Just recently I found this paper Learn Physics by Programming in Haskell which gives a very interesting Haskell oriented discussion on vectors and their use in mechanical problems. Reading it has inspired me to create my own implementation of vectors, bore you all senseless with it and see where it goes. What I’d like to do is:
- Create a module for Vectors. How far this will go I’m not sure yet but it will most likely be over several posts and might get as far as vector calculus. I’ll also try to look at Monoids, Monads etc. from a vector perspective and indeed see if there are meaningful vector oriented implementations of these classes. I don’t know yet. Almost all my Haskell posts are learning experiences for me!
- Apply the vector module to simple mechanical and geometric problems.
- Create some sort of visualisation of the solutions to these problems. Here the Haskell Gloss package is a prime contender.
- And it would be instructive to use QuickCheck too!
All nice and open ended and I’m writing the code more or less as I write the blog, so a lot may change! The code is in Github.
So here’s my first pass at a vector module.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 |
{-# LANGUAGE OverloadedStrings #-} {-# OPTIONS -Wall -fno-warn-type-defaults #-} module Vectors ( Scalar , Vector , (^+^) , (^-^) , (^*) , (*^) , (><) , (>.<) , (^/) , i , j , k , origin , normalise , neg , mag ) where type Scalar = Double type XYZ = (Scalar, Scalar, Scalar) newtype Vector = V { xyz :: XYZ} deriving (Show) -- Define operators using functions defines later infixl 6 ^+^ (^+^) :: Vector -> Vector -> Vector (^+^) = vAdd infixl 6 ^-^ (^-^) :: Vector -> Vector -> Vector (^-^) = vSub infixl 7 *^ (*^) :: Scalar -> Vector -> Vector (*^) = sMul infixl 7 ^* (^*) :: Vector -> Scalar -> Vector (^*) = sMul' infixl 7 ^/ (^/) :: Scalar -> Vector -> Vector (^/) = sDiv infixl 7 >.< (>.<) :: Vector -> Vector -> Scalar (>.<) = dot infixl 7 >< (><) :: Vector -> Vector -> Vector (><) = cross -- The 'traditional' i, j, k i, j, k, origin :: Vector i = V (1, 0, 0) j = V (0, 1, 0) k = V (0, 0, 1) origin = V (0, 0, 0) -- Normalise normalise :: Vector -> Vector normalise (V (0, 0, 0)) = V (0, 0, 0) normalise v = sDiv (mag v) v -- negate vector neg :: Vector -> Vector neg = mapVec ((-1) *) -- Magnitude of a vector mag :: Vector -> Scalar mag = sqrt . sumVec . mapVec (^2) -- ---------------------------------------------------------- -- Map a function over the internal tuple mapVec :: (Scalar -> Scalar) -> Vector -> Vector mapVec f (V (x, y, z)) = V (f x, f y, f z) -- Just sum the tuple entries sumVec :: Vector -> Scalar sumVec (V (x, y, z)) = x + y + z -- Apply a scalar function to x,y,z components of each vector zipWithVec :: (Scalar -> Scalar -> Scalar) -> Vector -> Vector -> Vector zipWithVec f (V (x, y, z)) (V (x', y', z')) = V (f x x', f y y', f z z') -- These are fairly obvious. The functions defined above -- make these functions simpler and cleaner. -- Just adding/subtracting the corresponding components, -- so zipWithVec using (+) or (-) vAdd :: Vector -> Vector -> Vector vAdd = zipWithVec (+) vSub :: Vector -> Vector -> Vector vSub = zipWithVec (-) -- Scalar then Vector - just map (*) over tuple entries sMul :: Scalar -> Vector -> Vector sMul s = mapVec (* s) -- or Vector then Scalar sMul' :: Vector -> Scalar -> Vector sMul' = flip sMul sDiv :: Scalar -> Vector -> Vector sDiv s = mapVec (/ s) -- dot product dot :: Vector -> Vector -> Scalar dot x = sumVec . zipWithVec (*) x -- Cross product - is there a neater way? -- note also that b×a=−a×b and a×a=0 - will use QuickCheck cross :: Vector -> Vector -> Vector cross (V (u1, u2, u3)) (V (v1, v2, v3)) = V (u2 * v3 - u3 * v2, u3 * v1 - u1 * v3, u1 * v2 - u2 * v1) |
The Vector itself is defined as a Haskell newtype. Internally the Vector has a three-tuple of Double (with type synonym Scalar). I went for a three-tuple rather than a list as all our vectors will be at most three dimensional and a tuple seemed a good fit. Maybe a list would be ‘better’ as it might simplify slightly the high-level functions mapVec, sumVec and zipWithVec that I’ve defined. The purpose of mapVec, sumVec and zipWithVec is to allow quite simple derivations of vector operations for addition, multiplication etc. along with the dot and cross product. (Having said all that I prefer the notation (x, y, z) to [x, y, z] as it looks more like how it would be written!)
Note also that there’s a number of custom operators, e.g. ‘^+^’, these allow us to write
a ^+^ b rather than
vAdd a b or
a `vAdd` b
We can load this into ghci and try a few examples…
Here we just declare a vector a with components (1,2,3)
|
1 2 3 4 |
ghci-> a = V (1,2,3) ghci-> a V {xyz = (1.0,2.0,3.0)} ghci-> |
The magnitude and negation (reversal) of vector a.
|
1 2 3 4 |
ghci-> mag a 3.7416573867739413 ghci-> neg a V {xyz = (-1.0,-2.0,-3.0)} |
This shows subtraction of vectors and the cross product.
|
1 2 3 4 5 6 7 8 9 |
ghci-> a = V (4,6,-1) ghci-> b = V (1,7,12) ghci-> ghci-> a ^-^ b V {xyz = (3.0,-1.0,-13.0)} ghci-> b ^-^ a V {xyz = (-3.0,1.0,13.0)} ghci-> a >< b V {xyz = (79.0,-49.0,22.0)} |
The dot product and its type.
|
1 2 3 4 |
ghci-> a >.< b 34.0 ghci-> :t a >.< b a >.< b :: Scalar |
In part two we will look at rendering Vectors to the screen to help visualise further work.
Thanks for reading!