packages feed

vect (empty) → 0.4.0

raw patch · 31 files changed

+4470/−0 lines, 31 filesdep +OpenGLdep +basedep +randombuild-type:Customsetup-changed

Dependencies added: OpenGL, base, random

Files

+ Data/Vect.hs view
@@ -0,0 +1,5 @@++-- | Importing this module is equivalent to importing "Data.Vect.Float".+module Data.Vect ( module Data.Vect.Float ) where+import Data.Vect.Float+
+ Data/Vect/Double.hs view
@@ -0,0 +1,23 @@+{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-}++module Data.Vect.Flt+  ( module Data.Vect.Flt.Base+  , module Data.Vect.Flt.Interpolate+  , module Data.Vect.Flt.Util.Dim2+  , module Data.Vect.Flt.Util.Dim3+  , module Data.Vect.Flt.Util.Projective+#ifdef VECT_OPENGL        +  , module Data.Vect.Flt.OpenGL       +#endif+  ) where++import Data.Vect.Flt.Base+import Data.Vect.Flt.Interpolate++import Data.Vect.Flt.Util.Dim2+import Data.Vect.Flt.Util.Dim3+import Data.Vect.Flt.Util.Projective++#ifdef VECT_OPENGL         +import Data.Vect.Flt.OpenGL       +#endif
+ Data/Vect/Double/Base.hs view
@@ -0,0 +1,735 @@+{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-}++module Data.Vect.Flt.Base where++import Control.Monad+import System.Random  +import Foreign++-- class declarations++class AbelianGroup g where+  (&+) :: g -> g -> g+  (&-) :: g -> g -> g+  neg  :: g -> g+  zero :: g++infixl 6 &++infixl 6 &- ++vecSum :: AbelianGroup g => [g] -> g+vecSum l = foldl (&+) zero l ++class (AbelianGroup r) => +      Ring r where+  (.*.) :: r -> r -> r+  one   :: r++infixl 7 .*. ++ringProduct :: Ring r => [r] -> r+ringProduct l = foldl (.*.) one l++class LeftModule r m where+  lmul :: r -> m -> m+  (*.) :: r -> m -> m+  (*.) = lmul++class RightModule m r where+  rmul :: m -> r -> m+  (.*) :: m -> r -> m+  (.*) = rmul++-- I'm not really sure about this.. may actually degrade the performance in some cases?  +{- RULES+"matrix multiplication left"   forall m n x.  (n .*. m) *. x = n *. (m *. x)  +"matrix multiplication right"  forall m n x.  x .* (m .*. n) = (x .* m) .* n+  -}++infixr 7 *.+infixl 7 .*++class AbelianGroup v => Vector v where+  mapVec    :: (Flt -> Flt) -> v -> v+  scalarMul :: Flt -> v -> v+  (*&) ::      Flt -> v -> v +  (&*) ::      v -> Flt -> v +  (*&) s v = scalarMul s v+  (&*) v s = scalarMul s v++infixr 7 *&+infixl 7 &*++{-# RULES+"scalar multiplication left"   forall s t x.  t *& (s *& x) = (t*s) *& x +"scalar multiplication right"  forall s t x.  (x &* s) &* t = x &* (s*t)  +  #-}++class DotProd v where+  (&.) :: v -> v -> Flt+  norm    :: v -> Flt+  normsqr :: v -> Flt+  len     :: v -> Flt+  lensqr  :: v -> Flt+  len = norm+  lensqr = normsqr+  dotprod :: v -> v -> Flt+  normsqr v = (v &. v)  +  norm = sqrt.lensqr+  dotprod = (&.)++infix 7 &.++{-# RULES+"len/square 1"   forall x.  (len x)*(len x) = lensqr x+"len/square 2"   forall x.  (len x)^2 = lensqr x+"norm/square 1"  forall x.  (norm x)*(norm x) = normsqr x+"norm/square 2"  forall x.  (norm x)^2 = normsqr x+  #-}++normalize :: (Vector v, DotProd v) => v -> v+normalize v = scalarMul (1.0/(len v)) v++distance :: (Vector v, DotProd v) => v -> v -> Flt+distance x y = norm (x &- y)++-- | the angle between two vectors+angle :: (Vector v, DotProd v) => v -> v -> Flt +angle x y = acos $ (x &. y) / (norm x * norm y)++-- | the angle between two unit vectors+angle' {- ' CPP is sensitive to primes -} :: (Vector v, UnitVector v u, DotProd v) => u -> u -> Flt +angle' x y = acos (fromNormal x &. fromNormal y)++{-# RULES+"normalize is idempotent"  forall x. normalize (normalize x) = normalize x+  #-}++class (Vector v, DotProd v) => UnitVector v u | v->u, u->v  where+  mkNormal         :: v -> u       -- ^ normalizes the input+  toNormalUnsafe   :: v -> u       -- ^ does not normalize the input!+  fromNormal       :: u -> v+  fromNormalRadius :: Flt -> u -> v+  fromNormalRadius t n = t *& fromNormal n ++-- | projects the first vector onto the direction of the second (unit) vector+project' :: (Vector v, UnitVector v u, DotProd v) => v -> u -> v+project' what dir = projectUnsafe what (fromNormal dir)++-- | direction (second argument) is assumed to be a /unit/ vector!+projectUnsafe :: (Vector v, DotProd v) => v -> v -> v+projectUnsafe what dir = what &- dir &* (what &. dir)++project :: (Vector v, DotProd v) => v -> v -> v+project what dir = what &- dir &* ((what &. dir) / (dir &. dir))++-- | since unit vectors are not a group, we need a separate function.+flipNormal :: UnitVector v n => n -> n +flipNormal = toNormalUnsafe . neg . fromNormal ++class CrossProd v where+  crossprod :: v -> v -> v+  (&^)      :: v -> v -> v+  (&^) = crossprod+  +class Pointwise v where+  pointwise :: v -> v -> v+  (&!)      :: v -> v -> v+  (&!) = pointwise ++infix 7 &^+infix 7 &!++class HasCoordinates v x | v->x where+  _1 :: v -> x+  _2 :: v -> x+  _3 :: v -> x+  _4 :: v -> x      ++-- | conversion between vectors (and matrices) of different dimensions+class Extend u v where+  extendZero :: u -> v          -- ^ example: @extendZero (Vec2 5 6) = Vec4 5 6 0 0@+  extendWith :: Flt -> u -> v   -- ^ example: @extendWith 1 (Vec2 5 6) = Vec4 5 6 1 1@+  trim :: v -> u                -- ^ example: @trim (Vec4 5 6 7 8) = Vec2 5 6@++-- | makes a diagonal matrix from a vector+class Diagonal s t | t->s where+  diag :: s -> t++class Matrix m where+  transpose :: m -> m +  inverse :: m -> m+  idmtx :: m++{-# RULES+"transpose is an involution"  forall m. transpose (transpose m) = m+"inverse is an involution"    forall m. inverse (inverse m) = m+  #-}+  +-- | Outer product (could be unified with Diagonal?)+class Tensor t v | t->v where+  outer :: v -> v -> t+    +class Determinant m where+  det :: m -> Flt    +    +-- Vec / Mat datatypes+ +data Vec2 = Vec2 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt +  deriving (Read,Show)+data Vec3 = Vec3 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt +  deriving (Read,Show)+data Vec4 = Vec4 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt +  deriving (Read,Show)++-- | these are /row/ vectors +data Mat2 = Mat2 !Vec2 !Vec2              deriving (Read,Show)+data Mat3 = Mat3 !Vec3 !Vec3 !Vec3        deriving (Read,Show)+data Mat4 = Mat4 !Vec4 !Vec4 !Vec4 !Vec4  deriving (Read,Show)++-- | The assumption when dealing with these is always that they are of unit length.+-- Also, interpolation works differently.+newtype Normal2 = Normal2 Vec2 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) +newtype Normal3 = Normal3 Vec3 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable,CrossProd) +newtype Normal4 = Normal4 Vec4 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) ++mkVec2 :: (Flt,Flt) -> Vec2+mkVec3 :: (Flt,Flt,Flt) -> Vec3+mkVec4 :: (Flt,Flt,Flt,Flt) -> Vec4++mkVec2 (x,y)     = Vec2 x y +mkVec3 (x,y,z)   = Vec3 x y z+mkVec4 (x,y,z,w) = Vec4 x y z w++-- Unit vectors+  +instance UnitVector Vec2 Normal2 where+  mkNormal v = Normal2 (normalize v)+  fromNormal (Normal2 v) = v +  toNormalUnsafe = Normal2++instance UnitVector Vec3 Normal3 where+  mkNormal v = Normal3 (normalize v)+  fromNormal (Normal3 v) = v +  toNormalUnsafe = Normal3++instance UnitVector Vec4 Normal4 where+  mkNormal v = Normal4 (normalize v)+  fromNormal (Normal4 v) = v +  toNormalUnsafe = Normal4++rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)+rndUnit g = +  if d > 0.01+    then ( v &* (1.0/d) , h )+    else rndUnit h+  where+    (v,h) = random g+    d = norm v+    +instance Random Normal2 where+  random g = let (v,h) = rndUnit g in (Normal2 v, h)  +  randomR _ = random++instance Random Normal3 where+  random g = let (v,h) = rndUnit g in (Normal3 v, h)  +  randomR _ = random++instance Random Normal4 where+  random g = let (v,h) = rndUnit g in (Normal4 v, h)  +  randomR _ = random++{-+instance Storable Normal2 where+  alignment _ = alignment (undefined::Vec2)+  sizeOf    _ = sizeOf    (undefined::Vec2)+  peek p = liftM (\v -> Normal2 v) (peek $ castPtr p)+  poke p (Normal2 v) = poke (castPtr p) v  + +instance Storable Normal3 where+  alignment _ = alignment (undefined::Vec3)+  sizeOf    _ = sizeOf    (undefined::Vec3)+  peek p = liftM (\v -> Normal3 v) (peek $ castPtr p)+  poke p (Normal3 v) = poke (castPtr p) v  ++instance Storable Normal4 where+  alignment _ = alignment (undefined::Vec4)+  sizeOf    _ = sizeOf    (undefined::Vec4)+  peek p = liftM (\v -> Normal4 v) (peek $ castPtr p)+  poke p (Normal4 v) = poke (castPtr p) v  +-}++-- Vec2 instances++instance HasCoordinates Vec2 Flt where+  _1 (Vec2 x _) = x+  _2 (Vec2 _ y) = y+  _3 _ = error "has only 2 coordinates"+  _4 _ = error "has only 2 coordinates"++instance AbelianGroup Vec2 where+  (&+) (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1+x2) (y1+y2) +  (&-) (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1-x2) (y1-y2)+  neg  (Vec2 x y)                = Vec2 (-x) (-y)+  zero = Vec2 0 0+  +instance Vector Vec2 where+  scalarMul s (Vec2 x y) = Vec2 (s*x) (s*y)+  mapVec    f (Vec2 x y) = Vec2 (f x) (f y)+  +instance DotProd Vec2 where+  (&.) (Vec2 x1 y1) (Vec2 x2 y2) = x1*x2 + y1*y2++instance Pointwise Vec2 where+  pointwise (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1*x2) (y1*y2)++instance Determinant (Vec2,Vec2) where+  det (Vec2 x1 y1 , Vec2 x2 y2) = x1*y2 - x2*y1  ++{-     +instance Show Vec2 where+  show (Vec2 x y) = "( " ++ show x ++ " , " ++ show y ++ " )"+-}++instance Random Vec2 where+  random = randomR (Vec2 (-1) (-1),Vec2 1 1)+  randomR (Vec2 a b, Vec2 c d) gen = +    let (x,gen1) = randomR (a,c) gen+        (y,gen2) = randomR (b,d) gen1+    in (Vec2 x y, gen2)+     +instance Storable Vec2 where+  sizeOf    _ = 2 * sizeOf (undefined::Flt)+  alignment _ = sizeOf (undefined::Flt)+  +  peek q = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    x <- peek        p +    y <- peekByteOff p k+    return (Vec2 x y)+    +  poke q (Vec2 x y) = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    poke        p   x+    pokeByteOff p k y+               +-- Mat2 instances++instance HasCoordinates Mat2 Vec2 where+  _1 (Mat2 x _) = x+  _2 (Mat2 _ y) = y+  _3 _ = error "has only 2 coordinates"+  _4 _ = error "has only 2 coordinates"++instance Matrix Mat2 where+  transpose (Mat2 row1 row2) = +    Mat2 (Vec2 (_1 row1) (_1 row2)) +         (Vec2 (_2 row1) (_2 row2)) +  idmtx = Mat2 (Vec2 1 0) (Vec2 0 1)+  inverse (Mat2 (Vec2 a b) (Vec2 c d)) = +    Mat2 (Vec2 (d*r) (-b*r)) (Vec2 (-c*r) (a*r)) +    where r = 1.0 / (a*d - b*c)++instance AbelianGroup Mat2 where+  (&+) (Mat2 r1 r2) (Mat2 s1 s2) = Mat2 (r1 &+ s1) (r2 &+ s2)+  (&-) (Mat2 r1 r2) (Mat2 s1 s2) = Mat2 (r1 &- s1) (r2 &- s2)+  neg  (Mat2 r1 r2)              = Mat2 (neg r1) (neg r2)  +  zero = Mat2 zero zero   -- (zero::Vec2) (zero::Vec2)++instance Vector Mat2 where+  scalarMul s (Mat2 r1 r2) = Mat2 (g r1) (g r2) where g = scalarMul s+  mapVec    f (Mat2 r1 r2) = Mat2 (g r1) (g r2) where g = mapVec f++instance Ring Mat2 where+  (.*.) (Mat2 r1 r2) n = +    let (Mat2 c1 c2) = transpose n+    in Mat2 (Vec2 (r1 &. c1) (r1 &. c2))+            (Vec2 (r2 &. c1) (r2 &. c2))+  one = idmtx ++instance LeftModule Mat2 Vec2 where+  lmul (Mat2 row1 row2) v = Vec2 (row1 &. v) (row2 &. v) +  +instance RightModule Vec2 Mat2 where+  rmul v mt = lmul (transpose mt) v++instance Diagonal Vec2 Mat2 where+  diag (Vec2 x y) = Mat2 (Vec2 x 0) (Vec2 0 y)++instance Tensor Mat2 Vec2 where+  outer (Vec2 a b) (Vec2 x y) = Mat2+    (Vec2 (a*x) (a*y))+    (Vec2 (b*x) (b*y))+{-+  outer v w = +    let full = Mat2 (Vec2 1 1) (Vec2 1 1)+    in  (diag v) .*. full .*. (diag w)+-}++instance Determinant Mat2 where+  det (Mat2 (Vec2 a b) (Vec2 c d)) = a*d - b*c ++{-+instance Show Mat2 where+  show (Mat2 r1 r2) = show r1 ++ "\n" ++ show r2+-}++instance Storable Mat2 where+  sizeOf    _ = 2 * sizeOf (undefined::Vec2)+  alignment _ = alignment  (undefined::Vec2)+  +  peek q = do+    let p = castPtr q :: Ptr Vec2+        k = sizeOf (undefined::Vec2)+    r1 <- peek        p +    r2 <- peekByteOff p k+    return (Mat2 r1 r2)+    +  poke q (Mat2 r1 r2) = do+    let p = castPtr q :: Ptr Vec2+        k = sizeOf (undefined::Vec2)+    poke        p   r1+    pokeByteOff p k r2++-- Vec3 instances++instance HasCoordinates Vec3 Flt where+  _1 (Vec3 x _ _) = x+  _2 (Vec3 _ y _) = y+  _3 (Vec3 _ _ z) = z+  _4 _ = error "has only 3 coordinates"++instance AbelianGroup Vec3 where+  (&+) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1+x2) (y1+y2) (z1+z2) +  (&-) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1-x2) (y1-y2) (z1-z2) +  neg  (Vec3 x y z)                    = Vec3 (-x) (-y) (-z)+  zero = Vec3 0 0 0+  +instance Vector Vec3 where+  scalarMul s (Vec3 x y z) = Vec3 (s*x) (s*y) (s*z)+  mapVec    f (Vec3 x y z) = Vec3 (f x) (f y) (f z)++instance DotProd Vec3 where+  (&.) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = x1*x2 + y1*y2 + z1*z2++instance Pointwise Vec3 where+  pointwise (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1*x2) (y1*y2) (z1*z2)++{-+instance Show Vec3 where+  show (Vec3 x y z) = "( " ++ show x ++ " , " ++ show y ++ " , " ++ show z ++ " )"+-}++instance Random Vec3 where+  random = randomR (Vec3 (-1) (-1) (-1),Vec3 1 1 1)+  randomR (Vec3 a b c, Vec3 d e f) gen = +    let (x,gen1) = randomR (a,d) gen+        (y,gen2) = randomR (b,e) gen1+        (z,gen3) = randomR (c,f) gen2  +    in (Vec3 x y z, gen3)+      +instance CrossProd Vec3 where+  crossprod (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (y1*z2-y2*z1) (z1*x2-z2*x1) (x1*y2-x2*y1) ++instance Determinant (Vec3,Vec3,Vec3) where+  det (u,v,w) = u &. (v &^ w)  + +instance Storable Vec3 where+  sizeOf    _ = 3 * sizeOf (undefined::Flt)+  alignment _ = sizeOf (undefined::Flt)+  +  peek q = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    x <- peek        p +    y <- peekByteOff p (k  )+    z <- peekByteOff p (k+k)+    return (Vec3 x y z)+    +  poke q (Vec3 x y z) = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    poke        p       x+    pokeByteOff p (k  ) y+    pokeByteOff p (k+k) z+   +-- Mat3 instances++instance HasCoordinates Mat3 Vec3 where+  _1 (Mat3 x _ _) = x+  _2 (Mat3 _ y _) = y+  _3 (Mat3 _ _ z) = z+  _4 _ = error "has only 3 coordinates"  ++instance Matrix Mat3 where++  transpose (Mat3 row1 row2 row3) = +    Mat3 (Vec3 (_1 row1) (_1 row2) (_1 row3)) +         (Vec3 (_2 row1) (_2 row2) (_2 row3)) +         (Vec3 (_3 row1) (_3 row2) (_3 row3)) +         +  idmtx = Mat3 (Vec3 1 0 0) (Vec3 0 1 0) (Vec3 0 0 1)+  +  inverse (Mat3 (Vec3 a b c) (Vec3 e f g) (Vec3 i j k)) = +    Mat3 (Vec3 (d11*r) (d21*r) (d31*r))  +         (Vec3 (d12*r) (d22*r) (d32*r))  +         (Vec3 (d13*r) (d23*r) (d33*r))  +    where+      r = 1.0 / ( a*d11 + b*d12 + c*d13 )++      d11 = f*k - g*j+      d12 = g*i - e*k+      d13 = e*j - f*i++      d31 = b*g - c*f+      d32 = c*e - a*g+      d33 = a*f - b*e++      d21 = c*j - b*k +      d22 = a*k - c*i +      d23 = b*i - a*j ++instance AbelianGroup Mat3 where+  (&+) (Mat3 r1 r2 r3) (Mat3 s1 s2 s3) = Mat3 (r1 &+ s1) (r2 &+ s2) (r3 &+ s3)+  (&-) (Mat3 r1 r2 r3) (Mat3 s1 s2 s3) = Mat3 (r1 &- s1) (r2 &- s2) (r3 &- s3)+  neg  (Mat3 r1 r2 r3)                 = Mat3 (neg r1) (neg r2) (neg r3) +  zero = Mat3 zero zero zero   -- (zero::Vec3) (zero::Vec3) (zero::Vec3)++instance Vector Mat3 where+  scalarMul s (Mat3 r1 r2 r3) = Mat3 (g r1) (g r2) (g r3) where g = scalarMul s+  mapVec    f (Mat3 r1 r2 r3) = Mat3 (g r1) (g r2) (g r3) where g = mapVec f++instance Ring Mat3 where+  (.*.) (Mat3 r1 r2 r3) n = +    let (Mat3 c1 c2 c3) = transpose n+    in Mat3 (Vec3 (r1 &. c1) (r1 &. c2) (r1 &. c3))+            (Vec3 (r2 &. c1) (r2 &. c2) (r2 &. c3))+            (Vec3 (r3 &. c1) (r3 &. c2) (r3 &. c3))+  one = idmtx ++instance LeftModule Mat3 Vec3 where+  lmul (Mat3 row1 row2 row3) v = Vec3 (row1 &. v) (row2 &. v) (row3 &. v)+  +instance RightModule Vec3 Mat3 where+  rmul v mt = lmul (transpose mt) v++instance Diagonal Vec3 Mat3 where+  diag (Vec3 x y z) = Mat3 (Vec3 x 0 0) (Vec3 0 y 0) (Vec3 0 0 z)++instance Tensor Mat3 Vec3 where+  outer (Vec3 a b c) (Vec3 x y z) = Mat3+    (Vec3 (a*x) (a*y) (a*z))+    (Vec3 (b*x) (b*y) (b*z))+    (Vec3 (c*x) (c*y) (c*z))+{-+  outer v w = +    let full = Mat3 (Vec3 1 1 1) (Vec3 1 1 1) (Vec3 1 1 1)+    in  (diag v) .*. full .*. (diag w)+-}++instance Determinant Mat3 where+  det (Mat3 r1 r2 r3) = det (r1,r2,r3)++{-+instance Show Mat3 where+  show (Mat3 r1 r2 r3) = show r1 ++ "\n" ++ show r2 ++ "\n" ++ show r3+-}++instance Storable Mat3 where+  sizeOf    _ = 3 * sizeOf (undefined::Vec3)+  alignment _ = alignment  (undefined::Vec3)+  +  peek q = do+    let p = castPtr q :: Ptr Vec3+        k = sizeOf (undefined::Vec3)+    r1 <- peek        p +    r2 <- peekByteOff p (k  )+    r3 <- peekByteOff p (k+k)+    return (Mat3 r1 r2 r3)+    +  poke q (Mat3 r1 r2 r3) = do+    let p = castPtr q :: Ptr Vec3+        k = sizeOf (undefined::Vec3)+    poke        p       r1+    pokeByteOff p (k  ) r2+    pokeByteOff p (k+k) r3++-- Vec4 instances++instance HasCoordinates Vec4 Flt where+  _1 (Vec4 x _ _ _) = x+  _2 (Vec4 _ y _ _) = y+  _3 (Vec4 _ _ z _) = z+  _4 (Vec4 _ _ _ w) = w++instance AbelianGroup Vec4 where+  (&+) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1+x2) (y1+y2) (z1+z2) (w1+w2)+  (&-) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1-x2) (y1-y2) (z1-z2) (w1-w2)+  neg  (Vec4 x y z w)                        = Vec4 (-x) (-y) (-z) (-w)+  zero = Vec4 0 0 0 0+  +instance Vector Vec4 where+  scalarMul s (Vec4 x y z w) = Vec4 (s*x) (s*y) (s*z) (s*w)+  mapVec    f (Vec4 x y z w) = Vec4 (f x) (f y) (f z) (f w)++instance DotProd Vec4 where+  (&.) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = x1*x2 + y1*y2 + z1*z2 + w1*w2++instance Pointwise Vec4 where+  pointwise (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1*x2) (y1*y2) (z1*z2) (w1*w2)++{-+instance Show Vec4 where+  show (Vec4 x y z w) = "( " ++ show x ++ " , " ++ show y ++ " , " ++ show z ++ " , " ++ show w ++ " )"+-}++instance Random Vec4 where+  random = randomR (Vec4 (-1) (-1) (-1) (-1),Vec4 1 1 1 1)+  randomR (Vec4 a b c d, Vec4 e f g h) gen = +    let (x,gen1) = randomR (a,e) gen+        (y,gen2) = randomR (b,f) gen1+        (z,gen3) = randomR (c,g) gen2  +        (w,gen4) = randomR (d,h) gen3  +    in (Vec4 x y z w, gen4)+           +instance Storable Vec4 where+  sizeOf    _ = 4 * sizeOf (undefined::Flt)+  alignment _ = sizeOf (undefined::Flt)+  +  peek q = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    x <- peek        p +    y <- peekByteOff p (k  )+    z <- peekByteOff p (k+k)+    w <- peekByteOff p (3*k)+    return (Vec4 x y z w)+    +  poke q (Vec4 x y z w) = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    poke        p       x+    pokeByteOff p (k  ) y+    pokeByteOff p (k+k) z+    pokeByteOff p (3*k) w++-- Mat4 instances++instance HasCoordinates Mat4 Vec4 where+  _1 (Mat4 x _ _ _) = x+  _2 (Mat4 _ y _ _) = y+  _3 (Mat4 _ _ z _) = z+  _4 (Mat4 _ _ _ w) = w++instance Matrix Mat4 where+  transpose (Mat4 row1 row2 row3 row4) = +    Mat4 (Vec4 (_1 row1) (_1 row2) (_1 row3) (_1 row4)) +         (Vec4 (_2 row1) (_2 row2) (_2 row3) (_2 row4)) +         (Vec4 (_3 row1) (_3 row2) (_3 row3) (_3 row4)) +         (Vec4 (_4 row1) (_4 row2) (_4 row3) (_4 row4)) +  idmtx = Mat4 (Vec4 1 0 0 0) (Vec4 0 1 0 0) (Vec4 0 0 1 0) (Vec4 0 0 0 1)+  inverse = error "inverse/Mat4: not implemented yet"++instance AbelianGroup Mat4 where+  (&+) (Mat4 r1 r2 r3 r4) (Mat4 s1 s2 s3 s4) = Mat4 (r1 &+ s1) (r2 &+ s2) (r3 &+ s3) (r4 &+ s4)+  (&-) (Mat4 r1 r2 r3 r4) (Mat4 s1 s2 s3 s4) = Mat4 (r1 &- s1) (r2 &- s2) (r3 &- s3) (r4 &- s4)+  neg  (Mat4 r1 r2 r3 r4)                    = Mat4 (neg r1) (neg r2) (neg r3) (neg r4) +  zero = Mat4 zero zero zero zero+  +instance Vector Mat4 where+  scalarMul s (Mat4 r1 r2 r3 r4) = Mat4 (g r1) (g r2) (g r3) (g r4) where g = scalarMul s+  mapVec    f (Mat4 r1 r2 r3 r4) = Mat4 (g r1) (g r2) (g r3) (g r4) where g = mapVec f++instance Ring Mat4 where+  (.*.) (Mat4 r1 r2 r3 r4) n = +    let (Mat4 c1 c2 c3 c4) = transpose n+    in Mat4 (Vec4 (r1 &. c1) (r1 &. c2) (r1 &. c3) (r1 &. c4))+            (Vec4 (r2 &. c1) (r2 &. c2) (r2 &. c3) (r2 &. c4))+            (Vec4 (r3 &. c1) (r3 &. c2) (r3 &. c3) (r3 &. c4))+            (Vec4 (r4 &. c1) (r4 &. c2) (r4 &. c3) (r4 &. c4))+  one = idmtx ++instance LeftModule Mat4 Vec4 where+  lmul (Mat4 row1 row2 row3 row4) v = Vec4 (row1 &. v) (row2 &. v) (row3 &. v) (row4 &. v)+  +instance RightModule Vec4 Mat4 where+  rmul v mt = lmul (transpose mt) v++instance Diagonal Vec4 Mat4 where+  diag (Vec4 x y z w) = Mat4 (Vec4 x 0 0 0) (Vec4 0 y 0 0) (Vec4 0 0 z 0) (Vec4 0 0 0 w)++instance Tensor Mat4 Vec4 where+  outer (Vec4 a b c d) (Vec4 x y z w) = Mat4+    (Vec4 (a*x) (a*y) (a*z) (a*w))+    (Vec4 (b*x) (b*y) (b*z) (b*w))+    (Vec4 (c*x) (c*y) (c*z) (c*w))+    (Vec4 (d*x) (d*y) (d*z) (d*w))+{-+  outer v w = +    let full = Mat4 (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1)+    in  (diag v) .*. full .*. (diag w)+-}++--instance Determinant Mat4 where+--  det (Mat4 r1 r2 r3 r4)++{-+instance Show Mat4 where+  show (Mat4 r1 r2 r3 r4) = show r1 ++ "\n" ++ show r2 ++ "\n" ++ show r3 ++ "\n" ++ show r4+-}++instance Storable Mat4 where+  sizeOf    _ = 4 * sizeOf (undefined::Vec4)+  alignment _ = alignment  (undefined::Vec4)+  +  peek q = do+    let p = castPtr q :: Ptr Vec4+        k = sizeOf (undefined::Vec4)+    r1 <- peek        p +    r2 <- peekByteOff p (k  )+    r3 <- peekByteOff p (k+k)+    r4 <- peekByteOff p (3*k)+    return (Mat4 r1 r2 r3 r4)+    +  poke q (Mat4 r1 r2 r3 r4) = do+    let p = castPtr q :: Ptr Vec4+        k = sizeOf (undefined::Vec4)+    poke        p       r1+    pokeByteOff p (k  ) r2+    pokeByteOff p (k+k) r3+    pokeByteOff p (3*k) r4++-- Extend instances++instance Extend Vec2 Vec3 where+  extendZero   (Vec2 x y) = Vec3 x y 0+  extendWith t (Vec2 x y) = Vec3 x y t+  trim (Vec3 x y _)       = Vec2 x y++instance Extend Vec2 Vec4 where+  extendZero   (Vec2 x y) = Vec4 x y 0 0+  extendWith t (Vec2 x y) = Vec4 x y t t+  trim (Vec4 x y _ _)     = Vec2 x y ++instance Extend Vec3 Vec4 where+  extendZero   (Vec3 x y z) = Vec4 x y z 0+  extendWith t (Vec3 x y z) = Vec4 x y z t+  trim (Vec4 x y z _)       = Vec3 x y z++instance Extend Mat2 Mat3 where+  extendZero (Mat2 p q) = Mat3 (extendZero p) (extendZero q) zero+  extendWith _ _ = error "extendWith is meaningless for matrices"+  trim (Mat3 p q _) = Mat2 (trim p) (trim q)++instance Extend Mat2 Mat4 where+  extendZero (Mat2 p q) = Mat4 (extendZero p) (extendZero q) zero zero+  extendWith _ _ = error "extendWith is meaningless for matrices"+  trim (Mat4 p q _ _) = Mat2 (trim p) (trim q)++instance Extend Mat3 Mat4 where+  extendZero (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) zero+  extendWith _ _ = error "extendWith is meaningless for matrices"+  trim (Mat4 p q r _) = Mat3 (trim p) (trim q) (trim r)+
+ Data/Vect/Double/GramSchmidt.hs view
@@ -0,0 +1,135 @@+{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-}++-- | Gram-Schmidt orthogonalization.+-- This module is not re-exported by "Data.Vect".++module Data.Vect.Flt.GramSchmidt +  ( GramSchmidt(..)+  )+  where++import Data.Vect.Flt.Base++-------------------------------------------------------++liftPair :: (a -> b) -> (a,a) -> (b,b)+liftPair f (x,y) = (f x, f y)++liftTriple :: (a -> b) -> (a,a,a) -> (b,b,b)+liftTriple f (x,y,z) = (f x, f y, f z)++liftQuadruple :: (a -> b) -> (a,a,a,a) -> (b,b,b,b)+liftQuadruple f (x,y,z,w) = (f x, f y, f z, f w)++-------------------------------------------------------+    +-- | produces orthogonal\/orthonormal vectors from a set of vectors    +class GramSchmidt a where+  gramSchmidt          :: a -> a   -- ^ does not normalize the vectors!+  gramSchmidtNormalize :: a -> a   -- ^ normalizes the vectors.++{-# RULES+"gramSchmidt is idempotent"  forall a. gramSchmidt (gramSchmidt a) = gramSchmidt a +"gramSchmidtNormalize is idempotent"  forall a. gramSchmidtNormalize (gramSchmidtNormalize a) = gramSchmidtNormalize a +  #-}++-------------------------------------------------------++instance GramSchmidt (Vec2,Vec2) where+  gramSchmidt = gramSchmidtPair+  gramSchmidtNormalize = gramSchmidtNormalizePair+  +instance GramSchmidt (Vec3,Vec3) where+  gramSchmidt = gramSchmidtPair+  gramSchmidtNormalize = gramSchmidtNormalizePair+  +instance GramSchmidt (Vec4,Vec4) where+  gramSchmidt = gramSchmidtPair+  gramSchmidtNormalize = gramSchmidtNormalizePair++----------++instance GramSchmidt (Normal2,Normal2) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal2!"+  gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal++instance GramSchmidt (Normal3,Normal3) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal3!"+  gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal++instance GramSchmidt (Normal4,Normal4) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal4!"+  gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal++----------+  +gramSchmidtPair :: (Vector v, DotProd v) => (v,v) -> (v,v)+gramSchmidtPair (u,v) = (u',v') where +  u' = u+  v' = project v u'     +  +gramSchmidtNormalizePair :: (Vector v, DotProd v) => (v,v) -> (v,v)+gramSchmidtNormalizePair (u,v) = (u',v') where+  u' = normalize u +  v' = normalize $ projectUnsafe v u'     ++----------++instance GramSchmidt (Vec3,Vec3,Vec3) where+  gramSchmidt = gramSchmidtTriple+  gramSchmidtNormalize = gramSchmidtNormalizeTriple+     +instance GramSchmidt (Vec4,Vec4,Vec4) where+  gramSchmidt = gramSchmidtTriple+  gramSchmidtNormalize = gramSchmidtNormalizeTriple++instance GramSchmidt (Normal3,Normal3,Normal3) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal3!"+  gramSchmidtNormalize = liftTriple toNormalUnsafe . gramSchmidtNormalizeTriple . liftTriple fromNormal++instance GramSchmidt (Normal4,Normal4,Normal4) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal4!"+  gramSchmidtNormalize = liftTriple toNormalUnsafe . gramSchmidtNormalizeTriple . liftTriple fromNormal++----------++gramSchmidtTriple :: (Vector v, DotProd v) => (v,v,v) -> (v,v,v)+gramSchmidtTriple (u,v,w) = (u',v',w') where +  u' = u+  v' = project v u'     +  w' = project (project w u') v' +  +gramSchmidtNormalizeTriple :: (Vector v, DotProd v) => (v,v,v) -> (v,v,v)+gramSchmidtNormalizeTriple (u,v,w) = (u',v',w') where+  u' = normalize $ u +  v' = normalize $ projectUnsafe v u'     +  w' = normalize $ projectUnsafe (projectUnsafe w u') v'     ++----------++instance GramSchmidt (Vec4,Vec4,Vec4,Vec4) where+  gramSchmidt          = gramSchmidtQuadruple+  gramSchmidtNormalize = gramSchmidtNormalizeQuadruple ++instance GramSchmidt (Normal4,Normal4,Normal4,Normal4) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal4!"+  gramSchmidtNormalize = liftQuadruple toNormalUnsafe . gramSchmidtNormalizeQuadruple . liftQuadruple fromNormal++----------+  +gramSchmidtQuadruple :: (Vector v, DotProd v) => (v,v,v,v) -> (v,v,v,v)+gramSchmidtQuadruple (u,v,w,z) = (u',v',w',z') where +  u' = u+  v' = project v u'     +  w' = project (project w u') v' +  z' = project (project (project z u') v') w'++gramSchmidtNormalizeQuadruple :: (Vector v, DotProd v) => (v,v,v,v) -> (v,v,v,v)+gramSchmidtNormalizeQuadruple (u,v,w,z) = (u',v',w',z') where+  u' = normalize $ u+  v' = normalize $ projectUnsafe v u'     +  w' = normalize $ projectUnsafe (projectUnsafe w u') v' +  z' = normalize $ projectUnsafe (projectUnsafe (projectUnsafe z u') v') w'+  +----------+
+ Data/Vect/Double/Interpolate.hs view
@@ -0,0 +1,43 @@+{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-}++-- | Interpolation of vectors. +-- Note: we interpolate unit vectors differently from ordinary vectors.++module Data.Vect.Flt.Interpolate where++import Data.Vect.Flt.Base+import Data.Vect.Flt.Util.Dim2 (sinCos',angle2')+import Data.Vect.Flt.Util.Dim3 (rotate3')++class Interpolate v where+  interpolate :: Flt -> v -> v -> v+  +instance Interpolate Flt where+  interpolate t x y = x + t*(y-x)++instance Interpolate Vec2 where interpolate t x y = x &+ t *& (y &- x)+instance Interpolate Vec3 where interpolate t x y = x &+ t *& (y &- x)+instance Interpolate Vec4 where interpolate t x y = x &+ t *& (y &- x)++instance Interpolate Normal2 where+  interpolate t nx ny = sinCos' $ ax + t*adiff where+    ax = angle2' nx+    ay = angle2' ny+    adiff = helper (ay - ax)+    helper d +      | d < -pi   = d + twopi+      | d >  pi   = d - twopi+      | otherwise = d+    twopi = 2*pi+    +instance Interpolate Normal3 where +  interpolate t nx ny = +    if maxAngle < 0.001  -- more or less ad-hoc critical angle+      then mkNormal $ interpolate t x y+      else toNormalUnsafe $ rotate3' (t*maxAngle) (mkNormal axis) x where+    x = fromNormal nx+    y = fromNormal ny+    axis = (x &^ y)+    maxAngle = acos (x &. y)+        +
+ Data/Vect/Double/OpenGL.hs view
@@ -0,0 +1,183 @@+{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-}++-- TODO: the pointer versions of these functions should be really implemented +-- via the pointer versions of the original opengl functions...++-- | OpenGL support, inclduing 'vertex', 'texCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.+ +module Data.Vect.Flt.OpenGL where++import Control.Monad+import Data.Vect.Flt.Base+import qualified Graphics.Rendering.OpenGL as GL++import Foreign++import Graphics.Rendering.OpenGL hiding (Normal3,rotate,translate,scale)++-------------------------------------------------------++{-# SPECIALISE radianToDegrees :: Float -> Float #-}+{-# SPECIALISE radianToDegrees :: Double -> Double #-}+radianToDegrees :: RealFrac a => a -> a+radianToDegrees x = x * 57.295779513082322++{-# SPECIALIZE degreesToRadian :: Float  -> Float  #-}+{-# SPECIALIZE degreesToRadian :: Double -> Double #-}+degreesToRadian :: Floating a => a -> a+degreesToRadian x = x * 1.7453292519943295e-2++-- | The angle is in radians. (WARNING: OpenGL uses degrees!)+rotate :: Flt -> Vec3 -> IO ()+rotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z)++translate :: Vec3 -> IO ()+translate (Vec3 x y z) = GL.translate (Vector3 x y z)++scale3 :: Vec3 -> IO ()+scale3 (Vec3 x y z) = GL.scale x y z++scale :: Flt -> IO ()+scale x = GL.scale x x x++-------------------------------------------------------++-- Vertex instances++instance GL.Vertex Vec2 where+  vertex (Vec2 x y) = GL.vertex (GL.Vertex2 x y)+  vertexv p = peek p >>= vertex +  +instance GL.Vertex Vec3 where+  vertex (Vec3 x y z) = GL.vertex (GL.Vertex3 x y z)+  vertexv p = peek p >>= vertex   +  +instance GL.Vertex Vec4 where+  vertex (Vec4 x y z w) = GL.vertex (GL.Vertex4 x y z w)+  vertexv p = peek p >>= vertex   ++-------------------------------------------------------++-- the Normal instance+-- note that there is no Normal2\/Normal4 in the OpenGL binding++instance GL.Normal Normal3 where+  normal (Normal3 (Vec3 x y z)) = GL.normal (GL.Normal3 x y z)+  normalv p = peek p >>= normal ++-------------------------------------------------------++-- Color instances+  +instance GL.Color Vec3 where+  color (Vec3 r g b) = GL.color (GL.Color3 r g b)+  colorv p = peek p >>= color++instance GL.Color Vec4 where+  color (Vec4 r g b a) = GL.color (GL.Color4 r g b a)+  colorv p = peek p >>= color++instance GL.SecondaryColor Vec3 where+  secondaryColor (Vec3 r g b) = GL.secondaryColor (GL.Color3 r g b)+  secondaryColorv p = peek p >>= secondaryColor++{-+-- there is no such thing?+instance GL.SecondaryColor Vec4 where+  secondaryColor (Vec4 r g b a) = GL.secondaryColor (GL.Color4 r g b a)+  secondaryColorv p = peek p >>= secondaryColor+-}++-------------------------------------------------------++-- TexCoord instances++instance GL.TexCoord Vec2 where+  texCoord (Vec2 u v) = GL.texCoord (GL.TexCoord2 u v)+  texCoordv p = peek p >>= texCoord+  multiTexCoord unit (Vec2 u v) = GL.multiTexCoord unit (GL.TexCoord2 u v)+  multiTexCoordv unit p = peek p >>= multiTexCoord unit++instance GL.TexCoord Vec3 where+  texCoord (Vec3 u v w) = GL.texCoord (GL.TexCoord3 u v w)+  texCoordv p = peek p >>= texCoord+  multiTexCoord unit (Vec3 u v w) = GL.multiTexCoord unit (GL.TexCoord3 u v w)+  multiTexCoordv unit p = peek p >>= multiTexCoord unit++instance GL.TexCoord Vec4 where+  texCoord (Vec4 u v w z) = GL.texCoord (GL.TexCoord4 u v w z)+  texCoordv p = peek p >>= texCoord+  multiTexCoord unit (Vec4 u v w z) = GL.multiTexCoord unit (GL.TexCoord4 u v w z)+  multiTexCoordv unit p = peek p >>= multiTexCoord unit++-------------------------------------------------------+    +-- Vertex Attributes (experimental)++class VertexAttrib' a where+  vertexAttrib :: GL.AttribLocation -> a -> IO ()+  +instance VertexAttrib' {- ' CPP is sensitive to primes -} Flt where+  vertexAttrib loc x = GL.vertexAttrib1 loc x++instance VertexAttrib' Vec2 where+  vertexAttrib loc (Vec2 x y) = GL.vertexAttrib2 loc x y++instance VertexAttrib' Vec3 where+  vertexAttrib loc (Vec3 x y z) = GL.vertexAttrib3 loc x y z ++instance VertexAttrib' Vec4 where+  vertexAttrib loc (Vec4 x y z w) = GL.vertexAttrib4 loc x y z w ++instance VertexAttrib' Normal2 where+  vertexAttrib loc (Normal2 (Vec2 x y)) = GL.vertexAttrib2 loc x y++instance VertexAttrib' Normal3 where+  vertexAttrib loc (Normal3 (Vec3 x y z)) = GL.vertexAttrib3 loc x y z++instance VertexAttrib' Normal4 where+  vertexAttrib loc (Normal4 (Vec4 x y z w)) = GL.vertexAttrib4 loc x y z w++-------------------------------------------------------++-- Uniform (again, experimental)++-- (note that the uniform location code in the OpenGL 2.2.1.1 is broken; +-- a work-around is to put a zero character at the end of uniform names)++{-+toFloat :: Flt -> Float+toFloat = realToFrac++fromFloat :: Float -> Flt+fromFloat = realToFrac+-}++-- Uniforms are always floats...+#ifdef VECT_Float++instance GL.Uniform Flt where+  uniform loc = GL.makeStateVar getter setter where+    getter = liftM (\(GL.Index1 x) -> x) $ get (uniform loc)+    setter x = ($=) (uniform loc) (Index1 x) +  uniformv loc cnt ptr = uniformv loc cnt (castPtr ptr :: Ptr (Index1 Flt))++instance GL.Uniform Vec2 where+  uniform loc = GL.makeStateVar getter setter where+    getter = liftM (\(GL.Vertex2 x y) -> Vec2 x y) $ get (uniform loc)+    setter (Vec2 x y) = ($=) (uniform loc) (Vertex2 x y) +  uniformv loc cnt ptr = uniformv loc (2*cnt) (castPtr ptr :: Ptr Flt)++instance GL.Uniform Vec3 where+  uniform loc = GL.makeStateVar getter setter where+    getter = liftM (\(GL.Vertex3 x y z) -> Vec3 x y z) $ get (uniform loc)+    setter (Vec3 x y z) = ($=) (uniform loc) (Vertex3 x y z) +  uniformv loc cnt ptr = uniformv loc (3*cnt) (castPtr ptr :: Ptr Flt)++instance GL.Uniform Vec4 where+  uniform loc = GL.makeStateVar getter setter where+    getter = liftM (\(GL.Vertex4 x y z w) -> Vec4 x y z w) $ get (uniform loc)+    setter (Vec4 x y z w) = ($=) (uniform loc) (Vertex4 x y z w) +  uniformv loc cnt ptr = uniformv loc (4*cnt) (castPtr ptr :: Ptr Flt)+    +#endif
+ Data/Vect/Double/Util/Dim2.hs view
@@ -0,0 +1,65 @@+{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-}++module Data.Vect.Flt.Util.Dim2 where++import Data.Vect.Flt.Base++-- |example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@.+structVec2 :: [Flt] -> [Vec2]+structVec2 [] = []+structVec2 (x:y:ls) = (Vec2 x y):(structVec2 ls) +structVec2 _ = error "structVec2"++destructVec2 :: [Vec2] -> [Flt]+destructVec2 [] = []+destructVec2 ((Vec2 x y):ls) = x:y:(destructVec2 ls)  ++det2 :: Vec2 -> Vec2 -> Flt+det2 u v = det (u,v)++vec2X :: Vec2+vec2Y :: Vec2++vec2X = Vec2 1 0 +vec2Y = Vec2 0 1 ++translate2X :: Flt -> Vec2 -> Vec2+translate2Y :: Flt -> Vec2 -> Vec2++translate2X t (Vec2 x y) = Vec2 (x+t) y +translate2Y t (Vec2 x y) = Vec2 x (y+t) ++-- | unit vector with given angle relative to the positive X axis (in the positive direction, that is, CCW).+-- A more precise name would be @cosSin@, but that sounds bad :)+sinCos :: Flt -> Vec2+sinCos a = Vec2 (cos a) (sin a)++sinCos' {- ' CPP is sensitive to primes -} :: Flt -> Normal2+sinCos' = toNormalUnsafe . sinCos++sinCosRadius :: Flt    -- ^ angle (in radians)+             -> Flt    -- ^ radius+             -> Vec2+sinCosRadius a r = Vec2 (r * cos a) (r * sin a)++-- | The angle relative to the positive X axis+angle2 :: Vec2 -> Flt+angle2 (Vec2 x y) = atan2 y x++angle2' {- ' CPP is sensitive to primes -} :: Normal2 -> Flt+angle2' = angle2 . fromNormal++-- |Rotation matrix by a given angle (in radians), counterclockwise.+rotMatrix2 :: Flt -> Mat2+rotMatrix2 a = Mat2 (Vec2 c s) (Vec2 (-s) c) where c = cos a; s = sin a++rotate2 :: Flt -> Vec2 -> Vec2+rotate2 a v = v .* (rotMatrix2 a) ++-- |Rotates counterclockwise by 90 degrees.+rotateCCW :: Vec2 -> Vec2+rotateCCW (Vec2 x y) = Vec2 (-y) x++-- |Rotates clockwise by 90 degrees.+rotateCW :: Vec2 -> Vec2+rotateCW (Vec2 x y) = Vec2 y (-x)
+ Data/Vect/Double/Util/Dim3.hs view
@@ -0,0 +1,76 @@+{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-}++module Data.Vect.Flt.Util.Dim3 where++import Data.Vect.Flt.Base++structVec3 :: [Flt] -> [Vec3]+structVec3 [] = []+structVec3 (x:y:z:ls) = (Vec3 x y z):(structVec3 ls) +structVec3 _ = error "structVec3"++destructVec3 :: [Vec3] -> [Flt]+destructVec3 [] = []+destructVec3 ((Vec3 x y z):ls) = x:y:z:(destructVec3 ls)  ++det3 :: Vec3 -> Vec3 -> Vec3 -> Flt+det3 u v w = det (u,v,w)++translate3X :: Flt -> Vec3 -> Vec3+translate3Y :: Flt -> Vec3 -> Vec3+translate3Z :: Flt -> Vec3 -> Vec3++translate3X t (Vec3 x y z) = Vec3 (x+t) y z +translate3Y t (Vec3 x y z) = Vec3 x (y+t) z +translate3Z t (Vec3 x y z) = Vec3 x y (z+t) ++vec3X :: Vec3+vec3Y :: Vec3+vec3Z :: Vec3++vec3X = Vec3 1 0 0+vec3Y = Vec3 0 1 0+vec3Z = Vec3 0 0 1++rotMatrixZ :: Flt -> Mat3+rotMatrixY :: Flt -> Mat3+rotMatrixX :: Flt -> Mat3++-- These are intended for multiplication on the /right/.+-- Should be consistent with the rotation around an arbitrary axis +-- (eg, @rotMatrixY a == rotate3 a vec3Y@)+rotMatrixZ a = Mat3 (Vec3 c s 0) (Vec3 (-s) c 0) (Vec3 0 0 1) where c = cos a; s = sin a+rotMatrixY a = Mat3 (Vec3 c 0 (-s)) (Vec3 0 1 0) (Vec3 s 0 c) where c = cos a; s = sin a+rotMatrixX a = Mat3 (Vec3 1 0 0) (Vec3 0 c s) (Vec3 0 (-s) c) where c = cos a; s = sin a++rotate3' :: {- ' CPP is sensitive to primes -} Flt       -- ^ angle (in radians)+         -> Normal3   -- ^ axis (should be a /unit/ vector!) +         -> Vec3      -- ^ vector+         -> Vec3      -- ^ result+rotate3' angle axis v = v .* (rotMatrix3' axis angle)++rotate3 :: Flt    -- ^ angle (in radians)+        -> Vec3   -- ^ axis (arbitrary nonzero vector)+        -> Vec3   -- ^ vector+        -> Vec3   -- ^ result+rotate3 angle axis v = v .* (rotMatrix3 axis angle)+      +-- |Rotation around an arbitrary 3D vector. The resulting 3x3 matrix is intended for multiplication on the /right/. +rotMatrix3 :: Vec3 -> Flt -> Mat3+rotMatrix3 v a = rotMatrix3' (mkNormal v) a++-- |Rotation around an arbitrary 3D /unit/ vector. The resulting 3x3 matrix is intended for multiplication on the /right/. +rotMatrix3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Mat3+rotMatrix3' (Normal3 v) a = +  let c = cos a+      s = sin a+      m1 = scalarMul (1-c) (outer v v)+      x = _1 v+      y = _2 v+      z = _3 v+      m2 = Mat3 (Vec3   c    ( s*z) (-s*y))+                (Vec3 (-s*z)   c    ( s*x))+                (Vec3 ( s*y) (-s*x)   c   )+  in (m1 &+ m2)++
+ Data/Vect/Double/Util/Dim4.hs view
@@ -0,0 +1,93 @@+{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-}++-- | Rotation around an arbitrary plane in four dimensions, and other miscellanea.+-- Not very useful for most people, and not re-exported by "Data.Vect".++module Data.Vect.Flt.Util.Dim4 where++import Data.Vect.Flt.Base+import Data.Vect.Flt.GramSchmidt++structVec4 :: [Flt] -> [Vec4]+structVec4 [] = []+structVec4 (x:y:z:w:ls) = (Vec4 x y z w):(structVec4 ls) +structVec4 _ = error "structVec4"++destructVec4 :: [Vec4] -> [Flt]+destructVec4 [] = []+destructVec4 ((Vec4 x y z w):ls) = x:y:z:w:(destructVec4 ls)  ++--det4 :: Vec4 -> Vec4 -> Vec4 -> Vec4 -> Flt+--det4 u v w z = det (u,v,w,z)++translate4X :: Flt -> Vec4 -> Vec4+translate4Y :: Flt -> Vec4 -> Vec4+translate4Z :: Flt -> Vec4 -> Vec4+translate4W :: Flt -> Vec4 -> Vec4++translate4X t (Vec4 x y z w) = Vec4 (x+t) y z w +translate4Y t (Vec4 x y z w) = Vec4 x (y+t) z w +translate4Z t (Vec4 x y z w) = Vec4 x y (z+t) w+translate4W t (Vec4 x y z w) = Vec4 x y z (w+t) ++vec4X :: Vec4+vec4Y :: Vec4+vec4Z :: Vec4+vec4W :: Vec4++vec4X = Vec4 1 0 0 0+vec4Y = Vec4 0 1 0 0+vec4Z = Vec4 0 0 1 0+vec4W = Vec4 0 0 0 1++---------------------------------------------------------------------------++-- |If @(x,y,u,v)@ is an orthonormal system, then (written in pseudo-code)+-- @biVector4 (x,y) = plusMinus (reverse $ biVector4 (u,v))@.+-- This is a helper function for the 4 dimensional rotation code.+-- If @(x,y,z,p,q,r) = biVector4 a b@, then the corresponding antisymmetric tensor is+--+-- > [  0  r  q  p ]+-- > [ -r  0  z -y ]+-- > [ -q -z  0  x ]+-- > [ -p  y -x  0 ]+biVector4 :: Vec4 -> Vec4 -> (Flt,Flt,Flt,Flt,Flt,Flt)+biVector4 (Vec4 x y z w) (Vec4 a b c d) = +  ( x*b-y*a , x*c-z*a , x*d-w*a , y*c-z*b , -y*d+w*b , z*d-w*c )++-- | the corresponding antisymmetric tensor+biVector4AsTensor :: Vec4 -> Vec4 -> Mat4+biVector4AsTensor v w = +  Mat4 ( Vec4   0  ( r) ( q) ( p) )+       ( Vec4 (-r)   0  ( z) (-y) )+       ( Vec4 (-q) (-z)   0  ( x) )+       ( Vec4 (-p) ( y) (-x)   0  )+  where +    (x,y,z,p,q,r) = biVector4 v w++-- | We assume that the axes are normalized and /orthogonal/ to each other!+rotate4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Vec4 -> Vec4+rotate4' angle axes v = v .* (rotMatrix4' angle axes)++-- | We assume only that the axes are independent vectors.+rotate4 :: Flt -> (Vec4,Vec4) -> Vec4 -> Vec4+rotate4 angle axes v = v .* (rotMatrix4 angle axes)++-- | Rotation matrix around a plane specified by two normalized and /orthogonal/ vectors.+-- Intended for multiplication on the /right/!+rotMatrix4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Mat4+rotMatrix4' angle (Normal4 v, Normal4 w) = m1 &+ (s *& m2) &+ m3 +  where+    c = cos angle ; s = sin angle+    m1 = scalarMul (1-c) ( outer v v  &+  outer w w )+    m2 = biVector4AsTensor v w+    m3 = diag (Vec4 c c c c)++-- | We assume only that the axes are independent vectors.+rotMatrix4 :: Flt -> (Vec4,Vec4) -> Mat4  +rotMatrix4 angle axes = +  rotMatrix4' angle $ liftPair toNormalUnsafe $ gramSchmidtNormalize axes +  where +    liftPair f (x,y) = (f x, f y)+    +
+ Data/Vect/Double/Util/Projective.hs view
@@ -0,0 +1,88 @@+{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-}++-- | Classic 4x4 projective matrices. Our convention is that they are intended for multiplication on+-- the /right/, that is, they are of the form+--+-- >     _____+-- > [  |     |  0  ]+-- > [  | 3x3 |  0  ]+-- > [  |_____|  0  ]+-- > [  p  q  r  1  ]+--+-- Please note that by default, OpenGL stores the matrices (in memory) by columns, while we +-- store them by rows; but OpenGL also use the opposite convention (so the OpenGL projective matrices +-- are intended for multiplication on the /left/). So in effect, they are the same when stored in the memory,+-- say with @poke :: Ptr Mat4 -> Mat4 -> IO ()@.++module Data.Vect.Flt.Util.Projective where++import Data.Vect.Flt.Base+import Data.Vect.Flt.Util.Dim3++import qualified Data.Vect.Flt.Util.Dim4 as Dim4++class ExtendProjective v e | v->e where+  extendProj     :: v -> e+  extendProjWith :: Flt -> v -> e+  extendProj = extendProjWith 1+  +instance ExtendProjective Vec2 Vec4 where+  extendProj       (Vec2 x y) = Vec4 x y 0 1+  extendProjWith w (Vec2 x y) = Vec4 x y 0 w+  +instance ExtendProjective Vec3 Vec4 where+  extendProj       (Vec3 x y z) = Vec4 x y z 1+  extendProjWith w (Vec3 x y z) = Vec4 x y z w++instance ExtendProjective Vec4 Vec4 where+  extendProj = id+  extendProjWith w (Vec4 x y z w') = let s = w/w' in Vec4 (s*x) (s*y) (s*z) w++instance ExtendProjective Mat2 Mat4 where+  extendProj       (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 1)+  extendProjWith w (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 w)++instance ExtendProjective Mat3 Mat4 where+  extendProj       (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 1)+  extendProjWith w (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 w)++rotMatrixProj :: Flt -> Normal3 -> Mat4+rotMatrixProj angle axis = extendProj $ rotMatrix3' axis angle++rotMatrixProj' :: {- ' CPP is sensitive to primes -} Flt -> Vec3 -> Mat4+rotMatrixProj' angle axis = extendProj $ rotMatrix3 axis angle++translMatrixProj :: Vec3 -> Mat4+translMatrixProj v = Mat4 Dim4.vec4X Dim4.vec4Y Dim4.vec4Z (extendProj v)++-- | we assume that the bottom-right corner is 1.+translWithProj :: Vec3 -> Mat4 -> Mat4+translWithProj v mat@(Mat4 r1 r2 r3 r4) = Mat4 r1 r2 r3 (extendProjWith 0 v &+ r4)++scaleMatrixProj :: Vec3 -> Mat4+scaleMatrixProj v = diag $ extendProj v++scaleMatrixUniformProj :: Flt -> Mat4+scaleMatrixUniformProj s = diag (Vec4 s s s 1)++class ProjectiveAction v where+  actProj :: v -> Mat4 -> v+ +instance ProjectiveAction Vec3 where+  actProj v m = trim $ (extendProj v) .* m ++instance ProjectiveAction Vec4 where+  actProj v m = v .* m ++-- | When acting on unit vectors, we ignore the translation part.+instance ProjectiveAction Normal3 where+  actProj (Normal3 v) m = Normal3 (v .* (trim m :: Mat3))++-- | Inverts a projective 4x4 matrix, assuming that the top-left 3x3 part is /orthogonal/,+-- and the bottom-right corner is 1.+invertProj :: Mat4 -> Mat4+invertProj mat@(Mat4 u v w t) = +  translWithProj t' $ extendProj $ transpose $ (trim mat :: Mat3)+  where+    t' = Vec3 (- u &. t) (- v &. t) (- w &. t)+
+ Data/Vect/Float.hs view
@@ -0,0 +1,23 @@+{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-}++module Data.Vect.Flt+  ( module Data.Vect.Flt.Base+  , module Data.Vect.Flt.Interpolate+  , module Data.Vect.Flt.Util.Dim2+  , module Data.Vect.Flt.Util.Dim3+  , module Data.Vect.Flt.Util.Projective+#ifdef VECT_OPENGL        +  , module Data.Vect.Flt.OpenGL       +#endif+  ) where++import Data.Vect.Flt.Base+import Data.Vect.Flt.Interpolate++import Data.Vect.Flt.Util.Dim2+import Data.Vect.Flt.Util.Dim3+import Data.Vect.Flt.Util.Projective++#ifdef VECT_OPENGL         +import Data.Vect.Flt.OpenGL       +#endif
+ Data/Vect/Float/Base.hs view
@@ -0,0 +1,735 @@+{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-}++module Data.Vect.Flt.Base where++import Control.Monad+import System.Random  +import Foreign++-- class declarations++class AbelianGroup g where+  (&+) :: g -> g -> g+  (&-) :: g -> g -> g+  neg  :: g -> g+  zero :: g++infixl 6 &++infixl 6 &- ++vecSum :: AbelianGroup g => [g] -> g+vecSum l = foldl (&+) zero l ++class (AbelianGroup r) => +      Ring r where+  (.*.) :: r -> r -> r+  one   :: r++infixl 7 .*. ++ringProduct :: Ring r => [r] -> r+ringProduct l = foldl (.*.) one l++class LeftModule r m where+  lmul :: r -> m -> m+  (*.) :: r -> m -> m+  (*.) = lmul++class RightModule m r where+  rmul :: m -> r -> m+  (.*) :: m -> r -> m+  (.*) = rmul++-- I'm not really sure about this.. may actually degrade the performance in some cases?  +{- RULES+"matrix multiplication left"   forall m n x.  (n .*. m) *. x = n *. (m *. x)  +"matrix multiplication right"  forall m n x.  x .* (m .*. n) = (x .* m) .* n+  -}++infixr 7 *.+infixl 7 .*++class AbelianGroup v => Vector v where+  mapVec    :: (Flt -> Flt) -> v -> v+  scalarMul :: Flt -> v -> v+  (*&) ::      Flt -> v -> v +  (&*) ::      v -> Flt -> v +  (*&) s v = scalarMul s v+  (&*) v s = scalarMul s v++infixr 7 *&+infixl 7 &*++{-# RULES+"scalar multiplication left"   forall s t x.  t *& (s *& x) = (t*s) *& x +"scalar multiplication right"  forall s t x.  (x &* s) &* t = x &* (s*t)  +  #-}++class DotProd v where+  (&.) :: v -> v -> Flt+  norm    :: v -> Flt+  normsqr :: v -> Flt+  len     :: v -> Flt+  lensqr  :: v -> Flt+  len = norm+  lensqr = normsqr+  dotprod :: v -> v -> Flt+  normsqr v = (v &. v)  +  norm = sqrt.lensqr+  dotprod = (&.)++infix 7 &.++{-# RULES+"len/square 1"   forall x.  (len x)*(len x) = lensqr x+"len/square 2"   forall x.  (len x)^2 = lensqr x+"norm/square 1"  forall x.  (norm x)*(norm x) = normsqr x+"norm/square 2"  forall x.  (norm x)^2 = normsqr x+  #-}++normalize :: (Vector v, DotProd v) => v -> v+normalize v = scalarMul (1.0/(len v)) v++distance :: (Vector v, DotProd v) => v -> v -> Flt+distance x y = norm (x &- y)++-- | the angle between two vectors+angle :: (Vector v, DotProd v) => v -> v -> Flt +angle x y = acos $ (x &. y) / (norm x * norm y)++-- | the angle between two unit vectors+angle' {- ' CPP is sensitive to primes -} :: (Vector v, UnitVector v u, DotProd v) => u -> u -> Flt +angle' x y = acos (fromNormal x &. fromNormal y)++{-# RULES+"normalize is idempotent"  forall x. normalize (normalize x) = normalize x+  #-}++class (Vector v, DotProd v) => UnitVector v u | v->u, u->v  where+  mkNormal         :: v -> u       -- ^ normalizes the input+  toNormalUnsafe   :: v -> u       -- ^ does not normalize the input!+  fromNormal       :: u -> v+  fromNormalRadius :: Flt -> u -> v+  fromNormalRadius t n = t *& fromNormal n ++-- | projects the first vector onto the direction of the second (unit) vector+project' :: (Vector v, UnitVector v u, DotProd v) => v -> u -> v+project' what dir = projectUnsafe what (fromNormal dir)++-- | direction (second argument) is assumed to be a /unit/ vector!+projectUnsafe :: (Vector v, DotProd v) => v -> v -> v+projectUnsafe what dir = what &- dir &* (what &. dir)++project :: (Vector v, DotProd v) => v -> v -> v+project what dir = what &- dir &* ((what &. dir) / (dir &. dir))++-- | since unit vectors are not a group, we need a separate function.+flipNormal :: UnitVector v n => n -> n +flipNormal = toNormalUnsafe . neg . fromNormal ++class CrossProd v where+  crossprod :: v -> v -> v+  (&^)      :: v -> v -> v+  (&^) = crossprod+  +class Pointwise v where+  pointwise :: v -> v -> v+  (&!)      :: v -> v -> v+  (&!) = pointwise ++infix 7 &^+infix 7 &!++class HasCoordinates v x | v->x where+  _1 :: v -> x+  _2 :: v -> x+  _3 :: v -> x+  _4 :: v -> x      ++-- | conversion between vectors (and matrices) of different dimensions+class Extend u v where+  extendZero :: u -> v          -- ^ example: @extendZero (Vec2 5 6) = Vec4 5 6 0 0@+  extendWith :: Flt -> u -> v   -- ^ example: @extendWith 1 (Vec2 5 6) = Vec4 5 6 1 1@+  trim :: v -> u                -- ^ example: @trim (Vec4 5 6 7 8) = Vec2 5 6@++-- | makes a diagonal matrix from a vector+class Diagonal s t | t->s where+  diag :: s -> t++class Matrix m where+  transpose :: m -> m +  inverse :: m -> m+  idmtx :: m++{-# RULES+"transpose is an involution"  forall m. transpose (transpose m) = m+"inverse is an involution"    forall m. inverse (inverse m) = m+  #-}+  +-- | Outer product (could be unified with Diagonal?)+class Tensor t v | t->v where+  outer :: v -> v -> t+    +class Determinant m where+  det :: m -> Flt    +    +-- Vec / Mat datatypes+ +data Vec2 = Vec2 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt +  deriving (Read,Show)+data Vec3 = Vec3 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt +  deriving (Read,Show)+data Vec4 = Vec4 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt +  deriving (Read,Show)++-- | these are /row/ vectors +data Mat2 = Mat2 !Vec2 !Vec2              deriving (Read,Show)+data Mat3 = Mat3 !Vec3 !Vec3 !Vec3        deriving (Read,Show)+data Mat4 = Mat4 !Vec4 !Vec4 !Vec4 !Vec4  deriving (Read,Show)++-- | The assumption when dealing with these is always that they are of unit length.+-- Also, interpolation works differently.+newtype Normal2 = Normal2 Vec2 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) +newtype Normal3 = Normal3 Vec3 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable,CrossProd) +newtype Normal4 = Normal4 Vec4 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) ++mkVec2 :: (Flt,Flt) -> Vec2+mkVec3 :: (Flt,Flt,Flt) -> Vec3+mkVec4 :: (Flt,Flt,Flt,Flt) -> Vec4++mkVec2 (x,y)     = Vec2 x y +mkVec3 (x,y,z)   = Vec3 x y z+mkVec4 (x,y,z,w) = Vec4 x y z w++-- Unit vectors+  +instance UnitVector Vec2 Normal2 where+  mkNormal v = Normal2 (normalize v)+  fromNormal (Normal2 v) = v +  toNormalUnsafe = Normal2++instance UnitVector Vec3 Normal3 where+  mkNormal v = Normal3 (normalize v)+  fromNormal (Normal3 v) = v +  toNormalUnsafe = Normal3++instance UnitVector Vec4 Normal4 where+  mkNormal v = Normal4 (normalize v)+  fromNormal (Normal4 v) = v +  toNormalUnsafe = Normal4++rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)+rndUnit g = +  if d > 0.01+    then ( v &* (1.0/d) , h )+    else rndUnit h+  where+    (v,h) = random g+    d = norm v+    +instance Random Normal2 where+  random g = let (v,h) = rndUnit g in (Normal2 v, h)  +  randomR _ = random++instance Random Normal3 where+  random g = let (v,h) = rndUnit g in (Normal3 v, h)  +  randomR _ = random++instance Random Normal4 where+  random g = let (v,h) = rndUnit g in (Normal4 v, h)  +  randomR _ = random++{-+instance Storable Normal2 where+  alignment _ = alignment (undefined::Vec2)+  sizeOf    _ = sizeOf    (undefined::Vec2)+  peek p = liftM (\v -> Normal2 v) (peek $ castPtr p)+  poke p (Normal2 v) = poke (castPtr p) v  + +instance Storable Normal3 where+  alignment _ = alignment (undefined::Vec3)+  sizeOf    _ = sizeOf    (undefined::Vec3)+  peek p = liftM (\v -> Normal3 v) (peek $ castPtr p)+  poke p (Normal3 v) = poke (castPtr p) v  ++instance Storable Normal4 where+  alignment _ = alignment (undefined::Vec4)+  sizeOf    _ = sizeOf    (undefined::Vec4)+  peek p = liftM (\v -> Normal4 v) (peek $ castPtr p)+  poke p (Normal4 v) = poke (castPtr p) v  +-}++-- Vec2 instances++instance HasCoordinates Vec2 Flt where+  _1 (Vec2 x _) = x+  _2 (Vec2 _ y) = y+  _3 _ = error "has only 2 coordinates"+  _4 _ = error "has only 2 coordinates"++instance AbelianGroup Vec2 where+  (&+) (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1+x2) (y1+y2) +  (&-) (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1-x2) (y1-y2)+  neg  (Vec2 x y)                = Vec2 (-x) (-y)+  zero = Vec2 0 0+  +instance Vector Vec2 where+  scalarMul s (Vec2 x y) = Vec2 (s*x) (s*y)+  mapVec    f (Vec2 x y) = Vec2 (f x) (f y)+  +instance DotProd Vec2 where+  (&.) (Vec2 x1 y1) (Vec2 x2 y2) = x1*x2 + y1*y2++instance Pointwise Vec2 where+  pointwise (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1*x2) (y1*y2)++instance Determinant (Vec2,Vec2) where+  det (Vec2 x1 y1 , Vec2 x2 y2) = x1*y2 - x2*y1  ++{-     +instance Show Vec2 where+  show (Vec2 x y) = "( " ++ show x ++ " , " ++ show y ++ " )"+-}++instance Random Vec2 where+  random = randomR (Vec2 (-1) (-1),Vec2 1 1)+  randomR (Vec2 a b, Vec2 c d) gen = +    let (x,gen1) = randomR (a,c) gen+        (y,gen2) = randomR (b,d) gen1+    in (Vec2 x y, gen2)+     +instance Storable Vec2 where+  sizeOf    _ = 2 * sizeOf (undefined::Flt)+  alignment _ = sizeOf (undefined::Flt)+  +  peek q = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    x <- peek        p +    y <- peekByteOff p k+    return (Vec2 x y)+    +  poke q (Vec2 x y) = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    poke        p   x+    pokeByteOff p k y+               +-- Mat2 instances++instance HasCoordinates Mat2 Vec2 where+  _1 (Mat2 x _) = x+  _2 (Mat2 _ y) = y+  _3 _ = error "has only 2 coordinates"+  _4 _ = error "has only 2 coordinates"++instance Matrix Mat2 where+  transpose (Mat2 row1 row2) = +    Mat2 (Vec2 (_1 row1) (_1 row2)) +         (Vec2 (_2 row1) (_2 row2)) +  idmtx = Mat2 (Vec2 1 0) (Vec2 0 1)+  inverse (Mat2 (Vec2 a b) (Vec2 c d)) = +    Mat2 (Vec2 (d*r) (-b*r)) (Vec2 (-c*r) (a*r)) +    where r = 1.0 / (a*d - b*c)++instance AbelianGroup Mat2 where+  (&+) (Mat2 r1 r2) (Mat2 s1 s2) = Mat2 (r1 &+ s1) (r2 &+ s2)+  (&-) (Mat2 r1 r2) (Mat2 s1 s2) = Mat2 (r1 &- s1) (r2 &- s2)+  neg  (Mat2 r1 r2)              = Mat2 (neg r1) (neg r2)  +  zero = Mat2 zero zero   -- (zero::Vec2) (zero::Vec2)++instance Vector Mat2 where+  scalarMul s (Mat2 r1 r2) = Mat2 (g r1) (g r2) where g = scalarMul s+  mapVec    f (Mat2 r1 r2) = Mat2 (g r1) (g r2) where g = mapVec f++instance Ring Mat2 where+  (.*.) (Mat2 r1 r2) n = +    let (Mat2 c1 c2) = transpose n+    in Mat2 (Vec2 (r1 &. c1) (r1 &. c2))+            (Vec2 (r2 &. c1) (r2 &. c2))+  one = idmtx ++instance LeftModule Mat2 Vec2 where+  lmul (Mat2 row1 row2) v = Vec2 (row1 &. v) (row2 &. v) +  +instance RightModule Vec2 Mat2 where+  rmul v mt = lmul (transpose mt) v++instance Diagonal Vec2 Mat2 where+  diag (Vec2 x y) = Mat2 (Vec2 x 0) (Vec2 0 y)++instance Tensor Mat2 Vec2 where+  outer (Vec2 a b) (Vec2 x y) = Mat2+    (Vec2 (a*x) (a*y))+    (Vec2 (b*x) (b*y))+{-+  outer v w = +    let full = Mat2 (Vec2 1 1) (Vec2 1 1)+    in  (diag v) .*. full .*. (diag w)+-}++instance Determinant Mat2 where+  det (Mat2 (Vec2 a b) (Vec2 c d)) = a*d - b*c ++{-+instance Show Mat2 where+  show (Mat2 r1 r2) = show r1 ++ "\n" ++ show r2+-}++instance Storable Mat2 where+  sizeOf    _ = 2 * sizeOf (undefined::Vec2)+  alignment _ = alignment  (undefined::Vec2)+  +  peek q = do+    let p = castPtr q :: Ptr Vec2+        k = sizeOf (undefined::Vec2)+    r1 <- peek        p +    r2 <- peekByteOff p k+    return (Mat2 r1 r2)+    +  poke q (Mat2 r1 r2) = do+    let p = castPtr q :: Ptr Vec2+        k = sizeOf (undefined::Vec2)+    poke        p   r1+    pokeByteOff p k r2++-- Vec3 instances++instance HasCoordinates Vec3 Flt where+  _1 (Vec3 x _ _) = x+  _2 (Vec3 _ y _) = y+  _3 (Vec3 _ _ z) = z+  _4 _ = error "has only 3 coordinates"++instance AbelianGroup Vec3 where+  (&+) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1+x2) (y1+y2) (z1+z2) +  (&-) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1-x2) (y1-y2) (z1-z2) +  neg  (Vec3 x y z)                    = Vec3 (-x) (-y) (-z)+  zero = Vec3 0 0 0+  +instance Vector Vec3 where+  scalarMul s (Vec3 x y z) = Vec3 (s*x) (s*y) (s*z)+  mapVec    f (Vec3 x y z) = Vec3 (f x) (f y) (f z)++instance DotProd Vec3 where+  (&.) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = x1*x2 + y1*y2 + z1*z2++instance Pointwise Vec3 where+  pointwise (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1*x2) (y1*y2) (z1*z2)++{-+instance Show Vec3 where+  show (Vec3 x y z) = "( " ++ show x ++ " , " ++ show y ++ " , " ++ show z ++ " )"+-}++instance Random Vec3 where+  random = randomR (Vec3 (-1) (-1) (-1),Vec3 1 1 1)+  randomR (Vec3 a b c, Vec3 d e f) gen = +    let (x,gen1) = randomR (a,d) gen+        (y,gen2) = randomR (b,e) gen1+        (z,gen3) = randomR (c,f) gen2  +    in (Vec3 x y z, gen3)+      +instance CrossProd Vec3 where+  crossprod (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (y1*z2-y2*z1) (z1*x2-z2*x1) (x1*y2-x2*y1) ++instance Determinant (Vec3,Vec3,Vec3) where+  det (u,v,w) = u &. (v &^ w)  + +instance Storable Vec3 where+  sizeOf    _ = 3 * sizeOf (undefined::Flt)+  alignment _ = sizeOf (undefined::Flt)+  +  peek q = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    x <- peek        p +    y <- peekByteOff p (k  )+    z <- peekByteOff p (k+k)+    return (Vec3 x y z)+    +  poke q (Vec3 x y z) = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    poke        p       x+    pokeByteOff p (k  ) y+    pokeByteOff p (k+k) z+   +-- Mat3 instances++instance HasCoordinates Mat3 Vec3 where+  _1 (Mat3 x _ _) = x+  _2 (Mat3 _ y _) = y+  _3 (Mat3 _ _ z) = z+  _4 _ = error "has only 3 coordinates"  ++instance Matrix Mat3 where++  transpose (Mat3 row1 row2 row3) = +    Mat3 (Vec3 (_1 row1) (_1 row2) (_1 row3)) +         (Vec3 (_2 row1) (_2 row2) (_2 row3)) +         (Vec3 (_3 row1) (_3 row2) (_3 row3)) +         +  idmtx = Mat3 (Vec3 1 0 0) (Vec3 0 1 0) (Vec3 0 0 1)+  +  inverse (Mat3 (Vec3 a b c) (Vec3 e f g) (Vec3 i j k)) = +    Mat3 (Vec3 (d11*r) (d21*r) (d31*r))  +         (Vec3 (d12*r) (d22*r) (d32*r))  +         (Vec3 (d13*r) (d23*r) (d33*r))  +    where+      r = 1.0 / ( a*d11 + b*d12 + c*d13 )++      d11 = f*k - g*j+      d12 = g*i - e*k+      d13 = e*j - f*i++      d31 = b*g - c*f+      d32 = c*e - a*g+      d33 = a*f - b*e++      d21 = c*j - b*k +      d22 = a*k - c*i +      d23 = b*i - a*j ++instance AbelianGroup Mat3 where+  (&+) (Mat3 r1 r2 r3) (Mat3 s1 s2 s3) = Mat3 (r1 &+ s1) (r2 &+ s2) (r3 &+ s3)+  (&-) (Mat3 r1 r2 r3) (Mat3 s1 s2 s3) = Mat3 (r1 &- s1) (r2 &- s2) (r3 &- s3)+  neg  (Mat3 r1 r2 r3)                 = Mat3 (neg r1) (neg r2) (neg r3) +  zero = Mat3 zero zero zero   -- (zero::Vec3) (zero::Vec3) (zero::Vec3)++instance Vector Mat3 where+  scalarMul s (Mat3 r1 r2 r3) = Mat3 (g r1) (g r2) (g r3) where g = scalarMul s+  mapVec    f (Mat3 r1 r2 r3) = Mat3 (g r1) (g r2) (g r3) where g = mapVec f++instance Ring Mat3 where+  (.*.) (Mat3 r1 r2 r3) n = +    let (Mat3 c1 c2 c3) = transpose n+    in Mat3 (Vec3 (r1 &. c1) (r1 &. c2) (r1 &. c3))+            (Vec3 (r2 &. c1) (r2 &. c2) (r2 &. c3))+            (Vec3 (r3 &. c1) (r3 &. c2) (r3 &. c3))+  one = idmtx ++instance LeftModule Mat3 Vec3 where+  lmul (Mat3 row1 row2 row3) v = Vec3 (row1 &. v) (row2 &. v) (row3 &. v)+  +instance RightModule Vec3 Mat3 where+  rmul v mt = lmul (transpose mt) v++instance Diagonal Vec3 Mat3 where+  diag (Vec3 x y z) = Mat3 (Vec3 x 0 0) (Vec3 0 y 0) (Vec3 0 0 z)++instance Tensor Mat3 Vec3 where+  outer (Vec3 a b c) (Vec3 x y z) = Mat3+    (Vec3 (a*x) (a*y) (a*z))+    (Vec3 (b*x) (b*y) (b*z))+    (Vec3 (c*x) (c*y) (c*z))+{-+  outer v w = +    let full = Mat3 (Vec3 1 1 1) (Vec3 1 1 1) (Vec3 1 1 1)+    in  (diag v) .*. full .*. (diag w)+-}++instance Determinant Mat3 where+  det (Mat3 r1 r2 r3) = det (r1,r2,r3)++{-+instance Show Mat3 where+  show (Mat3 r1 r2 r3) = show r1 ++ "\n" ++ show r2 ++ "\n" ++ show r3+-}++instance Storable Mat3 where+  sizeOf    _ = 3 * sizeOf (undefined::Vec3)+  alignment _ = alignment  (undefined::Vec3)+  +  peek q = do+    let p = castPtr q :: Ptr Vec3+        k = sizeOf (undefined::Vec3)+    r1 <- peek        p +    r2 <- peekByteOff p (k  )+    r3 <- peekByteOff p (k+k)+    return (Mat3 r1 r2 r3)+    +  poke q (Mat3 r1 r2 r3) = do+    let p = castPtr q :: Ptr Vec3+        k = sizeOf (undefined::Vec3)+    poke        p       r1+    pokeByteOff p (k  ) r2+    pokeByteOff p (k+k) r3++-- Vec4 instances++instance HasCoordinates Vec4 Flt where+  _1 (Vec4 x _ _ _) = x+  _2 (Vec4 _ y _ _) = y+  _3 (Vec4 _ _ z _) = z+  _4 (Vec4 _ _ _ w) = w++instance AbelianGroup Vec4 where+  (&+) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1+x2) (y1+y2) (z1+z2) (w1+w2)+  (&-) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1-x2) (y1-y2) (z1-z2) (w1-w2)+  neg  (Vec4 x y z w)                        = Vec4 (-x) (-y) (-z) (-w)+  zero = Vec4 0 0 0 0+  +instance Vector Vec4 where+  scalarMul s (Vec4 x y z w) = Vec4 (s*x) (s*y) (s*z) (s*w)+  mapVec    f (Vec4 x y z w) = Vec4 (f x) (f y) (f z) (f w)++instance DotProd Vec4 where+  (&.) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = x1*x2 + y1*y2 + z1*z2 + w1*w2++instance Pointwise Vec4 where+  pointwise (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1*x2) (y1*y2) (z1*z2) (w1*w2)++{-+instance Show Vec4 where+  show (Vec4 x y z w) = "( " ++ show x ++ " , " ++ show y ++ " , " ++ show z ++ " , " ++ show w ++ " )"+-}++instance Random Vec4 where+  random = randomR (Vec4 (-1) (-1) (-1) (-1),Vec4 1 1 1 1)+  randomR (Vec4 a b c d, Vec4 e f g h) gen = +    let (x,gen1) = randomR (a,e) gen+        (y,gen2) = randomR (b,f) gen1+        (z,gen3) = randomR (c,g) gen2  +        (w,gen4) = randomR (d,h) gen3  +    in (Vec4 x y z w, gen4)+           +instance Storable Vec4 where+  sizeOf    _ = 4 * sizeOf (undefined::Flt)+  alignment _ = sizeOf (undefined::Flt)+  +  peek q = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    x <- peek        p +    y <- peekByteOff p (k  )+    z <- peekByteOff p (k+k)+    w <- peekByteOff p (3*k)+    return (Vec4 x y z w)+    +  poke q (Vec4 x y z w) = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    poke        p       x+    pokeByteOff p (k  ) y+    pokeByteOff p (k+k) z+    pokeByteOff p (3*k) w++-- Mat4 instances++instance HasCoordinates Mat4 Vec4 where+  _1 (Mat4 x _ _ _) = x+  _2 (Mat4 _ y _ _) = y+  _3 (Mat4 _ _ z _) = z+  _4 (Mat4 _ _ _ w) = w++instance Matrix Mat4 where+  transpose (Mat4 row1 row2 row3 row4) = +    Mat4 (Vec4 (_1 row1) (_1 row2) (_1 row3) (_1 row4)) +         (Vec4 (_2 row1) (_2 row2) (_2 row3) (_2 row4)) +         (Vec4 (_3 row1) (_3 row2) (_3 row3) (_3 row4)) +         (Vec4 (_4 row1) (_4 row2) (_4 row3) (_4 row4)) +  idmtx = Mat4 (Vec4 1 0 0 0) (Vec4 0 1 0 0) (Vec4 0 0 1 0) (Vec4 0 0 0 1)+  inverse = error "inverse/Mat4: not implemented yet"++instance AbelianGroup Mat4 where+  (&+) (Mat4 r1 r2 r3 r4) (Mat4 s1 s2 s3 s4) = Mat4 (r1 &+ s1) (r2 &+ s2) (r3 &+ s3) (r4 &+ s4)+  (&-) (Mat4 r1 r2 r3 r4) (Mat4 s1 s2 s3 s4) = Mat4 (r1 &- s1) (r2 &- s2) (r3 &- s3) (r4 &- s4)+  neg  (Mat4 r1 r2 r3 r4)                    = Mat4 (neg r1) (neg r2) (neg r3) (neg r4) +  zero = Mat4 zero zero zero zero+  +instance Vector Mat4 where+  scalarMul s (Mat4 r1 r2 r3 r4) = Mat4 (g r1) (g r2) (g r3) (g r4) where g = scalarMul s+  mapVec    f (Mat4 r1 r2 r3 r4) = Mat4 (g r1) (g r2) (g r3) (g r4) where g = mapVec f++instance Ring Mat4 where+  (.*.) (Mat4 r1 r2 r3 r4) n = +    let (Mat4 c1 c2 c3 c4) = transpose n+    in Mat4 (Vec4 (r1 &. c1) (r1 &. c2) (r1 &. c3) (r1 &. c4))+            (Vec4 (r2 &. c1) (r2 &. c2) (r2 &. c3) (r2 &. c4))+            (Vec4 (r3 &. c1) (r3 &. c2) (r3 &. c3) (r3 &. c4))+            (Vec4 (r4 &. c1) (r4 &. c2) (r4 &. c3) (r4 &. c4))+  one = idmtx ++instance LeftModule Mat4 Vec4 where+  lmul (Mat4 row1 row2 row3 row4) v = Vec4 (row1 &. v) (row2 &. v) (row3 &. v) (row4 &. v)+  +instance RightModule Vec4 Mat4 where+  rmul v mt = lmul (transpose mt) v++instance Diagonal Vec4 Mat4 where+  diag (Vec4 x y z w) = Mat4 (Vec4 x 0 0 0) (Vec4 0 y 0 0) (Vec4 0 0 z 0) (Vec4 0 0 0 w)++instance Tensor Mat4 Vec4 where+  outer (Vec4 a b c d) (Vec4 x y z w) = Mat4+    (Vec4 (a*x) (a*y) (a*z) (a*w))+    (Vec4 (b*x) (b*y) (b*z) (b*w))+    (Vec4 (c*x) (c*y) (c*z) (c*w))+    (Vec4 (d*x) (d*y) (d*z) (d*w))+{-+  outer v w = +    let full = Mat4 (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1)+    in  (diag v) .*. full .*. (diag w)+-}++--instance Determinant Mat4 where+--  det (Mat4 r1 r2 r3 r4)++{-+instance Show Mat4 where+  show (Mat4 r1 r2 r3 r4) = show r1 ++ "\n" ++ show r2 ++ "\n" ++ show r3 ++ "\n" ++ show r4+-}++instance Storable Mat4 where+  sizeOf    _ = 4 * sizeOf (undefined::Vec4)+  alignment _ = alignment  (undefined::Vec4)+  +  peek q = do+    let p = castPtr q :: Ptr Vec4+        k = sizeOf (undefined::Vec4)+    r1 <- peek        p +    r2 <- peekByteOff p (k  )+    r3 <- peekByteOff p (k+k)+    r4 <- peekByteOff p (3*k)+    return (Mat4 r1 r2 r3 r4)+    +  poke q (Mat4 r1 r2 r3 r4) = do+    let p = castPtr q :: Ptr Vec4+        k = sizeOf (undefined::Vec4)+    poke        p       r1+    pokeByteOff p (k  ) r2+    pokeByteOff p (k+k) r3+    pokeByteOff p (3*k) r4++-- Extend instances++instance Extend Vec2 Vec3 where+  extendZero   (Vec2 x y) = Vec3 x y 0+  extendWith t (Vec2 x y) = Vec3 x y t+  trim (Vec3 x y _)       = Vec2 x y++instance Extend Vec2 Vec4 where+  extendZero   (Vec2 x y) = Vec4 x y 0 0+  extendWith t (Vec2 x y) = Vec4 x y t t+  trim (Vec4 x y _ _)     = Vec2 x y ++instance Extend Vec3 Vec4 where+  extendZero   (Vec3 x y z) = Vec4 x y z 0+  extendWith t (Vec3 x y z) = Vec4 x y z t+  trim (Vec4 x y z _)       = Vec3 x y z++instance Extend Mat2 Mat3 where+  extendZero (Mat2 p q) = Mat3 (extendZero p) (extendZero q) zero+  extendWith _ _ = error "extendWith is meaningless for matrices"+  trim (Mat3 p q _) = Mat2 (trim p) (trim q)++instance Extend Mat2 Mat4 where+  extendZero (Mat2 p q) = Mat4 (extendZero p) (extendZero q) zero zero+  extendWith _ _ = error "extendWith is meaningless for matrices"+  trim (Mat4 p q _ _) = Mat2 (trim p) (trim q)++instance Extend Mat3 Mat4 where+  extendZero (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) zero+  extendWith _ _ = error "extendWith is meaningless for matrices"+  trim (Mat4 p q r _) = Mat3 (trim p) (trim q) (trim r)+
+ Data/Vect/Float/GramSchmidt.hs view
@@ -0,0 +1,135 @@+{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-}++-- | Gram-Schmidt orthogonalization.+-- This module is not re-exported by "Data.Vect".++module Data.Vect.Flt.GramSchmidt +  ( GramSchmidt(..)+  )+  where++import Data.Vect.Flt.Base++-------------------------------------------------------++liftPair :: (a -> b) -> (a,a) -> (b,b)+liftPair f (x,y) = (f x, f y)++liftTriple :: (a -> b) -> (a,a,a) -> (b,b,b)+liftTriple f (x,y,z) = (f x, f y, f z)++liftQuadruple :: (a -> b) -> (a,a,a,a) -> (b,b,b,b)+liftQuadruple f (x,y,z,w) = (f x, f y, f z, f w)++-------------------------------------------------------+    +-- | produces orthogonal\/orthonormal vectors from a set of vectors    +class GramSchmidt a where+  gramSchmidt          :: a -> a   -- ^ does not normalize the vectors!+  gramSchmidtNormalize :: a -> a   -- ^ normalizes the vectors.++{-# RULES+"gramSchmidt is idempotent"  forall a. gramSchmidt (gramSchmidt a) = gramSchmidt a +"gramSchmidtNormalize is idempotent"  forall a. gramSchmidtNormalize (gramSchmidtNormalize a) = gramSchmidtNormalize a +  #-}++-------------------------------------------------------++instance GramSchmidt (Vec2,Vec2) where+  gramSchmidt = gramSchmidtPair+  gramSchmidtNormalize = gramSchmidtNormalizePair+  +instance GramSchmidt (Vec3,Vec3) where+  gramSchmidt = gramSchmidtPair+  gramSchmidtNormalize = gramSchmidtNormalizePair+  +instance GramSchmidt (Vec4,Vec4) where+  gramSchmidt = gramSchmidtPair+  gramSchmidtNormalize = gramSchmidtNormalizePair++----------++instance GramSchmidt (Normal2,Normal2) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal2!"+  gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal++instance GramSchmidt (Normal3,Normal3) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal3!"+  gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal++instance GramSchmidt (Normal4,Normal4) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal4!"+  gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal++----------+  +gramSchmidtPair :: (Vector v, DotProd v) => (v,v) -> (v,v)+gramSchmidtPair (u,v) = (u',v') where +  u' = u+  v' = project v u'     +  +gramSchmidtNormalizePair :: (Vector v, DotProd v) => (v,v) -> (v,v)+gramSchmidtNormalizePair (u,v) = (u',v') where+  u' = normalize u +  v' = normalize $ projectUnsafe v u'     ++----------++instance GramSchmidt (Vec3,Vec3,Vec3) where+  gramSchmidt = gramSchmidtTriple+  gramSchmidtNormalize = gramSchmidtNormalizeTriple+     +instance GramSchmidt (Vec4,Vec4,Vec4) where+  gramSchmidt = gramSchmidtTriple+  gramSchmidtNormalize = gramSchmidtNormalizeTriple++instance GramSchmidt (Normal3,Normal3,Normal3) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal3!"+  gramSchmidtNormalize = liftTriple toNormalUnsafe . gramSchmidtNormalizeTriple . liftTriple fromNormal++instance GramSchmidt (Normal4,Normal4,Normal4) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal4!"+  gramSchmidtNormalize = liftTriple toNormalUnsafe . gramSchmidtNormalizeTriple . liftTriple fromNormal++----------++gramSchmidtTriple :: (Vector v, DotProd v) => (v,v,v) -> (v,v,v)+gramSchmidtTriple (u,v,w) = (u',v',w') where +  u' = u+  v' = project v u'     +  w' = project (project w u') v' +  +gramSchmidtNormalizeTriple :: (Vector v, DotProd v) => (v,v,v) -> (v,v,v)+gramSchmidtNormalizeTriple (u,v,w) = (u',v',w') where+  u' = normalize $ u +  v' = normalize $ projectUnsafe v u'     +  w' = normalize $ projectUnsafe (projectUnsafe w u') v'     ++----------++instance GramSchmidt (Vec4,Vec4,Vec4,Vec4) where+  gramSchmidt          = gramSchmidtQuadruple+  gramSchmidtNormalize = gramSchmidtNormalizeQuadruple ++instance GramSchmidt (Normal4,Normal4,Normal4,Normal4) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal4!"+  gramSchmidtNormalize = liftQuadruple toNormalUnsafe . gramSchmidtNormalizeQuadruple . liftQuadruple fromNormal++----------+  +gramSchmidtQuadruple :: (Vector v, DotProd v) => (v,v,v,v) -> (v,v,v,v)+gramSchmidtQuadruple (u,v,w,z) = (u',v',w',z') where +  u' = u+  v' = project v u'     +  w' = project (project w u') v' +  z' = project (project (project z u') v') w'++gramSchmidtNormalizeQuadruple :: (Vector v, DotProd v) => (v,v,v,v) -> (v,v,v,v)+gramSchmidtNormalizeQuadruple (u,v,w,z) = (u',v',w',z') where+  u' = normalize $ u+  v' = normalize $ projectUnsafe v u'     +  w' = normalize $ projectUnsafe (projectUnsafe w u') v' +  z' = normalize $ projectUnsafe (projectUnsafe (projectUnsafe z u') v') w'+  +----------+
+ Data/Vect/Float/Interpolate.hs view
@@ -0,0 +1,43 @@+{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-}++-- | Interpolation of vectors. +-- Note: we interpolate unit vectors differently from ordinary vectors.++module Data.Vect.Flt.Interpolate where++import Data.Vect.Flt.Base+import Data.Vect.Flt.Util.Dim2 (sinCos',angle2')+import Data.Vect.Flt.Util.Dim3 (rotate3')++class Interpolate v where+  interpolate :: Flt -> v -> v -> v+  +instance Interpolate Flt where+  interpolate t x y = x + t*(y-x)++instance Interpolate Vec2 where interpolate t x y = x &+ t *& (y &- x)+instance Interpolate Vec3 where interpolate t x y = x &+ t *& (y &- x)+instance Interpolate Vec4 where interpolate t x y = x &+ t *& (y &- x)++instance Interpolate Normal2 where+  interpolate t nx ny = sinCos' $ ax + t*adiff where+    ax = angle2' nx+    ay = angle2' ny+    adiff = helper (ay - ax)+    helper d +      | d < -pi   = d + twopi+      | d >  pi   = d - twopi+      | otherwise = d+    twopi = 2*pi+    +instance Interpolate Normal3 where +  interpolate t nx ny = +    if maxAngle < 0.001  -- more or less ad-hoc critical angle+      then mkNormal $ interpolate t x y+      else toNormalUnsafe $ rotate3' (t*maxAngle) (mkNormal axis) x where+    x = fromNormal nx+    y = fromNormal ny+    axis = (x &^ y)+    maxAngle = acos (x &. y)+        +
+ Data/Vect/Float/OpenGL.hs view
@@ -0,0 +1,183 @@+{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-}++-- TODO: the pointer versions of these functions should be really implemented +-- via the pointer versions of the original opengl functions...++-- | OpenGL support, inclduing 'vertex', 'texCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.+ +module Data.Vect.Flt.OpenGL where++import Control.Monad+import Data.Vect.Flt.Base+import qualified Graphics.Rendering.OpenGL as GL++import Foreign++import Graphics.Rendering.OpenGL hiding (Normal3,rotate,translate,scale)++-------------------------------------------------------++{-# SPECIALISE radianToDegrees :: Float -> Float #-}+{-# SPECIALISE radianToDegrees :: Double -> Double #-}+radianToDegrees :: RealFrac a => a -> a+radianToDegrees x = x * 57.295779513082322++{-# SPECIALIZE degreesToRadian :: Float  -> Float  #-}+{-# SPECIALIZE degreesToRadian :: Double -> Double #-}+degreesToRadian :: Floating a => a -> a+degreesToRadian x = x * 1.7453292519943295e-2++-- | The angle is in radians. (WARNING: OpenGL uses degrees!)+rotate :: Flt -> Vec3 -> IO ()+rotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z)++translate :: Vec3 -> IO ()+translate (Vec3 x y z) = GL.translate (Vector3 x y z)++scale3 :: Vec3 -> IO ()+scale3 (Vec3 x y z) = GL.scale x y z++scale :: Flt -> IO ()+scale x = GL.scale x x x++-------------------------------------------------------++-- Vertex instances++instance GL.Vertex Vec2 where+  vertex (Vec2 x y) = GL.vertex (GL.Vertex2 x y)+  vertexv p = peek p >>= vertex +  +instance GL.Vertex Vec3 where+  vertex (Vec3 x y z) = GL.vertex (GL.Vertex3 x y z)+  vertexv p = peek p >>= vertex   +  +instance GL.Vertex Vec4 where+  vertex (Vec4 x y z w) = GL.vertex (GL.Vertex4 x y z w)+  vertexv p = peek p >>= vertex   ++-------------------------------------------------------++-- the Normal instance+-- note that there is no Normal2\/Normal4 in the OpenGL binding++instance GL.Normal Normal3 where+  normal (Normal3 (Vec3 x y z)) = GL.normal (GL.Normal3 x y z)+  normalv p = peek p >>= normal ++-------------------------------------------------------++-- Color instances+  +instance GL.Color Vec3 where+  color (Vec3 r g b) = GL.color (GL.Color3 r g b)+  colorv p = peek p >>= color++instance GL.Color Vec4 where+  color (Vec4 r g b a) = GL.color (GL.Color4 r g b a)+  colorv p = peek p >>= color++instance GL.SecondaryColor Vec3 where+  secondaryColor (Vec3 r g b) = GL.secondaryColor (GL.Color3 r g b)+  secondaryColorv p = peek p >>= secondaryColor++{-+-- there is no such thing?+instance GL.SecondaryColor Vec4 where+  secondaryColor (Vec4 r g b a) = GL.secondaryColor (GL.Color4 r g b a)+  secondaryColorv p = peek p >>= secondaryColor+-}++-------------------------------------------------------++-- TexCoord instances++instance GL.TexCoord Vec2 where+  texCoord (Vec2 u v) = GL.texCoord (GL.TexCoord2 u v)+  texCoordv p = peek p >>= texCoord+  multiTexCoord unit (Vec2 u v) = GL.multiTexCoord unit (GL.TexCoord2 u v)+  multiTexCoordv unit p = peek p >>= multiTexCoord unit++instance GL.TexCoord Vec3 where+  texCoord (Vec3 u v w) = GL.texCoord (GL.TexCoord3 u v w)+  texCoordv p = peek p >>= texCoord+  multiTexCoord unit (Vec3 u v w) = GL.multiTexCoord unit (GL.TexCoord3 u v w)+  multiTexCoordv unit p = peek p >>= multiTexCoord unit++instance GL.TexCoord Vec4 where+  texCoord (Vec4 u v w z) = GL.texCoord (GL.TexCoord4 u v w z)+  texCoordv p = peek p >>= texCoord+  multiTexCoord unit (Vec4 u v w z) = GL.multiTexCoord unit (GL.TexCoord4 u v w z)+  multiTexCoordv unit p = peek p >>= multiTexCoord unit++-------------------------------------------------------+    +-- Vertex Attributes (experimental)++class VertexAttrib' a where+  vertexAttrib :: GL.AttribLocation -> a -> IO ()+  +instance VertexAttrib' {- ' CPP is sensitive to primes -} Flt where+  vertexAttrib loc x = GL.vertexAttrib1 loc x++instance VertexAttrib' Vec2 where+  vertexAttrib loc (Vec2 x y) = GL.vertexAttrib2 loc x y++instance VertexAttrib' Vec3 where+  vertexAttrib loc (Vec3 x y z) = GL.vertexAttrib3 loc x y z ++instance VertexAttrib' Vec4 where+  vertexAttrib loc (Vec4 x y z w) = GL.vertexAttrib4 loc x y z w ++instance VertexAttrib' Normal2 where+  vertexAttrib loc (Normal2 (Vec2 x y)) = GL.vertexAttrib2 loc x y++instance VertexAttrib' Normal3 where+  vertexAttrib loc (Normal3 (Vec3 x y z)) = GL.vertexAttrib3 loc x y z++instance VertexAttrib' Normal4 where+  vertexAttrib loc (Normal4 (Vec4 x y z w)) = GL.vertexAttrib4 loc x y z w++-------------------------------------------------------++-- Uniform (again, experimental)++-- (note that the uniform location code in the OpenGL 2.2.1.1 is broken; +-- a work-around is to put a zero character at the end of uniform names)++{-+toFloat :: Flt -> Float+toFloat = realToFrac++fromFloat :: Float -> Flt+fromFloat = realToFrac+-}++-- Uniforms are always floats...+#ifdef VECT_Float++instance GL.Uniform Flt where+  uniform loc = GL.makeStateVar getter setter where+    getter = liftM (\(GL.Index1 x) -> x) $ get (uniform loc)+    setter x = ($=) (uniform loc) (Index1 x) +  uniformv loc cnt ptr = uniformv loc cnt (castPtr ptr :: Ptr (Index1 Flt))++instance GL.Uniform Vec2 where+  uniform loc = GL.makeStateVar getter setter where+    getter = liftM (\(GL.Vertex2 x y) -> Vec2 x y) $ get (uniform loc)+    setter (Vec2 x y) = ($=) (uniform loc) (Vertex2 x y) +  uniformv loc cnt ptr = uniformv loc (2*cnt) (castPtr ptr :: Ptr Flt)++instance GL.Uniform Vec3 where+  uniform loc = GL.makeStateVar getter setter where+    getter = liftM (\(GL.Vertex3 x y z) -> Vec3 x y z) $ get (uniform loc)+    setter (Vec3 x y z) = ($=) (uniform loc) (Vertex3 x y z) +  uniformv loc cnt ptr = uniformv loc (3*cnt) (castPtr ptr :: Ptr Flt)++instance GL.Uniform Vec4 where+  uniform loc = GL.makeStateVar getter setter where+    getter = liftM (\(GL.Vertex4 x y z w) -> Vec4 x y z w) $ get (uniform loc)+    setter (Vec4 x y z w) = ($=) (uniform loc) (Vertex4 x y z w) +  uniformv loc cnt ptr = uniformv loc (4*cnt) (castPtr ptr :: Ptr Flt)+    +#endif
+ Data/Vect/Float/Util/Dim2.hs view
@@ -0,0 +1,65 @@+{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-}++module Data.Vect.Flt.Util.Dim2 where++import Data.Vect.Flt.Base++-- |example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@.+structVec2 :: [Flt] -> [Vec2]+structVec2 [] = []+structVec2 (x:y:ls) = (Vec2 x y):(structVec2 ls) +structVec2 _ = error "structVec2"++destructVec2 :: [Vec2] -> [Flt]+destructVec2 [] = []+destructVec2 ((Vec2 x y):ls) = x:y:(destructVec2 ls)  ++det2 :: Vec2 -> Vec2 -> Flt+det2 u v = det (u,v)++vec2X :: Vec2+vec2Y :: Vec2++vec2X = Vec2 1 0 +vec2Y = Vec2 0 1 ++translate2X :: Flt -> Vec2 -> Vec2+translate2Y :: Flt -> Vec2 -> Vec2++translate2X t (Vec2 x y) = Vec2 (x+t) y +translate2Y t (Vec2 x y) = Vec2 x (y+t) ++-- | unit vector with given angle relative to the positive X axis (in the positive direction, that is, CCW).+-- A more precise name would be @cosSin@, but that sounds bad :)+sinCos :: Flt -> Vec2+sinCos a = Vec2 (cos a) (sin a)++sinCos' {- ' CPP is sensitive to primes -} :: Flt -> Normal2+sinCos' = toNormalUnsafe . sinCos++sinCosRadius :: Flt    -- ^ angle (in radians)+             -> Flt    -- ^ radius+             -> Vec2+sinCosRadius a r = Vec2 (r * cos a) (r * sin a)++-- | The angle relative to the positive X axis+angle2 :: Vec2 -> Flt+angle2 (Vec2 x y) = atan2 y x++angle2' {- ' CPP is sensitive to primes -} :: Normal2 -> Flt+angle2' = angle2 . fromNormal++-- |Rotation matrix by a given angle (in radians), counterclockwise.+rotMatrix2 :: Flt -> Mat2+rotMatrix2 a = Mat2 (Vec2 c s) (Vec2 (-s) c) where c = cos a; s = sin a++rotate2 :: Flt -> Vec2 -> Vec2+rotate2 a v = v .* (rotMatrix2 a) ++-- |Rotates counterclockwise by 90 degrees.+rotateCCW :: Vec2 -> Vec2+rotateCCW (Vec2 x y) = Vec2 (-y) x++-- |Rotates clockwise by 90 degrees.+rotateCW :: Vec2 -> Vec2+rotateCW (Vec2 x y) = Vec2 y (-x)
+ Data/Vect/Float/Util/Dim3.hs view
@@ -0,0 +1,76 @@+{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-}++module Data.Vect.Flt.Util.Dim3 where++import Data.Vect.Flt.Base++structVec3 :: [Flt] -> [Vec3]+structVec3 [] = []+structVec3 (x:y:z:ls) = (Vec3 x y z):(structVec3 ls) +structVec3 _ = error "structVec3"++destructVec3 :: [Vec3] -> [Flt]+destructVec3 [] = []+destructVec3 ((Vec3 x y z):ls) = x:y:z:(destructVec3 ls)  ++det3 :: Vec3 -> Vec3 -> Vec3 -> Flt+det3 u v w = det (u,v,w)++translate3X :: Flt -> Vec3 -> Vec3+translate3Y :: Flt -> Vec3 -> Vec3+translate3Z :: Flt -> Vec3 -> Vec3++translate3X t (Vec3 x y z) = Vec3 (x+t) y z +translate3Y t (Vec3 x y z) = Vec3 x (y+t) z +translate3Z t (Vec3 x y z) = Vec3 x y (z+t) ++vec3X :: Vec3+vec3Y :: Vec3+vec3Z :: Vec3++vec3X = Vec3 1 0 0+vec3Y = Vec3 0 1 0+vec3Z = Vec3 0 0 1++rotMatrixZ :: Flt -> Mat3+rotMatrixY :: Flt -> Mat3+rotMatrixX :: Flt -> Mat3++-- These are intended for multiplication on the /right/.+-- Should be consistent with the rotation around an arbitrary axis +-- (eg, @rotMatrixY a == rotate3 a vec3Y@)+rotMatrixZ a = Mat3 (Vec3 c s 0) (Vec3 (-s) c 0) (Vec3 0 0 1) where c = cos a; s = sin a+rotMatrixY a = Mat3 (Vec3 c 0 (-s)) (Vec3 0 1 0) (Vec3 s 0 c) where c = cos a; s = sin a+rotMatrixX a = Mat3 (Vec3 1 0 0) (Vec3 0 c s) (Vec3 0 (-s) c) where c = cos a; s = sin a++rotate3' :: {- ' CPP is sensitive to primes -} Flt       -- ^ angle (in radians)+         -> Normal3   -- ^ axis (should be a /unit/ vector!) +         -> Vec3      -- ^ vector+         -> Vec3      -- ^ result+rotate3' angle axis v = v .* (rotMatrix3' axis angle)++rotate3 :: Flt    -- ^ angle (in radians)+        -> Vec3   -- ^ axis (arbitrary nonzero vector)+        -> Vec3   -- ^ vector+        -> Vec3   -- ^ result+rotate3 angle axis v = v .* (rotMatrix3 axis angle)+      +-- |Rotation around an arbitrary 3D vector. The resulting 3x3 matrix is intended for multiplication on the /right/. +rotMatrix3 :: Vec3 -> Flt -> Mat3+rotMatrix3 v a = rotMatrix3' (mkNormal v) a++-- |Rotation around an arbitrary 3D /unit/ vector. The resulting 3x3 matrix is intended for multiplication on the /right/. +rotMatrix3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Mat3+rotMatrix3' (Normal3 v) a = +  let c = cos a+      s = sin a+      m1 = scalarMul (1-c) (outer v v)+      x = _1 v+      y = _2 v+      z = _3 v+      m2 = Mat3 (Vec3   c    ( s*z) (-s*y))+                (Vec3 (-s*z)   c    ( s*x))+                (Vec3 ( s*y) (-s*x)   c   )+  in (m1 &+ m2)++
+ Data/Vect/Float/Util/Dim4.hs view
@@ -0,0 +1,93 @@+{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-}++-- | Rotation around an arbitrary plane in four dimensions, and other miscellanea.+-- Not very useful for most people, and not re-exported by "Data.Vect".++module Data.Vect.Flt.Util.Dim4 where++import Data.Vect.Flt.Base+import Data.Vect.Flt.GramSchmidt++structVec4 :: [Flt] -> [Vec4]+structVec4 [] = []+structVec4 (x:y:z:w:ls) = (Vec4 x y z w):(structVec4 ls) +structVec4 _ = error "structVec4"++destructVec4 :: [Vec4] -> [Flt]+destructVec4 [] = []+destructVec4 ((Vec4 x y z w):ls) = x:y:z:w:(destructVec4 ls)  ++--det4 :: Vec4 -> Vec4 -> Vec4 -> Vec4 -> Flt+--det4 u v w z = det (u,v,w,z)++translate4X :: Flt -> Vec4 -> Vec4+translate4Y :: Flt -> Vec4 -> Vec4+translate4Z :: Flt -> Vec4 -> Vec4+translate4W :: Flt -> Vec4 -> Vec4++translate4X t (Vec4 x y z w) = Vec4 (x+t) y z w +translate4Y t (Vec4 x y z w) = Vec4 x (y+t) z w +translate4Z t (Vec4 x y z w) = Vec4 x y (z+t) w+translate4W t (Vec4 x y z w) = Vec4 x y z (w+t) ++vec4X :: Vec4+vec4Y :: Vec4+vec4Z :: Vec4+vec4W :: Vec4++vec4X = Vec4 1 0 0 0+vec4Y = Vec4 0 1 0 0+vec4Z = Vec4 0 0 1 0+vec4W = Vec4 0 0 0 1++---------------------------------------------------------------------------++-- |If @(x,y,u,v)@ is an orthonormal system, then (written in pseudo-code)+-- @biVector4 (x,y) = plusMinus (reverse $ biVector4 (u,v))@.+-- This is a helper function for the 4 dimensional rotation code.+-- If @(x,y,z,p,q,r) = biVector4 a b@, then the corresponding antisymmetric tensor is+--+-- > [  0  r  q  p ]+-- > [ -r  0  z -y ]+-- > [ -q -z  0  x ]+-- > [ -p  y -x  0 ]+biVector4 :: Vec4 -> Vec4 -> (Flt,Flt,Flt,Flt,Flt,Flt)+biVector4 (Vec4 x y z w) (Vec4 a b c d) = +  ( x*b-y*a , x*c-z*a , x*d-w*a , y*c-z*b , -y*d+w*b , z*d-w*c )++-- | the corresponding antisymmetric tensor+biVector4AsTensor :: Vec4 -> Vec4 -> Mat4+biVector4AsTensor v w = +  Mat4 ( Vec4   0  ( r) ( q) ( p) )+       ( Vec4 (-r)   0  ( z) (-y) )+       ( Vec4 (-q) (-z)   0  ( x) )+       ( Vec4 (-p) ( y) (-x)   0  )+  where +    (x,y,z,p,q,r) = biVector4 v w++-- | We assume that the axes are normalized and /orthogonal/ to each other!+rotate4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Vec4 -> Vec4+rotate4' angle axes v = v .* (rotMatrix4' angle axes)++-- | We assume only that the axes are independent vectors.+rotate4 :: Flt -> (Vec4,Vec4) -> Vec4 -> Vec4+rotate4 angle axes v = v .* (rotMatrix4 angle axes)++-- | Rotation matrix around a plane specified by two normalized and /orthogonal/ vectors.+-- Intended for multiplication on the /right/!+rotMatrix4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Mat4+rotMatrix4' angle (Normal4 v, Normal4 w) = m1 &+ (s *& m2) &+ m3 +  where+    c = cos angle ; s = sin angle+    m1 = scalarMul (1-c) ( outer v v  &+  outer w w )+    m2 = biVector4AsTensor v w+    m3 = diag (Vec4 c c c c)++-- | We assume only that the axes are independent vectors.+rotMatrix4 :: Flt -> (Vec4,Vec4) -> Mat4  +rotMatrix4 angle axes = +  rotMatrix4' angle $ liftPair toNormalUnsafe $ gramSchmidtNormalize axes +  where +    liftPair f (x,y) = (f x, f y)+    +
+ Data/Vect/Float/Util/Projective.hs view
@@ -0,0 +1,88 @@+{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-}++-- | Classic 4x4 projective matrices. Our convention is that they are intended for multiplication on+-- the /right/, that is, they are of the form+--+-- >     _____+-- > [  |     |  0  ]+-- > [  | 3x3 |  0  ]+-- > [  |_____|  0  ]+-- > [  p  q  r  1  ]+--+-- Please note that by default, OpenGL stores the matrices (in memory) by columns, while we +-- store them by rows; but OpenGL also use the opposite convention (so the OpenGL projective matrices +-- are intended for multiplication on the /left/). So in effect, they are the same when stored in the memory,+-- say with @poke :: Ptr Mat4 -> Mat4 -> IO ()@.++module Data.Vect.Flt.Util.Projective where++import Data.Vect.Flt.Base+import Data.Vect.Flt.Util.Dim3++import qualified Data.Vect.Flt.Util.Dim4 as Dim4++class ExtendProjective v e | v->e where+  extendProj     :: v -> e+  extendProjWith :: Flt -> v -> e+  extendProj = extendProjWith 1+  +instance ExtendProjective Vec2 Vec4 where+  extendProj       (Vec2 x y) = Vec4 x y 0 1+  extendProjWith w (Vec2 x y) = Vec4 x y 0 w+  +instance ExtendProjective Vec3 Vec4 where+  extendProj       (Vec3 x y z) = Vec4 x y z 1+  extendProjWith w (Vec3 x y z) = Vec4 x y z w++instance ExtendProjective Vec4 Vec4 where+  extendProj = id+  extendProjWith w (Vec4 x y z w') = let s = w/w' in Vec4 (s*x) (s*y) (s*z) w++instance ExtendProjective Mat2 Mat4 where+  extendProj       (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 1)+  extendProjWith w (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 w)++instance ExtendProjective Mat3 Mat4 where+  extendProj       (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 1)+  extendProjWith w (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 w)++rotMatrixProj :: Flt -> Normal3 -> Mat4+rotMatrixProj angle axis = extendProj $ rotMatrix3' axis angle++rotMatrixProj' :: {- ' CPP is sensitive to primes -} Flt -> Vec3 -> Mat4+rotMatrixProj' angle axis = extendProj $ rotMatrix3 axis angle++translMatrixProj :: Vec3 -> Mat4+translMatrixProj v = Mat4 Dim4.vec4X Dim4.vec4Y Dim4.vec4Z (extendProj v)++-- | we assume that the bottom-right corner is 1.+translWithProj :: Vec3 -> Mat4 -> Mat4+translWithProj v mat@(Mat4 r1 r2 r3 r4) = Mat4 r1 r2 r3 (extendProjWith 0 v &+ r4)++scaleMatrixProj :: Vec3 -> Mat4+scaleMatrixProj v = diag $ extendProj v++scaleMatrixUniformProj :: Flt -> Mat4+scaleMatrixUniformProj s = diag (Vec4 s s s 1)++class ProjectiveAction v where+  actProj :: v -> Mat4 -> v+ +instance ProjectiveAction Vec3 where+  actProj v m = trim $ (extendProj v) .* m ++instance ProjectiveAction Vec4 where+  actProj v m = v .* m ++-- | When acting on unit vectors, we ignore the translation part.+instance ProjectiveAction Normal3 where+  actProj (Normal3 v) m = Normal3 (v .* (trim m :: Mat3))++-- | Inverts a projective 4x4 matrix, assuming that the top-left 3x3 part is /orthogonal/,+-- and the bottom-right corner is 1.+invertProj :: Mat4 -> Mat4+invertProj mat@(Mat4 u v w t) = +  translWithProj t' $ extendProj $ transpose $ (trim mat :: Mat3)+  where+    t' = Vec3 (- u &. t) (- v &. t) (- w &. t)+
+ LICENSE view
@@ -0,0 +1,29 @@+Copyright (c) 2008, Balazs Komuves+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.+ +- Neither names of the copyright holders nor the names of the contributors+may be used to endorse or promote products derived from this software without+specific prior written permission. ++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.lhs view
@@ -0,0 +1,44 @@+#! /usr/bin/env runhaskell+>+> import Control.Monad+> import Distribution.Simple+> import Distribution.PackageDescription+> import System.IO+> import System.Directory+>+> copyFileWithPrefix src tgt prefix = +>   readFile src >>= \txt -> writeFile tgt (prefix ++ txt)+>+> copyFiles srcdir tgtdir prefix = do+>   files <- getDirectoryContents srcdir+>   forM_ files $ \fname -> do+>     let src = srcdir ++ fname+>         tgt = tgtdir ++ fname+>     doesFileExist src >>= \b -> when b $ copyFileWithPrefix src tgt prefix+>     doesDirectoryExist src>>= \b -> when ( b && fname /= "." && fname /= ".." ) $ do+>       createDirectoryIfMissing False tgt+>       copyFiles (src ++ "/") (tgt ++ "/") prefix+>+> thePrefix flt = "{-# OPTIONS_GHC -DFlt=" ++ flt ++ " -DVECT_" ++ flt ++ " #-}\n"+>+> myPreBuildHook args buildflags = do+>   createDirectoryIfMissing False "Data/Vect/Float"+>   createDirectoryIfMissing False "Data/Vect/Double"+>   copyFileWithPrefix "src/flt.hs" "Data/Vect/Float.hs"  (thePrefix "Float")+>   copyFileWithPrefix "src/flt.hs" "Data/Vect/Double.hs" (thePrefix "Double")+>   copyFiles "src/flt/" "Data/Vect/Float/"  (thePrefix "Float")+>   copyFiles "src/flt/" "Data/Vect/Double/" (thePrefix "Double")+>   return $ emptyHookedBuildInfo  +>+> myPostCleanHook args cleanflags pdep mlocalbuildinfo = do+>   removeDirectoryRecursive "Data/Vect/Float"+>   removeDirectoryRecursive "Data/Vect/Double"+>+> myUserHooks = simpleUserHooks +>   { preBuild = myPreBuildHook +>   , postClean = myPostCleanHook+>   }+>+> main = do+>   defaultMainWithHooks myUserHooks+>
+ src/flt.hs view
@@ -0,0 +1,22 @@++module Data.Vect.Flt+  ( module Data.Vect.Flt.Base+  , module Data.Vect.Flt.Interpolate+  , module Data.Vect.Flt.Util.Dim2+  , module Data.Vect.Flt.Util.Dim3+  , module Data.Vect.Flt.Util.Projective+#ifdef VECT_OPENGL        +  , module Data.Vect.Flt.OpenGL       +#endif+  ) where++import Data.Vect.Flt.Base+import Data.Vect.Flt.Interpolate++import Data.Vect.Flt.Util.Dim2+import Data.Vect.Flt.Util.Dim3+import Data.Vect.Flt.Util.Projective++#ifdef VECT_OPENGL         +import Data.Vect.Flt.OpenGL       +#endif
+ src/flt/Base.hs view
@@ -0,0 +1,734 @@++module Data.Vect.Flt.Base where++import Control.Monad+import System.Random  +import Foreign++-- class declarations++class AbelianGroup g where+  (&+) :: g -> g -> g+  (&-) :: g -> g -> g+  neg  :: g -> g+  zero :: g++infixl 6 &++infixl 6 &- ++vecSum :: AbelianGroup g => [g] -> g+vecSum l = foldl (&+) zero l ++class (AbelianGroup r) => +      Ring r where+  (.*.) :: r -> r -> r+  one   :: r++infixl 7 .*. ++ringProduct :: Ring r => [r] -> r+ringProduct l = foldl (.*.) one l++class LeftModule r m where+  lmul :: r -> m -> m+  (*.) :: r -> m -> m+  (*.) = lmul++class RightModule m r where+  rmul :: m -> r -> m+  (.*) :: m -> r -> m+  (.*) = rmul++-- I'm not really sure about this.. may actually degrade the performance in some cases?  +{- RULES+"matrix multiplication left"   forall m n x.  (n .*. m) *. x = n *. (m *. x)  +"matrix multiplication right"  forall m n x.  x .* (m .*. n) = (x .* m) .* n+  -}++infixr 7 *.+infixl 7 .*++class AbelianGroup v => Vector v where+  mapVec    :: (Flt -> Flt) -> v -> v+  scalarMul :: Flt -> v -> v+  (*&) ::      Flt -> v -> v +  (&*) ::      v -> Flt -> v +  (*&) s v = scalarMul s v+  (&*) v s = scalarMul s v++infixr 7 *&+infixl 7 &*++{-# RULES+"scalar multiplication left"   forall s t x.  t *& (s *& x) = (t*s) *& x +"scalar multiplication right"  forall s t x.  (x &* s) &* t = x &* (s*t)  +  #-}++class DotProd v where+  (&.) :: v -> v -> Flt+  norm    :: v -> Flt+  normsqr :: v -> Flt+  len     :: v -> Flt+  lensqr  :: v -> Flt+  len = norm+  lensqr = normsqr+  dotprod :: v -> v -> Flt+  normsqr v = (v &. v)  +  norm = sqrt.lensqr+  dotprod = (&.)++infix 7 &.++{-# RULES+"len/square 1"   forall x.  (len x)*(len x) = lensqr x+"len/square 2"   forall x.  (len x)^2 = lensqr x+"norm/square 1"  forall x.  (norm x)*(norm x) = normsqr x+"norm/square 2"  forall x.  (norm x)^2 = normsqr x+  #-}++normalize :: (Vector v, DotProd v) => v -> v+normalize v = scalarMul (1.0/(len v)) v++distance :: (Vector v, DotProd v) => v -> v -> Flt+distance x y = norm (x &- y)++-- | the angle between two vectors+angle :: (Vector v, DotProd v) => v -> v -> Flt +angle x y = acos $ (x &. y) / (norm x * norm y)++-- | the angle between two unit vectors+angle' {- ' CPP is sensitive to primes -} :: (Vector v, UnitVector v u, DotProd v) => u -> u -> Flt +angle' x y = acos (fromNormal x &. fromNormal y)++{-# RULES+"normalize is idempotent"  forall x. normalize (normalize x) = normalize x+  #-}++class (Vector v, DotProd v) => UnitVector v u | v->u, u->v  where+  mkNormal         :: v -> u       -- ^ normalizes the input+  toNormalUnsafe   :: v -> u       -- ^ does not normalize the input!+  fromNormal       :: u -> v+  fromNormalRadius :: Flt -> u -> v+  fromNormalRadius t n = t *& fromNormal n ++-- | projects the first vector onto the direction of the second (unit) vector+project' :: (Vector v, UnitVector v u, DotProd v) => v -> u -> v+project' what dir = projectUnsafe what (fromNormal dir)++-- | direction (second argument) is assumed to be a /unit/ vector!+projectUnsafe :: (Vector v, DotProd v) => v -> v -> v+projectUnsafe what dir = what &- dir &* (what &. dir)++project :: (Vector v, DotProd v) => v -> v -> v+project what dir = what &- dir &* ((what &. dir) / (dir &. dir))++-- | since unit vectors are not a group, we need a separate function.+flipNormal :: UnitVector v n => n -> n +flipNormal = toNormalUnsafe . neg . fromNormal ++class CrossProd v where+  crossprod :: v -> v -> v+  (&^)      :: v -> v -> v+  (&^) = crossprod+  +class Pointwise v where+  pointwise :: v -> v -> v+  (&!)      :: v -> v -> v+  (&!) = pointwise ++infix 7 &^+infix 7 &!++class HasCoordinates v x | v->x where+  _1 :: v -> x+  _2 :: v -> x+  _3 :: v -> x+  _4 :: v -> x      ++-- | conversion between vectors (and matrices) of different dimensions+class Extend u v where+  extendZero :: u -> v          -- ^ example: @extendZero (Vec2 5 6) = Vec4 5 6 0 0@+  extendWith :: Flt -> u -> v   -- ^ example: @extendWith 1 (Vec2 5 6) = Vec4 5 6 1 1@+  trim :: v -> u                -- ^ example: @trim (Vec4 5 6 7 8) = Vec2 5 6@++-- | makes a diagonal matrix from a vector+class Diagonal s t | t->s where+  diag :: s -> t++class Matrix m where+  transpose :: m -> m +  inverse :: m -> m+  idmtx :: m++{-# RULES+"transpose is an involution"  forall m. transpose (transpose m) = m+"inverse is an involution"    forall m. inverse (inverse m) = m+  #-}+  +-- | Outer product (could be unified with Diagonal?)+class Tensor t v | t->v where+  outer :: v -> v -> t+    +class Determinant m where+  det :: m -> Flt    +    +-- Vec / Mat datatypes+ +data Vec2 = Vec2 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt +  deriving (Read,Show)+data Vec3 = Vec3 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt +  deriving (Read,Show)+data Vec4 = Vec4 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt +  deriving (Read,Show)++-- | these are /row/ vectors +data Mat2 = Mat2 !Vec2 !Vec2              deriving (Read,Show)+data Mat3 = Mat3 !Vec3 !Vec3 !Vec3        deriving (Read,Show)+data Mat4 = Mat4 !Vec4 !Vec4 !Vec4 !Vec4  deriving (Read,Show)++-- | The assumption when dealing with these is always that they are of unit length.+-- Also, interpolation works differently.+newtype Normal2 = Normal2 Vec2 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) +newtype Normal3 = Normal3 Vec3 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable,CrossProd) +newtype Normal4 = Normal4 Vec4 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) ++mkVec2 :: (Flt,Flt) -> Vec2+mkVec3 :: (Flt,Flt,Flt) -> Vec3+mkVec4 :: (Flt,Flt,Flt,Flt) -> Vec4++mkVec2 (x,y)     = Vec2 x y +mkVec3 (x,y,z)   = Vec3 x y z+mkVec4 (x,y,z,w) = Vec4 x y z w++-- Unit vectors+  +instance UnitVector Vec2 Normal2 where+  mkNormal v = Normal2 (normalize v)+  fromNormal (Normal2 v) = v +  toNormalUnsafe = Normal2++instance UnitVector Vec3 Normal3 where+  mkNormal v = Normal3 (normalize v)+  fromNormal (Normal3 v) = v +  toNormalUnsafe = Normal3++instance UnitVector Vec4 Normal4 where+  mkNormal v = Normal4 (normalize v)+  fromNormal (Normal4 v) = v +  toNormalUnsafe = Normal4++rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)+rndUnit g = +  if d > 0.01+    then ( v &* (1.0/d) , h )+    else rndUnit h+  where+    (v,h) = random g+    d = norm v+    +instance Random Normal2 where+  random g = let (v,h) = rndUnit g in (Normal2 v, h)  +  randomR _ = random++instance Random Normal3 where+  random g = let (v,h) = rndUnit g in (Normal3 v, h)  +  randomR _ = random++instance Random Normal4 where+  random g = let (v,h) = rndUnit g in (Normal4 v, h)  +  randomR _ = random++{-+instance Storable Normal2 where+  alignment _ = alignment (undefined::Vec2)+  sizeOf    _ = sizeOf    (undefined::Vec2)+  peek p = liftM (\v -> Normal2 v) (peek $ castPtr p)+  poke p (Normal2 v) = poke (castPtr p) v  + +instance Storable Normal3 where+  alignment _ = alignment (undefined::Vec3)+  sizeOf    _ = sizeOf    (undefined::Vec3)+  peek p = liftM (\v -> Normal3 v) (peek $ castPtr p)+  poke p (Normal3 v) = poke (castPtr p) v  ++instance Storable Normal4 where+  alignment _ = alignment (undefined::Vec4)+  sizeOf    _ = sizeOf    (undefined::Vec4)+  peek p = liftM (\v -> Normal4 v) (peek $ castPtr p)+  poke p (Normal4 v) = poke (castPtr p) v  +-}++-- Vec2 instances++instance HasCoordinates Vec2 Flt where+  _1 (Vec2 x _) = x+  _2 (Vec2 _ y) = y+  _3 _ = error "has only 2 coordinates"+  _4 _ = error "has only 2 coordinates"++instance AbelianGroup Vec2 where+  (&+) (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1+x2) (y1+y2) +  (&-) (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1-x2) (y1-y2)+  neg  (Vec2 x y)                = Vec2 (-x) (-y)+  zero = Vec2 0 0+  +instance Vector Vec2 where+  scalarMul s (Vec2 x y) = Vec2 (s*x) (s*y)+  mapVec    f (Vec2 x y) = Vec2 (f x) (f y)+  +instance DotProd Vec2 where+  (&.) (Vec2 x1 y1) (Vec2 x2 y2) = x1*x2 + y1*y2++instance Pointwise Vec2 where+  pointwise (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1*x2) (y1*y2)++instance Determinant (Vec2,Vec2) where+  det (Vec2 x1 y1 , Vec2 x2 y2) = x1*y2 - x2*y1  ++{-     +instance Show Vec2 where+  show (Vec2 x y) = "( " ++ show x ++ " , " ++ show y ++ " )"+-}++instance Random Vec2 where+  random = randomR (Vec2 (-1) (-1),Vec2 1 1)+  randomR (Vec2 a b, Vec2 c d) gen = +    let (x,gen1) = randomR (a,c) gen+        (y,gen2) = randomR (b,d) gen1+    in (Vec2 x y, gen2)+     +instance Storable Vec2 where+  sizeOf    _ = 2 * sizeOf (undefined::Flt)+  alignment _ = sizeOf (undefined::Flt)+  +  peek q = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    x <- peek        p +    y <- peekByteOff p k+    return (Vec2 x y)+    +  poke q (Vec2 x y) = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    poke        p   x+    pokeByteOff p k y+               +-- Mat2 instances++instance HasCoordinates Mat2 Vec2 where+  _1 (Mat2 x _) = x+  _2 (Mat2 _ y) = y+  _3 _ = error "has only 2 coordinates"+  _4 _ = error "has only 2 coordinates"++instance Matrix Mat2 where+  transpose (Mat2 row1 row2) = +    Mat2 (Vec2 (_1 row1) (_1 row2)) +         (Vec2 (_2 row1) (_2 row2)) +  idmtx = Mat2 (Vec2 1 0) (Vec2 0 1)+  inverse (Mat2 (Vec2 a b) (Vec2 c d)) = +    Mat2 (Vec2 (d*r) (-b*r)) (Vec2 (-c*r) (a*r)) +    where r = 1.0 / (a*d - b*c)++instance AbelianGroup Mat2 where+  (&+) (Mat2 r1 r2) (Mat2 s1 s2) = Mat2 (r1 &+ s1) (r2 &+ s2)+  (&-) (Mat2 r1 r2) (Mat2 s1 s2) = Mat2 (r1 &- s1) (r2 &- s2)+  neg  (Mat2 r1 r2)              = Mat2 (neg r1) (neg r2)  +  zero = Mat2 zero zero   -- (zero::Vec2) (zero::Vec2)++instance Vector Mat2 where+  scalarMul s (Mat2 r1 r2) = Mat2 (g r1) (g r2) where g = scalarMul s+  mapVec    f (Mat2 r1 r2) = Mat2 (g r1) (g r2) where g = mapVec f++instance Ring Mat2 where+  (.*.) (Mat2 r1 r2) n = +    let (Mat2 c1 c2) = transpose n+    in Mat2 (Vec2 (r1 &. c1) (r1 &. c2))+            (Vec2 (r2 &. c1) (r2 &. c2))+  one = idmtx ++instance LeftModule Mat2 Vec2 where+  lmul (Mat2 row1 row2) v = Vec2 (row1 &. v) (row2 &. v) +  +instance RightModule Vec2 Mat2 where+  rmul v mt = lmul (transpose mt) v++instance Diagonal Vec2 Mat2 where+  diag (Vec2 x y) = Mat2 (Vec2 x 0) (Vec2 0 y)++instance Tensor Mat2 Vec2 where+  outer (Vec2 a b) (Vec2 x y) = Mat2+    (Vec2 (a*x) (a*y))+    (Vec2 (b*x) (b*y))+{-+  outer v w = +    let full = Mat2 (Vec2 1 1) (Vec2 1 1)+    in  (diag v) .*. full .*. (diag w)+-}++instance Determinant Mat2 where+  det (Mat2 (Vec2 a b) (Vec2 c d)) = a*d - b*c ++{-+instance Show Mat2 where+  show (Mat2 r1 r2) = show r1 ++ "\n" ++ show r2+-}++instance Storable Mat2 where+  sizeOf    _ = 2 * sizeOf (undefined::Vec2)+  alignment _ = alignment  (undefined::Vec2)+  +  peek q = do+    let p = castPtr q :: Ptr Vec2+        k = sizeOf (undefined::Vec2)+    r1 <- peek        p +    r2 <- peekByteOff p k+    return (Mat2 r1 r2)+    +  poke q (Mat2 r1 r2) = do+    let p = castPtr q :: Ptr Vec2+        k = sizeOf (undefined::Vec2)+    poke        p   r1+    pokeByteOff p k r2++-- Vec3 instances++instance HasCoordinates Vec3 Flt where+  _1 (Vec3 x _ _) = x+  _2 (Vec3 _ y _) = y+  _3 (Vec3 _ _ z) = z+  _4 _ = error "has only 3 coordinates"++instance AbelianGroup Vec3 where+  (&+) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1+x2) (y1+y2) (z1+z2) +  (&-) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1-x2) (y1-y2) (z1-z2) +  neg  (Vec3 x y z)                    = Vec3 (-x) (-y) (-z)+  zero = Vec3 0 0 0+  +instance Vector Vec3 where+  scalarMul s (Vec3 x y z) = Vec3 (s*x) (s*y) (s*z)+  mapVec    f (Vec3 x y z) = Vec3 (f x) (f y) (f z)++instance DotProd Vec3 where+  (&.) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = x1*x2 + y1*y2 + z1*z2++instance Pointwise Vec3 where+  pointwise (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1*x2) (y1*y2) (z1*z2)++{-+instance Show Vec3 where+  show (Vec3 x y z) = "( " ++ show x ++ " , " ++ show y ++ " , " ++ show z ++ " )"+-}++instance Random Vec3 where+  random = randomR (Vec3 (-1) (-1) (-1),Vec3 1 1 1)+  randomR (Vec3 a b c, Vec3 d e f) gen = +    let (x,gen1) = randomR (a,d) gen+        (y,gen2) = randomR (b,e) gen1+        (z,gen3) = randomR (c,f) gen2  +    in (Vec3 x y z, gen3)+      +instance CrossProd Vec3 where+  crossprod (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (y1*z2-y2*z1) (z1*x2-z2*x1) (x1*y2-x2*y1) ++instance Determinant (Vec3,Vec3,Vec3) where+  det (u,v,w) = u &. (v &^ w)  + +instance Storable Vec3 where+  sizeOf    _ = 3 * sizeOf (undefined::Flt)+  alignment _ = sizeOf (undefined::Flt)+  +  peek q = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    x <- peek        p +    y <- peekByteOff p (k  )+    z <- peekByteOff p (k+k)+    return (Vec3 x y z)+    +  poke q (Vec3 x y z) = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    poke        p       x+    pokeByteOff p (k  ) y+    pokeByteOff p (k+k) z+   +-- Mat3 instances++instance HasCoordinates Mat3 Vec3 where+  _1 (Mat3 x _ _) = x+  _2 (Mat3 _ y _) = y+  _3 (Mat3 _ _ z) = z+  _4 _ = error "has only 3 coordinates"  ++instance Matrix Mat3 where++  transpose (Mat3 row1 row2 row3) = +    Mat3 (Vec3 (_1 row1) (_1 row2) (_1 row3)) +         (Vec3 (_2 row1) (_2 row2) (_2 row3)) +         (Vec3 (_3 row1) (_3 row2) (_3 row3)) +         +  idmtx = Mat3 (Vec3 1 0 0) (Vec3 0 1 0) (Vec3 0 0 1)+  +  inverse (Mat3 (Vec3 a b c) (Vec3 e f g) (Vec3 i j k)) = +    Mat3 (Vec3 (d11*r) (d21*r) (d31*r))  +         (Vec3 (d12*r) (d22*r) (d32*r))  +         (Vec3 (d13*r) (d23*r) (d33*r))  +    where+      r = 1.0 / ( a*d11 + b*d12 + c*d13 )++      d11 = f*k - g*j+      d12 = g*i - e*k+      d13 = e*j - f*i++      d31 = b*g - c*f+      d32 = c*e - a*g+      d33 = a*f - b*e++      d21 = c*j - b*k +      d22 = a*k - c*i +      d23 = b*i - a*j ++instance AbelianGroup Mat3 where+  (&+) (Mat3 r1 r2 r3) (Mat3 s1 s2 s3) = Mat3 (r1 &+ s1) (r2 &+ s2) (r3 &+ s3)+  (&-) (Mat3 r1 r2 r3) (Mat3 s1 s2 s3) = Mat3 (r1 &- s1) (r2 &- s2) (r3 &- s3)+  neg  (Mat3 r1 r2 r3)                 = Mat3 (neg r1) (neg r2) (neg r3) +  zero = Mat3 zero zero zero   -- (zero::Vec3) (zero::Vec3) (zero::Vec3)++instance Vector Mat3 where+  scalarMul s (Mat3 r1 r2 r3) = Mat3 (g r1) (g r2) (g r3) where g = scalarMul s+  mapVec    f (Mat3 r1 r2 r3) = Mat3 (g r1) (g r2) (g r3) where g = mapVec f++instance Ring Mat3 where+  (.*.) (Mat3 r1 r2 r3) n = +    let (Mat3 c1 c2 c3) = transpose n+    in Mat3 (Vec3 (r1 &. c1) (r1 &. c2) (r1 &. c3))+            (Vec3 (r2 &. c1) (r2 &. c2) (r2 &. c3))+            (Vec3 (r3 &. c1) (r3 &. c2) (r3 &. c3))+  one = idmtx ++instance LeftModule Mat3 Vec3 where+  lmul (Mat3 row1 row2 row3) v = Vec3 (row1 &. v) (row2 &. v) (row3 &. v)+  +instance RightModule Vec3 Mat3 where+  rmul v mt = lmul (transpose mt) v++instance Diagonal Vec3 Mat3 where+  diag (Vec3 x y z) = Mat3 (Vec3 x 0 0) (Vec3 0 y 0) (Vec3 0 0 z)++instance Tensor Mat3 Vec3 where+  outer (Vec3 a b c) (Vec3 x y z) = Mat3+    (Vec3 (a*x) (a*y) (a*z))+    (Vec3 (b*x) (b*y) (b*z))+    (Vec3 (c*x) (c*y) (c*z))+{-+  outer v w = +    let full = Mat3 (Vec3 1 1 1) (Vec3 1 1 1) (Vec3 1 1 1)+    in  (diag v) .*. full .*. (diag w)+-}++instance Determinant Mat3 where+  det (Mat3 r1 r2 r3) = det (r1,r2,r3)++{-+instance Show Mat3 where+  show (Mat3 r1 r2 r3) = show r1 ++ "\n" ++ show r2 ++ "\n" ++ show r3+-}++instance Storable Mat3 where+  sizeOf    _ = 3 * sizeOf (undefined::Vec3)+  alignment _ = alignment  (undefined::Vec3)+  +  peek q = do+    let p = castPtr q :: Ptr Vec3+        k = sizeOf (undefined::Vec3)+    r1 <- peek        p +    r2 <- peekByteOff p (k  )+    r3 <- peekByteOff p (k+k)+    return (Mat3 r1 r2 r3)+    +  poke q (Mat3 r1 r2 r3) = do+    let p = castPtr q :: Ptr Vec3+        k = sizeOf (undefined::Vec3)+    poke        p       r1+    pokeByteOff p (k  ) r2+    pokeByteOff p (k+k) r3++-- Vec4 instances++instance HasCoordinates Vec4 Flt where+  _1 (Vec4 x _ _ _) = x+  _2 (Vec4 _ y _ _) = y+  _3 (Vec4 _ _ z _) = z+  _4 (Vec4 _ _ _ w) = w++instance AbelianGroup Vec4 where+  (&+) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1+x2) (y1+y2) (z1+z2) (w1+w2)+  (&-) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1-x2) (y1-y2) (z1-z2) (w1-w2)+  neg  (Vec4 x y z w)                        = Vec4 (-x) (-y) (-z) (-w)+  zero = Vec4 0 0 0 0+  +instance Vector Vec4 where+  scalarMul s (Vec4 x y z w) = Vec4 (s*x) (s*y) (s*z) (s*w)+  mapVec    f (Vec4 x y z w) = Vec4 (f x) (f y) (f z) (f w)++instance DotProd Vec4 where+  (&.) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = x1*x2 + y1*y2 + z1*z2 + w1*w2++instance Pointwise Vec4 where+  pointwise (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1*x2) (y1*y2) (z1*z2) (w1*w2)++{-+instance Show Vec4 where+  show (Vec4 x y z w) = "( " ++ show x ++ " , " ++ show y ++ " , " ++ show z ++ " , " ++ show w ++ " )"+-}++instance Random Vec4 where+  random = randomR (Vec4 (-1) (-1) (-1) (-1),Vec4 1 1 1 1)+  randomR (Vec4 a b c d, Vec4 e f g h) gen = +    let (x,gen1) = randomR (a,e) gen+        (y,gen2) = randomR (b,f) gen1+        (z,gen3) = randomR (c,g) gen2  +        (w,gen4) = randomR (d,h) gen3  +    in (Vec4 x y z w, gen4)+           +instance Storable Vec4 where+  sizeOf    _ = 4 * sizeOf (undefined::Flt)+  alignment _ = sizeOf (undefined::Flt)+  +  peek q = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    x <- peek        p +    y <- peekByteOff p (k  )+    z <- peekByteOff p (k+k)+    w <- peekByteOff p (3*k)+    return (Vec4 x y z w)+    +  poke q (Vec4 x y z w) = do+    let p = castPtr q :: Ptr Flt+        k = sizeOf (undefined::Flt)+    poke        p       x+    pokeByteOff p (k  ) y+    pokeByteOff p (k+k) z+    pokeByteOff p (3*k) w++-- Mat4 instances++instance HasCoordinates Mat4 Vec4 where+  _1 (Mat4 x _ _ _) = x+  _2 (Mat4 _ y _ _) = y+  _3 (Mat4 _ _ z _) = z+  _4 (Mat4 _ _ _ w) = w++instance Matrix Mat4 where+  transpose (Mat4 row1 row2 row3 row4) = +    Mat4 (Vec4 (_1 row1) (_1 row2) (_1 row3) (_1 row4)) +         (Vec4 (_2 row1) (_2 row2) (_2 row3) (_2 row4)) +         (Vec4 (_3 row1) (_3 row2) (_3 row3) (_3 row4)) +         (Vec4 (_4 row1) (_4 row2) (_4 row3) (_4 row4)) +  idmtx = Mat4 (Vec4 1 0 0 0) (Vec4 0 1 0 0) (Vec4 0 0 1 0) (Vec4 0 0 0 1)+  inverse = error "inverse/Mat4: not implemented yet"++instance AbelianGroup Mat4 where+  (&+) (Mat4 r1 r2 r3 r4) (Mat4 s1 s2 s3 s4) = Mat4 (r1 &+ s1) (r2 &+ s2) (r3 &+ s3) (r4 &+ s4)+  (&-) (Mat4 r1 r2 r3 r4) (Mat4 s1 s2 s3 s4) = Mat4 (r1 &- s1) (r2 &- s2) (r3 &- s3) (r4 &- s4)+  neg  (Mat4 r1 r2 r3 r4)                    = Mat4 (neg r1) (neg r2) (neg r3) (neg r4) +  zero = Mat4 zero zero zero zero+  +instance Vector Mat4 where+  scalarMul s (Mat4 r1 r2 r3 r4) = Mat4 (g r1) (g r2) (g r3) (g r4) where g = scalarMul s+  mapVec    f (Mat4 r1 r2 r3 r4) = Mat4 (g r1) (g r2) (g r3) (g r4) where g = mapVec f++instance Ring Mat4 where+  (.*.) (Mat4 r1 r2 r3 r4) n = +    let (Mat4 c1 c2 c3 c4) = transpose n+    in Mat4 (Vec4 (r1 &. c1) (r1 &. c2) (r1 &. c3) (r1 &. c4))+            (Vec4 (r2 &. c1) (r2 &. c2) (r2 &. c3) (r2 &. c4))+            (Vec4 (r3 &. c1) (r3 &. c2) (r3 &. c3) (r3 &. c4))+            (Vec4 (r4 &. c1) (r4 &. c2) (r4 &. c3) (r4 &. c4))+  one = idmtx ++instance LeftModule Mat4 Vec4 where+  lmul (Mat4 row1 row2 row3 row4) v = Vec4 (row1 &. v) (row2 &. v) (row3 &. v) (row4 &. v)+  +instance RightModule Vec4 Mat4 where+  rmul v mt = lmul (transpose mt) v++instance Diagonal Vec4 Mat4 where+  diag (Vec4 x y z w) = Mat4 (Vec4 x 0 0 0) (Vec4 0 y 0 0) (Vec4 0 0 z 0) (Vec4 0 0 0 w)++instance Tensor Mat4 Vec4 where+  outer (Vec4 a b c d) (Vec4 x y z w) = Mat4+    (Vec4 (a*x) (a*y) (a*z) (a*w))+    (Vec4 (b*x) (b*y) (b*z) (b*w))+    (Vec4 (c*x) (c*y) (c*z) (c*w))+    (Vec4 (d*x) (d*y) (d*z) (d*w))+{-+  outer v w = +    let full = Mat4 (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1)+    in  (diag v) .*. full .*. (diag w)+-}++--instance Determinant Mat4 where+--  det (Mat4 r1 r2 r3 r4)++{-+instance Show Mat4 where+  show (Mat4 r1 r2 r3 r4) = show r1 ++ "\n" ++ show r2 ++ "\n" ++ show r3 ++ "\n" ++ show r4+-}++instance Storable Mat4 where+  sizeOf    _ = 4 * sizeOf (undefined::Vec4)+  alignment _ = alignment  (undefined::Vec4)+  +  peek q = do+    let p = castPtr q :: Ptr Vec4+        k = sizeOf (undefined::Vec4)+    r1 <- peek        p +    r2 <- peekByteOff p (k  )+    r3 <- peekByteOff p (k+k)+    r4 <- peekByteOff p (3*k)+    return (Mat4 r1 r2 r3 r4)+    +  poke q (Mat4 r1 r2 r3 r4) = do+    let p = castPtr q :: Ptr Vec4+        k = sizeOf (undefined::Vec4)+    poke        p       r1+    pokeByteOff p (k  ) r2+    pokeByteOff p (k+k) r3+    pokeByteOff p (3*k) r4++-- Extend instances++instance Extend Vec2 Vec3 where+  extendZero   (Vec2 x y) = Vec3 x y 0+  extendWith t (Vec2 x y) = Vec3 x y t+  trim (Vec3 x y _)       = Vec2 x y++instance Extend Vec2 Vec4 where+  extendZero   (Vec2 x y) = Vec4 x y 0 0+  extendWith t (Vec2 x y) = Vec4 x y t t+  trim (Vec4 x y _ _)     = Vec2 x y ++instance Extend Vec3 Vec4 where+  extendZero   (Vec3 x y z) = Vec4 x y z 0+  extendWith t (Vec3 x y z) = Vec4 x y z t+  trim (Vec4 x y z _)       = Vec3 x y z++instance Extend Mat2 Mat3 where+  extendZero (Mat2 p q) = Mat3 (extendZero p) (extendZero q) zero+  extendWith _ _ = error "extendWith is meaningless for matrices"+  trim (Mat3 p q _) = Mat2 (trim p) (trim q)++instance Extend Mat2 Mat4 where+  extendZero (Mat2 p q) = Mat4 (extendZero p) (extendZero q) zero zero+  extendWith _ _ = error "extendWith is meaningless for matrices"+  trim (Mat4 p q _ _) = Mat2 (trim p) (trim q)++instance Extend Mat3 Mat4 where+  extendZero (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) zero+  extendWith _ _ = error "extendWith is meaningless for matrices"+  trim (Mat4 p q r _) = Mat3 (trim p) (trim q) (trim r)+
+ src/flt/GramSchmidt.hs view
@@ -0,0 +1,134 @@++-- | Gram-Schmidt orthogonalization.+-- This module is not re-exported by "Data.Vect".++module Data.Vect.Flt.GramSchmidt +  ( GramSchmidt(..)+  )+  where++import Data.Vect.Flt.Base++-------------------------------------------------------++liftPair :: (a -> b) -> (a,a) -> (b,b)+liftPair f (x,y) = (f x, f y)++liftTriple :: (a -> b) -> (a,a,a) -> (b,b,b)+liftTriple f (x,y,z) = (f x, f y, f z)++liftQuadruple :: (a -> b) -> (a,a,a,a) -> (b,b,b,b)+liftQuadruple f (x,y,z,w) = (f x, f y, f z, f w)++-------------------------------------------------------+    +-- | produces orthogonal\/orthonormal vectors from a set of vectors    +class GramSchmidt a where+  gramSchmidt          :: a -> a   -- ^ does not normalize the vectors!+  gramSchmidtNormalize :: a -> a   -- ^ normalizes the vectors.++{-# RULES+"gramSchmidt is idempotent"  forall a. gramSchmidt (gramSchmidt a) = gramSchmidt a +"gramSchmidtNormalize is idempotent"  forall a. gramSchmidtNormalize (gramSchmidtNormalize a) = gramSchmidtNormalize a +  #-}++-------------------------------------------------------++instance GramSchmidt (Vec2,Vec2) where+  gramSchmidt = gramSchmidtPair+  gramSchmidtNormalize = gramSchmidtNormalizePair+  +instance GramSchmidt (Vec3,Vec3) where+  gramSchmidt = gramSchmidtPair+  gramSchmidtNormalize = gramSchmidtNormalizePair+  +instance GramSchmidt (Vec4,Vec4) where+  gramSchmidt = gramSchmidtPair+  gramSchmidtNormalize = gramSchmidtNormalizePair++----------++instance GramSchmidt (Normal2,Normal2) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal2!"+  gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal++instance GramSchmidt (Normal3,Normal3) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal3!"+  gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal++instance GramSchmidt (Normal4,Normal4) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal4!"+  gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal++----------+  +gramSchmidtPair :: (Vector v, DotProd v) => (v,v) -> (v,v)+gramSchmidtPair (u,v) = (u',v') where +  u' = u+  v' = project v u'     +  +gramSchmidtNormalizePair :: (Vector v, DotProd v) => (v,v) -> (v,v)+gramSchmidtNormalizePair (u,v) = (u',v') where+  u' = normalize u +  v' = normalize $ projectUnsafe v u'     ++----------++instance GramSchmidt (Vec3,Vec3,Vec3) where+  gramSchmidt = gramSchmidtTriple+  gramSchmidtNormalize = gramSchmidtNormalizeTriple+     +instance GramSchmidt (Vec4,Vec4,Vec4) where+  gramSchmidt = gramSchmidtTriple+  gramSchmidtNormalize = gramSchmidtNormalizeTriple++instance GramSchmidt (Normal3,Normal3,Normal3) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal3!"+  gramSchmidtNormalize = liftTriple toNormalUnsafe . gramSchmidtNormalizeTriple . liftTriple fromNormal++instance GramSchmidt (Normal4,Normal4,Normal4) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal4!"+  gramSchmidtNormalize = liftTriple toNormalUnsafe . gramSchmidtNormalizeTriple . liftTriple fromNormal++----------++gramSchmidtTriple :: (Vector v, DotProd v) => (v,v,v) -> (v,v,v)+gramSchmidtTriple (u,v,w) = (u',v',w') where +  u' = u+  v' = project v u'     +  w' = project (project w u') v' +  +gramSchmidtNormalizeTriple :: (Vector v, DotProd v) => (v,v,v) -> (v,v,v)+gramSchmidtNormalizeTriple (u,v,w) = (u',v',w') where+  u' = normalize $ u +  v' = normalize $ projectUnsafe v u'     +  w' = normalize $ projectUnsafe (projectUnsafe w u') v'     ++----------++instance GramSchmidt (Vec4,Vec4,Vec4,Vec4) where+  gramSchmidt          = gramSchmidtQuadruple+  gramSchmidtNormalize = gramSchmidtNormalizeQuadruple ++instance GramSchmidt (Normal4,Normal4,Normal4,Normal4) where+  gramSchmidt          = error "use 'gramSchmidtNormalize' for Normal4!"+  gramSchmidtNormalize = liftQuadruple toNormalUnsafe . gramSchmidtNormalizeQuadruple . liftQuadruple fromNormal++----------+  +gramSchmidtQuadruple :: (Vector v, DotProd v) => (v,v,v,v) -> (v,v,v,v)+gramSchmidtQuadruple (u,v,w,z) = (u',v',w',z') where +  u' = u+  v' = project v u'     +  w' = project (project w u') v' +  z' = project (project (project z u') v') w'++gramSchmidtNormalizeQuadruple :: (Vector v, DotProd v) => (v,v,v,v) -> (v,v,v,v)+gramSchmidtNormalizeQuadruple (u,v,w,z) = (u',v',w',z') where+  u' = normalize $ u+  v' = normalize $ projectUnsafe v u'     +  w' = normalize $ projectUnsafe (projectUnsafe w u') v' +  z' = normalize $ projectUnsafe (projectUnsafe (projectUnsafe z u') v') w'+  +----------+
+ src/flt/Interpolate.hs view
@@ -0,0 +1,42 @@++-- | Interpolation of vectors. +-- Note: we interpolate unit vectors differently from ordinary vectors.++module Data.Vect.Flt.Interpolate where++import Data.Vect.Flt.Base+import Data.Vect.Flt.Util.Dim2 (sinCos',angle2')+import Data.Vect.Flt.Util.Dim3 (rotate3')++class Interpolate v where+  interpolate :: Flt -> v -> v -> v+  +instance Interpolate Flt where+  interpolate t x y = x + t*(y-x)++instance Interpolate Vec2 where interpolate t x y = x &+ t *& (y &- x)+instance Interpolate Vec3 where interpolate t x y = x &+ t *& (y &- x)+instance Interpolate Vec4 where interpolate t x y = x &+ t *& (y &- x)++instance Interpolate Normal2 where+  interpolate t nx ny = sinCos' $ ax + t*adiff where+    ax = angle2' nx+    ay = angle2' ny+    adiff = helper (ay - ax)+    helper d +      | d < -pi   = d + twopi+      | d >  pi   = d - twopi+      | otherwise = d+    twopi = 2*pi+    +instance Interpolate Normal3 where +  interpolate t nx ny = +    if maxAngle < 0.001  -- more or less ad-hoc critical angle+      then mkNormal $ interpolate t x y+      else toNormalUnsafe $ rotate3' (t*maxAngle) (mkNormal axis) x where+    x = fromNormal nx+    y = fromNormal ny+    axis = (x &^ y)+    maxAngle = acos (x &. y)+        +
+ src/flt/OpenGL.hs view
@@ -0,0 +1,182 @@++-- TODO: the pointer versions of these functions should be really implemented +-- via the pointer versions of the original opengl functions...++-- | OpenGL support, inclduing 'vertex', 'texCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.+ +module Data.Vect.Flt.OpenGL where++import Control.Monad+import Data.Vect.Flt.Base+import qualified Graphics.Rendering.OpenGL as GL++import Foreign++import Graphics.Rendering.OpenGL hiding (Normal3,rotate,translate,scale)++-------------------------------------------------------++{-# SPECIALISE radianToDegrees :: Float -> Float #-}+{-# SPECIALISE radianToDegrees :: Double -> Double #-}+radianToDegrees :: RealFrac a => a -> a+radianToDegrees x = x * 57.295779513082322++{-# SPECIALIZE degreesToRadian :: Float  -> Float  #-}+{-# SPECIALIZE degreesToRadian :: Double -> Double #-}+degreesToRadian :: Floating a => a -> a+degreesToRadian x = x * 1.7453292519943295e-2++-- | The angle is in radians. (WARNING: OpenGL uses degrees!)+rotate :: Flt -> Vec3 -> IO ()+rotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z)++translate :: Vec3 -> IO ()+translate (Vec3 x y z) = GL.translate (Vector3 x y z)++scale3 :: Vec3 -> IO ()+scale3 (Vec3 x y z) = GL.scale x y z++scale :: Flt -> IO ()+scale x = GL.scale x x x++-------------------------------------------------------++-- Vertex instances++instance GL.Vertex Vec2 where+  vertex (Vec2 x y) = GL.vertex (GL.Vertex2 x y)+  vertexv p = peek p >>= vertex +  +instance GL.Vertex Vec3 where+  vertex (Vec3 x y z) = GL.vertex (GL.Vertex3 x y z)+  vertexv p = peek p >>= vertex   +  +instance GL.Vertex Vec4 where+  vertex (Vec4 x y z w) = GL.vertex (GL.Vertex4 x y z w)+  vertexv p = peek p >>= vertex   ++-------------------------------------------------------++-- the Normal instance+-- note that there is no Normal2\/Normal4 in the OpenGL binding++instance GL.Normal Normal3 where+  normal (Normal3 (Vec3 x y z)) = GL.normal (GL.Normal3 x y z)+  normalv p = peek p >>= normal ++-------------------------------------------------------++-- Color instances+  +instance GL.Color Vec3 where+  color (Vec3 r g b) = GL.color (GL.Color3 r g b)+  colorv p = peek p >>= color++instance GL.Color Vec4 where+  color (Vec4 r g b a) = GL.color (GL.Color4 r g b a)+  colorv p = peek p >>= color++instance GL.SecondaryColor Vec3 where+  secondaryColor (Vec3 r g b) = GL.secondaryColor (GL.Color3 r g b)+  secondaryColorv p = peek p >>= secondaryColor++{-+-- there is no such thing?+instance GL.SecondaryColor Vec4 where+  secondaryColor (Vec4 r g b a) = GL.secondaryColor (GL.Color4 r g b a)+  secondaryColorv p = peek p >>= secondaryColor+-}++-------------------------------------------------------++-- TexCoord instances++instance GL.TexCoord Vec2 where+  texCoord (Vec2 u v) = GL.texCoord (GL.TexCoord2 u v)+  texCoordv p = peek p >>= texCoord+  multiTexCoord unit (Vec2 u v) = GL.multiTexCoord unit (GL.TexCoord2 u v)+  multiTexCoordv unit p = peek p >>= multiTexCoord unit++instance GL.TexCoord Vec3 where+  texCoord (Vec3 u v w) = GL.texCoord (GL.TexCoord3 u v w)+  texCoordv p = peek p >>= texCoord+  multiTexCoord unit (Vec3 u v w) = GL.multiTexCoord unit (GL.TexCoord3 u v w)+  multiTexCoordv unit p = peek p >>= multiTexCoord unit++instance GL.TexCoord Vec4 where+  texCoord (Vec4 u v w z) = GL.texCoord (GL.TexCoord4 u v w z)+  texCoordv p = peek p >>= texCoord+  multiTexCoord unit (Vec4 u v w z) = GL.multiTexCoord unit (GL.TexCoord4 u v w z)+  multiTexCoordv unit p = peek p >>= multiTexCoord unit++-------------------------------------------------------+    +-- Vertex Attributes (experimental)++class VertexAttrib' a where+  vertexAttrib :: GL.AttribLocation -> a -> IO ()+  +instance VertexAttrib' {- ' CPP is sensitive to primes -} Flt where+  vertexAttrib loc x = GL.vertexAttrib1 loc x++instance VertexAttrib' Vec2 where+  vertexAttrib loc (Vec2 x y) = GL.vertexAttrib2 loc x y++instance VertexAttrib' Vec3 where+  vertexAttrib loc (Vec3 x y z) = GL.vertexAttrib3 loc x y z ++instance VertexAttrib' Vec4 where+  vertexAttrib loc (Vec4 x y z w) = GL.vertexAttrib4 loc x y z w ++instance VertexAttrib' Normal2 where+  vertexAttrib loc (Normal2 (Vec2 x y)) = GL.vertexAttrib2 loc x y++instance VertexAttrib' Normal3 where+  vertexAttrib loc (Normal3 (Vec3 x y z)) = GL.vertexAttrib3 loc x y z++instance VertexAttrib' Normal4 where+  vertexAttrib loc (Normal4 (Vec4 x y z w)) = GL.vertexAttrib4 loc x y z w++-------------------------------------------------------++-- Uniform (again, experimental)++-- (note that the uniform location code in the OpenGL 2.2.1.1 is broken; +-- a work-around is to put a zero character at the end of uniform names)++{-+toFloat :: Flt -> Float+toFloat = realToFrac++fromFloat :: Float -> Flt+fromFloat = realToFrac+-}++-- Uniforms are always floats...+#ifdef VECT_Float++instance GL.Uniform Flt where+  uniform loc = GL.makeStateVar getter setter where+    getter = liftM (\(GL.Index1 x) -> x) $ get (uniform loc)+    setter x = ($=) (uniform loc) (Index1 x) +  uniformv loc cnt ptr = uniformv loc cnt (castPtr ptr :: Ptr (Index1 Flt))++instance GL.Uniform Vec2 where+  uniform loc = GL.makeStateVar getter setter where+    getter = liftM (\(GL.Vertex2 x y) -> Vec2 x y) $ get (uniform loc)+    setter (Vec2 x y) = ($=) (uniform loc) (Vertex2 x y) +  uniformv loc cnt ptr = uniformv loc (2*cnt) (castPtr ptr :: Ptr Flt)++instance GL.Uniform Vec3 where+  uniform loc = GL.makeStateVar getter setter where+    getter = liftM (\(GL.Vertex3 x y z) -> Vec3 x y z) $ get (uniform loc)+    setter (Vec3 x y z) = ($=) (uniform loc) (Vertex3 x y z) +  uniformv loc cnt ptr = uniformv loc (3*cnt) (castPtr ptr :: Ptr Flt)++instance GL.Uniform Vec4 where+  uniform loc = GL.makeStateVar getter setter where+    getter = liftM (\(GL.Vertex4 x y z w) -> Vec4 x y z w) $ get (uniform loc)+    setter (Vec4 x y z w) = ($=) (uniform loc) (Vertex4 x y z w) +  uniformv loc cnt ptr = uniformv loc (4*cnt) (castPtr ptr :: Ptr Flt)+    +#endif
+ src/flt/Util/Dim2.hs view
@@ -0,0 +1,64 @@++module Data.Vect.Flt.Util.Dim2 where++import Data.Vect.Flt.Base++-- |example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@.+structVec2 :: [Flt] -> [Vec2]+structVec2 [] = []+structVec2 (x:y:ls) = (Vec2 x y):(structVec2 ls) +structVec2 _ = error "structVec2"++destructVec2 :: [Vec2] -> [Flt]+destructVec2 [] = []+destructVec2 ((Vec2 x y):ls) = x:y:(destructVec2 ls)  ++det2 :: Vec2 -> Vec2 -> Flt+det2 u v = det (u,v)++vec2X :: Vec2+vec2Y :: Vec2++vec2X = Vec2 1 0 +vec2Y = Vec2 0 1 ++translate2X :: Flt -> Vec2 -> Vec2+translate2Y :: Flt -> Vec2 -> Vec2++translate2X t (Vec2 x y) = Vec2 (x+t) y +translate2Y t (Vec2 x y) = Vec2 x (y+t) ++-- | unit vector with given angle relative to the positive X axis (in the positive direction, that is, CCW).+-- A more precise name would be @cosSin@, but that sounds bad :)+sinCos :: Flt -> Vec2+sinCos a = Vec2 (cos a) (sin a)++sinCos' {- ' CPP is sensitive to primes -} :: Flt -> Normal2+sinCos' = toNormalUnsafe . sinCos++sinCosRadius :: Flt    -- ^ angle (in radians)+             -> Flt    -- ^ radius+             -> Vec2+sinCosRadius a r = Vec2 (r * cos a) (r * sin a)++-- | The angle relative to the positive X axis+angle2 :: Vec2 -> Flt+angle2 (Vec2 x y) = atan2 y x++angle2' {- ' CPP is sensitive to primes -} :: Normal2 -> Flt+angle2' = angle2 . fromNormal++-- |Rotation matrix by a given angle (in radians), counterclockwise.+rotMatrix2 :: Flt -> Mat2+rotMatrix2 a = Mat2 (Vec2 c s) (Vec2 (-s) c) where c = cos a; s = sin a++rotate2 :: Flt -> Vec2 -> Vec2+rotate2 a v = v .* (rotMatrix2 a) ++-- |Rotates counterclockwise by 90 degrees.+rotateCCW :: Vec2 -> Vec2+rotateCCW (Vec2 x y) = Vec2 (-y) x++-- |Rotates clockwise by 90 degrees.+rotateCW :: Vec2 -> Vec2+rotateCW (Vec2 x y) = Vec2 y (-x)
+ src/flt/Util/Dim3.hs view
@@ -0,0 +1,75 @@++module Data.Vect.Flt.Util.Dim3 where++import Data.Vect.Flt.Base++structVec3 :: [Flt] -> [Vec3]+structVec3 [] = []+structVec3 (x:y:z:ls) = (Vec3 x y z):(structVec3 ls) +structVec3 _ = error "structVec3"++destructVec3 :: [Vec3] -> [Flt]+destructVec3 [] = []+destructVec3 ((Vec3 x y z):ls) = x:y:z:(destructVec3 ls)  ++det3 :: Vec3 -> Vec3 -> Vec3 -> Flt+det3 u v w = det (u,v,w)++translate3X :: Flt -> Vec3 -> Vec3+translate3Y :: Flt -> Vec3 -> Vec3+translate3Z :: Flt -> Vec3 -> Vec3++translate3X t (Vec3 x y z) = Vec3 (x+t) y z +translate3Y t (Vec3 x y z) = Vec3 x (y+t) z +translate3Z t (Vec3 x y z) = Vec3 x y (z+t) ++vec3X :: Vec3+vec3Y :: Vec3+vec3Z :: Vec3++vec3X = Vec3 1 0 0+vec3Y = Vec3 0 1 0+vec3Z = Vec3 0 0 1++rotMatrixZ :: Flt -> Mat3+rotMatrixY :: Flt -> Mat3+rotMatrixX :: Flt -> Mat3++-- These are intended for multiplication on the /right/.+-- Should be consistent with the rotation around an arbitrary axis +-- (eg, @rotMatrixY a == rotate3 a vec3Y@)+rotMatrixZ a = Mat3 (Vec3 c s 0) (Vec3 (-s) c 0) (Vec3 0 0 1) where c = cos a; s = sin a+rotMatrixY a = Mat3 (Vec3 c 0 (-s)) (Vec3 0 1 0) (Vec3 s 0 c) where c = cos a; s = sin a+rotMatrixX a = Mat3 (Vec3 1 0 0) (Vec3 0 c s) (Vec3 0 (-s) c) where c = cos a; s = sin a++rotate3' :: {- ' CPP is sensitive to primes -} Flt       -- ^ angle (in radians)+         -> Normal3   -- ^ axis (should be a /unit/ vector!) +         -> Vec3      -- ^ vector+         -> Vec3      -- ^ result+rotate3' angle axis v = v .* (rotMatrix3' axis angle)++rotate3 :: Flt    -- ^ angle (in radians)+        -> Vec3   -- ^ axis (arbitrary nonzero vector)+        -> Vec3   -- ^ vector+        -> Vec3   -- ^ result+rotate3 angle axis v = v .* (rotMatrix3 axis angle)+      +-- |Rotation around an arbitrary 3D vector. The resulting 3x3 matrix is intended for multiplication on the /right/. +rotMatrix3 :: Vec3 -> Flt -> Mat3+rotMatrix3 v a = rotMatrix3' (mkNormal v) a++-- |Rotation around an arbitrary 3D /unit/ vector. The resulting 3x3 matrix is intended for multiplication on the /right/. +rotMatrix3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Mat3+rotMatrix3' (Normal3 v) a = +  let c = cos a+      s = sin a+      m1 = scalarMul (1-c) (outer v v)+      x = _1 v+      y = _2 v+      z = _3 v+      m2 = Mat3 (Vec3   c    ( s*z) (-s*y))+                (Vec3 (-s*z)   c    ( s*x))+                (Vec3 ( s*y) (-s*x)   c   )+  in (m1 &+ m2)++
+ src/flt/Util/Dim4.hs view
@@ -0,0 +1,92 @@++-- | Rotation around an arbitrary plane in four dimensions, and other miscellanea.+-- Not very useful for most people, and not re-exported by "Data.Vect".++module Data.Vect.Flt.Util.Dim4 where++import Data.Vect.Flt.Base+import Data.Vect.Flt.GramSchmidt++structVec4 :: [Flt] -> [Vec4]+structVec4 [] = []+structVec4 (x:y:z:w:ls) = (Vec4 x y z w):(structVec4 ls) +structVec4 _ = error "structVec4"++destructVec4 :: [Vec4] -> [Flt]+destructVec4 [] = []+destructVec4 ((Vec4 x y z w):ls) = x:y:z:w:(destructVec4 ls)  ++--det4 :: Vec4 -> Vec4 -> Vec4 -> Vec4 -> Flt+--det4 u v w z = det (u,v,w,z)++translate4X :: Flt -> Vec4 -> Vec4+translate4Y :: Flt -> Vec4 -> Vec4+translate4Z :: Flt -> Vec4 -> Vec4+translate4W :: Flt -> Vec4 -> Vec4++translate4X t (Vec4 x y z w) = Vec4 (x+t) y z w +translate4Y t (Vec4 x y z w) = Vec4 x (y+t) z w +translate4Z t (Vec4 x y z w) = Vec4 x y (z+t) w+translate4W t (Vec4 x y z w) = Vec4 x y z (w+t) ++vec4X :: Vec4+vec4Y :: Vec4+vec4Z :: Vec4+vec4W :: Vec4++vec4X = Vec4 1 0 0 0+vec4Y = Vec4 0 1 0 0+vec4Z = Vec4 0 0 1 0+vec4W = Vec4 0 0 0 1++---------------------------------------------------------------------------++-- |If @(x,y,u,v)@ is an orthonormal system, then (written in pseudo-code)+-- @biVector4 (x,y) = plusMinus (reverse $ biVector4 (u,v))@.+-- This is a helper function for the 4 dimensional rotation code.+-- If @(x,y,z,p,q,r) = biVector4 a b@, then the corresponding antisymmetric tensor is+--+-- > [  0  r  q  p ]+-- > [ -r  0  z -y ]+-- > [ -q -z  0  x ]+-- > [ -p  y -x  0 ]+biVector4 :: Vec4 -> Vec4 -> (Flt,Flt,Flt,Flt,Flt,Flt)+biVector4 (Vec4 x y z w) (Vec4 a b c d) = +  ( x*b-y*a , x*c-z*a , x*d-w*a , y*c-z*b , -y*d+w*b , z*d-w*c )++-- | the corresponding antisymmetric tensor+biVector4AsTensor :: Vec4 -> Vec4 -> Mat4+biVector4AsTensor v w = +  Mat4 ( Vec4   0  ( r) ( q) ( p) )+       ( Vec4 (-r)   0  ( z) (-y) )+       ( Vec4 (-q) (-z)   0  ( x) )+       ( Vec4 (-p) ( y) (-x)   0  )+  where +    (x,y,z,p,q,r) = biVector4 v w++-- | We assume that the axes are normalized and /orthogonal/ to each other!+rotate4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Vec4 -> Vec4+rotate4' angle axes v = v .* (rotMatrix4' angle axes)++-- | We assume only that the axes are independent vectors.+rotate4 :: Flt -> (Vec4,Vec4) -> Vec4 -> Vec4+rotate4 angle axes v = v .* (rotMatrix4 angle axes)++-- | Rotation matrix around a plane specified by two normalized and /orthogonal/ vectors.+-- Intended for multiplication on the /right/!+rotMatrix4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Mat4+rotMatrix4' angle (Normal4 v, Normal4 w) = m1 &+ (s *& m2) &+ m3 +  where+    c = cos angle ; s = sin angle+    m1 = scalarMul (1-c) ( outer v v  &+  outer w w )+    m2 = biVector4AsTensor v w+    m3 = diag (Vec4 c c c c)++-- | We assume only that the axes are independent vectors.+rotMatrix4 :: Flt -> (Vec4,Vec4) -> Mat4  +rotMatrix4 angle axes = +  rotMatrix4' angle $ liftPair toNormalUnsafe $ gramSchmidtNormalize axes +  where +    liftPair f (x,y) = (f x, f y)+    +
+ src/flt/Util/Projective.hs view
@@ -0,0 +1,87 @@++-- | Classic 4x4 projective matrices. Our convention is that they are intended for multiplication on+-- the /right/, that is, they are of the form+--+-- >     _____+-- > [  |     |  0  ]+-- > [  | 3x3 |  0  ]+-- > [  |_____|  0  ]+-- > [  p  q  r  1  ]+--+-- Please note that by default, OpenGL stores the matrices (in memory) by columns, while we +-- store them by rows; but OpenGL also use the opposite convention (so the OpenGL projective matrices +-- are intended for multiplication on the /left/). So in effect, they are the same when stored in the memory,+-- say with @poke :: Ptr Mat4 -> Mat4 -> IO ()@.++module Data.Vect.Flt.Util.Projective where++import Data.Vect.Flt.Base+import Data.Vect.Flt.Util.Dim3++import qualified Data.Vect.Flt.Util.Dim4 as Dim4++class ExtendProjective v e | v->e where+  extendProj     :: v -> e+  extendProjWith :: Flt -> v -> e+  extendProj = extendProjWith 1+  +instance ExtendProjective Vec2 Vec4 where+  extendProj       (Vec2 x y) = Vec4 x y 0 1+  extendProjWith w (Vec2 x y) = Vec4 x y 0 w+  +instance ExtendProjective Vec3 Vec4 where+  extendProj       (Vec3 x y z) = Vec4 x y z 1+  extendProjWith w (Vec3 x y z) = Vec4 x y z w++instance ExtendProjective Vec4 Vec4 where+  extendProj = id+  extendProjWith w (Vec4 x y z w') = let s = w/w' in Vec4 (s*x) (s*y) (s*z) w++instance ExtendProjective Mat2 Mat4 where+  extendProj       (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 1)+  extendProjWith w (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 w)++instance ExtendProjective Mat3 Mat4 where+  extendProj       (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 1)+  extendProjWith w (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 w)++rotMatrixProj :: Flt -> Normal3 -> Mat4+rotMatrixProj angle axis = extendProj $ rotMatrix3' axis angle++rotMatrixProj' :: {- ' CPP is sensitive to primes -} Flt -> Vec3 -> Mat4+rotMatrixProj' angle axis = extendProj $ rotMatrix3 axis angle++translMatrixProj :: Vec3 -> Mat4+translMatrixProj v = Mat4 Dim4.vec4X Dim4.vec4Y Dim4.vec4Z (extendProj v)++-- | we assume that the bottom-right corner is 1.+translWithProj :: Vec3 -> Mat4 -> Mat4+translWithProj v mat@(Mat4 r1 r2 r3 r4) = Mat4 r1 r2 r3 (extendProjWith 0 v &+ r4)++scaleMatrixProj :: Vec3 -> Mat4+scaleMatrixProj v = diag $ extendProj v++scaleMatrixUniformProj :: Flt -> Mat4+scaleMatrixUniformProj s = diag (Vec4 s s s 1)++class ProjectiveAction v where+  actProj :: v -> Mat4 -> v+ +instance ProjectiveAction Vec3 where+  actProj v m = trim $ (extendProj v) .* m ++instance ProjectiveAction Vec4 where+  actProj v m = v .* m ++-- | When acting on unit vectors, we ignore the translation part.+instance ProjectiveAction Normal3 where+  actProj (Normal3 v) m = Normal3 (v .* (trim m :: Mat3))++-- | Inverts a projective 4x4 matrix, assuming that the top-left 3x3 part is /orthogonal/,+-- and the bottom-right corner is 1.+invertProj :: Mat4 -> Mat4+invertProj mat@(Mat4 u v w t) = +  translWithProj t' $ extendProj $ transpose $ (trim mat :: Mat3)+  where+    t' = Vec3 (- u &. t) (- v &. t) (- w &. t)+
+ vect.cabal view
@@ -0,0 +1,78 @@+Name:                vect+Version:             0.4.0+Synopsis:            A low-dimensional linear algebra library, tailored to computer graphics.+Description:         A low-dimensional (2, 3 and 4) linear algebra library, +                     with lots of useful functions. Intended usage is primarily +                     computer graphics (basic OpenGL support is included).+                     Projective 4 dimensional operations, as used in eg. +                     OpenGL, are also supported.+                     The base field is either Float or Double.+License:             BSD3+License-file:        LICENSE+Author:              Balazs Komuves+Copyright:           (c) 2008-2009 Balazs Komuves+Maintainer:          bkomuves (plus) hackage (at) gmail (dot) com+Homepage:            http://code.haskell.org/~bkomuves/+Stability:           Experimental+Category:            Graphics, Math+Tested-With:         GHC == 6.10.1+Cabal-Version:       >= 1.6+Build-Type:          Custom++extra-source-files:  src/flt.hs,+                     src/flt/Base.hs,+                     src/flt/Interpolate.hs,+                     src/flt/OpenGL.hs,+                     src/flt/GramSchmidt.hs,+                     src/flt/Util/Dim2.hs,+                     src/flt/Util/Dim3.hs,+                     src/flt/Util/Dim4.hs,+                     src/flt/Util/Projective.hs+                     +Flag splitBase+  Description: Choose the new smaller, split-up base package.++Flag OpenGL+  Description: Compile with OpenGL support +  Default:     True  +  +Library+  if flag(splitBase)+    Build-Depends:       base >= 3, random +  else+    Build-Depends:       base <  3++  if flag(OpenGL)+    Build-Depends:       OpenGL+    cpp-options:         -DVECT_OPENGL+    Exposed-Modules:     Data.Vect.Float.OpenGL,                         +                         Data.Vect.Double.OpenGL                         ++  Exposed-Modules:     Data.Vect, +  +                       Data.Vect.Float,+                       Data.Vect.Float.Base, +                       Data.Vect.Float.Interpolate,+                       Data.Vect.Float.GramSchmidt,+                       Data.Vect.Float.Util.Dim2, +                       Data.Vect.Float.Util.Dim3, +                       Data.Vect.Float.Util.Dim4,+                       Data.Vect.Float.Util.Projective,++                       Data.Vect.Double,+                       Data.Vect.Double.Base, +                       Data.Vect.Double.Interpolate,+                       Data.Vect.Double.GramSchmidt,+                       Data.Vect.Double.Util.Dim2, +                       Data.Vect.Double.Util.Dim3, +                       Data.Vect.Double.Util.Dim4,+                       Data.Vect.Double.Util.Projective++  Hs-Source-Dirs:      .+  Extensions:          ForeignFunctionInterface, CPP,+                       MultiParamTypeClasses, FunctionalDependencies,+                       FlexibleInstances, TypeSynonymInstances,+                       GeneralizedNewtypeDeriving++  ghc-options:         -Wall+