gasp 1.2.0.0 → 1.3.0.0

raw patch · 3 files changed

+326/−152 lines, 3 filesdep +QuickCheckPVP ok

version bump matches the API change (PVP)

Dependencies added: QuickCheck

API changes (from Hackage documentation)

Files

• Algebra/Classes.hs +224/−91
• Algebra/Linear.hs +100/−60
• gasp.cabal +2/−1

@@ -1,8 +1,15 @@+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE AllowAmbiguousTypes #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-}+{-# LANGUAGE GADTs #-} {-# LANGUAGE MultiParamTypeClasses, ConstraintKinds, FlexibleContexts, FlexibleInstances, DeriveGeneric #-} module Algebra.Classes where -import Prelude as Algebra.Classes (Int,Integer,Float,Double, Foldable (..), (==), Monoid(..), Ord(..)-                                  ,Real(..), Enum(..), snd, Rational, Functor(..), Eq(..), Bool(..), Semigroup(..))+import Prelude (Int,Integer,Float,Double, (==), Monoid(..), Ord(..), Foldable,+                foldMap, Char,+                 Real(..), Enum(..), snd, Rational, Functor(..), Eq(..), Bool(..), Semigroup(..), Show(..), uncurry)+ import qualified Prelude import qualified Data.Ratio import qualified Data.Map.Strict as M@@ -10,8 +17,13 @@ import Foreign.C import Data.Word import Data.Binary+import Data.Complex import GHC.Generics+import Test.QuickCheck+import Control.Applicative +-- import Data.Functor.Utils ((#.))+ infixl 6 - infixl 6 + @@ -20,6 +32,8 @@ infixl 7 / infixl 7 `mod` infixl 7 `div`+infixr 8 ^+infixr 8 ^+  type Natural = Integer @@ -48,25 +62,62 @@ instance Additive a => Multiplicative (Exponential a) where   Exponential a * Exponential b = Exponential (a + b)   one = Exponential zero-  Exponential a ^ n = Exponential (times n a)+  Exponential a ^+ n = Exponential (times n a)  instance Group a => Division (Exponential a) where   recip (Exponential a) = Exponential (negate a)   Exponential a / Exponential b = Exponential (a - b) +timesDefault :: (Additive a2, Prelude.Integral a1) => a1 -> a2 -> a2+timesDefault n0 = if n0 < 0 then Prelude.error "Algebra.Classes.times: negative number of times" else go n0+    where go 0 _ = zero+          go n x = if r == 0 then y + y else x + y + y+            where (m,r) = n `Prelude.divMod` 2+                  y = go m x+ -- | Additive monoid class Additive a where   (+) :: a -> a -> a   zero :: a   times :: Natural -> a -> a-  times 0 _ = zero-  times n x = if r == 0 then y + y else x + y + y-    where (m,r) = n `Prelude.divMod` 2-          y = times m x+  times = timesDefault -add :: (Foldable t, Additive a) => t a -> a-add xs = fromSum (foldMap Sum xs)+class (Arbitrary a, Show a) => TestEqual a where+  (=.=) :: a -> a -> Property +infix 0 =.=++instance Multiplicative Property where+  one = property True+  (*) = (.&&.)++nameLaw :: Testable prop => Prelude.String -> prop -> Property+nameLaw x p = label x (counterexample x p)++law_zero_plus :: forall a. (Additive a, TestEqual a) => a -> Property+law_zero_plus n = nameLaw "zero/plus" (zero + n =.= n)++law_plus_zero :: (Additive a, TestEqual a) => a -> Property+law_plus_zero n = nameLaw "plus/zero" (n + zero =.= n)++law_plus_assoc :: (Additive a, TestEqual a) => a -> a -> a -> Property+law_plus_assoc m n o = nameLaw "plus/assoc" (n + (m + o) =.= (n + m) + o)++law_times :: (TestEqual a, Additive a) => Positive Integer -> a -> Property+law_times (Positive m) n = nameLaw "times" (times m n =.= timesDefault m n)++laws_additive :: forall a. (Additive a, TestEqual a) => Property+laws_additive = product [property (law_zero_plus @a)+                        ,property (law_plus_zero @a)+                        ,property (law_plus_assoc @a)+                        ,property (law_times @a)]++instance TestEqual Int where (=.=) = (===)+++sum :: (Foldable t, Additive a) => t a -> a+sum xs = fromSum (foldMap Sum xs)+ instance Additive Integer where   (+) = (Prelude.+)   zero = 0@@ -115,6 +166,10 @@ class Additive r => DecidableZero r where   isZero :: r -> Bool +law_decidable_zero :: forall a. (DecidableZero a, TestEqual a) => Property+law_decidable_zero = property (isZero (zero @a))++ instance DecidableZero Integer where   isZero = (== 0) instance DecidableZero CInt where@@ -136,6 +191,13 @@  class Additive a => AbelianAdditive a   -- just a law.++law_plus_comm :: (TestEqual a, Additive a) => a -> a -> Property+law_plus_comm m n = nameLaw "plus/comm" (m + n =.= n + m)++laws_abelian_additive :: forall a. (Group a, TestEqual a) => Property+laws_abelian_additive = laws_additive @a .&&. product [property (law_plus_comm @a)]+ instance AbelianAdditive Integer instance AbelianAdditive CInt instance AbelianAdditive Int@@ -143,6 +205,9 @@ instance AbelianAdditive Float instance (Ord k,AbelianAdditive v) => AbelianAdditive (Map k v) +multDefault :: Group a => Natural -> a -> a+multDefault n x = if n < 0 then negate (times (negate n) x) else times n x+ class Additive a => Group a where   {-# MINIMAL (negate | (-)) #-}   (-) :: a -> a -> a@@ -150,8 +215,23 @@   negate :: a -> a   negate b = zero - b   mult :: Integer -> a -> a-  mult n x = if n < 0 then negate (times (negate n) x) else times n x+  mult = multDefault +law_negate_minus :: (TestEqual a, Group a) => a -> a -> Property+law_negate_minus m n = nameLaw "minus/negate" (m + negate n =.= m - n)++law_mult :: (TestEqual a, Group a) => Integer -> a -> Property+law_mult m n = nameLaw "mult" (mult m n =.= multDefault m n)+++laws_group :: forall a. (Group a, TestEqual a) => Property+laws_group = laws_additive @a .&&. product [property (law_negate_minus @a)+                                           ,property (law_mult @a)]++laws_abelian_group :: forall a. (Group a, TestEqual a) => Property+laws_abelian_group = laws_group @a .&&. product [property (law_plus_comm @a)]++ instance Group Integer where   (-) = (Prelude.-)   negate = Prelude.negate@@ -194,6 +274,36 @@ class (AbelianAdditive a, PreRing scalar) => Module scalar a where   (*^) :: scalar -> a -> a +law_module_zero :: forall s a. (Module s a, TestEqual a) => s -> Property+law_module_zero s = nameLaw "module/zero" (s *^ zero =.= zero @a)++law_module_one :: forall s a. (Module s a, TestEqual a) => a -> Property+law_module_one x = nameLaw "module/one" ((one @s) *^ x =.= x)++law_module_sum :: forall s a. (Module s a, TestEqual a) => s -> a -> a -> Property+law_module_sum s x y = nameLaw "module/distr/left" (s *^ (x + y) =.= s*^x + s *^ y)++law_module_sum_left :: forall s a. (Module s a, TestEqual a) => s -> s -> a -> Property+law_module_sum_left s t x = nameLaw "module/distr/right" ((s + t) *^ x =.= s*^x + t *^ x)++law_module_mul :: forall s a. (Module s a, TestEqual a) => s -> s -> a -> Property+law_module_mul s t x = nameLaw "module/mul/assoc" ((s * t) *^ x =.= s *^ t *^ x)++laws_module :: forall s a. (Module s a, TestEqual a, Arbitrary s, Show s) => Property+laws_module = laws_additive @a .&&. product [property (law_module_zero @s @a)+                                            ,property (law_module_one @s @a)+                                            ,property (law_module_sum @s @a)+                                            ,property (law_module_sum_left @s @a)+                                            ,property (law_module_mul @s @a)+                                            ]++-- Comparision of maps with absence of a key equivalent to zero value.+instance (Ord x, Show x, Arbitrary x,TestEqual a,Additive a) => TestEqual (Map x a) where+  x =.= y = product (uncurry (=.=) <$> M.unionWith collapse ((,zero) <$> x) ((zero,) <$> y))+    where collapse :: (a,b) -> (c,d) -> (a,d)+          collapse (a,_) (_,b) = (a,b)++ instance Module Integer Integer where   (*^) = (*) @@ -209,64 +319,67 @@ instance Module Float Float where   (*^) = (*) -instance (Ord k, Module v v) => Module v (Map k v) where-  s *^ m = fmap (s *) m+instance (Ord k, Module a b) => Module a (Map k b) where+  s *^ m = fmap (s *^) m  -- | Multiplicative monoid class Multiplicative a where   (*) :: a -> a -> a   one :: a-  (^) :: a -> Natural -> a+  (^+) :: a -> Natural -> a -  (^) _ 0 = one-  (^) x n = if r == 0 then y * y else x * y * y-    where (m,r) = n `Prelude.divMod` 2-          y = (^) y m+  x0 ^+ n0 = if n0 < 0 then Prelude.error "Algebra.Classes.^: negative exponent" else go x0 n0+    where go _ 0 = one+          go x n = if r == 0 then y * y else x * y * y+            where (m,r) = n `Prelude.divMod` 2+                  y = go x m -multiply :: (Multiplicative a, Foldable f) => f a -> a-multiply xs = fromProduct (foldMap Product xs) +product :: (Multiplicative a, Foldable f) => f a -> a+product xs = fromProduct (foldMap Product xs)+ instance Multiplicative Integer where   (*) = (Prelude.*)   one = 1-  (^) = (Prelude.^)+  (^+) = (Prelude.^)  instance Multiplicative CInt where   (*) = (Prelude.*)   one = 1-  (^) = (Prelude.^)+  (^+) = (Prelude.^)  instance Multiplicative Word32 where   (*) = (Prelude.*)   one = 1-  (^) = (Prelude.^)+  (^+) = (Prelude.^)  instance Multiplicative Word16 where   (*) = (Prelude.*)   one = 1-  (^) = (Prelude.^)+  (^+) = (Prelude.^)  instance Multiplicative Word8 where   (*) = (Prelude.*)   one = 1-  (^) = (Prelude.^)+  (^+) = (Prelude.^)  instance Multiplicative Int where   (*) = (Prelude.*)   one = 1-  (^) = (Prelude.^)+  (^+) = (Prelude.^)  instance Multiplicative Double where   (*) = (Prelude.*)   one = 1-  (^) = (Prelude.^)+  (^+) = (Prelude.^)  instance Multiplicative Float where   (*) = (Prelude.*)   one = 1-  (^) = (Prelude.^)+  (^+) = (Prelude.^)  + type SemiRing a = (Multiplicative a, AbelianAdditive a) type PreRing a = (SemiRing a, Group a) @@ -300,11 +413,19 @@   (/) :: a -> a -> a   x / y           =  x * recip y +  (^) :: a -> Integer -> a+  b ^ n | n < 0 = recip b ^+ negate n+        | True  = b ^+ n+ instance Division Double where   (/) = (Prelude./)+  recip = Prelude.recip+  (^) = (Prelude.^^)  instance Division Float where   (/) = (Prelude./)+  recip = Prelude.recip+  (^) = (Prelude.^^)  class (Ring a, Division a) => Field a where   fromRational :: Rational -> a@@ -317,7 +438,6 @@ instance Field Float where   fromRational = Prelude.fromRational -type VectorSpace scalar a = (Field scalar, Module scalar a)  class Ring a => EuclideanDomain a where     {-# MINIMAL (stdUnit | normalize) , (divMod | (div , mod)) #-}@@ -370,19 +490,7 @@     rem       =  Prelude.rem     toInteger = Prelude.toInteger -{--Note: the following is not quite what we intuitively want, because -class Field a => AlgebraicallyClosed  a where-  sqrt :: a -> (a,a)--AlgebraicallyClosed numbers have two square roots.---}--data Ratio a = !a :% !a  deriving (Eq)-type MyRational = Ratio Integer- gcd             :: (Integral a) => a -> a -> a {-# NOINLINE [1] gcd #-} gcd x y         =  gcd' (stdAssociate x) (stdAssociate y)@@ -391,50 +499,12 @@    gcd' a 0  =  a    gcd' a b  =  gcd' b (a `rem` b) -{---- | 'reduce' is a subsidiary function used only in this module.--- It normalises a ratio by dividing both numerator and denominator by--- their greatest common divisor.-reduce :: (Eq a, Integral a) => a -> a -> Ratio a-{-# SPECIALISE reduce :: Integer -> Integer -> MyRational #-}-reduce _ 0              =  error "reduce: division by zero"-reduce x y              =  (x `quot` d) :% (y `quot` d)-                           where d = gcd x y+{- -} -(%) :: Integral a => a -> a -> Ratio a-x % y                   =  reduce (x * stdUnit y) (stdAssociate y)  -instance Integral a => AbelianAdditive (Ratio a) where--instance Integral a => Additive (Ratio a) where-    (x:%y) + (x':%y')   =  reduce (x*y' + x'*y) (y*y')-    zero = 0-    times n (x :% y) = reduce (times n x) y--instance Integral a => Multiplicative (Ratio a) where-    (x:%y) * (x':%y')   =  reduce (x * x') (y * y')-    one = 1 :% 1--instance Integral a => Group (Ratio a) where-    (x:%y) - (x':%y')   =  reduce (x*y' - x'*y) (y*y')-    negate (x:%y)       =  (-x) :% y--instance Integral a => EuclideanDomain (Ratio a) where-    stdAssociate (x:%y)          =  stdAssociate x :% y-    stdUnit (x:%_)       =  stdUnit x :% 1--instance Integral a => Ring (Ratio a) where-    fromInteger x       =  fromInteger x :% 1--instance Integral a => Division (Ratio a) where-  recip (x:%y) = y:%x--instance Integral a => Field (Ratio a) where----}-+--------------------------+-- Ratio instances instance Prelude.Integral a => Additive (Data.Ratio.Ratio a) where   zero = Prelude.fromInteger 0   (+) = (Prelude.+)@@ -448,31 +518,94 @@ instance Prelude.Integral a => Multiplicative (Data.Ratio.Ratio a) where   one = Prelude.fromInteger 1   (*) = (Prelude.*)+  (^+) = (Prelude.^)  instance Prelude.Integral a => Division (Data.Ratio.Ratio a) where   recip = Prelude.recip   (/) = (Prelude./)-+  (^) = (Prelude.^^) instance Prelude.Integral a => Module (Data.Ratio.Ratio a) (Data.Ratio.Ratio a) where   (*^) = (*)- instance Prelude.Integral a => Ring (Data.Ratio.Ratio a) where   fromInteger = Prelude.fromInteger- instance Prelude.Integral a => Field (Data.Ratio.Ratio a) where   fromRational = Prelude.fromRational ++----------------------+-- Complex instances+instance Module Rational Double where+    r *^ d = fromRational r * d+instance Additive a => Additive (Complex a) where+    (x:+y) + (x':+y')   =  (x+x') :+ (y+y')+    zero = zero :+ zero+instance Ring a => Multiplicative (Complex a) where+    (x:+y) * (x':+y')   =  (x*x'-y*y') :+ (x*y'+y*x')+    one = one :+ zero+instance Group a => Group  (Complex a) where+    (x:+y) - (x':+y')   =  (x-x') :+ (y-y')+    negate (x:+y)       =  negate x :+ negate y+instance AbelianAdditive a => AbelianAdditive (Complex a)+instance Ring a => Module (Complex a) (Complex a) where+  (*^) = (*)+instance Ring a => Module a (Complex a) where+  s *^ (x :+ y) =  (s *^ x :+ s *^ y)+instance Ring a => Ring (Complex a) where+    fromInteger n  =  fromInteger n :+ zero++instance  Field a => Division (Complex a)  where+    {-# SPECIALISE instance Division (Complex Double) #-}+    (x:+y) / (x':+y')   =  (x*x'+y*y') / d :+ (y*x'-x*y') / d+      where d   = x'*x' + y'*y'++instance Field a => Field (Complex a) where+    fromRational a =  fromRational a :+ zero++{-data Expr a where+  Embed :: a -> Expr a+  Add :: Expr a -> Expr a -> Expr a+  Mul :: Expr a -> Expr a -> Expr a+  Zero :: Expr a+  One :: Expr a+  deriving (Prelude.Show)+++instance Additive (Expr a) where+  zero = Zero+  Zero + x = x+  x + Zero = x+  x + y = Add x y++instance Multiplicative (Expr a) where+  one = One+  One * x = x+  x * One = x+  x * y = Mul x y+-}++-- Syntax+ ifThenElse :: Bool -> t -> t -> t ifThenElse True a _ = a ifThenElse False _ a = a  -data InitialAdditive = InitialAdditive :+ InitialAdditive | Zero-  deriving (Prelude.Show) -instance Additive InitialAdditive where-  zero = Zero-  (+) = (:+)+-- >>> times 5 (Embed "x")+-- Add (Add (Embed "x") (Add (Embed "x") (Embed "x"))) (Add (Embed "x") (Embed "x")) -instance Module Rational Double where-  r *^ d = fromRational r * d++-- >>> (Embed "x")+-- Zero+++{-+Note: the following is not quite what we intuitively want, because++class Field a => AlgebraicallyClosed  a where+  sqrt :: a -> (a,a)++AlgebraicallyClosed numbers have two square roots.++-}+

