packages feed

linear (empty) → 0.2

raw patch · 21 files changed

+1178/−0 lines, 21 filesdep +basedep +directorydep +distributivesetup-changed

Dependencies added: base, directory, distributive, doctest, filepath, lens

Files

+ .ghci view
@@ -0,0 +1,1 @@+:set -isrc -idist/build/autogen -optP-include -optPdist/build/autogen/cabal_macros.h
+ .gitignore view
@@ -0,0 +1,5 @@+dist+docs+wiki+TAGS+tags
+ .travis.yml view
@@ -0,0 +1,8 @@+language: haskell+# Uncomment the next 4 lines whenever hackage is down.+# before_install:+#   - mkdir -p ~/.cabal+#   - cp config ~/.cabal/config+#   - cabal update+notifications:+  irc: "irc.freenode.org#haskell-lens"
+ .vim.custom view
@@ -0,0 +1,21 @@+" Add the following to your .vimrc to automatically load this on startup+" if filereadable(".vim.custom")+"     so .vim.custom+" endif++function StripTrailingWhitespace()+  let myline=line(".")+  let mycolumn = col(".")+  silent %s/  *$//+  call cursor(myline, mycolumn)+endfunction++syntax on+set tags=TAGS;/+set listchars=tab:‗‗,trail:‗+set list++map <F2> :exec ":!hasktags -x -c --ignore src"<CR><CR>++au BufWritePre *.hs,*.markdown silent! cal StripTrailingWhitespace()+au BufWritePost *.hs silent! :exec ":!hasktags -x -c --ignore src"
+ CHANGELOG.markdown view
@@ -0,0 +1,3 @@+0.2+---+* Initial hackage release
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright 2011-12 Edward Kmett++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++1. Redistributions of source code must retain the above copyright+   notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright+   notice, this list of conditions and the following disclaimer in the+   documentation and/or other materials provided with the distribution.++3. Neither the name of the author nor the names of his contributors+   may be used to endorse or promote products derived from this software+   without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``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 AUTHORS 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.
+ README.markdown view
@@ -0,0 +1,13 @@+linear+======++[![Build Status](https://secure.travis-ci.org/ekmett/linear.png?branch=master)](http://travis-ci.org/ekmett/linear)++Contact Information+-------------------++Contributions and bug reports are welcome!++Please feel free to contact me through github or on the #haskell IRC channel on irc.freenode.net.++-Edward Kmett
+ Setup.lhs view
@@ -0,0 +1,7 @@+#!/usr/bin/runhaskell+> module Main (main) where++> import Distribution.Simple++> main :: IO ()+> main = defaultMain
+ config view
@@ -0,0 +1,16 @@+-- This provides a custom ~/.cabal/config file for use when hackage is down that should work on unix+--+-- This is particularly useful for travis-ci to get it to stop complaining+-- about a broken build when everything is still correct on our end.+--+-- This uses Luite Stegeman's mirror of hackage provided by his 'hdiff' site instead+--+-- To enable this, uncomment the before_script in .travis.yml++remote-repo: hdiff.luite.com:http://hdiff.luite.com/packages/archive+remote-repo-cache: ~/.cabal/packages+world-file: ~/.cabal/world+build-summary: ~/.cabal/logs/build.log+remote-build-reporting: anonymous+install-dirs user+install-dirs global
+ linear.cabal view
@@ -0,0 +1,62 @@+name:          linear+category:      Math, Algebra+version:       0.2+license:       BSD3+cabal-version: >= 1.8+license-file:  LICENSE+author:        Edward A. Kmett+maintainer:    Edward A. Kmett <ekmett@gmail.com>+stability:     provisional+homepage:      http://github.com/ekmett/linear/+bug-reports:   http://github.com/ekmett/linear/issues+copyright:     Copyright (C) 2012 Edward A. Kmett+synopsis:      Linear Algebra+description:   Types and combinators for low-dimension-count linear algebra on free vectors spaces+build-type:    Simple+tested-with:   GHC == 7.4.1+extra-source-files:+  .travis.yml+  .ghci+  .gitignore+  .vim.custom+  config+  README.markdown+  CHANGELOG.markdown++source-repository head+  type: git+  location: git://github.com/ekmett/linear.git++library+  build-depends:+    base == 4.*,+    distributive >= 0.2.2 && < 0.3,+    lens == 2.9.*++  exposed-modules:+    Linear.Conjugate+    Linear.Epsilon+    Linear.Matrix+    Linear.Metric+    Linear.Plucker+    Linear.Quaternion+    Linear.V2+    Linear.V3+    Linear.V4+    Linear.Vector++  ghc-options: -Wall -fwarn-tabs -O2 -fdicts-cheap -funbox-strict-fields+  hs-source-dirs: src++-- Verify the results of the examples+test-suite doctests+  type:    exitcode-stdio-1.0+  main-is: doctests.hs+  build-depends:+    base == 4.*,+    directory >= 1.0 && < 1.2,+    doctest >= 0.8 && <= 0.9,+    filepath >= 1.3 && < 1.4+  ghc-options: -Wall -Werror -threaded+  hs-source-dirs: tests+
+ src/Linear/Conjugate.hs view
@@ -0,0 +1,18 @@+module Linear.Conjugate+  ( Conjugate(..)+  ) where++import Data.Complex hiding (conjugate)++-- | An involutive ring+class Num a => Conjugate a where+  -- | Conjugate a value. This defaults to the trivial involution.+  conjugate :: a -> a+  conjugate = id++instance Conjugate Double+instance Conjugate Float+instance (Conjugate a, RealFloat a) => Conjugate (Complex a) where+  {-# SPECIALIZE instance Conjugate (Complex Float) #-}+  {-# SPECIALIZE instance Conjugate (Complex Double) #-}+  conjugate (a :+ b) = conjugate a :+ negate b
+ src/Linear/Epsilon.hs view
@@ -0,0 +1,26 @@+-----------------------------------------------------------------------------+-- |+-- Module      :  Linear.Epsilon+-- Copyright   :  (C) 2012 Edward Kmett+-- License     :  BSD-style (see the file LICENSE)+-- Maintainer  :  Edward Kmett <ekmett@gmail.com>+-- Stability   :  provisional+-- Portability :  portable+--+-----------------------------------------------------------------------------+module Linear.Epsilon+  ( Epsilon(..)+  ) where++-- | Provides a fairly subjective test to see if a quantity is near zero.+class Num a => Epsilon a where+  -- | Determine if a quantity is near zero.+  nearZero :: a -> Bool++-- | @'abs' a '<=' 1e-6@+instance Epsilon Float where+  nearZero a = abs a <= 1e-6++-- | @'abs' a '<=' 1e-12@+instance Epsilon Double where+  nearZero a = abs a <= 1e-12
+ src/Linear/Matrix.hs view
@@ -0,0 +1,95 @@+module Linear.Matrix+  ( (!*!), (!*) , (*!)+  , adjoint+  , M33, M44, M43, m33_to_m44, m43_to_m44+  , eye3, eye4+  , trace+  , translation+  , fromQuaternion+  , mkTransformation+  ) where++import Control.Applicative+import Control.Lens+import Data.Distributive+import Data.Foldable as Foldable+import Linear.Quaternion+import Linear.V3+import Linear.V4+import Linear.Metric+import Linear.Conjugate++infixl 7 !*!+-- | matrix product+(!*!) :: (Functor m, Foldable r, Applicative r, Distributive n, Num a) => m (r a) -> r (n a) -> m (n a)+f !*! g = fmap (\r -> Foldable.foldr (+) 0 . liftA2 (*) r <$> g') f+  where g' = distribute g++-- | matrix * column vector+infixl 7 *!+(!*) :: (Functor m, Metric r, Num a) => m (r a) -> r a -> m a+m !* v = dot v <$> m++infixl 7 !*+-- | row vector * matrix+(*!) :: (Metric r, Distributive n, Num a) => r a -> r (n a) -> n a+f *! g = dot f <$> distribute g++-- | hermitian conjugate or conjugate transpose+adjoint :: (Functor m, Distributive n, Conjugate a) => m (n a) -> n (m a)+adjoint = collect (fmap conjugate)+{-# INLINE adjoint #-}++-- | Compute the trace of a matrix+trace :: (Monad f, Foldable f, Num a) => f (f a) -> a+trace m = Foldable.sum (m >>= id)+{-# INLINE trace #-}++-- | Matrices use a row-major representation.+type M33 a = V3 (V3 a)+type M44 a = V4 (V4 a)+type M43 a = V4 (V3 a)++-- | Build a rotation matrix from a 'Quaternion'.+fromQuaternion :: Num a => Quaternion a -> M33 a+fromQuaternion (Quaternion a (V3 b c d)) =+  V3 (V3 (a*a+b*b-c*c-d*d) (2*b*c-2*a*d) (2*b*d+2*a*c))+     (V3 (2*b*c+2*a*d) (a*a-b*b+c*c-d*d) (2*c*d-2*a*b))+     (V3 (2*b*d-2*a*c) (2*c*d+2*a*b) (a*a-b*b-c*c+d*d))++mkTransformationMat :: Num a => M33 a -> V3 a -> M44 a+mkTransformationMat (V3 r1 r2 r3) (V3 tx ty tz) =+  V4 (snoc3 r1 tx) (snoc3 r2 ty) (snoc3 r3 tz) (set _w 1 0)+  where snoc3 (V3 x y z) w = V4 x y z w++-- |Build a transformation matrix from a rotation expressed as a+-- 'Quaternion' and a translation vector.+mkTransformation :: Num a => Quaternion a -> V3 a -> M44 a+mkTransformation = mkTransformationMat . fromQuaternion++m43_to_m44 :: Num a => M43 a -> M44 a+m43_to_m44+  (V4 (V3 a b c)+      (V3 d e f)+      (V3 g h i)+      (V3 j k l)) =+  (V4 (V4 a b c 0)+      (V4 d e f 0)+      (V4 g h i 0)+      (V4 j k l 1))++m33_to_m44 :: Num a => M33 a -> M44 a+m33_to_m44 (V3 r1 r2 r3) = V4 (vector r1) (vector r2) (vector r3) (point 0)++-- |3x3 identity matrix.+eye3 :: Num a => M33 a+eye3 = V3 (set _x 1 0) (set _y 1 0) (set _z 1 0)++-- |4x4 identity matrix.+eye4 :: Num a => M44 a+eye4 = V4 (set _x 1 0) (set _y 1 0) (set _z 1 0) (set _w 1 0)++-- |Extract the translation vector (first three entries of the last+-- column) from a 3x4 or 4x4 matrix+translation :: (R3 t, R4 v, Functor f, Functor t) => (V3 a -> f (V3 a)) -> t (v a) -> f (t a)+translation = (. fmap (^._w)) . _xyz
+ src/Linear/Metric.hs view
@@ -0,0 +1,40 @@+module Linear.Metric+  ( Metric(..), normalize+  ) where++import Control.Applicative+import Linear.Epsilon++-- | A free inner product/metric space+class Applicative f => Metric f where+  -- | Compute the inner product of two vectors or (equivalently)+  -- convert a vector @f a@ into a covector @f a -> a@.+  dot :: Num a => f a -> f a -> a++  -- | Compute the squared norm. The name quadrance arises from+  -- Norman J. Wildberger's rational trigonometry.+  quadrance :: Num a => f a -> a+  quadrance v = dot v v++  -- | Compute the quadrance of the difference+  qd :: Num a => f a -> f a -> a+  qd f g = quadrance (liftA2 (-) f g)++  -- | Compute the distance between two vectors in a metric space+  distance :: Floating a => f a -> f a -> a+  distance f g = norm (liftA2 (-) f g)++  -- | Compute the norm of a vector in a metric space+  norm :: Floating a => f a -> a+  norm v = sqrt (dot v v)++  -- | Convert a non-zero vector to unit vector.+  signorm :: Floating a => f a -> f a+  signorm v = fmap (/m) v where+    m = norm v++-- | Normalize a 'Metric' functor to have unit 'norm'. This function+-- does not change the functor if its 'norm' is 0 or 1.+normalize :: (Floating a, Metric f, Epsilon a) => f a -> f a+normalize v = if nearZero l || nearZero (1-l) then v else fmap (/sqrt l) v+  where l = quadrance v
+ src/Linear/Plucker.hs view
@@ -0,0 +1,114 @@+module Linear.Plucker+  ( Plucker(..)+  , squaredError+  , isotropic+  , (><)+  , plucker+  , intersects+  ) where++import Control.Applicative+import Data.Distributive+import Data.Foldable as Foldable+import Data.Monoid+import Data.Traversable+import Linear.Epsilon+import Linear.Metric+import Control.Lens+import Linear.V4++-- | Plücker coordinates for lines in a 3-dimensional space.+data Plucker a = Plucker a a a a a a deriving (Eq,Ord,Show,Read)++instance Functor Plucker where+  fmap g (Plucker a b c d e f) = Plucker (g a) (g b) (g c) (g d) (g e) (g f)++instance Applicative Plucker where+  pure a = Plucker a a a a a a+  Plucker a b c d e f <*> Plucker g h i j k l =+    Plucker (a g) (b h) (c i) (d j) (e k) (f l)++instance Monad Plucker where+  return a = Plucker a a a a a a+  (>>=) = bindRep++instance Distributive Plucker where+  distribute = distributeRep++instance Representable Plucker where+  rep f = Plucker (f p01) (f p02) (f p03) (f p23) (f p31) (f p12)++instance Foldable Plucker where+  foldMap g (Plucker a b c d e f) =+    g a `mappend` g b `mappend` g c `mappend` g d `mappend` g e `mappend` g f++instance Traversable Plucker where+  traverse g (Plucker a b c d e f) =+    Plucker <$> g a <*> g b <*> g c <*> g d <*> g e <*> g f++instance Num a => Num (Plucker a) where+  (+) = liftA2 (+)+  (*) = liftA2 (*)+  negate = fmap negate+  abs = fmap abs+  signum = fmap signum+  fromInteger = pure . fromInteger++instance Fractional a => Fractional (Plucker a) where+  recip = fmap recip+  (/) = liftA2 (/)+  fromRational = pure . fromRational++-- | Given a pair of points represented by homogeneous coordinates generate Plücker coordinates+-- for the line through them.+plucker :: Num a => V4 a -> V4 a -> Plucker a+plucker (V4 a b c d)+        (V4 e f g h) =+  Plucker (a*f-b*e)+          (a*g-c*e)+          (a*d-h*e)+          (c*h-d*g)+          (d*f-b*h)+          (b*g-c*f)++-- | These elements form a basis for the Plücker space, or the Grassmanian manifold @Gr(2,V4)@.+p01, p02, p03, p23, p31, p12 :: Functor f => (a -> f a) -> Plucker a -> f (Plucker a)+p01 g (Plucker a b c d e f) = (\a' -> Plucker a' b c d e f) <$> g a+p02 g (Plucker a b c d e f) = (\b' -> Plucker a b' c d e f) <$> g b+p03 g (Plucker a b c d e f) = (\c' -> Plucker a b c' d e f) <$> g c+p23 g (Plucker a b c d e f) = (\d' -> Plucker a b c d' e f) <$> g d+p31 g (Plucker a b c d e f) = (\e' -> Plucker a b c d e' f) <$> g e+p12 g (Plucker a b c d e f) = Plucker a b c d e <$> g f+{-# INLINE p01 #-}+{-# INLINE p02 #-}+{-# INLINE p03 #-}+{-# INLINE p23 #-}+{-# INLINE p31 #-}+{-# INLINE p12 #-}++-- | Valid Plücker coordinates @p@ will have @'squaredError' p '==' 0@+--+-- That said, floating point makes a mockery of this claim, so you may want to use 'nearZero'.+squaredError :: (Eq a, Num a) => Plucker a -> a+squaredError v = v >< v++-- | This isn't th actual metric because this bilinear form gives rise to an isotropic quadratic space+infixl 5 ><+(><) :: Num a => Plucker a -> Plucker a -> a+Plucker a b c d e f >< Plucker g h i j k l = a*g+b*h+c*i-d*j-e*k-f*l++-- | Checks if the line is near-isotropic (isotropic vectors in this quadratic space represent lines in real 3d space)+isotropic :: Epsilon a => Plucker a -> Bool+isotropic a = nearZero (a >< a)++-- | Checks if the two vectors intersect (or nearly intersect)+intersects :: Epsilon a => Plucker a -> Plucker a -> Bool+intersects a b = nearZero (a >< b)++instance Metric Plucker where+  dot (Plucker a b c d e f) (Plucker g h i j k l) = a*g+b*h+c*i+d*j+e*k+f*l++instance Epsilon a => Epsilon (Plucker a) where+  nearZero = nearZero . quadrance++-- TODO: drag some more stuff out of my thesis
+ src/Linear/Quaternion.hs view
@@ -0,0 +1,337 @@+{-# LANGUAGE DeriveDataTypeable, PatternGuards, ScopedTypeVariables #-}+module Linear.Quaternion+  ( Quaternion(..)+  , Complicated(..)+  , Hamiltonian(..)+  , slerp+  , asinq+  , acosq+  , atanq+  , asinhq+  , acoshq+  , atanhq+  , absi+  , pow+  , rotate+  , axisAngle+  ) where+import Control.Applicative+import Control.Lens+import Data.Complex (Complex((:+)))+import Data.Data+import Data.Distributive+import Data.Foldable+import qualified Data.Foldable as F+import Data.Monoid+import Foreign.Ptr (castPtr, plusPtr)+import Foreign.Storable (Storable(..))+import Linear.Epsilon+import Linear.Conjugate+import Linear.Metric+import Linear.V3+import Linear.Vector+import Prelude hiding (any)++data Quaternion a = Quaternion a {-# UNPACK #-}!(V3 a)+                    deriving (Eq,Ord,Read,Show,Data,Typeable)++instance Functor Quaternion where+  fmap f (Quaternion e v) = Quaternion (f e) (fmap f v)+  a <$ _ = Quaternion a (V3 a a a)++instance Applicative Quaternion where+  pure a = Quaternion a (pure a)+  Quaternion f fv <*> Quaternion a v = Quaternion (f a) (fv <*> v)++instance Monad Quaternion where+  return = pure+  (>>=) = bindRep -- the diagonal of a sedenion is super useful!++instance Representable Quaternion where+  rep f = Quaternion (f _e) (V3 (f _i) (f _j) (f _k))++instance Foldable Quaternion where+  foldMap f (Quaternion e v) = f e `mappend` foldMap f v+  foldr f z (Quaternion e v) = f e (F.foldr f z v)++instance Traversable Quaternion where+  traverse f (Quaternion e v) = Quaternion <$> f e <*> traverse f v++instance forall a. Storable a => Storable (Quaternion a) where+  sizeOf _ = 4 * sizeOf (undefined::a)+  alignment _ = alignment (undefined::a)+  poke ptr (Quaternion e v) = poke (castPtr ptr) e >>+                              poke (castPtr (ptr `plusPtr` sz)) v+    where sz = sizeOf (undefined::a)+  peek ptr = Quaternion <$> peek (castPtr ptr)+                        <*> peek (castPtr (ptr `plusPtr` sz))+    where sz = sizeOf (undefined::a)++instance RealFloat a => Num (Quaternion a) where+  {-# SPECIALIZE instance Num (Quaternion Float) #-}+  {-# SPECIALIZE instance Num (Quaternion Double) #-}+  (+) = liftA2 (+)+  (-) = liftA2 (-)+  negate = fmap negate+  Quaternion s1 v1 * Quaternion s2 v2 = Quaternion (s1*s2 - (v1 `dot` v2)) $+                                        (v1 `cross` v2) + s1*^v2 + s2*^v1+  fromInteger x = Quaternion (fromInteger x) 0+  abs z = Quaternion (norm z) 0+  signum q@(Quaternion e (V3 i j k))+    | m == 0.0 = q+    | not (isInfinite m || isNaN m) = q ^/ sqrt m+    | any isNaN q = qNaN+    | not (ii || ij || ik) = Quaternion 1 (V3 0 0 0)+    | not (ie || ij || ik) = Quaternion 0 (V3 1 0 0)+    | not (ie || ii || ik) = Quaternion 0 (V3 0 1 0)+    | not (ie || ii || ij) = Quaternion 0 (V3 0 0 1)+    | otherwise = qNaN+    where+      m = quadrance q+      ie = isInfinite e+      ii = isInfinite i+      ij = isInfinite j+      ik = isInfinite k++  -- abs    = error "Quaternion.abs: use norm"+  -- signum = error "Quaternion.signum: use signorm"++qNaN :: RealFloat a => Quaternion a+qNaN = Quaternion fNaN (V3 fNaN fNaN fNaN) where fNaN = 0/0++-- {-# RULES "abs/norm" abs x = Quaternion (norm x) 0 #-}+-- {-# RULES "signum/signorm" signum = signorm #-}++-- this will attempt to rewrite calls to abs to use norm intead when it is available.++instance RealFloat a => Fractional (Quaternion a) where+  {-# SPECIALIZE instance Fractional (Quaternion Float) #-}+  {-# SPECIALIZE instance Fractional (Quaternion Double) #-}+  Quaternion q0 (V3 q1 q2 q3) / Quaternion r0 (V3 r1 r2 r3) =+    Quaternion (r0*q0+r1*q1+r2*q2+r3*q3)+               (V3 (r0*q1-r1*q0-r2*q3+r3*q2)+                   (r0*q2+r1*q3-r2*q0-r3*q1)+                   (r0*q3-r1*q2+r2*q1-r3*q0))+               ^/ (r0*r0 + r1*r1 + r2*r2 + r3*r3)+  recip q = q ^/ quadrance q+  fromRational x = Quaternion (fromRational x) 0++instance Metric Quaternion where+  Quaternion e v `dot` Quaternion e' v' = e*e' + (v `dot` v')++class Complicated t where+  _e :: Functor f => (a -> f a) -> t a -> f (t a)+  _i :: Functor f => (a -> f a) -> t a -> f (t a)++instance Complicated Complex where+  _e f (a :+ b) = (:+ b) <$> f a+  _i f (a :+ b) = (a :+) <$> f b++instance Complicated Quaternion where+  _e f (Quaternion a v) = (\a' -> Quaternion a' v) <$> f a+  _i f (Quaternion a v) = Quaternion a <$> traverseOf _x f v+  --_i f (Quaternion a (V3 b c d)) = (\b' -> Quaternion a (V3 b' c d)) <$> f b++class Complicated t => Hamiltonian t where+  _j :: Functor f => (a -> f a) -> t a -> f (t a)+  _k :: Functor f => (a -> f a) -> t a -> f (t a)+  _ijk :: Functor f => (V3 a -> f (V3 a)) -> t a -> f (t a)++instance Hamiltonian Quaternion where+  _j f (Quaternion a v) = Quaternion a <$> traverseOf _y f v+  _k f (Quaternion a v) = Quaternion a <$> traverseOf _z f v+  -- _j f (Quaternion a (V3 b c d)) = (\c' -> Quaternion a (V3 b c' d)) <$> f c+  -- _k f (Quaternion a (V3 b c d)) = Quaternion a . V3 b c <$> f d++  _ijk f (Quaternion a v) = Quaternion a <$> f v++instance Distributive Quaternion where+  distribute = distributeRep++instance (Conjugate a, RealFloat a) => Conjugate (Quaternion a) where+  conjugate (Quaternion e v) = Quaternion (conjugate e) (negate v)++reimagine :: RealFloat a => a -> a -> Quaternion a -> Quaternion a+reimagine r s (Quaternion _ v)+  | isNaN s || isInfinite s = let aux 0 = 0+                                  aux x = s * x+                              in Quaternion r (aux <$> v)+  | otherwise = Quaternion r (v^*s)++-- | quadrance of the imaginary component+qi :: Num a => Quaternion a -> a+qi (Quaternion _ v) = quadrance v++-- | norm of the imaginary component+absi :: Floating a => Quaternion a -> a+absi = sqrt . qi++-- | raise a 'Quaternion' to a scalar power+pow :: RealFloat a => Quaternion a -> a -> Quaternion a+pow q t = exp (t *^ log q)++-- ehh..+instance RealFloat a => Floating (Quaternion a) where+  {-# SPECIALIZE instance Floating (Quaternion Float) #-}+  {-# SPECIALIZE instance Floating (Quaternion Double) #-}+  pi = Quaternion pi 0+  exp q@(Quaternion e v)+    | qiq == 0 = Quaternion (exp e) v+    | ai <- sqrt qiq, ee <- exp e = reimagine (ee * cos ai) (ee * (sin ai / ai)) q+    where qiq = qi q+  log q@(Quaternion e v@(V3 _i j k))+    | qiq == 0 = if e >= 0+                 then Quaternion (log e) v+                 else Quaternion (log (negate e)) (V3 pi j k) -- mmm, pi+    | ai <- sqrt qiq, m <- sqrt (e*e + qiq) = reimagine (log m) (atan2 m e / ai) q+    where qiq = qi q+  x ** y = exp (y * log x)+  sqrt q@(Quaternion e v)+    | m   == 0 = q+    | qiq == 0 = if e > 0+                 then Quaternion (sqrt e) 0+                 else Quaternion 0 (V3 (sqrt (negate e)) 0 0)+    | im <- sqrt (0.5*(m-e)) / sqrt qiq = Quaternion (0.5*(m+e)) (v^*im)+    where qiq = qi q+          m = sqrt (e*e + qiq)+  cos q@(Quaternion e v)+    | qiq == 0 = Quaternion (cos e) v+    | ai <- sqrt qiq = reimagine (cos e * cosh ai) (- sin e * (sinh ai / ai)) q+    where qiq = qi q+  sin q@(Quaternion e v)+    | qiq == 0 = Quaternion (sin e) v+    | ai <- sqrt qiq = reimagine (sin e * cosh ai) (cos e * (sinh ai / ai)) q+    where qiq = qi q+  tan q@(Quaternion e v)+    | qiq == 0 = Quaternion (tan e) v+    | ai <- sqrt qiq, ce <- cos e, sai <- sinh ai, d <- ce*ce + sai*sai =+      reimagine (ce * sin e / d) (cosh ai * (sai / ai) / d) q+    where qiq = qi q+  sinh q@(Quaternion e v)+    | qiq == 0 = Quaternion (sinh e) v+    | ai <- sqrt qiq = reimagine (sinh e * cos ai) (cosh e * (sin ai / ai)) q+    where qiq = qi q+  cosh q@(Quaternion e v)+    | qiq == 0 = Quaternion (cosh e) v+    | ai <- sqrt qiq = reimagine (cosh e * cos ai) ((sinh e * sin ai) / ai) q+    where qiq = qi q+  tanh q@(Quaternion e v)+    | qiq == 0 = Quaternion (tanh e) v+    | ai <- sqrt qiq, se <- sinh e, cai <- cos ai, d <- se*se + cai*cai =+      reimagine ((cosh e * se) / d) ((cai * (sin ai / ai)) / d) q+    where qiq = qi q++  asin q = cut asin q+  acos q = cut acos q+  atan q = cut atan q++  asinh q = cut asinh q+  acosh q = cut acosh q+  atanh q = cut atanh q++-- | Helper for calculating with specific branch cuts+cut :: RealFloat a => (Complex a -> Complex a) -> Quaternion a -> Quaternion a+cut f q@(Quaternion e v)+  | qiq == 0 = Quaternion a (_x.~b$v)+  | otherwise = reimagine a (b / ai) q+  where qiq = qi q+        ai = sqrt qiq+        a :+ b = f (e :+ ai)++-- | Helper for calculating with specific branch cuts+cutWith :: RealFloat a => Complex a -> Quaternion a -> Quaternion a+cutWith (r :+ im) q@(Quaternion e v)+  | e /= 0 || qiq == 0 || isNaN qiq || isInfinite qiq = error "bad cut"+  | s <- im / sqrt qiq = Quaternion r (v^*s)+  where qiq = qi q++-- | 'asin' with a specified branch cut.+asinq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a+asinq q@(Quaternion e _) u+  | qiq /= 0.0 || e >= -1 && e <= 1 = asin q+  | otherwise = cutWith (asin (e :+ sqrt qiq)) u+  where qiq = qi q++-- | 'acos' with a specified branch cut.+acosq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a+acosq q@(Quaternion e _) u+  | qiq /= 0.0 || e >= -1 && e <= 1 = acos q+  | otherwise = cutWith (acos (e :+ sqrt qiq)) u+  where qiq = qi q++-- | 'atan' with a specified branch cut.+atanq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a+atanq q@(Quaternion e _) u+  | e /= 0.0 || qiq >= -1 && qiq <= 1 = atan q+  | otherwise = cutWith (atan (e :+ sqrt qiq)) u+  where qiq = qi q++-- | 'asinh' with a specified branch cut.+asinhq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a+asinhq q@(Quaternion e _) u+  | e /= 0.0 || qiq >= -1 && qiq <= 1 = asinh q+  | otherwise = cutWith (asinh (e :+ sqrt qiq)) u+  where qiq = qi q++-- | 'acosh' with a specified branch cut.+acoshq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a+acoshq q@(Quaternion e _) u+  | qiq /= 0.0 || e >= 1 = asinh q+  | otherwise = cutWith (acosh (e :+ sqrt qiq)) u+  where qiq = qi q++-- | 'atanh' with a specified branch cut.+atanhq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a+atanhq q@(Quaternion e _) u+  | qiq /= 0.0 || e > -1 && e < 1 = atanh q+  | otherwise = cutWith (atanh (e :+ sqrt qiq)) u+  where qiq = qi q++-- | Spherical linear interpolation between two quaternions.+slerp :: RealFloat a => Quaternion a -> Quaternion a -> a -> Quaternion a+slerp q p t+  | 1.0 - cosphi < 1e-8 = q+  | phi <- acos cosphi, r <- recip (sin phi)+  = (sin ((1-t)*phi)*r *^ q ^+^ f (sin (t*phi)*r) *^ p) ^/ sin phi+  where+   dqp = dot q p+   (cosphi, f) = if dqp < 0 then (-dqp, negate) else (dqp, id)+{-# SPECIALIZE slerp :: Quaternion Float -> Quaternion Float -> Float -> Quaternion Float #-}+{-# SPECIALIZE slerp :: Quaternion Double -> Quaternion Double -> Double -> Quaternion Double #-}++--slerp :: RealFloat a => Quaternion a -> Quaternion a -> a -> Quaternion a+--slerp q0 q1 = let q10 = q1 / q0 in \t -> pow q10 t * q0++-- | Apply a rotation to a vector.+rotate :: (Conjugate a, RealFloat a) => Quaternion a -> V3 a -> V3 a+rotate q v = (q * Quaternion 0 v * conjugate q)^._ijk++{-+rotate :: Num a => Quaternion a -> V3 a -> V3 a+rotate (Quaternion a' b c d) (V3 x y z) = V3+  (2*((t8+t10)*x+(t6- t4)*y+(t3+t7)*z)+x)+  (2*((t4+ t6)*y+(t5+t10)*y+(t9-t2)*z)+y)+  (2*((t7- t3)*z+(t2+ t9)*z+(t5+t8)*z)+z)+  where+    a = -a'+    t2 = a*b+    t3 = a*c+    t4 = a*d+    t5 = -b*b+    t6 = b*c+    t7 = b*d+    t8 = -c*c+    t9 = c*d+    t10 = -d*d+-}+{-# SPECIALIZE rotate :: Quaternion Float -> V3 Float -> V3 Float #-}+{-# SPECIALIZE rotate :: Quaternion Double -> V3 Double -> V3 Double #-}++instance (RealFloat a, Epsilon a) => Epsilon (Quaternion a) where+  nearZero = nearZero . quadrance++-- | @'axisAngle' axis theta@ builds a 'Quaternion' representing a+-- rotation of @theta@ radians about @axis@.+axisAngle :: (Epsilon a, Floating a) => V3 a -> a -> Quaternion a+axisAngle axis theta = normalize $ Quaternion (cos half) $ (sin half) *^ axis+  where half = theta / 2
+ src/Linear/V2.hs view
@@ -0,0 +1,92 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE ScopedTypeVariables #-}+-- {-# OPTIONS_GHC -fno-warn-name-shadowing #-}+module Linear.V2+  ( V2(..)+  , R2(..)+  , perp+  ) where++import Control.Applicative+import Control.Lens+import Data.Data+import Data.Distributive+import Data.Foldable+import Data.Monoid+import Foreign.Ptr (castPtr)+import Foreign.Storable (Storable(..))+import Linear.Metric+import Linear.Epsilon++-- | A 2-dimensional vector+data V2 a = V2 a a deriving (Eq,Ord,Show,Read,Data,Typeable)++instance Functor V2 where+  fmap f (V2 a b) = V2 (f a) (f b)++instance Foldable V2 where+  foldMap f (V2 a b) = f a `mappend` f b++instance Traversable V2 where+  traverse f (V2 a b) = V2 <$> f a <*> f b++instance Applicative V2 where+  pure a = V2 a a+  V2 a b <*> V2 d e = V2 (a d) (b e)++instance Monad V2 where+  return a = V2 a a+  (>>=) = bindRep++instance Num a => Num (V2 a) where+  (+) = liftA2 (+)+  (*) = liftA2 (*)+  negate = fmap negate+  abs = fmap abs+  signum = fmap signum+  fromInteger = pure . fromInteger++instance Fractional a => Fractional (V2 a) where+  recip = fmap recip+  (/) = liftA2 (/)+  fromRational = pure . fromRational++instance Metric V2 where+  dot (V2 a b) (V2 c d) = a * c + b * d++-- | A space that distinguishes 2 orthogonal basis vectors '_x' and '_y', but may have more.+class R2 t where+  _x :: Functor f => (a -> f a) -> t a -> f (t a)+  _x = _xy._x++  _y :: Functor f => (a -> f a) -> t a -> f (t a)+  _y = _xy._y++  _xy :: Functor f => (V2 a -> f (V2 a)) -> t a -> f (t a)++instance R2 V2 where+  _x f (V2 a b) = (`V2` b) <$> f a+  _y f (V2 a b) = (V2 a) <$> f b+  _xy = id++instance Representable V2 where+  rep f = V2 (f _x) (f _y)++instance Distributive V2 where+  distribute f = V2 (fmap (^._x) f) (fmap (^._y) f)++-- | the counter-clockwise perpendicular vector+perp :: Num a => V2 a -> V2 a+perp (V2 a b) = V2 (negate b) a++instance Epsilon a => Epsilon (V2 a) where+  nearZero = nearZero . quadrance++instance forall a. Storable a => Storable (V2 a) where+  sizeOf _ = 2 * sizeOf (undefined::a)+  alignment _ = alignment (undefined::a)+  poke ptr (V2 x y) = poke ptr' x >> pokeElemOff ptr' 1 y+    where ptr' = castPtr ptr+  peek ptr = V2 <$> peek ptr' <*> peekElemOff ptr' 1+    where ptr' = castPtr ptr
+ src/Linear/V3.hs view
@@ -0,0 +1,98 @@+{-# LANGUAGE DeriveDataTypeable, ScopedTypeVariables #-}+module Linear.V3+  ( V3(..)+  , cross, triple+  , R2(..)+  , R3(..)+  ) where++import Control.Applicative+import Control.Lens+import Data.Data+import Data.Distributive+import Data.Foldable+import Data.Monoid+import Foreign.Ptr (castPtr)+import Foreign.Storable (Storable(..))+import Linear.Epsilon+import Linear.Metric+import Linear.V2++-- | A 3-dimensional vector+data V3 a = V3 a a a deriving (Eq,Ord,Show,Read,Data,Typeable)++instance Functor V3 where+  fmap f (V3 a b c) = V3 (f a) (f b) (f c)++instance Foldable V3 where+  foldMap f (V3 a b c) = f a `mappend` f b `mappend` f c++instance Traversable V3 where+  traverse f (V3 a b c) = V3 <$> f a <*> f b <*> f c++instance Applicative V3 where+  pure a = V3 a a a+  V3 a b c <*> V3 d e f = V3 (a d) (b e) (c f)++instance Monad V3 where+  return a = V3 a a a+  (>>=) = bindRep++instance Num a => Num (V3 a) where+  (+) = liftA2 (+)+  (*) = liftA2 (*)+  negate = fmap negate+  abs = fmap abs+  signum = fmap signum+  fromInteger = pure . fromInteger++instance Fractional a => Fractional (V3 a) where+  recip = fmap recip+  (/) = liftA2 (/)+  fromRational = pure . fromRational++instance Metric V3 where+  dot (V3 a b c) (V3 d e f) = a * d + b * e + c * f+  {-# INLINABLE dot #-}++instance Distributive V3 where+  distribute f = V3 (fmap (^._x) f) (fmap (^._y) f) (fmap (^._z) f)++-- | A space that distinguishes 3 orthogonal basis vectors: '_x', '_y', and '_z'. (It may have more)+class R2 t => R3 t where+  _z :: Functor f => (a -> f a) -> t a -> f (t a)+  _xyz :: Functor f => (V3 a -> f (V3 a)) -> t a -> f (t a)++instance R2 V3 where+  _x f (V3 a b c) = (\a' -> V3 a' b c) <$> f a+  _y f (V3 a b c) = (\b' -> V3 a b' c) <$> f b+  _xy f (V3 a b c) = (\(V2 a' b') -> V3 a' b' c) <$> f (V2 a b)++instance R3 V3 where+  _z f (V3 a b c) = V3 a b <$> f c+  _xyz = id++instance Representable V3 where+  rep f = V3 (f _x) (f _y) (f _z)++instance forall a. Storable a => Storable (V3 a) where+  sizeOf _ = 3 * sizeOf (undefined::a)+  alignment _ = alignment (undefined::a)+  poke ptr (V3 x y z) = do poke ptr' x+                           pokeElemOff ptr' 1 y+                           pokeElemOff ptr' 2 z+    where ptr' = castPtr ptr+  peek ptr = V3 <$> peek ptr' <*> peekElemOff ptr' 1 <*> peekElemOff ptr' 2+    where ptr' = castPtr ptr++-- | cross product+cross :: Num a => V3 a -> V3 a -> V3 a+cross (V3 a b c) (V3 d e f) = V3 (b*f-c*e) (c*d-a*f) (a*e-b*d)+{-# INLINABLE cross #-}++-- | scalar triple product+triple :: Num a => V3 a -> V3 a -> V3 a -> a+triple a b c = dot a (cross b c)++instance Epsilon a => Epsilon (V3 a) where+  nearZero = nearZero . quadrance
+ src/Linear/V4.hs view
@@ -0,0 +1,104 @@+{-# LANGUAGE DeriveDataTypeable, ScopedTypeVariables #-}+module Linear.V4+  ( V4(..)+  , vector, point+  , R2(..)+  , R3(..)+  , R4(..)+  ) where++import Control.Applicative+import Control.Lens+import Data.Data+import Data.Distributive+import Data.Foldable+import Data.Monoid+import Foreign.Ptr (castPtr)+import Foreign.Storable (Storable(..))+import Linear.Epsilon+import Linear.Metric+import Linear.V2+import Linear.V3++-- | A 4-dimensional vector.+data V4 a = V4 a a a a deriving (Eq,Ord,Show,Read,Data,Typeable)++instance Functor V4 where+  fmap f (V4 a b c d) = V4 (f a) (f b) (f c) (f d)++instance Foldable V4 where+  foldMap f (V4 a b c d) = f a `mappend` f b `mappend` f c `mappend` f d++instance Traversable V4 where+  traverse f (V4 a b c d) = V4 <$> f a <*> f b <*> f c <*> f d++instance Applicative V4 where+  pure a = V4 a a a a+  V4 a b c d <*> V4 e f g h = V4 (a e) (b f) (c g) (d h)++instance Monad V4 where+  return a = V4 a a a a+  (>>=) = bindRep++instance Num a => Num (V4 a) where+  (+) = liftA2 (+)+  (*) = liftA2 (*)+  negate = fmap negate+  abs = fmap abs+  signum = fmap signum+  fromInteger = pure . fromInteger++instance Fractional a => Fractional (V4 a) where+  recip = fmap recip+  (/) = liftA2 (/)+  fromRational = pure . fromRational++instance Metric V4 where+  dot (V4 a b c d) (V4 e f g h) = a * e + b * f + c * g + d * h++instance Distributive V4 where+  distribute f = V4 (fmap (^._x) f) (fmap (^._y) f) (fmap (^._z) f) (fmap (^._w) f)++-- | A space that distinguishes orthogonal basis vectors '_x', '_y', '_z', '_w'. (It may have more.)+class R3 t => R4 t where+  _w :: Functor f => (a -> f a) -> t a -> f (t a)+  _xyzw :: Functor f => (V4 a -> f (V4 a)) -> t a -> f (t a)++instance R2 V4 where+  _x f (V4 a b c d) = (\a' -> V4 a' b c d) <$> f a+  _y f (V4 a b c d) = (\b' -> V4 a b' c d) <$> f b+  _xy f (V4 a b c d) = (\(V2 a' b') -> V4 a' b' c d) <$> f (V2 a b)++instance R3 V4 where+  _z f (V4 a b c d) = (\c' -> V4 a b c' d) <$> f c+  _xyz f (V4 a b c d) = (\(V3 a' b' c') -> V4 a' b' c' d) <$> f (V3 a b c)++instance R4 V4 where+  _w f (V4 a b c d) = V4 a b c <$> f d+  _xyzw = id++instance Representable V4 where+  rep f = V4 (f _x) (f _y) (f _z) (f _w)++instance forall a. Storable a => Storable (V4 a) where+  sizeOf _ = 4 * sizeOf (undefined::a)+  alignment _ = alignment (undefined::a)+  poke ptr (V4 x y z w) = do poke ptr' x+                             pokeElemOff ptr' 1 y+                             pokeElemOff ptr' 2 z+                             pokeElemOff ptr' 3 w+    where ptr' = castPtr ptr+  peek ptr = V4 <$> peek ptr' <*> peekElemOff ptr' 1+                <*> peekElemOff ptr' 2 <*> peekElemOff ptr' 3+    where ptr' = castPtr ptr++-- | Convert a 3-dimensional affine vector into a 4-dimensional homogeneous vector.+vector :: Num a => V3 a -> V4 a+vector (V3 a b c) = V4 a b c 0++-- | Convert a 3-dimensional affine point into a 4-dimensional homogeneous vector.+point :: Num a => V3 a -> V4 a+point (V3 a b c) = V4 a b c 1++instance Epsilon a => Epsilon (V4 a) where+  nearZero = nearZero . quadrance
+ src/Linear/Vector.hs view
@@ -0,0 +1,60 @@+-----------------------------------------------------------------------------+-- |+-- Module      :  Linear.Epsilon+-- Copyright   :  (C) 2012 Edward Kmett+-- License     :  BSD-style (see the file LICENSE)+-- Maintainer  :  Edward Kmett <ekmett@gmail.com>+-- Stability   :  provisional+-- Portability :  portable+--+-- Operations on free vector spaces.+-----------------------------------------------------------------------------+module Linear.Vector+  ( (^+^)+  , gnegate+  , (^-^)+  , (^*)+  , (*^)+  , (^/)+  , lerp+  ) where++import Control.Applicative++infixl 6 ^+^, ^-^+infixl 7 ^*, *^, ^/++-- | Compute the sum of two vectors+(^+^) :: (Applicative f, Num a) => f a -> f a -> f a+(^+^) = liftA2 (+)+{-# INLINE (^+^) #-}++-- | Compute the negation of a vector+gnegate :: (Functor f, Num a) => f a -> f a+gnegate = fmap negate+{-# INLINE gnegate #-}++-- | Compute the difference between two vectors+(^-^) :: (Applicative f, Num a) => f a -> f a -> f a+(^-^) = liftA2 (-)+{-# INLINE (^-^) #-}++-- | Compute the left scalar product+(*^) :: (Functor f, Num a) => a -> f a -> f a+(*^) a = fmap (a*)+{-# INLINE (*^) #-}++-- | Compute the right scalar product+(^*) :: (Functor f, Num a) => f a -> a -> f a+f ^* a = fmap (*a) f+{-# INLINE (^*) #-}++-- | Compute division by a scalar on the right.+(^/) :: (Functor f, Fractional a) => f a -> a -> f a+f ^/ a = fmap (/a) f+{-# INLINE (^/) #-}++-- | Linearly interpolate between two vectors.+lerp :: (Applicative f, Num a) => a -> f a -> f a -> f a+lerp alpha u v = alpha *^ u ^+^ (1 - alpha) *^ v+{-# INLINE lerp #-}
+ tests/doctests.hs view
@@ -0,0 +1,28 @@+module Main where++import Test.DocTest+import System.Directory+import System.FilePath+import Control.Applicative+import Control.Monad+import Data.List++main :: IO ()+main = getSources >>= \sources -> doctest $+    "-isrc"+  : "-idist/build/autogen"+  : "-optP-include"+  : "-optPdist/build/autogen/cabal_macros.h"+  : sources++getSources :: IO [FilePath]+getSources = filter (isSuffixOf ".hs") <$> go "src"+  where+    go dir = do+      (dirs, files) <- getFilesAndDirectories dir+      (files ++) . concat <$> mapM go dirs++getFilesAndDirectories :: FilePath -> IO ([FilePath], [FilePath])+getFilesAndDirectories dir = do+  c <- map (dir </>) . filter (`notElem` ["..", "."]) <$> getDirectoryContents dir+  (,) <$> filterM doesDirectoryExist c <*> filterM doesFileExist c