HaskellForMaths 0.2.2 → 0.3.1

raw patch · 29 files changed

+2291/−56 lines, 29 filesdep ~basePVP ok

version bump matches the API change (PVP)

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• HaskellForMaths.cabal +13/−4
• Math/Algebra/Commutative/MPoly.hs +1/−0
• Math/Algebra/Group/PermutationGroup.hs +4/−2
• Math/Algebra/Group/Subquotients.hs +20/−5
• Math/Algebra/NonCommutative/NCPoly.hs +10/−4
• Math/Algebras/AffinePlane.hs +90/−0
• Math/Algebras/Commutative.hs +160/−0
• Math/Algebras/GroupAlgebra.hs +54/−0
• Math/Algebras/LaurentPoly.hs +81/−0
• Math/Algebras/Matrix.hs +91/−0
• Math/Algebras/NonCommutative.hs +149/−0
• Math/Algebras/Quaternions.hs +58/−0
• Math/Algebras/Structures.hs +160/−0
• Math/Algebras/TensorAlgebra.hs +54/−0
• Math/Algebras/TensorProduct.hs +47/−0
• Math/Algebras/VectorSpace.hs +130/−0
• Math/Combinatorics/Design.hs +44/−12
• Math/Combinatorics/Hypergraph.hs +31/−12
• Math/Combinatorics/LatinSquares.hs +8/−1
• Math/Combinatorics/StronglyRegularGraph.hs +2/−1
• Math/Projects/ChevalleyGroup/Classical.hs +13/−5
• Math/Projects/ChevalleyGroup/Exceptional.hs +15/−10
• Math/QuantumAlgebra/OrientedTangle.hs +210/−0
• Math/QuantumAlgebra/QuantumPlane.hs +236/−0
• Math/QuantumAlgebra/Tangle.hs +177/−0
• Math/QuantumAlgebra/TensorCategory.hs +128/−0
• Math/Test/TAlgebras/TGroupAlgebra.hs +55/−0
• Math/Test/TAlgebras/TQuaternions.hs +39/−0
• Math/Test/TAlgebras/TStructures.hs +211/−0

@@ -1,5 +1,5 @@    Name:                HaskellForMaths
-   Version:             0.2.2
+   Version:             0.3.1
    Category:            Math
    Description:         A library of maths code in the areas of combinatorics, group theory, commutative algebra, and non-commutative algebra. The library is mainly intended for educational purposes, but does have efficient implementations of several fundamental algorithms.
    Synopsis:            Combinatorics, group theory, commutative algebra, non-commutative algebra
@@ -7,7 +7,7 @@    License-file:        license.txt
    Author:              David Amos
    Maintainer:          haskellformaths-at-gmail-dot-com
-   Homepage:            http://www.polyomino.f2s.com/haskellformathsv2/HaskellForMathsv2.html
+   Homepage:            http://haskellformaths.blogspot.com/
    Build-Type:          Simple
    Cabal-Version:       >=1.2
 
@@ -22,9 +22,12 @@         Math/Test/TRootSystem.hs,
         Math/Test/TSubquotients.hs,
         Math/Test/TestAll.hs
+        Math/Test/TAlgebras/TStructures.hs
+        Math/Test/TAlgebras/TQuaternions.hs
+        Math/Test/TAlgebras/TGroupAlgebra.hs
 
    Library
-     Build-Depends:     base >=3 && < 4, containers, array, random, QuickCheck
+     Build-Depends:     base >= 2 && < 5, containers, array, random, QuickCheck
      Exposed-modules:
         Math.Algebra.LinearAlgebra,
         Math.Algebra.Commutative.Monomial,
@@ -34,6 +37,10 @@         Math.Algebra.Group.RandomSchreierSims, Math.Algebra.Group.Subquotients,
         Math.Algebra.Group.StringRewriting, Math.Algebra.Group.CayleyGraph,
         Math.Algebra.NonCommutative.NCPoly, Math.Algebra.NonCommutative.GSBasis, Math.Algebra.NonCommutative.TensorAlgebra,
+        Math.Algebras.AffinePlane, Math.Algebras.Commutative, Math.Algebras.GroupAlgebra,
+        Math.Algebras.LaurentPoly, Math.Algebras.Matrix, Math.Algebras.NonCommutative,
+        Math.Algebras.Quaternions, Math.Algebras.Structures, Math.Algebras.TensorAlgebra,
+        Math.Algebras.TensorProduct, Math.Algebras.VectorSpace,
         Math.Combinatorics.Graph, Math.Combinatorics.GraphAuts, Math.Combinatorics.StronglyRegularGraph,
         Math.Combinatorics.Design, Math.Combinatorics.FiniteGeometry, Math.Combinatorics.Hypergraph,
         Math.Combinatorics.LatinSquares,
@@ -42,6 +49,8 @@         Math.Projects.Rubik, Math.Projects.MiniquaternionGeometry,
         Math.Projects.ChevalleyGroup.Classical, Math.Projects.ChevalleyGroup.Exceptional,
         Math.Projects.KnotTheory.Braid,
-        Math.Projects.KnotTheory.LaurentMPoly, Math.Projects.KnotTheory.TemperleyLieb, Math.Projects.KnotTheory.IwahoriHecke
+        Math.Projects.KnotTheory.LaurentMPoly, Math.Projects.KnotTheory.TemperleyLieb, Math.Projects.KnotTheory.IwahoriHecke,
+        Math.QuantumAlgebra.OrientedTangle,
+        Math.QuantumAlgebra.QuantumPlane, Math.QuantumAlgebra.Tangle, Math.QuantumAlgebra.TensorCategory
 
      ghc-options:       -w