@@ -21,16 +21,36 @@ {-# LANGUAGE TypeOperators #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE ViewPatterns #-}+{-# LANGUAGE RebindableSyntax #-}  module Algebra.Linear where -import Algebra.Classes-import Prelude (cos,sin,Floating(..),Functor(..),Show(..),Eq(..),Int,fst,flip,($),Ord)+import Algebra.Classes hiding ((*<))+import Prelude (cos,sin,Floating(..),Functor(..),Show(..),Eq(..),Int,fst,($),Ord,Double) import Control.Applicative-import Data.Foldable+import Data.Foldable hiding (sum,product) import Data.Traversable import Control.Monad.State import Algebra.Category++infixr 7 *<++type VectorSpace scalar a = (Field scalar, Module scalar a, Group a)+-- Because of the existence of bases, vector spaces can always be made representable (Traversable, Applicative) functors.+-- So we'd be better off using the following definition:++-- | Representation of vector as traversable functor+type VectorR v = (Applicative v,Traversable v)+-- ... but this is missing the link with *^ for module.  We should be+-- able to add forall s. PreRing s => Module s (v s), but GHC does not+-- like it. (In fact, QuantifiedConstraints is very buggy in ghc 8.6)++class VectorR v => InnerProdSpace v where+  inner :: Field s => v s -> v s -> s++--------------------------------------------------------------+-- Construction of finite vectors+ data VZero a = VZero deriving (Functor,Foldable,Traversable,Show,Eq,Ord) instance Applicative VZero where   pure _ = VZero@@ -53,7 +73,10 @@ pattern V3' :: forall a. a -> a -> a -> V3' a pattern V3' x y z = VNext (V2' x y) z --- | Make a Euclidean vector out of a traversable functor+--------------------------------------------+-- Euclidean spaces with a (inner product)++-- | Make a Euclidean vector out of a traversable functor. (The p) newtype Euclid f a = Euclid {fromEuclid :: f a} deriving (Functor,Foldable,Traversable,Show,Eq,Ord,Applicative)  type V3 = Euclid V3'@@ -78,49 +101,38 @@ pureMat :: (Applicative v, Applicative w) => s -> Mat s v w pureMat x = Mat (pure (pure x)) -(>*<) :: (Applicative v, Applicative w) =>-               Mat (a -> s) v w -> Mat a v w -> Mat s v w-Mat f >*< Mat x = Mat (((<*>) <$> f) <*> x)--(>$<) :: (Applicative v, Applicative w) =>-               (a -> s) -> Mat a v w -> Mat s v w-f >$< x = pureMat f >*< x- instance (Applicative f,Applicative g,Additive a) => Additive (Mat a f g) where   zero = pureMat zero-  x + y =  (+) >$< x >*< y+  x + y = matFlat ((+) <$> flatMat x <*> flatMat y) instance (Applicative f,Applicative g,AbelianAdditive a) => AbelianAdditive (Mat a f g) where instance (Applicative f,Applicative g,Group a) => Group (Mat a f g) where-  negate x = negate >$< x-  x - y = (-) >$< x >*< y+  negate x = matFlat (negate <$> flatMat x)+  x - y = matFlat ((-) <$> flatMat x <*> flatMat y)  instance (Applicative f, Applicative g,Module s a) => Module s (Mat a f g) where   s *^ Mat t = Mat (((s*^) <$>) <$> t) -class VectorSpace (Scalar v) v => InnerProdSpace v where-  type Scalar v-  dotProd :: v -> v -> Scalar v +-- | Hadamard product (⊙) :: Applicative v => Multiplicative s => v s -> v s -> v s x ⊙ y = (*) <$> x <*> y -instance (Ring a, Field a, Applicative f, Foldable f) => InnerProdSpace (Euclid f a) where-  type Scalar (Euclid f a) = a-  dotProd x y = add (x ⊙ y)+instance (VectorR f) => InnerProdSpace (Euclid f) where+  inner x y = sum (x ⊙ y) -- fixme -(·) :: InnerProdSpace v => v -> v -> Scalar v-(·) = dotProd+(·) :: Field s => InnerProdSpace v => v s -> v s -> s+(·) = inner -sqNorm :: InnerProdSpace v => v -> Scalar v-sqNorm x = dotProd x x+sqNorm :: Field s => InnerProdSpace v => v s -> s+sqNorm x = inner x x -norm :: InnerProdSpace v => Floating (Scalar v) => v  -> Scalar v+norm :: Field s => InnerProdSpace v => Floating s => v s  -> s norm = sqrt . sqNorm -normalize :: Floating (Scalar v) => InnerProdSpace v => v -> v+normalize :: (VectorSpace s (v s)) => Floating s => InnerProdSpace v => v s -> v s normalize v = recip (norm v) *^ v --- | Cross product https://en.wikipedia.org/wiki/Cross_product+-- | Cross product in 3 dimensions https://en.wikipedia.org/wiki/Cross_product (×) :: Ring a => V3 a -> V3 a -> V3 a (V3 a1 a2 a3) × (V3 b1 b2 b3) = V3 (a2*b3 - a3*b2)  (negate (a1*b3 - a3*b1)) (a1*b2 - a2*b1) @@ -129,52 +141,77 @@   where increment = do x <- get; put (x+1); return x  type SqMat v s = Mat s v v-newtype Mat s w v = Mat {fromMat :: v (w s)} deriving Show +-- | Matrix type. (w s) is a column. (v s) is a row.+newtype Mat s w v = Mat {fromMat :: w (v s)} deriving Show++-- | View of the matrix as a composition of functors.+newtype Flat w v s = Flat {fromFlat :: w (v s)} deriving (Show,Functor,Foldable)+flatMat :: Mat s w v -> Flat w v s+flatMat (Mat x) = (Flat x)+matFlat :: Flat w v s -> Mat s w v+matFlat (Flat x) = (Mat x)++instance (Applicative w, Applicative v) => Applicative (Flat w v) where+  pure x = Flat (pure (pure x))+  Flat f <*> Flat a = Flat (((<*>) <$> f) <*> a)++ instance Ring s => Category (Mat s) where-  type Con v = (Applicative v, Traversable v)+  type Con v = VectorR v   (.) = matMul   id = identity -type Mat3x3 s = SqMat V3' s-type Mat2x2 s = SqMat V2' s+type Mat3x3 s = SqMat V3 s+type Mat2x2 s = SqMat V2 s -pattern Mat2x2 :: forall s. s -> s -> s -> s -> Mat s V2' V2'-pattern Mat2x2 a b c d = Mat (V2' (V2' a b) (V2' c d))+pattern Mat2x2 :: forall s. s -> s -> s -> s -> Mat s V2 V2+pattern Mat2x2 a b c d = Mat (V2 (V2 a c)+                                 (V2 b d)) -pattern Mat3x3 :: forall s. s -> s -> s -> s -> s -> s -> s -> s -> s -> Mat s V3' V3'-pattern Mat3x3 a b c d e f g h i = Mat (V3' (V3' a b c) (V3' d e f) (V3' g h i))+pattern Mat3x3 :: forall s. s -> s -> s -> s -> s -> s -> s -> s -> s -> Mat s V3 V3+pattern Mat3x3 a b c d e f g h i = Mat (V3 (V3 a d g)+                                           (V3 b e h)+                                           (V3 c f i)) -matVecMul :: (Foldable f1, Ring b, Applicative f1, Functor f2) => Mat b f1 f2 -> Euclid f1 b -> Euclid f2 b-matVecMul (Mat m) v = Euclid (euclideanDotProd v <$> (Euclid <$> m))-   where euclideanDotProd x y = add (Euclid x ⊙ Euclid y)+-- | Vector scaling. If Module a (f a), then (*^) must be the same as (*<).+(*<) :: (Functor f, Multiplicative b) => b -> f b -> f b+s *< v = (s*) <$> v -rotation2d :: Floating a => a -> Mat2x2 a-rotation2d θ = Mat $ V2' (V2' (cos θ) (-sin θ))-                         (V2' (sin θ)  (cos θ)) +(<+>) :: (Applicative f, Additive b) => f b -> f b -> f b+u <+> v = (+) <$> u <*> v+++matVecMul :: forall s v w. (Ring s, Foldable v,Applicative v,Applicative w) => Mat s v w -> v s -> w s+matVecMul (Mat m) x = foldr (<+>) (pure zero) ((*<) <$> x <*> m) -- If GHC gets fixed: use VectorR constraint instead of Applicative, and add instead of foldr.++rotation2d :: (Group a,Floating a) => a -> Mat2x2 a+rotation2d θ = transpose $ Mat $ V2 (V2 (cos θ) (-sin θ))+                                    (V2 (sin θ)  (cos θ))+ -- >>> rotation2d (pi/2) -- Mat {fromMat = V2' (V2' 6.123233995736766e-17 (-1.0)) (V2' 1.0 6.123233995736766e-17)}  crossProductMatrix :: Group a => V3 a -> Mat3x3 a-crossProductMatrix (V3 a1 a2 a3) = Mat (V3'  (V3' zero (negate a3) a2)-                                             (V3' a3 zero (negate a1))-                                             (V3' (negate a2) a1 zero))+crossProductMatrix (V3 a1 a2 a3) = Mat3x3 zero  (-a3) a2+                                          a3    zero  (-a1)+                                          (-a2) a1    zero  -- | Tensor product (⊗) :: (Applicative v, Applicative w, Multiplicative s)-    => Euclid w s -> Euclid v s -> Mat s w v-Euclid v1 ⊗ Euclid v2 = flip (tensorWith (*)) v1 v2+    => w s -> v s -> Mat s w v+v1 ⊗ v2 = tensorWith (*) v2 v1  tensorWith :: (Applicative v, Applicative w)            => (s -> t -> u) -> w s -> v t -> Mat u v w-tensorWith f v1 v2 = flip f >$< Mat (pure v2) >*< Mat (pure <$> v1)+tensorWith f v1 v2 = matFlat (f <$> Flat (pure v1) <*> Flat (pure <$> v2))  identity :: Traversable v => Ring s => Applicative v => SqMat v s identity = tensorWith (\x y -> if x == y then one else zero) index index -diagonal :: Traversable v => Ring s => Applicative v => Euclid v s -> SqMat v s-diagonal (Euclid v) = tensorWith (\x (y,a) -> if x == y then a else zero) index ((,) <$> index <*> v)+diagonal :: Traversable v => Ring s => Applicative v => v s -> SqMat v s+diagonal v = tensorWith (\x (y,a) -> if x == y then a else zero) index ((,) <$> index <*> v)  -- | 3d rotation around given axis rotation3d :: Ring a => Floating a => a -> V3 a -> Mat3x3 a@@ -182,27 +219,30 @@                  sin θ *^ crossProductMatrix u +                  (1 - cos θ) *^ (u ⊗ u) +-- | 3d rotation mapping the direction of 'from' to that of 'to' rotationFromTo :: (Floating a, Module a a,Field a)                => V3 a -> V3 a -> Mat3x3 a rotationFromTo from to = c *^ identity + s *^ crossProductMatrix v + (1-c) *^ (v ⊗ v)   where y = to         x = from-        v = x × y-        c = dotProd x y-        s = norm v+        v = x × y -- axis of rotation+        c = inner x y -- cos of angle+        s = norm v -- sin of angle -transpose :: Applicative f => Traversable g => Mat a f g -> Mat a g f-transpose = Mat . sequenceA . fromMat+-- >>> let u = (V3 (1::Double) 0 0); v = (V3 0 1 1); in (rotationFromTo u v) `matVecMul` u+-- Euclid {fromEuclid = VNext (VNext (VNext VZero 0.0) 1.4142135623730951) 1.4142135623730951} -matMul' :: (Traversable u, Ring s, Applicative w, Applicative v, Applicative u) => Mat s v u -> Mat s u w -> Mat s v w-matMul' (transpose -> Mat y) (Mat x)  = tensorWith (\a b -> add (a ⊙ b)) x y+transpose :: Applicative g => Traversable f => Mat a f g -> Mat a g f+transpose = Mat . sequenceA . fromMat  matMul :: (Traversable u, Ring s, Applicative w, Applicative v, Applicative u) => Mat s u w -> Mat s v u -> Mat s v w-matMul = flip matMul'+matMul a (Mat b) = Mat (matVecMul a <$> b)  --- >>> let t1 = rotation2d (1::Double) in matMul t1 (transpose t1)--- Mat {fromMat = V2' (V2' 1.0 0.0) (V2' 0.0 1.0)}+-- >>> let t1 = rotation2d (1::Double) in matMul (transpose t1) t1+-- Mat {fromMat = VNext (VNext VZero (VNext (VNext VZero 1.0) 0.0)) (VNext (VNext VZero 0.0) 1.0)}++  -- The group of Orthogonal matrices, using "Multiplicative" for respecting conventions a bit better newtype OrthoMat v s = OrthoMat (SqMat v s)

@@ -1,5 +1,5 @@ name:           gasp-version:        1.2.0.0+version:        1.3.0.0 category:       Algebra synopsis:       A framework of algebraic classes description:@@ -22,6 +22,7 @@   build-depends: containers   build-depends: binary   build-depends: mtl+  build-depends: QuickCheck    exposed-modules:        Algebra.Classes