@@ -2,6 +2,7 @@ 
 {-# OPTIONS_GHC -fglasgow-exts #-}
 
+-- |A module providing a type for (commutative) multivariate polynomials, with support for various term orders.
 module Math.Algebra.Commutative.MPoly where
 
 import qualified Data.Map as M

@@ -8,12 +8,14 @@ 
 import Math.Common.ListSet (toListSet, union, (\\) ) -- a version of union which assumes the arguments are ascending sets (no repeated elements)
 
+infix 8 ^-, ~^
+
 rotateL (x:xs) = xs ++ [x]
 
 
 -- PERMUTATIONS
 
--- |Type for permutations, considered as group elements.
+-- |A type for permutations, considered as functions or actions which can be performed on an underlying set.
 newtype Permutation a = P (M.Map a a) deriving (Eq,Ord)
 
 fromPairs xys | isValid   = fromPairs' xys
@@ -41,7 +43,7 @@            Just y  -> y
            Nothing -> x -- if x `notElem` supp (P g), then x is not moved
 
--- |b -^ g returns the image of an edge or block b under the action of the permutation g
+-- |b -^ g returns the image of an edge or block b under the action of the permutation g.
 -- The dash is meant to be a mnemonic for edge or line or block.
 (-^) :: (Ord t) => [t] -> Permutation t -> [t]
 xs -^ g = L.sort [x .^ g | x <- xs]

@@ -58,6 +58,7 @@ -}  +isTransitive :: (Ord t) => [Permutation t] -> Bool isTransitive gs = length (orbits gs) == 1  @@ -75,8 +76,12 @@ -}  -- Seress p81--- A transitive constituent homomorphism is the restriction of G <= Sym(omega) to an orbit delta <= omega--- This function returns the kernel and the image+-- |Given a group gs and a transitive constituent ys, return the kernel and image of the transitive constituent homomorphism.+-- That is, suppose that gs acts on a set xs, and ys is a subset of xs on which gs acts transitively.+-- Then the transitive constituent homomorphism is the restriction of the action of gs to an action on the ys.+transitiveConstituentHomomorphism+  :: (Ord a, Show a) =>+     [Permutation a] -> [a] -> ([Permutation a], [Permutation a]) transitiveConstituentHomomorphism gs delta     | delta == closure delta [(.^ g) | g <- gs] -- delta is closed under action of gs, hence a union of orbits         = transitiveConstituentHomomorphism' gs delta@@ -112,14 +117,17 @@ -- Because the support of the permutations is not constrained to be [1..n], we have to use a map instead of an array -- This probably affects the complexity, but isn't a problem in practice --- Find all block systems+-- |Given a transitive group gs, find all non-trivial block systems. That is, if gs act on xs,+-- find all the ways that the xs can be divided into blocks, such that the gs also have a permutation action on the blocks+blockSystems :: (Ord t) => [Permutation t] -> [[[t]]] blockSystems gs     | isTransitive gs = toListSet $ filter (/= [x:xs]) $ map (minimalBlock gs) [ [x,x'] | x' <- xs ]     | otherwise = error "blockSystems: not transitive"     where x:xs = foldl union [] $ map supp gs  --- More efficient version if we have an sgs+-- |A more efficient version of blockSystems, if we have an sgs+blockSystemsSGS :: (Ord a) => [Permutation a] -> [[[a]]] blockSystemsSGS gs = toListSet $ filter (/= [x:xs]) $ map (minimalBlock gs) [ [x,x'] | x' <- rs ]     where x:xs = foldl union [] $ map supp gs           hs = filter (\g -> x < minsupp g) gs -- sgs for stabiliser Gx@@ -133,14 +141,21 @@ -- see Holt RandomStab function  +-- |A permutation group is primitive if it has no non-trivial block systems+isPrimitive :: (Ord t) => [Permutation t] -> Bool isPrimitive gs = null (blockSystems gs) +isPrimitiveSGS :: (Ord a) => [Permutation a] -> Bool isPrimitiveSGS gs = null (blockSystemsSGS gs)  -- There are other optimisations we haven't done -- see Holt p86 -+-- |Given a transitive group gs, and a block system for gs, return the kernel and image of the block homomorphism+-- (the homomorphism onto the action of gs on the blocks)+blockHomomorphism+  :: (Ord t, Show t) =>+     [Permutation t] -> [[t]] -> ([Permutation t], [Permutation [t]]) blockHomomorphism gs bs     | bs == closure bs [(-^ g) | g <- gs] -- bs is closed under action of gs         = blockHomomorphism' gs bs

@@ -1,5 +1,6 @@ -- Copyright (c) David Amos, 2008. All rights reserved.
 
+-- |A module providing a type for non-commutative polynomials.
 module Math.Algebra.NonCommutative.NCPoly where
 
 import Data.List as L
@@ -14,10 +15,13 @@     compare (M xs) (M ys) = compare (length xs,xs) (length ys,ys)
 -- Glex ordering
 
-instance Show v => Show (Monomial v) where
+instance (Eq v, Show v) => Show (Monomial v) where
     show (M xs) | null xs = "1"
-                | otherwise = concatMap show xs
--- !! we can do better - we should show "xx" as "x^2"
+                | otherwise = concatMap showPower (L.group xs)
+        where showPower [v] = showVar v
+              showPower vs@(v:_) = showVar v ++ "^" ++ show (length vs)
+              showVar v = filter (/= '"') (show v)
+-- Taken from NonComMonomial - why don't we just use it directly
 
 instance (Eq v, Show v) => Num (Monomial v) where
     M xs * M ys = M (xs ++ ys)
@@ -54,7 +58,6 @@                                  then "+(" ++ c:cs ++ ")"
                                  else if c == '-' then c:cs else '+':c:cs
 
-
 instance (Ord v, Show v, Num r) => Num (NPoly r v) where
     NP ts + NP us = NP (mergeTerms ts us)
     negate (NP ts) = NP $ map (\(m,c) -> (m,-c)) ts
@@ -97,6 +100,9 @@     show Y = "y"
     show Z = "z"
 
+-- |Create a non-commutative variable for use in forming non-commutative polynomials.
+-- For example, we could define x = var "x", y = var "y". Then x*y /= y*x.
+var :: (Num k) => v -> NPoly k v
 var v = NP [(M [v], 1)]
 
 x = var X :: NPoly Q Var

@@ -0,0 +1,90 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}++-- |A module defining the affine plane and its symmetries+module Math.Algebras.AffinePlane where++import Math.Algebra.Field.Base hiding (powers)+import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Structures+import Math.Algebras.Commutative+++data XY = X | Y deriving (Eq, Ord)++instance Show XY where show X = "x"; show Y = "y"++x = glexVar X :: GlexPoly Q XY+y = glexVar Y :: GlexPoly Q XY+++data ABCD = A | B | C | D deriving (Eq, Ord)++instance Show ABCD where show A = "a"; show B = "b"; show C = "c"; show D = "d"++a,b,c,d :: Monomial m => Vect Q (m ABCD)+a = var A+b = var B+c = var C+d = var D+++-- SL2++newtype SL2 v = SL2 (GlexMonomial v) deriving (Eq,Ord)++instance Show v => Show (SL2 v) where show (SL2 m) = show m++instance Algebra Q (SL2 ABCD) where -- to do this for Num k instead of Q we would need a,b,c,d defined for Num k+    unit 0 = zero -- V []+    unit x = V [(munit,x)] where munit = SL2 (Glex 0 [])+    mult x = x''' where+        x' = mult $ fmap (\(T (SL2 a) (SL2 b)) -> T a b) x -- perform the multiplication in GlexPoly+        x'' = x' %% [a*d-b*c-1] -- :: GlexPoly Q ABCD] -- quotient by ad-bc=1 in GlexPoly Q ABCD+        x''' = fmap SL2 x'' -- ie wrap the monomials up as SL2 again+        -- mmult (Glex si xis) (Glex sj yjs) = Glex (si+sj) $ addmerge xis yjs++sl2Var v = V [(SL2 (Glex 1 [(v,1)]), 1)] -- :: Vect Q (SL2 ABCD)+++-- For example:+-- > a*d :: Vect Q (SL2 ABCD)+-- bc+1++instance Monomial SL2 where+    var = sl2Var+    powers (SL2 (Glex _ xis)) = xis++++instance Coalgebra Q (SL2 ABCD) where+    counit x = case x `bind` cu of+               V [] -> 0+               V [(SL2 (Glex 0 []), c)] -> c+        where cu A = 1 :: Vect Q (SL2 ABCD)+              cu B = 0+              cu C = 0+              cu D = 1+    comult x = x `bind` cm+        where cm A = a `te` a + b `te` c+              cm B = a `te` b + b `te` d+              cm C = c `te` a + d `te` c+              cm D = c `te` b + d `te` d+-- In other words+-- counit (a b) = (1 0)+--        (c d)   (0 1)+-- comult (a b) = (a1 b1) `te` (a2 b2)+--        (c d)   (c1 d1)      (c2 d2)++instance Bialgebra Q (SL2 ABCD) where {}++instance HopfAlgebra Q (SL2 ABCD) where+    antipode x = x `bind` antipode'+        where antipode' A = d+              antipode' B = b+              antipode' C = c+              antipode' D = a+-- in the GL2 case we would need 1/det factor as well+

@@ -0,0 +1,160 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}++-- |A module defining the algebra of commutative polynomials over a field k+module Math.Algebras.Commutative where++import Math.Algebra.Field.Base hiding (powers)+import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Structures+++-- GLEX MONOMIALS++data GlexMonomial v = Glex Int [(v,Int)] deriving (Eq)+-- The initial Int is the degree of the monomial. Storing it speeds up equality tests and comparisons++-- type GlexMonomialS = GlexMonomial String++instance Ord v => Ord (GlexMonomial v) where+    compare (Glex si xis) (Glex sj yjs) = compare (-si, xis) (-sj, yjs)++instance Show v => Show (GlexMonomial v) where+    show (Glex _ []) = "1"+    show (Glex _ xis) = concatMap (\(x,i) -> if i==1 then showVar x else showVar x ++ "^" ++ show i) xis+        where showVar x = filter ( /= '"' ) (show x) -- in case v == String++{-+-- GlexMonomial is a functor and a monad+-- However, this isn't all that much use, and to make proper use of it we'd need a "nf" function+-- So leaving this commented out++-- map one basis to another+instance Functor GlexMonomial where+    fmap f (Glex si xis) = Glex si [(f x, i) | (x,i) <- xis]+-- Note that as we can't assume the Ord instance, we would need to call "nf" afterwards++-- GlexMonomial is the free commutative monoid, and hence a monad+instance Monad GlexMonomial where+    return x = Glex 1 [(x,1)]+    (Glex _ xis) >>= f = let parts = [(i, sj, yjs) | (x,i) <- xis, let Glex sj yjs = f x]+                         in Glex (sum [i*sj | (i,sj,_) <- parts])+                                 (concatMap (\(i,_,yjs)->map (\(y,j)->(y,i*j)) yjs) parts)+    -- this isn't really much use - it's variable substitution, but we're only allowed to substitute monomials for each var+-- Note that as we can't assume the Ord instance, we would need to call "nf" afterwards+-}++-- This is the monoid algebra for commutative monomials (which are the free commutative monoid)+instance (Num k, Ord v) => Algebra k (GlexMonomial v) where+    unit 0 = zero -- V []+    unit x = V [(munit,x)] where munit = Glex 0 []+    mult (V ts) = nf $ fmap (\(T a b) -> a `mmult` b) (V ts)+        where mmult (Glex si xis) (Glex sj yjs) = Glex (si+sj) $ addmerge xis yjs++{-+-- This is just the Set Coalgebra, so better to use a generic instance+-- Also, not used anywhere. Hence commented out+instance Num k => Coalgebra k (GlexMonomial v) where+    counit (V ts) = sum [x | (m,x) <- ts] -- trace+    comult = fmap (\m -> T m m) -- diagonal+    -- comult (V ts) = V [(T m m, x) | (m, x) <- ts]+-}+type GlexPoly k v = Vect k (GlexMonomial v)++++glexVar v = V [(Glex 1 [(v,1)], 1)]+++class Monomial m where+    var :: v -> Vect Q (m v)+    powers :: m v -> [(v,Int)]++-- |In effect, we have (Num k, Monomial m) => Monad (\v -> Vect k (m v)), with return = var, and (>>=) = bind.+-- However, we can't express this directly in Haskell, firstly because of the Ord b constraint,+-- secondly because Haskell doesn't support type functions.+bind :: (Monomial m, Num k, Ord b, Show b, Algebra k b) =>+     Vect k (m v) -> (v -> Vect k b) -> Vect k b+V ts `bind` f = sum [c `smultL` product [f x ^ i | (x,i) <- powers m] | (m, c) <- ts] +-- V ts `bind` f = sum [product [f x ^ i | (x,i) <- powers m] * unit c | (m, c) <- ts] +++instance Monomial GlexMonomial where+    var = glexVar+    powers (Glex _ xis) = xis+++-- DIVISION++lt (V (t:ts)) = t++class DivisionBasis b where+    dividesB :: b -> b -> Bool+    divB :: b -> b -> b++dividesT (b1,x1) (b2,x2) = dividesB b1 b2+divT (b1,x1) (b2,x2) = (divB b1 b2, x1/x2)++-- given f, gs, find as, r such that f = sum (zipWith (*) as gs) + r, with r not divisible by any g+quotRemMP f gs = quotRemMP' f (replicate n 0, 0) where+    n = length gs+    quotRemMP' 0 (us,r) = (us,r)+    quotRemMP' h (us,r) = divisionStep h (gs,[],us,r)+    divisionStep h (g:gs,us',u:us,r) =+        if lt g `dividesT` lt h+        then let t = V [lt h `divT` lt g]+                 h' = h - t*g+                 u' = u+t+             in quotRemMP' h' (reverse us' ++ u':us, r)+        else divisionStep h (gs,u:us',us,r)+    divisionStep h ([],us',[],r) =+        let (lth,h') = splitlt h+        in quotRemMP' h' (reverse us', r+lth)+    splitlt (V (t:ts)) = (V [t], V ts)++infixl 7 %%++(%%) :: (Fractional k, Ord b, Show b, Algebra k b, DivisionBasis b)+     => Vect k b -> [Vect k b] -> Vect k b+f %% gs = r where (_,r) = quotRemMP f gs+++instance Ord v => DivisionBasis (GlexMonomial v) where+    dividesB (Glex si xis) (Glex sj yjs) = si <= sj && dividesB' xis yjs where+        dividesB' ((x,i):xis) ((y,j):yjs) =+            case compare x y of+            LT -> False+            GT -> dividesB' ((x,i):xis) yjs+            EQ -> if i<=j then dividesB' xis yjs else False+        dividesB' [] _ = True+        dividesB' _ [] = False+    divB (Glex si xis) (Glex sj yjs) = Glex (si-sj) $ divB' xis yjs where+        divB' ((x,i):xis) ((y,j):yjs) =+            case compare x y of+            LT -> (x,i) : divB' xis ((y,j):yjs)+            EQ -> if i == j then divB' xis yjs else (x,i-j) : divB' xis yjs -- we don't bother to check i > j+            GT -> error "divB'" -- (y,-j) : divB' ((x,i):xis) yjs+        divB' xis [] = xis+        divB' [] yjs = error "divB'"++{-+-- Need to thread this through Maybe properly, so perhaps use do notation+divB2 (Glex si xis) (Glex sj yjs)+    | si < sj = Nothing+    | otherwise = case divB' xis yjs of+                  Nothing -> Nothing+                  Just zks -> Glex (si-sj) zks+    where divB' ((x,i):xis) ((y,j):yjs) =+              case compare x y of+              LT -> (x,i) : divB' xis ((y,j):yjs)+              EQ -> case compare i j of+                    LT -> Nothing+                    EQ -> divB' xis yjs+                    GT -> (x,i-j) : divB' xis yjs+              GT -> Nothing+-}+-- !! could change divB to return Maybe, and avoid need for dividesB++

@@ -0,0 +1,54 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE  MultiParamTypeClasses, FlexibleInstances #-}++module Math.Algebras.GroupAlgebra where++import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Structures++import Math.Algebra.Group.PermutationGroup hiding (action)++import Math.Algebra.Field.Base+++instance Mon (Permutation Int) where+    munit = 1+    mmult = (*)++type GroupAlgebra k = Vect k (Permutation Int)++-- Monoid Algebra instance+instance Num k => Algebra k (Permutation Int) where+    unit 0 = zero -- V []+    unit x = V [(munit,x)]+    mult = nf . fmap (\(T a b) -> a `mmult` b)++-- Set Coalgebra instance+-- instance SetCoalgebra (Permutation Int) where {}++instance Num k => Coalgebra k (Permutation Int) where+    counit (V ts) = sum [x | (m,x) <- ts] -- trace+    comult = fmap (\m -> T m m) -- diagonal++instance Num k => Bialgebra k (Permutation Int) where {}+-- should check that the algebra and coalgebra structures are compatible++instance (Num k) => HopfAlgebra k (Permutation Int) where+    antipode (V ts) = nf $ V [(g^-1,x) | (g,x) <- ts]++-- inject permutation into group algebra+ip :: [[Int]] -> GroupAlgebra Q+ip cs = return $ p cs+++instance Num k => Module k (Permutation Int) Int where+    action = nf . fmap (\(T g x) -> x .^ g)++-- use *. instead+-- r *> m = action (r `te` m)++++

@@ -0,0 +1,81 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE TypeSynonymInstances, FlexibleInstances, MultiParamTypeClasses #-}+++module Math.Algebras.LaurentPoly where+++import Math.Algebra.Field.Base hiding (powers)++import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Structures+import qualified Data.List as L++import Math.Algebras.Commutative -- for DivisionBasis and quotRemMP+++-- LAURENT MONOMIALS++data LaurentMonomial = LM Int [(String,Int)] deriving (Eq,Ord)+{-+instance Ord LaurentMonomial where+    compare (LM si xis) (LM sj yjs) = compare (-si, xis) (-sj, yjs)+-}+instance Show LaurentMonomial where+    show (LM 0 []) = "1"+    show (LM _ xis) = concatMap (\(x,i) -> if i==1 then x else x ++ "^" ++ show i) xis++instance Mon LaurentMonomial where+    munit = LM 0 []+    mmult (LM si xis) (LM sj yjs) = LM (si+sj) $ addmerge xis yjs++instance Num k => Algebra k LaurentMonomial where+    unit 0 = zero -- V []+    unit x = V [(munit,x)] +    mult (V ts) = nf $ fmap (\(T a b) -> a `mmult` b) (V ts)+    -- mult (V ts) = nf $ V [(a `mmult` b, x) | (T a b, x) <- ts]++{-+-- This is just the Set Coalgebra, so better to use a generic instance+-- Also, not used anywhere. Hence commented out+instance Num k => Coalgebra k LaurentMonomial where+    counit (V ts) = sum [x | (m,x) <- ts] -- trace+    comult = fmap (\m -> T m m)+-}++type LaurentPoly k = Vect k LaurentMonomial++lvar v = V [(LM 1 [(v,1)], 1)] :: LaurentPoly Q++instance Fractional k => Fractional (LaurentPoly k) where+    recip (V [(LM si xis,c)]) = V [(LM (-si) $ map (\(x,i)->(x,-i)) xis, recip c)]+    recip _ = error "LaurentPoly.recip: only defined for single terms"++q = lvar "q"+q' = 1/q+++{-+-- division doesn't terminate with the derived Ord instance+-- if we use the graded Ord instance instead, division doesn't continue into negative powers+-- so we get the negative powers as remained, even if they're divisible+instance DivisionBasis LaurentMonomial where+    dividesB (LM si xis) (LM sj yjs) = si <= sj && dividesB' xis yjs where+        dividesB' ((x,i):xis) ((y,j):yjs) =+            case compare x y of+            LT -> False+            GT -> dividesB' ((x,i):xis) yjs+            EQ -> if i<=j then dividesB' xis yjs else False+        dividesB' [] _ = True+        dividesB' _ [] = False+    divB (LM si xis) (LM sj yjs) = LM (si-sj) $ divB' xis yjs where+        divB' ((x,i):xis) ((y,j):yjs) =+            case compare x y of+            LT -> (x,i) : divB' xis ((y,j):yjs)+            EQ -> if i == j then divB' xis yjs else (x,i-j) : divB' xis yjs -- we don't bother to check i > j+            GT -> error "divB'" -- (y,-j) : divB' ((x,i):xis) yjs+        divB' xis [] = xis+        divB' [] yjs = error "divB'"+-}

@@ -0,0 +1,91 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE  MultiParamTypeClasses, FlexibleInstances #-}+++module Math.Algebras.Matrix where++import Math.Algebra.Field.Base+import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Structures+++-- Mat2++delta i j | i == j    = 1+          | otherwise = 0++data Mat2 = E2 Int Int deriving (Eq,Ord,Show)+-- E i j represents the elementary matrix with a 1 at the (i,j) position, and 0s elsewhere++instance Num k => Algebra k Mat2 where+    -- unit 0 = zero -- V []+    unit x = x `smultL` V [(E2 i i, 1) | i <- [1..2] ]+    -- mult ab = nf $ ab >>= mult' where+    mult = linear mult' where+        mult' (T (E2 i j) (E2 k l)) = delta j k `smultL` return (E2 i l)++-- In other words+-- unit x = x (1 0)+--            (0 1)+-- mult (a1 b1) `te` (a2 b2) = (a1 b1) * (a2 b2) = (a b)+--      (c1 d1)      (c2 d2)   (c1 d1)   (c2 d2)   (c d)++instance Num k => Module k Mat2 EBasis where+    -- action ax = nf $ ax >>= action' where+    action = linear action' where+        action' (T (E2 i j) (E k)) = delta j k `smultL` return (E i)++-- In other words+-- action (a b) `te` (x) = (ax+by)+--        (c d)      (y)   (cx+dy)++toMat2 [[a,b],[c,d]] = sum $ zipWith (\x e -> unit x * return e) [a,b,c,d] [E2 1 1, E2 1 2, E2 2 1, E2 2 2]+-- fromMat2++toEB2 [x,y] = foldl add zero $ zipWith (\x e -> x `smultL` return e) [x,y] [E 1, E 2]++toEB xs = foldl add zero $ zipWith (\x e -> x `smultL` return e) xs (map E [1..])++++data Mat2' = E2' Int Int deriving (Eq,Ord,Show)+-- E' i j represents the dual basis element corresponding to E i j++-- Kassel p42+instance Num k => Coalgebra k Mat2' where+    counit (V ts) = sum [xij * delta i j | (E2' i j, xij) <- ts]+    -- comult (V ts) = V $ concatMap (\(E2' i j,xij) -> [(T (E2' i k) (E2' k j), xij) | k <- [1..2]]) ts+    comult = linear (\(E2' i j) -> foldl (<+>) zero [return (T (E2' i k) (E2' k j)) | k <- [1..2]])+-- In other words+-- counit (a b) = (1 0)+--        (c d)   (0 1)+-- comult (a b) = (a1 b1) `te` (a2 b2)+--        (c d)   (c1 d1)      (c2 d2)+-- ??+-- ?? How does this act on Mat2?+-- ?? What is the relationship between this and SL2 ABCD, which it seems to resemble++++-- !! Now do the quickchecks++++data M3 = E3 Int Int deriving (Eq,Ord,Show)+-- E i j represents the elementary matrix with a 1 at the (i,j) position, and 0s elsewhere++instance Num k => Algebra k M3 where+    unit 0 = zero -- V []+    unit x = V [(E3 i i, x) | i <- [1..3] ]+    mult (V ts) = nf $ V $ map (\(T (E3 i j) (E3 k l), x) -> (E3 i l, delta j k * x)) ts++{-+-- Kassel p42+-- In this coalgebra instance, the E3 i j are to be interpreted as the dual basis, not the original basis+instance Num k => Coalgebra k M3 where+    counit (V ts) = sum [xij * delta i j | (E3 i j, xij) <- ts]+    comult (V ts) = V $ concatMap (\(E3 i j,xij) -> [(T (E3 i k) (E3 k j), xij) | k <- [1..3]]) ts+-- (is this order preserving?)+-}

@@ -0,0 +1,149 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}++-- |A module defining the algebra of non-commutative polynomials over a field k+module Math.Algebras.NonCommutative where++import Math.Algebra.Field.Base hiding (powers)+import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Structures+import qualified Data.List as L+++data NonComMonomial v = NCM Int [v] deriving (Eq)++instance Ord v => Ord (NonComMonomial v) where+    compare (NCM lx xs) (NCM ly ys) = compare (-lx, xs) (-ly, ys)+-- ie Glex ordering++instance (Eq v, Show v) => Show (NonComMonomial v) where+    show (NCM _ []) = "1"+    show (NCM _ vs) = concatMap showPower (L.group vs)+        where showPower [v] = showVar v+              showPower vs@(v:_) = showVar v ++ "^" ++ show (length vs)+              showVar v = filter (/= '"') (show v)++instance Mon (NonComMonomial v) where+    munit = NCM 0 []+    mmult (NCM i xs) (NCM j ys) = NCM (i+j) (xs++ys)++instance (Num k, Ord v) => Algebra k (NonComMonomial v) where+    unit 0 = zero -- V []+    unit x = V [(munit,x)]+    mult = nf . fmap (\(T a b) -> a `mmult` b)++{-+-- This is the monoid algebra for non-commutative monomials (which is the free monoid)+instance (Num k, Ord v) => Algebra k (NonComMonomial v) where+    unit 0 = zero -- V []+    unit x = V [(munit,x)] where munit = NCM 0 []+    mult (V ts) = nf $ fmap (\(T a b) -> a `mmult` b) (V ts)+        where mmult (NCM lu us) (NCM lv vs) = NCM (lu+lv) (us++vs)+    -- mult (V ts) = nf $ V [(a `mmult` b, x) | (T a b, x) <- ts]+-}++{-+-- This is just the Set Coalgebra, so better to use a generic instance+-- Also, not used anywhere. Hence commented out+instance Num k => Coalgebra k (NonComMonomial v) where+    counit (V ts) = sum [x | (m,x) <- ts] -- trace+    comult = fmap (\m -> T m m)+-}+++class Monomial m where+    var :: v -> Vect Q (m v)+    powers :: Eq v => m v -> [(v,Int)]++V ts `bind` f = sum [c `smultL` product [f x ^ i | (x,i) <- powers m] | (m, c) <- ts] +-- V ts `bind` f = sum [product [f x ^ i | (x,i) <- powers m] * unit c | (m, c) <- ts] ++instance Monomial NonComMonomial where+    var v = V [(NCM 1 [v],1)]+    powers (NCM _ vs) = map power (L.group vs)+        where power vs@(v:_) = (v,length vs)+++type NCPoly v = Vect Q (NonComMonomial v)++{-+x,y,z :: NCPoly String+x = var "x"+y = var "y"+z = var "z"+-}+++-- DIVISION++class DivisionBasis m where+    divM :: m -> m -> Maybe (m,m)+    -- divM a b tries to find l, r such that a = lbr+{-+    findOverlap :: m -> m -> Maybe (m,m,m)+    -- given two monomials f g, find if possible a,b,c with f=ab g=bc+-}++instance Eq v => DivisionBasis (NonComMonomial v) where+    divM (NCM _ a) (NCM _ b) = divM' [] a where+        divM' ls (r:rs) =+            if b `L.isPrefixOf` (r:rs)+            then Just (ncm $ reverse ls, ncm $ drop (length b) (r:rs))+            else divM' (r:ls) rs+        divM' _ [] = Nothing+{-+    findOverlap (NCM _ xs) (NCM _ ys) = findOverlap' [] xs ys where+        findOverlap' as [] cs = Nothing -- (reverse as, [], cs)+        findOverlap' as (b:bs) cs =+            if (b:bs) `L.isPrefixOf` cs+            then Just (ncm $ reverse as, ncm $ b:bs, ncm $ drop (length (b:bs)) cs)+            else findOverlap' (b:as) bs cs+-}+ncm xs = NCM (length xs) xs++lm (V ((m,c):ts)) = m+lc (V ((m,c):ts)) = c+lt (V (t:ts)) = V [t]++-- given f, gs, find ls, rs, f' such that f = sum (zipWith3 (*) ls gs rs) + f', with f' not divisible by any g+quotRemNP f gs | all (/=0) gs = quotRemNP' f (replicate n (0,0), 0)+               | otherwise = error "quotRemNP: division by zero"+    where+    n = length gs+    quotRemNP' 0 (lrs,f') = (lrs,f')+    quotRemNP' h (lrs,f') = divisionStep h (gs,[],lrs,f')+    divisionStep h (g:gs, lrs', (l,r):lrs, f') =+        case lm h `divM` lm g of+        Just (l',r') -> let l'' = V [(l',lc h / lc g)]+                            r'' = V [(r',1)]+                            h' = h - l'' * g * r''+                        in quotRemNP' h' (reverse lrs' ++ (l+l'',r+r''):lrs, f')+        Nothing -> divisionStep h (gs,(l,r):lrs',lrs,f')+    divisionStep h ([],lrs',[],f') =+        let lth = lt h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms+        in quotRemNP' (h-lth) (reverse lrs', f'+lth)++-- It is only marginally (5-10%) more space/time efficient not to track the (lazily unevaluated) factors+remNP f gs | all (/=0) gs = remNP' f 0+           | otherwise = error "remNP: division by zero"+    where+    n = length gs+    remNP' 0 f' = f'+    remNP' h f' = divisionStep h gs f'+    divisionStep h (g:gs) f' =+        case lm h `divM` lm g of+        Just (l',r') -> let l'' = V [(l',lc h / lc g)]+                            r'' = V [(r',1)]+                            h' = h - l'' * g * r''+                        in remNP' h' f'+        Nothing -> divisionStep h gs f'+    divisionStep h [] f' =+        let lth = lt h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms+        in remNP' (h-lth) (f'+lth)++infixl 7 %%+-- f %% gs = r where (_,r) = quotRemNP f gs+f %% gs = remNP f gs+

@@ -0,0 +1,58 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, NoMonomorphismRestriction #-}+++module Math.Algebras.Quaternions where++import Math.Algebra.Field.Base+import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Structures+++-- QUATERNIONS++data HBasis = One | I | J | K deriving (Eq,Ord)++type Quaternion k = Vect k HBasis++instance Show HBasis where+    show One = "1"+    show I = "i"+    show J = "j"+    show K = "k"++instance (Num k) => Algebra k HBasis where+    unit 0 = zero -- V []+    unit x = V [(One,x)]+    -- mult x = nf (x >>= m)+    mult = linear m+         where m (T One b) = return b+               m (T b One) = return b+               m (T I I) = unit (-1)+               m (T J J) = unit (-1)+               m (T K K) = unit (-1)+               m (T I J) = return K+               m (T J I) = -1 *> return K+               m (T J K) = return I+               m (T K J) = -1 *> return I+               m (T K I) = return J+               m (T I K) = -1 *> return J++i,j,k :: Num k => Quaternion k+i = return I+j = return J+k = return K++{-+-- Set coalgebra instance+instance Num k => Coalgebra k HBasis where+    counit (V ts) = sum [x | (m,x) <- ts] -- trace+    comult = fmap (\m -> T m m)           -- diagonal+-}++instance Num k => Coalgebra k HBasis where+    counit (V ts) = sum [x | (One,x) <- ts]+    comult = linear cm+        where cm m = if m == One then return (T m m) else return (T m One) <+> return (T One m)

@@ -0,0 +1,160 @@+-- Copyright (c) David Amos, 2010. All rights reserved.++{-# LANGUAGE MultiParamTypeClasses, NoMonomorphismRestriction #-}+{-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-}+{-# LANGUAGE IncoherentInstances #-}++module Math.Algebras.Structures where++import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+++-- MONOID++-- |Monoid+class Mon m where+    munit :: m+    mmult :: m -> m -> m+++-- ALGEBRAS, COALGEBRAS, BIALGEBRAS, HOPF ALGEBRAS++-- |"Vect k b is a k-algebra"+class Algebra k b where+    unit :: k -> Vect k b+    mult :: Vect k (Tensor b b) -> Vect k b++-- |"Vect k b is a k-coalgebra"+class Coalgebra k b where+    counit :: Vect k b -> k+    comult :: Vect k b -> Vect k (Tensor b b)+++-- |A bialgebra is an algebra which is also a coalgebra, subject to some compatibility conditions+class (Algebra k b, Coalgebra k b) => Bialgebra k b where {}++class Bialgebra k b => HopfAlgebra k b where+    antipode :: Vect k b -> Vect k b+++instance (Num k, Eq b, Ord b, Show b, Algebra k b) => Num (Vect k b) where+    x+y = add x y+    negate (V ts) = V $ map (\(b,x) -> (b, negate x)) ts+    x*y = mult (x `te` y)+    fromInteger n = unit (fromInteger n)+    abs _ = error "Prelude.Num.abs: inappropriate abstraction"+    signum _ = error "Prelude.Num.signum: inappropriate abstraction"+++-- This is the Frobenius form, provided some conditions are met+-- pairing = counit . mult++{-+-- A class to be used to declare that a type b should be given the set coalgebra structure+class SetCoalgebra b where {}++instance (Num k, SetCoalgebra b) => Coalgebra k b where+    counit (V ts) = sum [x | (m,x) <- ts] -- trace+    comult = fmap (\m -> T m m) -- diagonal+-}+++instance Num k => Algebra k () where+    unit 0 = zero -- V []+    unit x = V [( (),x)]+    mult (V [(T () (),x)]) = V [( (),x)]++instance Num k => Coalgebra k () where+    counit (V []) = 0+    counit (V [( (),x)]) = x+    comult (V [( (),x)]) = V [(T () (),x)]++type Trivial k = Vect k ()++unit' :: (Num k, Algebra k b) => Trivial k -> Vect k b+unit' = unit . unwrap where unwrap = counit :: Num k => Trivial k -> k++counit' :: (Num k, Coalgebra k b) => Vect k b -> Trivial k+counit' = wrap . counit where wrap = unit :: Num k => k -> Trivial k++-- unit' and counit' enable us to form tensors of these functions+++-- Kassel p32+instance (Num k, Ord a, Ord b, Algebra k a, Algebra k b) => Algebra k (Tensor a b) where+    unit 0 = V []+    unit x = x `smultL` (unit 1 `te` unit 1)+    -- mult x = nf $ x >>= m where+    mult = linear m where+        m (T (T a b) (T a' b')) = (mult $ return $ T a a') `te` (mult $ return $ T b b')++-- Kassel p42+instance (Num k, Ord a, Ord b, Coalgebra k a, Coalgebra k b) => Coalgebra k (Tensor a b) where+    counit = counit . (counit' `tf` counit')+    -- counit = counit . linear (\(T x y) -> counit' (return x) * counit' (return y))+    comult = assocL . (id `tf` assocR) . (id `tf` (twist `tf` id))+           . (id `tf` assocL) . assocR . (comult `tf` comult)+++newtype SetCoalgebra b = SC b deriving (Eq,Ord,Show)++instance Num k => Coalgebra k (SetCoalgebra b) where+    counit (V ts) = sum [x | (m,x) <- ts] -- trace+    comult = fmap (\m -> T m m)           -- diagonal+++newtype MonoidCoalgebra m = MC m deriving (Eq,Ord,Show)++instance (Num k, Ord m, Mon m) => Coalgebra k (MonoidCoalgebra m) where+    counit (V ts) = sum [if m == MC munit then x else 0 | (m,x) <- ts]+    comult = linear cm+        where cm m = if m == MC munit then return (T m m) else return (T m (MC munit)) <+> return (T (MC munit) m)+-- Brzezinski and Wisbauer, Corings and Comodules, p5++-- Both of the above can be used to define coalgebra structure on polynomial algebras+-- by using the definitions above on the generators (ie the indeterminates) and then extending multiplicatively+-- They are then guaranteed to be algebra morphisms?+++-- MODULES AND COMODULES++class Algebra k a => Module k a m where+    action :: Vect k (Tensor a m) -> Vect k m++r *. m = action (r `te` m)++class Coalgebra k c => Comodule k c n where+    coaction :: Vect k n -> Vect k (Tensor c n)+++instance Algebra k a => Module k a a where+    action = mult++instance Coalgebra k c => Comodule k c c where+    coaction = comult++-- module and comodule instances for tensor products++-- Kassel p57-8++instance (Num k, Ord a, Ord u, Ord v, Algebra k a, Module k a u, Module k a v)+         => Module k (Tensor a a) (Tensor u v) where+    -- action x = nf $ x >>= action'+    action = linear action'+        where action' (T (T a a') (T u v)) = (action $ return $ T a u) `te` (action $ return $ T a' v)++instance (Num k, Ord a, Ord u, Ord v, Bialgebra k a, Module k a u, Module k a v)+         => Module k a (Tensor u v) where+    -- action x = nf $ x >>= action'+    action = linear action'+        where action' (T a (T u v)) = action $ (comult $ return a) `te` (return $ T u v)+-- !! Overlapping instances+-- If a == Tensor b b, then we have overlapping instance with the previous definition+-- On the other hand, if a == Tensor u v, then we have overlapping instance with the earlier instance++-- Kassel p63+instance (Num k, Ord a, Ord m, Ord n, Bialgebra k a, Comodule k a m, Comodule k a n)+         => Comodule k a (Tensor m n) where+    coaction = (mult `tf` id) . twistm . (coaction `tf` coaction)+        where twistm x = nf $ fmap (\(T (T h m) (T h' n)) -> T (T h h') (T m n)) x

@@ -0,0 +1,54 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}++-- |A module defining the tensor algebra, symmetric algebra, and exterior (or alternating) algebra+module Math.Algebras.TensorAlgebra where++import qualified Data.List as L++import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Structures++import Math.Algebra.Field.Base+++data TensorAlgebra a = TA Int [a] deriving (Eq,Ord,Show)++instance Mon (TensorAlgebra a) where+    munit = TA 0 []+    mmult (TA i xs) (TA j ys) = TA (i+j) (xs++ys)++instance (Num k, Ord a) => Algebra k (TensorAlgebra a) where+    unit 0 = zero -- V []+    unit x = V [(munit,x)]+    mult = nf . fmap (\(T a b) -> a `mmult` b)+++data SymmetricAlgebra a = Sym Int [a] deriving (Eq,Ord,Show)++instance Ord a => Mon (SymmetricAlgebra a) where+    munit = Sym 0 []+    mmult (Sym i xs) (Sym j ys) = Sym (i+j) $ L.sort (xs++ys)++instance (Num k, Ord a) => Algebra k (SymmetricAlgebra a) where+    unit 0 = zero -- V []+    unit x = V [(munit,x)]+    mult = nf . fmap (\(T a b) -> a `mmult` b)+++data ExteriorAlgebra a = Ext Int [a] deriving (Eq,Ord,Show)++instance (Num k, Ord a) => Algebra k (ExteriorAlgebra a) where+    unit 0 = zero -- V []+    unit x = V [(Ext 0 [],x)]+    mult xy = nf $ xy >>= (\(T (Ext i xs) (Ext j ys)) -> signedMerge 1 (0,[]) (i,xs) (j,ys))+        where signedMerge s (k,zs) (i,x:xs) (j,y:ys) =+                  case compare x y of+                  EQ -> zero+                  LT -> signedMerge s (k+1,x:zs) (i-1,xs) (j,y:ys)+                  GT -> let s' = if even i then s else -s -- we had to commute y past x:xs, with i sign changes+                        in signedMerge s' (k+1,y:zs) (i,x:xs) (j-1,ys)+              signedMerge s (k,zs) (i,xs) (0,[]) = s *> (return $ Ext (k+i) $ reverse zs ++ xs)+              signedMerge s (k,zs) (0,[]) (j,ys) = s *> (return $ Ext (k+j) $ reverse zs ++ ys)

@@ -0,0 +1,47 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE NoMonomorphismRestriction #-}++-- |A module defining tensor products of vector spaces+module Math.Algebras.TensorProduct where++import Math.Algebras.VectorSpace+++data Tensor a b = T a b deriving (Eq, Ord, Show)+-- or T !a !b, forcing strictness, but not proven to be better+++-- |Tensor product of two elements+te :: Num k => Vect k a -> Vect k b -> Vect k (Tensor a b)+te (V us) (V vs) = V [(T ei ej, xi*xj) | (ei,xi) <- us, (ej,xj) <- vs]+-- preserves order - that is, if the inputs are correctly ordered, so is the output++-- Implicit assumption - f and g are linear+-- |Tensor product of two (linear) functions+tf :: (Num k, Ord a', Ord b') => (Vect k a -> Vect k a') -> (Vect k b -> Vect k b')+   -> Vect k (Tensor a b) -> Vect k (Tensor a' b')+tf f g (V ts) = sum [te (f $ V [(a, 1)]) (g $ V [(b, x)]) | (T a b, x) <- ts]+    where sum = foldl add zero -- (V [])+++-- tensor isomorphisms++-- in fact, this definition works for any Functor f, not just (Vect k)+assocL :: Vect k (Tensor u (Tensor v w)) -> Vect k (Tensor (Tensor u v) w)+assocL = fmap (\(T a (T b c)) -> T (T a b) c)++assocR :: Vect k (Tensor (Tensor u v) w) -> Vect k (Tensor u (Tensor v w))+assocR = fmap (\(T (T a b) c) -> T a (T b c))++inUnitL = fmap (\a -> T () a)++inUnitR = fmap (\a -> T a ())++outUnitL = fmap (\(T () a) -> a)++outUnitR = fmap (\(T a ()) -> a)++twist v = nf $ fmap (\(T a b) -> T b a) v+-- note the nf call, as f is not order-preserving+

@@ -0,0 +1,130 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE NoMonomorphismRestriction #-}++-- |A module defining the type and operations of free k-vector spaces over a basis b (for a field k)+module Math.Algebras.VectorSpace where++import qualified Data.List as L+import qualified Data.Set as S -- only needed for toSet++-- toSet = S.toList . S.fromList++infixr 7 *>+infixl 7 <*+infixl 6 <+>+++-- |Given a field type k (ie a Fractional instance), Vect k b is the type of the free k-vector space over the basis type b.+-- Elements of Vect k b consist of k-linear combinations of elements of b.+data Vect k b = V [(b,k)] deriving (Eq,Ord)++instance (Num k, Show b) => Show (Vect k b) where+    show (V []) = "0"+    show (V ts) = concatWithPlus $ map showTerm ts+        where showTerm (b,x) | show b == "1" = show x+                             | show x == "1" = show b+                             | show x == "-1" = "-" ++ show b+                             | otherwise = (if isAtomic (show x) then show x else "(" ++ show x ++ ")")+                                           -- (let (c:cs) = show x in+                                           -- if any (`elem` "+-") cs then "(" ++ show x ++ ")" else show x)+                                           ++ show b+              concatWithPlus (t1:t2:ts) = if head t2 == '-'+                                          then t1 ++ concatWithPlus (t2:ts)+                                          else t1 ++ '+' : concatWithPlus (t2:ts)+              concatWithPlus [t] = t+              isAtomic (c:cs) = isAtomic' cs+              isAtomic' ('^':'-':cs) = isAtomic' cs+              isAtomic' ('+':cs) = False+              isAtomic' ('-':cs) = False+              isAtomic' (c:cs) = isAtomic' cs+              isAtomic' [] = True++-- |The zero vector+zero :: Vect k b+zero = V []++-- |Addition of vectors+add :: (Ord b, Num k) => Vect k b -> Vect k b -> Vect k b+add (V ts) (V us) = V $ addmerge ts us++-- |Addition of vectors (same as add)+(<+>) :: (Ord b, Num k) => Vect k b -> Vect k b -> Vect k b+(<+>) = add++addmerge ((a,x):ts) ((b,y):us) =+    case compare a b of+    LT -> (a,x) : addmerge ts ((b,y):us)+    EQ -> if x+y == 0 then addmerge ts us else (a,x+y) : addmerge ts us+    GT -> (b,y) : addmerge ((a,x):ts) us+addmerge ts [] = ts+addmerge [] us = us++-- |Negation of vector+neg :: (Num k) => Vect k b -> Vect k b+neg (V ts) = V $ map (\(b,x) -> (b,-x)) ts++-- |Scalar multiplication (on the left)+smultL :: (Num k) => k -> Vect k b -> Vect k b+smultL 0 _ = zero -- V []+smultL k (V ts) = V [(ei,k*xi) | (ei,xi) <- ts]++-- |Same as smultL. Mnemonic is "multiply through (from the left)"+(*>) :: (Num k) => k -> Vect k b -> Vect k b+(*>) = smultL++-- |Scalar multiplication on the right+smultR :: (Num k) => Vect k b -> k -> Vect k b+smultR _ 0 = zero -- V []+smultR (V ts) k = V [(ei,xi*k) | (ei,xi) <- ts]++-- |Same as smultR. Mnemonic is "multiply through (from the right)"+(<*) :: (Num k) => Vect k b -> k -> Vect k b+(<*) = smultR++-- same as return+-- injection of basis elt into vector space+-- inject b = V [(b,1)]++-- same as fmap+-- liftFromBasis f (V ts) = V [(f b, x) | (b, x) <- ts]+-- if f is not order-preserving, then you need to call nf afterwards++-- |Convert an element of Vect k b into normal form. Normal form consists in having the basis elements in ascending order,+-- with no duplicates, and all coefficients non-zero+nf :: (Ord b, Num k) => Vect k b -> Vect k b+nf (V ts) = V $ nf' $ L.sortBy compareFst ts where+    nf' ((b1,x1):(b2,x2):ts) =+        case compare b1 b2 of+        LT -> if x1 == 0 then nf' ((b2,x2):ts) else (b1,x1) : nf' ((b2,x2):ts)+        EQ -> if x1+x2 == 0 then nf' ts else nf' ((b1,x1+x2):ts)+        GT -> error "nf': not pre-sorted"+    nf' [(b,x)] = if x == 0 then [] else [(b,x)]+    nf' [] = []+    compareFst (b1,x1) (b2,x2) = compare b1 b2+    -- compareFst = curry ( uncurry compare . (fst *** fst) )+++-- lift a function on the basis to a function on the vector space+instance Functor (Vect k) where+    fmap f (V ts) = V [(f b, x) | (b,x) <- ts]+-- Note that if f is not order-preserving, then we need to call "nf" afterwards++instance Num k => Monad (Vect k) where+    return a = V [(a,1)]+    V ts >>= f = V $ concat [ [(b,y*x) | let V us = f a, (b,y) <- us] | (a,x) <- ts]+    -- Note that as we can't assume Ord a in the Monad instance, we need to call "nf" afterwards++linear :: (Ord b, Num k) => (a -> Vect k b) -> Vect k a -> Vect k b+linear f v = nf $ v >>= f++newtype EBasis = E Int deriving (Eq,Ord)++instance Show EBasis where show (E i) = "e" ++ show i++e i = return (E i)+e1 = e 1+e2 = e 2+e3 = e 3++-- dual (E i) = E (-i)

@@ -1,4 +1,4 @@--- Copyright (c) David Amos, 2008. All rights reserved.
+-- Copyright (c) 2008, David Amos. All rights reserved.
 
 module Math.Combinatorics.Design where
 
@@ -12,7 +12,7 @@ import Math.Algebra.Field.Extension
 import Math.Algebra.Group.PermutationGroup hiding (elts, order, isMember)
 import Math.Algebra.Group.SchreierSims as SS
-import Math.Combinatorics.Graph as G hiding (to1n)
+import Math.Combinatorics.Graph as G hiding (to1n, incidenceMatrix)
 import Math.Combinatorics.GraphAuts (refine', isSingleton, graphAuts, incidenceAuts) -- , removeGens)
 import Math.Combinatorics.FiniteGeometry
 
@@ -89,9 +89,10 @@ isSquare d@(D xs bs) = is2Design d && length xs == length bs
 
 
--- incidence matrix of a design
--- (rows and columns indexed by blocks and points respectively)
--- (this follows Cameron & van Lint, though elsewhere in the literature it is sometimes the other way round)
+-- (We follow Cameron & van Lint.)
+-- |The incidence matrix of a design, with rows indexed by blocks and columns by points.
+-- (Note that in the literature, the opposite convention is sometimes used instead.)
+incidenceMatrix :: (Eq t) => Design t -> [[Int]]
 incidenceMatrix (D xs bs) = [ [if x `elem` b then 1 else 0 | x <- xs] | b <- bs]
 
 
@@ -109,13 +110,15 @@     vs = vertices graph
     es = edges graph
 
--- the affine plane AG(2,Fq) - a 2-(q^2,q,1) design
+-- |The affine plane AG(2,Fq), a 2-(q^2,q,1) design
+ag2 :: (FiniteField k, Ord k) => [k] -> Design [k]
 ag2 fq = design (points, lines) where
     points = ptsAG 2 fq
     lines = map line $ tail $ ptsPG 2 fq
     line [a,b,c] = [ [x,y] | [x,y] <- points, a*x+b*y+c==0 ]
 
--- the projective plane PG(2,Fq) - a square 2-(q^2+q+1,q+1,1) design
+-- |The projective plane PG(2,Fq), a square 2-(q^2+q+1,q+1,1) design
+pg2 :: (FiniteField k, Ord k) => [k] -> Design [k]
 pg2 fq = design (points, lines) where
     points = ptsPG 2 fq
     lines = L.sort $ map line points
@@ -168,18 +171,22 @@ -- NEW DESIGNS FROM OLD
 
 -- Dual of a design. Cameron & van Lint p11
+-- |The dual of a design
+dual :: (Ord t) => Design t -> Design [t]
 dual (D xs bs) = design (bs, map beta xs) where
     beta x = filter (x `elem`) bs
 
 -- Derived design relative to a point. Cameron & van Lint p11
 -- Derived design of a t-(v,k,lambda) is a t-1-(v-1,k-1,lambda) design.
+derivedDesign :: (Ord t) => Design t -> t -> Design t
 derivedDesign (D xs bs) p = design (xs L.\\ [p], [b L.\\ [p] | b <- bs, p `elem` b])
 
 -- Residual design relative to a point. Cameron & van Lint p13
 -- Point-residual of a t-(v,k,lambda) is a t-1-(v-1,k,mu).
+pointResidual :: (Ord t) => Design t -> t -> Design t
 pointResidual (D xs bs) p = design (xs L.\\ [p], [b | b <- bs, p `notElem` b])
 
--- Complementary design relative to a point. Cameron & van Lint p13
+-- Complementary design. Cameron & van Lint p13
 -- Complement of a t-(v,k,lambda) is a t-(v,v-k,mu).
 complementaryDesign (D xs bs) = design (xs, [xs L.\\ b | b <- bs])
 
@@ -187,6 +194,7 @@ -- This is only a design if (xs,bs) is a square design
 -- It may have repeated blocks - but if so, residuals of the complement will not
 -- Block-residual of a 2-(v,k,lambda) is a 2-(v-k,k-lambda,lambda).
+blockResidual :: (Ord t) => Design t -> [t] -> Design t
 blockResidual d@(D xs bs) b | isSquare d = design (xs L.\\ b, [b' L.\\ b | b' <- bs, b' /= b])
 
 
@@ -195,11 +203,15 @@ isDesignAut (D xs bs) g | supp g `isSubset` xs = all (`S.member` bs') [b -^ g | b <- bs]
     where bs' = S.fromList bs
 
+-- |The incidence graph of a design
+incidenceGraph :: (Ord a) => Design a -> Graph (Either a [a])
 incidenceGraph (D xs bs) = G vs es where -- graph (vs,es) where
     vs = L.sort $ map Left xs ++ map Right bs
     es = L.sort [ [Left x, Right b] | x <- xs, b <- bs, x `elem` b ]
 
 
+-- |Find a strong generating set for the automorphism group of a design
+designAuts :: (Ord t) => Design t -> [Permutation t]
 designAuts d = incidenceAuts $ incidenceGraph d
 
 -- We find design auts by finding graph auts of the incidence graph of the design
@@ -229,16 +241,21 @@ deltaM24 = p [[-1],[0],[1,18,4,2,6],[3],[5,21,20,10,7],[8,16,13,9,12],[11,19,22,14,17],[15]]
 -- this is t -> t^3 / 9 (for t a quadratic residue), t -> 9 t^3 (t a non-residue)
 
+-- |Generators for the Mathieu group M24, a finite simple group of order 244823040
+m24 :: [Permutation Integer]
 m24 = [alphaL2_23, betaL2_23, gammaL2_23, deltaM24]
 
+-- |A strong generating set for the Mathieu group M24, a finite simple group of order 244823040
+m24sgs :: [Permutation Integer]
 m24sgs = sgs m24
--- order 244823040
 
+-- |A strong generating set for the Mathieu group M23, a finite simple group of order 10200960
+m23sgs :: [Permutation Integer]
 m23sgs = filter (\g -> (-1).^g == -1) m24sgs
--- order 10200960
 
+-- |A strong generating set for the Mathieu group M22, a finite simple group of order 443520
+m22sgs :: [Permutation Integer]
 m22sgs = filter (\g -> 0.^g == 0) m23sgs
--- order 443520
 
 -- sgs uses the base implied by the Ord instance, which will be [-1,0,..]
 
@@ -248,16 +265,21 @@ octad = [0,1,2,3,4,7,10,12]
 -- Conway&Sloane p276 - this is a weight 8 codeword from Golay code G24
 
-
+-- |The Steiner system S(5,8,24), with 759 blocks, whose automorphism group is M24
+s_5_8_24 :: Design Integer
 s_5_8_24 = design ([-1..22], octad -^^ l2_23)
 -- S(5,8,24) constructed as the image of a single octad under the action of PSL(2,23)
 -- 759 blocks ( (24 `choose` 5) `div` (8 `choose` 5) )
 -- Automorphism group is M24
 
+-- |The Steiner system S(4,7,23), with 253 blocks, whose automorphism group is M23
+s_4_7_23 :: Design Integer
 s_4_7_23 = derivedDesign s_5_8_24 (-1)
 -- 253 blocks ( (23 `choose` 4) `div` (7 `choose` 4) )
 -- Automorphism group is M23
 
+-- |The Steiner system S(3,6,22), with 77 blocks, whose automorphism group is M22
+s_3_6_22 :: Design Integer
 s_3_6_22 = derivedDesign s_4_7_23 0
 -- 77 blocks
 -- Automorphism group is M22
@@ -293,20 +315,30 @@ -- the squares (quadratic residues) in F11
 -- http://en.wikipedia.org/wiki/Steiner_system
 
+-- |The Steiner system S(5,6,12), with 132 blocks, whose automorphism group is M12
+s_5_6_12 :: Design Integer
 s_5_6_12 = design ([-1..10], hexad -^^ l2_11)
 -- S(5,6,12) constructed as the image of a single hexad under the action of PSL(2,11)
 -- 132 blocks ( (12 `choose` 5) `div` (6 `choose` 5) )
 -- Automorphism group is M12
 
+-- |The Steiner system S(4,5,11), with 66 blocks, whose automorphism group is M11
+s_4_5_11 :: Design Integer
 s_4_5_11 = derivedDesign s_5_6_12 (-1)
 -- 66 blocks
 -- Automorphism group is M11
 
+-- |Generators for the Mathieu group M12, a finite simple group of order 95040
+m12 :: [Permutation Integer]
 m12 = [alphaL2_11, betaL2_11, gammaL2_11, deltaM12]
 
+-- |A strong generating set for the Mathieu group M12, a finite simple group of order 95040
+m12sgs :: [Permutation Integer]
 m12sgs = sgs m12
 -- order 95040
 
+-- |A strong generating set for the Mathieu group M11, a finite simple group of order 7920
+m11sgs :: [Permutation Integer]
 m11sgs = filter (\g -> (-1).^g == -1) m12sgs
 -- order 7920
 

@@ -1,5 +1,6 @@--- Copyright (c) David Amos, 2009. All rights reserved.
+-- Copyright (c) 2009, David Amos. All rights reserved.
 
+-- |A module defining a type for hypergraphs.
 module Math.Combinatorics.Hypergraph where
 
 import qualified Data.List as L
@@ -24,7 +25,8 @@ -- this still doesn't guarantee that all bs are subset of xs
 
 
--- uniform hypergraph - all blocks are same size
+-- |Is this hypergraph uniform - meaning that all blocks are of the same size
+isUniform :: (Ord a) => Hypergraph a -> Bool
 isUniform h@(H xs bs) = isSetSystem xs bs && same (map length bs)
 
 same (x:xs) = all (==x) xs
@@ -46,14 +48,15 @@ 
 -- INCIDENCE GRAPH
 
-data Incidence a = P a | B [a] deriving (Eq, Ord, Show)
+-- data Incidence a = P a | B [a] deriving (Eq, Ord, Show)
 
 -- compare Design, where we just use Left, Right
 
 -- Also called the Levi graph
+incidenceGraph :: (Ord a) => Hypergraph a -> Graph (Either a [a])
 incidenceGraph (H xs bs) = G vs es where
-    vs = map P xs ++ map B bs
-    es = L.sort [ [P x, B b] | b <- bs, x <- b]
+    vs = map Left xs ++ map Right bs
+    es = L.sort [ [Left x, Right b] | b <- bs, x <- b]
 
 
 -- INCIDENCE MATRIX
@@ -81,19 +84,24 @@ -- We can represent various incidence structures as hypergraphs,
 -- by identifying the lines with the sets of points that they contain
 
+isPartialLinearSpace :: (Ord a) => Hypergraph a -> Bool
 isPartialLinearSpace h@(H ps ls) =
     isSetSystem ps ls &&
     all ( (<=1) . length ) [filter (pair `isSubset`) ls | pair <- combinationsOf 2 ps]
     -- any two points are incident with at most one line
 
 -- Godsil & Royle, p79
+-- |Is this hypergraph a projective plane - meaning that any two lines meet in a unique point,
+-- and any two points lie on a unique line
+isProjectivePlane :: (Ord a) => Hypergraph a -> Bool
 isProjectivePlane h@(H ps ls) =
     isSetSystem ps ls &&
     all ( (==1) . length) [intersect l1 l2 | [l1,l2] <- combinationsOf 2 ls] && -- any two lines meet in a unique point
     all ( (==1) . length) [ filter ([p1,p2] `isSubset`) ls | [p1,p2] <- combinationsOf 2 ps] -- any two points lie in a unique line
 
--- a projective plane with a triangle
--- this is a weak non-degeneracy condition, which eliminates all points on the same line, or all lines through the same point
+-- |Is this hypergraph a projective plane with a triangle.
+-- This is a weak non-degeneracy condition, which eliminates all points on the same line, or all lines through the same point.
+isProjectivePlaneTri :: (Ord a) => Hypergraph a -> Bool
 isProjectivePlaneTri h@(H ps ls) =
     isProjectivePlane h && any triangle (combinationsOf 3 ps)
     where triangle t@[p1,p2,p3] =
@@ -101,8 +109,9 @@                    (not . null) [l | l <- ls, [p1,p3] `isSubset` l, p2 `notElem` l] &&
                    (not . null) [l | l <- ls, [p2,p3] `isSubset` l, p1 `notElem` l] 
 
--- a projective plane with a quadrangle
--- this is a stronger non-degeneracy condition
+-- |Is this hypergraph a projective plane with a quadrangle.
+-- This is a stronger non-degeneracy condition.
+isProjectivePlaneQuad :: (Ord a) => Hypergraph a -> Bool
 isProjectivePlaneQuad h@(H ps ls) =
     isProjectivePlane h && any quadrangle (combinationsOf 4 ps)
     where quadrangle q = all (not . collinear) (combinationsOf 3 q) -- no three points collinear
@@ -116,6 +125,7 @@ -- GENERALIZED QUADRANGLES
 
 -- Godsil & Royle p81
+isGeneralizedQuadrangle :: (Ord a) => Hypergraph a -> Bool
 isGeneralizedQuadrangle h@(H ps ls) =
     isPartialLinearSpace h &&
     all (\(l,p) -> unique [p' | p' <- l, collinear (pair p p')]) [(l,p) | l <- ls, p <- ps, p `notElem` l] &&
@@ -143,29 +153,36 @@ -- CONFIGURATIONS
 
 -- http://en.wikipedia.org/wiki/Projective_configuration
+-- |Is this hypergraph a (projective) configuration.
+isConfiguration :: (Ord a) => Hypergraph a -> Bool
 isConfiguration h@(H ps ls) =
     isUniform h && -- a set system, with each line incident with the same number of points
     same [length (filter (p `elem`) ls) | p <- ps] -- each point is incident with the same number of lines
 
 
-
-
+fanoPlane :: Hypergraph Integer
 fanoPlane = toHypergraph [1..7] [[1,2,4],[2,3,5],[3,4,6],[4,5,7],[5,6,1],[6,7,2],[7,1,3]]
 
+-- |The Heawood graph is the incidence graph of the Fano plane
+heawoodGraph :: Graph (Either Integer [Integer])
 heawoodGraph = incidenceGraph fanoPlane
 
 
+desarguesConfiguration :: Hypergraph [Integer]
 desarguesConfiguration = H xs bs where
     xs = combinationsOf 2 [1..5]
     bs = [ [x | x <- xs, x `isSubset` b] | b <- combinationsOf 3 [1..5] ]
 
+desarguesGraph :: Graph (Either [Integer] [[Integer]])
 desarguesGraph = incidenceGraph desarguesConfiguration
 
 
+pappusConfiguration :: Hypergraph Integer
 pappusConfiguration = H xs bs where
     xs = [1..9]
     bs = L.sort [ [1,2,3], [4,5,6], [7,8,9], [1,5,9], [1,6,8], [2,4,9], [3,4,8], [2,6,7], [3,5,7] ]
 
+pappusGraph :: Graph (Either Integer [Integer])
 pappusGraph = incidenceGraph pappusConfiguration
 
 
@@ -187,6 +204,7 @@     -- (so these are the projective planes over 7 points)
 -}
 -- Godsil & Royle p69
+coxeterGraph :: Graph [Integer]
 coxeterGraph = G vs es where
     g = p [[1..7]]
     vs = L.sort $ concatMap (orbitB [g]) [[1,2,4],[3,5,7],[3,6,7],[5,6,7]]
@@ -204,7 +222,8 @@                            d2 <- duads, d2 > d1, disjoint d1 d2,
                            d3 <- duads, d3 > d2, disjoint d1 d3, disjoint d2 d3 ]
 
--- Tutte 8-cage
+-- |The Tutte-Coxeter graph, also called the Tutte 8-cage
+tutteCoxeterGraph :: Graph (Either [Integer] [[Integer]])
 tutteCoxeterGraph = incidenceGraph $ H duads synthemes
 
 

@@ -18,6 +18,7 @@  -- LATIN SQUARES +findLatinSqs :: (Eq a) => [a] -> [[[a]]] findLatinSqs xs = findLatinSqs' 1 [xs] where     n = length xs     findLatinSqs' i rows@@ -28,6 +29,7 @@                                     | r <- rs, r `notElem` col]     findRows [] ls _ = [reverse ls] +isLatinSq :: (Ord a) => [[a]] -> Bool isLatinSq rows = all isOneOfEach rows && all isOneOfEach cols where     cols = L.transpose rows @@ -54,6 +56,8 @@  -- ORTHOGONAL AND MUTUALLY ORTHOGONAL LATINS SQUARES +-- |Are the two latin squares orthogonal?+isOrthogonal :: (Ord a, Ord b) => [[a]] -> [[b]] -> Bool isOrthogonal greeks latins = isOneOfEach pairs     where pairs = zip (concat greeks) (concat latins) @@ -65,10 +69,13 @@         else findMOLS' i ls rs     findMOLS' _ _ [] = [] +-- |Are the latin squares mutually orthogonal (ie each pair is orthogonal)?+isMOLS :: (Ord a) => [[[a]]] -> Bool isMOLS (greek:latins) = all (isOrthogonal greek) latins && isMOLS latins isMOLS [] = True --- MOLS from a projective plane+-- |MOLS from a projective plane+fromProjectivePlane :: (Ord k, Num k) => Design [k] -> [[[Int]]] fromProjectivePlane (D xs bs) = map toLS parallelClasses where     k = [x | [0,1,x] <- xs] -- the field we're working over     n = length k            -- the order of the projective plane

@@ -1,5 +1,6 @@--- Copyright (c) David Amos, 2008. All rights reserved.
+-- Copyright (c) 2008, David Amos. All rights reserved.
 
+-- |A module defining various strongly regular graphs, including the Clebsch, Hoffman-Singleton, Higman-Sims, and McLaughlin graphs
 module Math.Combinatorics.StronglyRegularGraph where
 
 import qualified Data.List as L

@@ -1,4 +1,4 @@--- Copyright (c) David Amos, 2008. All rights reserved.
+-- Copyright (c) 2008, David Amos. All rights reserved.
 
 module Math.Projects.ChevalleyGroup.Classical where
 
@@ -20,15 +20,18 @@ 
 -- LINEAR GROUPS
 
--- SL(n,Fq) is generated by elementary transvections
--- sl :: FiniteField k => Int -> k -> [[[k]]]
+-- |The special linear group SL(n,Fq), generated by elementary transvections, returned as matrices
+sl :: FiniteField k => Int -> [k] -> [[[k]]]
 sl n fq = [elemTransvection n (r,c) l | r <- [1..n], c <- [1..n], r /= c, l <- fq']
     where fq' = basisFq undefined -- tail fq
     -- Carter p68 - x_r(t1) x_r(t2) == x_r(t1+t2) - this is true in general, not just in this case
 
 elemTransvection n (r,c) l = fMatrix n (\i j -> if i == j then 1 else if (i,j) == (r,c) then l else 0)
 
--- PSL(n,Fq) == A(n,Fq) == SL(n,Fq)/Z
+-- |The projective special linear group PSL(n,Fq) == A(n,Fq) == SL(n,Fq)/Z,
+-- returned as permutations of the points of PG(n-1,Fq).
+-- This is a finite simple group provided n>2 or q>3.
+l :: (FiniteField k, Ord k) => Int -> [k] -> [Permutation [k]]
 l n fq = [fromPairs $ [(p, pnf (p <*>> m)) | p <- ps] | m <- sl n fq]
     where ps = ptsPG (n-1) fq
 
@@ -39,6 +42,8 @@ -- SYMPLECTIC GROUPS
 -- Carter p186 and 181-3
 
+-- |The symplectic group Sp(2n,Fq), returned as matrices
+sp2 :: FiniteField k => Int -> [k] -> [[[k]]]
 sp2 n fq =
     [_I <<+>> t *>> (e i j <<->> e (-j) (-i)) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
     [_I <<->> t *>> (e (-i) (-j) <<->> e j i) | i <- [1..n], j <- [i+1..n], t <- fq' ] ++
@@ -52,7 +57,10 @@         e i j = e' (if i > 0 then i else n-i) (if j > 0 then j else n-j)
         e' i j = fMatrix (2*n) (\k l -> if (k,l) == (i,j) then 1 else 0)
 
--- PSp2n(Fq) == Cn(Fq) == Sp2n(Fq)/Z
+-- |The projective symplectic group PSp(2n,Fq) == Cn(Fq) == Sp(2n,Fq)/Z,
+-- returned as permutations of the points of PG(2n-1,Fq).
+-- This is a finite simple group for n>1, except for PSp(4,F2).
+s2 :: (FiniteField k, Ord k) => Int -> [k] -> [Permutation [k]]
 s2 n fq = [fromPairs $ [(p, pnf (p <*>> m)) | p <- ps] | m <- sp2 n fq]
     where ps = ptsPG (2*n-1) fq
 

@@ -1,4 +1,4 @@--- Copyright (c) David Amos, 2008. All rights reserved.
+-- Copyright (c) 2008, David Amos. All rights reserved.
 
 module Math.Projects.ChevalleyGroup.Exceptional where
 
@@ -9,7 +9,8 @@ import Math.Algebra.LinearAlgebra
 
 import Math.Algebra.Group.PermutationGroup hiding (fromList)
-import Math.Algebra.Group.SchreierSims as SS
+-- import Math.Algebra.Group.SchreierSims as SS
+import Math.Algebra.Group.RandomSchreierSims as RSS
 
 import Math.Combinatorics.FiniteGeometry (ptsAG)
 -- import ClassicalChevalleyGroup (ptsAG)
@@ -24,13 +25,14 @@ 
 newtype Octonion k = O [(Int,k)] deriving (Eq, Ord)
 
-i0 = O [(0,1)] :: Octonion Q
-i1 = O [(1,1)] :: Octonion Q
-i2 = O [(2,1)] :: Octonion Q
-i3 = O [(3,1)] :: Octonion Q
-i4 = O [(4,1)] :: Octonion Q
-i5 = O [(5,1)] :: Octonion Q
-i6 = O [(6,1)] :: Octonion Q
+i0, i1, i2, i3, i4, i5, i6 :: Octonion Q
+i0 = O [(0,1)]
+i1 = O [(1,1)]
+i2 = O [(2,1)]
+i3 = O [(3,1)]
+i4 = O [(4,1)]
+i5 = O [(5,1)]
+i6 = O [(6,1)]
 
 fromList as = O $ filter ((/=0) . snd) $ zip [-1..6] as
 
@@ -155,8 +157,11 @@ beta3' = fromPairs [(x, x %^ beta3) | x <- unitImagOctonions f3]
 gamma3' = fromPairs [(x, x %^ gamma3) | x <- unitImagOctonions f3]
 
+-- |Generators for G2(3), a finite simple group of order 4245696,
+-- as a permutation group on the 702 unit imaginary octonions over F3
+g2_3 :: [Permutation (Octonion F3)]
+g2_3 = [alpha3', beta3', gamma3']
 -- These three together generate a group of order 4245696, which is therefore the whole of G2(3)
--- (But takes nearly 10 minutes to construct the BSGS in interpreter)
 
 -- Unit imaginary octonions form one orbit under the action of G2
 

@@ -0,0 +1,210 @@+-- Copyright (c) David Amos, 2010. All rights reserved.++{-# LANGUAGE TypeFamilies, EmptyDataDecls #-}+++module Math.QuantumAlgebra.OrientedTangle where++import Math.Algebra.Field.Base+import Math.Algebras.LaurentPoly -- hiding (lvar, q, q')++import Math.QuantumAlgebra.TensorCategory++import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Structures++-- import MathExperiments.Algebra.TAlgebra+++-- ORIENTED TANGLE CATEGORY++data Oriented = Plus | Minus deriving (Eq,Ord,Show)++data HorizDir = ToL | ToR deriving (Eq,Ord,Show)++data OrientedTangle++-- In GHCi 6.12.1, we appear to be limited to 8 value constructors within an associated data family+instance Category OrientedTangle where+    data Ob OrientedTangle = OT [Oriented] deriving (Eq,Ord,Show)+    data Ar OrientedTangle = IdT [Oriented]+                           | CapT HorizDir+                           | CupT HorizDir+                           | XPlus | XMinus+                           | SeqT [Ar OrientedTangle]+                           | ParT [Ar OrientedTangle]+                           deriving (Eq,Ord,Show)+    id_ (OT os) = IdT os+    source (IdT os) = OT os+    source (CapT _) = OT []+    source (CupT toR) = OT [Plus,Minus]+    source (CupT toL) = OT [Minus,Plus]+    source XPlus = OT [Plus,Plus]+    source XMinus = OT [Plus,Plus]+    source (ParT as) = OT $ concatMap ((\(OT os) -> os) . source) as+    source (SeqT as) = source (head as)+    target (IdT os) = OT os+    target (CapT toR) = OT [Minus,Plus]+    target (CapT toL) = OT [Plus,Minus]+    target (CupT _) = OT []+    target XPlus = OT [Plus,Plus]+    target XMinus = OT [Plus,Plus]+    target (ParT as) = OT $ concatMap ((\(OT os) -> os) . target) as+    target (SeqT as) = target (last as)+    a >>> b | target a == source b = SeqT [a,b]++instance TensorCategory OrientedTangle where+    tunit = OT []+    tob (OT as) (OT bs) = OT (as++bs)+    tar a b = ParT [a,b]+++++++idV = id+idV' = id++evalV  = \(T (E i) (E j)) -> if i + j == 0 then return () else zero+evalV' = \(T (E i) (E j)) -> if i + j == 0 then return () else zero++coevalV  m = foldl (<+>) zero [e i `te` e (-i) | i <- [1..m] ]+coevalV' m = foldl (<+>) zero [e (-i) `te` e i | i <- [1..m] ]++lambda m = q' ^ m -- q^-m++c m (T (E i) (E j)) = case compare i j of+                      EQ -> (lambda m * q) *> return (T (E i) (E i))+                      LT -> lambda m *> return (T (E j) (E i))+                      GT -> lambda m *> (return (T (E j) (E i)) <+> (q - q') *> return (T (E i) (E j)))++-- inverse of c+c' m (T (E i) (E j)) = case compare i j of+                       EQ -> (1/(lambda m * q)) *> return (T (E i) (E i))+                       LT -> (1/lambda m) *> (return (T (E j) (E i)) <+> (q'-q) *> return (T (E i) (E j)))+                       GT -> (1/lambda m) *> return (T (E j) (E i))++testcc' m v = nf $ v >>= c m >>= c' m++mu m (E i) = (1 / (lambda m * q ^ (2*i-1))) *> return (E i)++mu' m (E i) = (lambda m * q ^ (2*i-1)) *> return (E i)++-- The following are modified from Kassel. We compose diagrams downwards, whereas he composes them upwards.++capRL m = coevalV m++capLR m = do+    T i j <- coevalV' m+    k <- mu' m j+    return (T i k)++cupRL m = evalV++cupLR m (T i j) = do+    k <- mu m i+    evalV' (T k j)    +-- linear evalV' . (linear (mu' m) `tf` idV)++++xplus m = c m++xminus m = c' m++yplus m (T p q) = do+    T r s <- capRL m+    T t u <- xplus m (T q r)+    cupRL m (T p t)+    return (T u s)++yminus m (T p q) = do+    T r s <- capRL m+    T t u <- xminus m (T q r)+    cupRL m (T p t)+    return (T u s)++tplus m (T p q) = do+    T r s <- capLR m+    T t u <- xplus m (T s p)+    cupLR m (T u q)+    return (T r t)++tminus m (T p q) = do+    T r s <- capLR m+    T t u <- xminus m (T s p)+    cupLR m (T u q)+    return (T r t)++zplus m (T p q) = do+    T r u <- capLR m+    T s t <- capLR m+    T v w <- xplus m (T t u)+    cupLR m (T v q)+    cupLR m (T w p)+    return (T r s)++zminus m (T p q) = do+    T r u <- capLR m+    T s t <- capLR m+    T v w <- xminus m (T t u)+    cupLR m (T v q)+    cupLR m (T w p)+    return (T r s)++{-+Then we have for example the following:+> let v = e1 `te` e2 in nf $ v >>= xplus 2 >>= xminus 2+T e1 e2+> let v = e (-1) `te` e2 in nf $ v >>= yplus 2 >>= tminus 2+T e-1 e2+> let v = e (-1) `te` e (-2) in nf $ v >>= zplus 2 >>= zminus 2+T e-1 e-2++-}+++oloop m = nf $ do+    T a b <- capLR m+    cupRL m (T a b)++-- oriented trefoil+otrefoil m = nf $ do+    T p q <- capLR m+    T r s <- capLR m+    T t u <- tminus m (T q r)+    T v w <- zminus m (T p t)+    T x y <- xminus m (T u s)+    cupRL m (T w x)+    cupRL m (T v y)++-- oriented the other way+otrefoil' m = nf $ do+    T p q <- capRL m+    T r s <- capRL m+    T t u <- yminus m (T q r)+    T v w <- xminus m (T p t)+    T x y <- zminus m (T u s)+    cupLR m (T w x)+    cupLR m (T v y)+++{-+-- REPRESENTATIONS OF THE TANGLE CATEGORY IN VECTOR SPACE CATEGORY+-- But we need to convert the above code to use TensorAlgebra first++kauffman :: Ar Tangle -> TangleRep [Oriented] -> TangleRep [Oriented]+kauffman (IdT n) = id -- could be tf of n ids+kauffman CapT = linear cap+kauffman CupT = linear cup+kauffman OverT = linear over+kauffman UnderT = linear under+kauffman (SeqT fs) = foldl (>>>) id $ map kauffman fs+    where g >>> h = h . g+kauffman (ParT [f]) = kauffman f+kauffman (ParT (f:fs)) = tf m (kauffman f) (kauffman (ParT fs))+    where OT m = source f+          tf m f' fs' = linear (\xs -> let (ls,rs) = splitAt m xs in f' (return ls) * fs' (return rs) )+-}

@@ -0,0 +1,236 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE TypeSynonymInstances, FlexibleInstances, MultiParamTypeClasses, NoMonomorphismRestriction #-}++-- |A module defining the quantum plane and its symmetries+module Math.QuantumAlgebra.QuantumPlane where++-- Refs:+-- Kassel, Quantum Groups+-- Street, Quantum Groups++import Math.Algebra.Field.Base hiding (powers)++import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Structures+import Math.Algebras.LaurentPoly+import Math.Algebras.NonCommutative+import qualified Data.List as L++++qvar v = let V [(m,1)] = var v in V [(m,1 :: LaurentPoly Q)]+++a = qvar "a"+b = qvar "b"+c = qvar "c"+d = qvar "d"++detq = a*d-unit q'*b*c+++x = qvar "x"+y = qvar "y"+-- z = qvar "z"++u = qvar "u"+v = qvar "v"+++-- Quantum plane Aq20++aq20 = [y*x-unit q*x*y]+-- Kassel p72, Street p10++newtype Aq20 v = Aq20 (NonComMonomial v) deriving (Eq,Ord)++instance (Eq v, Show v) => Show (Aq20 v) where show (Aq20 m) = show m++instance Monomial Aq20 where+    var v = V [(Aq20 (NCM 1 [v]),1)]+    powers (Aq20 m) = powers m++instance Algebra (LaurentPoly Q) (Aq20 String) where+    unit 0 = zero -- V []+    unit x = V [(munit,x)] where munit = Aq20 (NCM 0 [])+    mult x = x''' where+        x' = mult $ fmap (\(T (Aq20 a) (Aq20 b)) -> T a b) x -- unwrap and multiply+        x'' = x' %% aq20 -- quotient by m2q relations while unwrapped+        x''' = fmap Aq20 x'' -- wrap the monomials up as Aq20 again+++-- Quantum superplane Aq02++aq02 = [u^2, v^2, u*v+unit q*v*u]+-- Street p10++newtype Aq02 v = Aq02 (NonComMonomial v) deriving (Eq,Ord)++instance (Eq v, Show v) => Show (Aq02 v) where show (Aq02 m) = show m++instance Monomial Aq02 where+    var v = V [(Aq02 (NCM 1 [v]),1)]+    powers (Aq02 m) = powers m++instance Algebra (LaurentPoly Q) (Aq02 String) where+    unit 0 = zero -- V []+    unit x = V [(munit,x)] where munit = Aq02 (NCM 0 [])+    mult x = x''' where+        x' = mult $ fmap (\(T (Aq02 a) (Aq02 b)) -> T a b) x -- unwrap and multiply+        x'' = x' %% aq02 -- quotient by m2q relations while unwrapped+        x''' = fmap Aq02 x'' -- wrap the monomials up as Aq02 again+++-- M2q++m2q = [a*b-unit q'*b*a, a*c-unit q'*c*a, c*d-unit q'*d*c, b*d-unit q'*d*b,+       b*c-c*b, a*d-d*a-unit (q'-q)*b*c]+-- Kassel p78, Street p9+-- I think this is already a Groebner basis++newtype M2q v = M2q (NonComMonomial v) deriving (Eq,Ord)++instance (Eq v, Show v) => Show (M2q v) where show (M2q m) = show m++instance Monomial M2q where+    var v = V [(M2q (NCM 1 [v]),1)]+    powers (M2q m) = powers m++instance Algebra (LaurentPoly Q) (M2q String) where+    unit 0 = zero -- V []+    unit x = V [(munit,x)] where munit = M2q (NCM 0 [])+    mult x = x''' where+        x' = mult $ fmap (\(T (M2q a) (M2q b)) -> T a b) x -- unwrap and multiply+        x'' = x' %% m2q -- quotient by m2q relations while unwrapped+        x''' = fmap M2q x'' -- wrap the monomials up as M2q again++-- Kassel p82-3+instance Coalgebra (LaurentPoly Q) (M2q String) where+    counit x = case x `bind` cu of+               V [] -> 0+               V [(M2q (NCM 0 []), c)] -> c+        where cu "a" = 1 :: Vect (LaurentPoly Q) (M2q String)+              cu "b" = 0+              cu "c" = 0+              cu "d" = 1+    comult x = x `bind` cm+        where cm "a" = a `te` a + b `te` c+              cm "b" = a `te` b + b `te` d+              cm "c" = c `te` a + d `te` c+              cm "d" = c `te` b + d `te` d++instance Bialgebra (LaurentPoly Q) (M2q String) where {}++{-+-- The following shows that the M2q relations are *sufficient*+-- for M2q to be symmetries of Aq20 and Aq02++> let x' = a*x+b*y :: Vect (LaurentPoly Q) (NonComMonomial String)+> let y' = c*x+d*y :: Vect (LaurentPoly Q) (NonComMonomial String)+> (y'*x'-unit q*x'*y') %% (m2q ++ aq20 ++ [s*t-t*s | s <- [a,b,c,d], t <- [x,y]])+0++> let u' = a*u+b*v :: Vect (LaurentPoly Q) (NonComMonomial String)+> let v' = c*u+d*v :: Vect (LaurentPoly Q) (NonComMonomial String)+> (u'^2) %% (m2q ++ aq02 ++ [s*t-t*s | s <- [a,b,c,d], t <- [u,v]])+0+> (v'^2) %% (m2q ++ aq02 ++ [s*t-t*s | s <- [a,b,c,d], t <- [u,v]])+0+> (u'*v'+unit q*v'*u') %% (m2q ++ aq02 ++ [s*t-t*s | s <- [a,b,c,d], t <- [u,v]])+0++-- To show that the M2q relations are necessary,+-- set the coefficients of x^2, yx, y^2, and vu == 0 in all of the following+> (y'*x'-unit q*x'*y') %% (aq20 ++ [p*q-q*p | p <- [a,b,c,d], q <- [x,y]])+-qx^2ac+x^2ca-yxad-qyxbc+q^-1yxcb+yxda-qy^2bd+y^2db+> (u'^2) %% (aq02 ++ [p*q-q*p | p <- [a,b,c,d], q <- [u,v]])+-qvuab+vuba+> (v'^2) %% (aq02 ++ [p*q-q*p | p <- [a,b,c,d], q <- [u,v]])+-qvucd+vudc+> (u'*v'+unit q*v'*u') %% (aq02 ++ [p*q-q*p | p <- [a,b,c,d], q <- [u,v]])+-qvuad+vubc-q^2vucb+qvuda++-- yx => -ad-qbc+q^-1cb+da == 0+-- vu => -qad+bc-q^2cb+qda == 0+-- qyx-vu => -q^2bc+cb-bc+q^2cb == 0 => bc == cb+-- Now substitute back into yx++-- We could probably have got gb to do this for us+-}++-- Kassel p85+instance Comodule (LaurentPoly Q) (M2q String) (Aq20 String) where+    coaction xy = xy `bind` ca where+        ca "x" = (a `te` x) + (b `te` y) -- we can use (+) instead of add since Aq20 is an algebra +        ca "y" = (c `te` x) + (d `te` y)+-- coaction (x) = (a b) `te` (x)+--          (y)   (c d)      (y)+++-- SL2q++sl2q = [a*b-unit q'*b*a, a*c-unit q'*c*a, c*d-unit q'*d*c, b*d-unit q'*d*b,+        b*c-c*b, a*d-d*a-unit (q'-q)*b*c,+        -unit q*c*b + d*a - 1] -- det q, but reduced+--        a*d-unit q'*b*c-1] -- det_q+-- We have to hand-reduce detq, or else call gb++++newtype SL2q v = SL2q (NonComMonomial v) deriving (Eq,Ord)++instance (Eq v, Show v) => Show (SL2q v) where show (SL2q m) = show m++instance Monomial SL2q where+    var v = V [(SL2q (NCM 1 [v]),1)]+    powers (SL2q m) = powers m++instance Algebra (LaurentPoly Q) (SL2q String) where+    unit 0 = zero -- V []+    unit x = V [(munit,x)] where munit = SL2q (NCM 0 [])+    mult x = x''' where+        x' = mult $ fmap (\(T (SL2q a) (SL2q b)) -> T a b) x -- unwrap and multiply+        x'' = x' %% sl2q -- quotient by sl2q relations while unwrapped+        x''' = fmap SL2q x'' -- wrap the monomials up as SL2q again++instance Coalgebra (LaurentPoly Q) (SL2q String) where+    counit x = case x `bind` cu of+               V [] -> 0+               V [(SL2q (NCM 0 []), c)] -> c+        where cu "a" = 1 :: Vect (LaurentPoly Q) (SL2q String)+              cu "b" = 0+              cu "c" = 0+              cu "d" = 1+    comult x = x `bind` cm+        where cm "a" = a `te` a + b `te` c+              cm "b" = a `te` b + b `te` d+              cm "c" = c `te` a + d `te` c+              cm "d" = c `te` b + d `te` d++instance Bialgebra (LaurentPoly Q) (SL2q String) where {}++-- Kassel p84+instance HopfAlgebra (LaurentPoly Q) (SL2q String) where+    antipode x = x `bind` antipode'+        where antipode' "a" = d+              antipode' "b" = - unit q * b+              antipode' "c" = - unit q' * c+              antipode' "d" = a+-- in the GL2q case we would need 1/detq factor as well++++-- !! The following probably needs to be rehoused in separate module at some point+-- YANG-BAXTER OPERATOR++-- This is a Yang-Baxter operator, but not the only possible such+-- Street, p93+yb x = nf $ x >>= yb' where+    yb' (T a b) = case compare a b of+                 GT -> return (T b a)+                 LT -> return (T b a) + unit (q-q') * return (T a b)+                 EQ -> unit q * return (T a a)++

@@ -0,0 +1,177 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE MultiParamTypeClasses, TypeFamilies, FlexibleInstances, EmptyDataDecls #-}++-- |A module defining the category of tangles, and representations into the category of vector spaces+-- (specifically, knot invariants).+module Math.QuantumAlgebra.Tangle where++-- import qualified Data.List as L++import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Structures++import Math.Algebra.Field.Base+import Math.Algebras.LaurentPoly++import Math.QuantumAlgebra.TensorCategory+++instance Mon [a] where+    munit = []+    mmult = (++)++-- type TensorAlgebra k a = Vect k [a]++instance (Num k, Ord a) => Algebra k [a] where+    unit 0 = zero -- V []+    unit x = V [(munit,x)]+    mult = nf . fmap (\(T a b) -> a `mmult` b)++-- Could make TensorAlgebra k a into an instance of Category, TensorCategory+    ++-- TANGLE CATEGORY+-- (Unoriented)++data Tangle++instance Category Tangle where+    data Ob Tangle = OT Int deriving (Eq,Ord,Show)+    data Ar Tangle = IdT Int+                   | CapT+                   | CupT+                   | OverT+                   | UnderT+--                   | SeqT (Ar Tangle) (Ar Tangle)+                   | SeqT [Ar Tangle]+--                   | ParT (Ar Tangle) (Ar Tangle)+                   | ParT [Ar Tangle]+                   deriving (Eq,Ord,Show)+    id_ (OT n) = IdT n+    source (IdT n) = OT n+    source CapT = OT 0+    source CupT = OT 2+    source OverT = OT 2+    source UnderT = OT 2+--    source (ParT a b) = OT (sa + sb) where OT sa = source a; OT sb = source b+    source (ParT as) = OT $ sum [sa | a <- as, let OT sa = source a]+--    source (SeqT a b) = source a+    source (SeqT as) = source (head as)+    target (IdT n) = OT n+    target CapT = OT 2+    target CupT = OT 0+    target OverT = OT 2+    target UnderT = OT 2+--    target (ParT a b) = OT (ta + tb) where OT ta = target a; OT tb = target b+    target (ParT as) = OT $ sum [ta | a <- as, let OT ta = target a]+--    target (SeqT a b) = target b+    target (SeqT as) = target (last as)+--    a >>> b | target a == source b = SeqT a b+    a >>> b | target a == source b = SeqT [a,b]++instance TensorCategory Tangle where+    tunit = OT 0+    tob (OT a) (OT b) = OT (a+b)+--    tar a b = ParT a b+    tar a b = ParT [a,b]++++-- KAUFFMAN BRACKET++data Oriented = Plus | Minus deriving (Eq,Ord,Show)++type TangleRep b = Vect (LaurentPoly Q) b+++-- adapted from http://blog.sigfpe.com/2008/10/untangling-with-continued-fractions.html+cap :: [Oriented] -> TangleRep [Oriented]+cap [] = return [Plus, Minus] <+> (-q^2) *> return [Minus, Plus]++cup :: [Oriented] -> TangleRep [Oriented]+cup [Plus, Minus] = (-q'^2) *> return []+cup [Minus, Plus] = return []+cup _ = zero++-- also called xminus+over :: [Oriented] -> TangleRep [Oriented]+over [u, v] = q  *> do {[] <- cup [u, v]; cap []}+          <+> q' *> return [u, v]++{-+-- if you expand "over" into terms, you find that it equals the following,+-- which strongly resembles c' below+over' (T i j) = case compare i j of+                EQ -> q' *> return (T i i)                                       -- ++ -> q' ++, -- -> q' -- +                LT -> q  *> return (T j i)                                       -- +- -> q -++                GT -> q  *> (return (T j i) <+> (q'^2 - q^2) *> return (T i j))  -- -+ -> q +- + (q'-q^3) -++-}+-- also called xplus+under :: [Oriented] -> TangleRep [Oriented]+under [u, v] = q' *> do {[] <- cup [u, v]; cap []}+           <+> q  *> return [u, v]++{-+-- if you expand "under" into terms, you find that it equals the following,+-- which strongly resembles c below+under' (T i j) = case compare i j of+                 EQ -> q  *> return (T i i)                                       -- ++ -> q ++, -- -> q -- +                 LT -> q' *> (return (T j i) <+> (q^2 - q'^2) *> return (T i j))  -- +- -> q' -+ + (q-q^-3) -++                 GT -> q' *> return (T j i)                                       -- -+ -> q' +-+-}+loop = nf $ do {[i, j] <- cap []; cup [i, j]}++trefoil = nf $ do+    [i, j] <- cap []+    [k, l] <- cap []+    [m, n] <- under [j, k]+    [p, q] <- over [i, m]+    [r, s] <- over [n, l]+    cup [p, s]+    cup [q, r]+++-- KAUFFMAN BRACKET AS A REPRESENTATION FROM TANGLE TO VECT++-- But this isn't quite the Kauffman bracket - we still need to divide by (-q^2-q^-2)+kauffman :: Ar Tangle -> TangleRep [Oriented] -> TangleRep [Oriented]+kauffman (IdT n) = id -- could be tf of n ids+kauffman CapT = linear cap+kauffman CupT = linear cup+kauffman OverT = linear over+kauffman UnderT = linear under+kauffman (SeqT fs) = foldl (>>>) id $ map kauffman fs+    where g >>> h = h . g+kauffman (ParT [f]) = kauffman f+kauffman (ParT (f:fs)) = tf m (kauffman f) (kauffman (ParT fs))+    where OT m = source f+          tf m f' fs' = linear (\xs -> let (ls,rs) = splitAt m xs in f' (return ls) * fs' (return rs) )+{-+kauffman (ParT f g) = tf m n (kauffman f) (kauffman g)+    where OT m = source f+          OT n = source g+          tf m n f' g' = linear (\xs -> let (ls,rs) = splitAt m xs in f' (return ls) * g' (return rs) )+-}++-- loopT = SeqT CapT CupT+loopT = SeqT [CapT, CupT]++{-+trefoilT = (ParT CapT CapT) `SeqT` (ParT (IdT 1) (ParT UnderT (IdT 1)))+    `SeqT` (ParT OverT OverT) `SeqT` (ParT (IdT 1) (ParT CupT (IdT 1))) `SeqT` CupT++trefoilT = ParT [CapT, CapT]+    `SeqT` ParT [IdT 1, UnderT, IdT 1]+    `SeqT` ParT [OverT, OverT]+    `SeqT` ParT [IdT 1, CupT, IdT 1]+    `SeqT` CupT+-}+trefoilT = SeqT [+    ParT [CapT, CapT],+    ParT [IdT 1, UnderT, IdT 1],+    ParT [OverT, OverT],+    ParT [IdT 1, CupT, IdT 1],+    CupT]+-- eg kauffman (trefoilT) (return [])

@@ -0,0 +1,128 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE TypeFamilies, EmptyDataDecls #-}++-- |A module defining classes and example instances of categories and tensor categories+module Math.QuantumAlgebra.TensorCategory where++import Math.Algebra.Group.PermutationGroup+++class Category c where+    data Ob c :: *+    data Ar c :: *+    id_ :: Ob c -> Ar c+    source, target :: Ar c -> Ob c+    (>>>) :: Ar c -> Ar c -> Ar c++-- Kassel p282+-- The following is actually definition of a strict tensor category+class Category c => TensorCategory c where+    tunit :: Ob c+    tob :: Ob c -> Ob c -> Ob c -- tensor product of objects+    tar :: Ar c -> Ar c -> Ar c -- tensor product of arrows++class TensorCategory c => StrictTensorCategory c where {}+-- we want to be able to declare some tensor categories as strict++class TensorCategory c => WeakTensorCategory c where+    assoc :: Ob c -> Ob c -> Ob c -- (u `tob` v) `tob` w -> u `tob` (v `tob` w)+    lunit :: Ob c -> Ob c         -- unit `tob` v -> v+    runit :: Ob c -> Ob c         -- v `tob` unit -> v++{-+instance (TensorCategory c, Eq (Ar c), Show (Ar c)) => Num (Ar c) where+    (*) = tar+-}++-- SYMMETRIC GROUPOID++data SymmetricGroupoid++instance Category SymmetricGroupoid where+    data Ob SymmetricGroupoid = OS Int deriving (Eq,Ord,Show)+    data Ar SymmetricGroupoid = AS Int (Permutation Int) deriving (Eq,Ord,Show)+    id_ (OS n) = AS n 1+    source (AS n _) = OS n+    target (AS n _) = OS n+    AS m g >>> AS n h | m == n = AS m (g*h)++instance TensorCategory SymmetricGroupoid where+    tunit = OS 0+    tob (OS m) (OS n) = OS (m+n)+    tar (AS m g) (AS n h) = AS (m+n) (g * h~^k)+        where k = p [[1..m+n]] ^ m+--    tar (AS m g) (AS n h) = AS (m+n) (fromPairs $ toPairs g ++ map (\(x,y)->(x+m,y+m)) (toPairs h))+++-- BRAID CATEGORY++data Braid++instance Category Braid where+    data Ob Braid = OB Int deriving (Eq,Ord,Show)+    data Ar Braid = AB Int [Int] deriving (Eq,Ord,Show)+    id_ (OB n) = AB n []+    source (AB n _) = OB n+    target (AB n _) = OB n+    AB m is >>> AB n js | m == n = AB m (is ++ js)++s n i | 0 < i && i < n = AB n [i]++instance TensorCategory Braid where+    tunit = OB 0+    tob (OB a) (OB b) = OB (a+b)+    tar (AB m is) (AB n js) = AB (m+n) (is ++ map (+m) js)+++++data Cob2+-- works very similar to Tangle category++instance Category Cob2 where+    data Ob Cob2 = O Int deriving (Eq,Ord,Show)+    data Ar Cob2 = Id Int+                 | Unit+                 | Mult+                 | Counit+                 | Comult+                 | Par (Ar Cob2) (Ar Cob2)+                 | Seq (Ar Cob2) (Ar Cob2)+                 deriving (Eq,Ord,Show)+    id_ (O n) = Id n+    source (Id n) = O n+    source Unit = O 0+    source Mult = O 2+    source Counit = O 1+    source Comult = O 1+    source (Par a b) = O (sa + sb) where O sa = source a; O sb = source b+    source (Seq a b) = source a+    target (Id n) = O n+    target Unit = O 1+    target Mult = O 1+    target Counit = O 0+    target Comult = O 2+    target (Par a b) = O (ta + tb) where O ta = target a; O tb = target b+    target (Seq a b) = target b+    a >>> b | target a == source b = Seq a b++instance TensorCategory Cob2 where+    tunit = O 0+    tob (O a) (O b) = O (a+b)+    tar a b = Par a b++-- rewrite a Cob2 so that it is a Seq of Pars+-- (this isn't necessarily going to help us towards a normal form - there may not even be one+rewrite (Par (Seq a1 a2) (Seq b1 b2)) =+    Seq (Par idSourceA b1')+        ( (Seq (Par idSourceA b2')+               (Seq (Par a1' idTargetB)+                    (Par a2' idTargetB) ) ) )+    where idSourceA = id_ (source a1)+          idTargetB = id_ (target b2)+          a1' = rewrite a1+          a2' = rewrite a2+          b1' = rewrite b1+          b2' = rewrite b2+rewrite x = x

@@ -0,0 +1,55 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-}++module Math.Test.TAlgebras.TGroupAlgebra where++import Test.QuickCheck++import Math.Algebra.Group.PermutationGroup+import Math.Test.TPermutationGroup -- for instance Arbitrary (Permutation Int)++import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Structures+import Math.Algebras.GroupAlgebra++import Math.Test.TAlgebras.TStructures++{-+instance Arbitrary (TensorAlgebra Integer) where+    arbitrary = do ts <- arbitrary :: Gen [([Int], Integer)]+                   return $ nf $ V ts+-}+instance Arbitrary (GroupAlgebra Integer) where+    arbitrary = do ts <- arbitrary :: Gen [(Permutation Int, Integer)]+                   return $ nf $ V ts+++{-+prop_Algebra_TensorAlgebra (k,x,y,z) = prop_Algebra (k,x,y,z)+    where types = (k,x,y,z) :: (Integer, TensorAlgebra Integer, TensorAlgebra Integer, TensorAlgebra Integer)++prop_Coalgebra_TensorAlgebra x = prop_Coalgebra x+    where types = x :: TensorAlgebra Integer+-}++prop_Algebra_GroupAlgebra (k,x,y,z) = prop_Algebra (k,x,y,z)+    where types = (k,x,y,z) :: (Integer, GroupAlgebra Integer, GroupAlgebra Integer, GroupAlgebra Integer)++-- have to split the 8-tuple into two 4-tuples to avoid having to write Arbitrary instance+prop_Algebra_Linear_GroupAlgebra ((k,l,m,n),(x,y,z,w)) = prop_Algebra_Linear (k,l,m,n,x,y,z,w)+    where types = (k,l,m,n,x,y,z,w) :: (Integer, Integer, Integer, Integer,+                   GroupAlgebra Integer, GroupAlgebra Integer, GroupAlgebra Integer, GroupAlgebra Integer)++prop_Coalgebra_GroupAlgebra x = prop_Coalgebra x+    where types = x :: GroupAlgebra Integer++prop_Coalgebra_Linear_GroupAlgebra (k,l,x,y) = prop_Coalgebra_Linear (k,l,x,y)+    where types = (k,l,x,y) :: (Integer, Integer, GroupAlgebra Integer, GroupAlgebra Integer)++prop_Bialgebra_GroupAlgebra (k,x,y) = prop_Bialgebra (k,x,y)+    where types = (k,x,y) :: (Integer, GroupAlgebra Integer, GroupAlgebra Integer)++prop_HopfAlgebra_GroupAlgebra x = prop_HopfAlgebra x+    where types = x :: GroupAlgebra Integer

@@ -0,0 +1,39 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-}++module Math.Test.TAlgebras.TQuaternions where++import Test.QuickCheck++import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Quaternions++import Math.Test.TAlgebras.TStructures++instance Arbitrary HBasis where+    arbitrary = elements [One,I,J,K]++instance Arbitrary (Quaternion Integer) where+    arbitrary = do ts <- arbitrary :: Gen [(HBasis, Integer)]+                   return $ nf $ V ts++++prop_Algebra_Quaternion (k,x,y,z) = prop_Algebra (k,x,y,z)+    where types = (k,x,y,z) :: (Integer, Quaternion Integer, Quaternion Integer, Quaternion Integer)++prop_Coalgebra_Quaternion x = prop_Coalgebra x+    where types = x :: Quaternion Integer++-- Fails - the algebra and coalgebra structures I've given are not compatible+prop_Bialgebra_Quaternion (k,x,y) = prop_Bialgebra (k,x,y)+    where types = (k,x,y) :: (Integer, Quaternion Integer, Quaternion Integer)++{-+prop_FrobeniusRelation_Quaternion (x,y) = prop_FrobeniusRelation (x,y)+    where types = (x,y) :: (Quaternion Integer, Quaternion Integer)+-- !! fails, because the counit we have given is not a Frobenius form+-}+

@@ -0,0 +1,211 @@+-- Copyright (c) 2010, David Amos. All rights reserved.++{-# LANGUAGE NoMonomorphismRestriction #-}+++module Math.Test.TAlgebras.TStructures where++-- import Test.QuickCheck+-- don't actually need, as we don't define any Arbitrary instances here++import Control.Arrow ( (>>>) ) -- actually you can get this from Category?+++import Math.Algebras.VectorSpace+import Math.Algebras.TensorProduct+import Math.Algebras.Structures -- what we're testing+-- import MathExperiments.Algebra.MonoidAlgebra+-- import MathExperiments.Algebra.Examples++{-+prop_VectorSpace (k,l,x,y,z) =+    smultL k (smultL l x) == smultL (k*l) x &&+    add x y == add y z &&+    add x (add y z) == add (add x y) z &&+    add x zero == x &&+    add zero x == x+-- !! check definition - have I forgotten anything - yes, additive inverses+-}++prop_Linear f (k,l,x,y) =+    f (add (smultL k x) (smultL l y)) == add (smultL k (f x)) (smultL l (f y))+-- now use this to show algebra and coalgebra ops are linear++-- in this version we supply z of the intended return type of f,+-- so that we can make sure we select the correct instance for f polymorphic in return type+prop_Linear' f (k,l,x,y,z) =+    f (add (smultL k x) (smultL l y)) `add` z == add (smultL k (f x)) (smultL l (f y)) `add` z+++-- prop_Bilinear could be defined in terms of prop_Linear over tensor product+-- if we had a way to convert a bilinear function to a tensor function++prop_Algebra_Linear ::+    (Num k, Ord b, Algebra k b) =>+    (k, k, k, k, Vect k b, Vect k b, Vect k b, Vect k b) -> Bool+prop_Algebra_Linear (k,l,m,n,x,y,z,w) =+--    (unit (k * m + l * n) :: Vect k b) == (add (smultL k (unit m)) (smultL l (unit n)) :: Vect k b) &&+    prop_Linear' unit' (k,l,wrap m, wrap n, x) &&+    prop_Linear mult (k,l, x `te` y, z `te` w)+    where wrap = (\c -> V [((),c)]) :: k -> Trivial k++prop_Coalgebra_Linear (k,l,x,y) =+    prop_Linear counit' (k,l,x,y) &&+    prop_Linear comult (k,l,x,y)+-- now need instances for GroupAlgebra etc+++-- ALGEBRAS++prop_Algebra (k,x,y,z) =+    mult (x `te` mult (y `te` z)) == mult (mult (x `te` y) `te` z)  && -- associativity+    smultL k x == mult (unit k `te` x)                             && -- left unit+    -- mult (k' `te` x) == (mult . (unit' `tf` id)) (k' `te` x)         && -- left unit+    smultR x k == mult (x `te` unit k)                        -- && -- right unit+    -- mult (x `te` k') == (mult . (id `tf` unit')) (x `te` k')            -- right unit+    where k' = V [( (),k)]+-- additionally, unit and mult must be linear++prop_Commutative (x,y) =+    let xy = x `te` y+    in (mult . twist) xy == mult xy+++-- COALGEBRAS++prop_Coalgebra x =+    ((comult `tf` id) . comult) x == (assocL . (id `tf` comult) . comult) x && -- coassociativity+    ((counit' `tf` id) . comult) x == V [((),1)] `te` x                     && -- left counit+    ((id `tf` counit') . comult) x == x `te` V [((),1)]                        -- right counit+-- additionally, counit and comult must be linear++prop_Cocommutative x =+    (twist . comult) x == comult x+++-- MORPHISMS++prop_AlgebraMorphism f (k,l,x,y) =+    prop_Linear f (k,l,x,y) &&+    -- (f . unit) k == unit k &&+    (f . mult) (x `te` y) == (mult . (f `tf` f)) (x `te` y) ++-- in this version we supply z of the intended return type of f,+-- so that we can make sure we select the correct instance for f polymorphic in return type+prop_AlgebraMorphism' f (k,l,x,y,z) =+    prop_Linear f (k,l,x,y) &&+    (f . unit) k + z == unit k + z &&+    (f . mult) (x `te` y) == (mult . (f `tf` f)) (x `te` y) ++prop_CoalgebraMorphism f x =+    -- prop_Linear f (k,l,x,y) &&+    (counit . f) x == counit x &&+    ( (f `tf` f) . comult) x == (comult . f) x+++-- BIALGEBRAS++prop_Bialgebra1 (x,y) =+    let xy = x `te` y in+    (comult . mult) xy ==+    ( (mult `tf` mult) .+      assocL . (id `tf` assocR) .+      (id `tf` (twist `tf` id)) .+      (id `tf` assocL) . assocR .+      (comult `tf` comult) ) xy++prop_Bialgebra2 (k,xy) =+    (comult . unit') k' + xy == ((unit' `tf` unit') . iso) k' + xy+    where iso = fmap (\ () -> T () () ) -- the isomorphism k ~= k tensor k+          k' = unit k :: Trivial Integer -- inject into the trivial algebra+-- the +xy is just to force the other expression to be of the right type++prop_Bialgebra3 (x,y) =+    (counit' . mult) xy == (iso . (counit' `tf` counit')) xy+    where xy = x `te` y+          iso = fmap (\(T () ()) -> ())++prop_Bialgebra4 (k,x) =+    id k == (counit . (\a -> a+x-x) . unit) k+-- so we are using the x just to force the intermediate value to be of the right type++prop_Bialgebra (k,x,y) =+    prop_Bialgebra1 (x,y) &&+    prop_Bialgebra2 (k,x `te` y) &&+    prop_Bialgebra3 (x,y) &&+    prop_Bialgebra4 (k,x)+++prop_HopfAlgebra x =+    (unit . counit) x == (mult . (antipode `tf` id) . comult) x &&+    (unit . counit) x == (mult . (id `tf` antipode) . comult) x++-- Street p87+-- we also require that f be invertible+prop_YangBaxter f (x,y,z) =+    ( (f `tf` id) >>> assocR >>> (id `tf` f) >>> assocL >>> (f `tf` id) >>> assocR ) xyz == +    ( assocR >>> (id `tf` f) >>> assocL >>> (f `tf` id) >>> assocR >>> (id `tf` f) ) xyz+    where xyz = ( (x `te` y) `te` z) ++++-- MODULES AND COMODULES++prop_Module_Linear (k,l,x,y) = prop_Linear action (k,l,x,y)++prop_Module_Assoc (r,s,m) =+    (action . (mult `tf` id)) ((r `te` s) `te` m) == (action . (id `tf` action)) (r `te` (s `te` m))+{-+prop_Module_Unit (k,m) =+    (action . (unit' `tf` id)) k' ==  +-}++++++++++++++frobeniusLeft1 = (id `tf` mult) . assocR . (comult `tf` id)++frobeniusLeft2 x = nf $ x >>= fl+    where fl (T i j) = do+              T k l <- comultM i+              m <- idM j+              p <- idM k+              q <- multM (T l m)+              return (T p q)++frobeniusMiddle1 = comult . mult++frobeniusMiddle2 x = nf $ x >>= fm+    where fm (T i j) = do+              k <- multM (T i j)+              T l m <- comultM k+              return (T l m)++prop_FrobeniusRelation (x,y) =+    let xy = x `te` y+    in frobeniusLeft1 xy == frobeniusMiddle1 xy++-- (inject == return)++multM = mult . return -- inject+comultM = comult . return -- inject+idM = id . return++-- can we do the same with unit, counit?+-- unit takes k as input, so isn't in the monad+-- counit gives k as output - what would we do with it+-- so perhaps we have to use unit' and counit'+++++