diff --git a/HaskellForMaths.cabal b/HaskellForMaths.cabal --- a/HaskellForMaths.cabal +++ b/HaskellForMaths.cabal @@ -1,5 +1,5 @@ Name: HaskellForMaths - Version: 0.3.4 + Version: 0.4.0 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 @@ -12,7 +12,6 @@ Cabal-Version: >=1.2 Extra-source-files: - Math/Test/TCommutativeAlgebra.hs, Math/Test/TDesign.hs, Math/Test/TField.hs, Math/Test/TFiniteGeometry.hs, @@ -34,13 +33,12 @@ Math/Test/TCombinatorics/TMatroid.hs Math/Test/TCommutativeAlgebra/TGroebnerBasis.hs Math/Test/TCore/TField.hs + Math/Test/TProjects/TMiniquaternionGeometry.hs Library - Build-Depends: base >= 2 && < 5, containers, array, random, QuickCheck + Build-Depends: base >= 2 && < 5, containers, array, random Exposed-modules: Math.Algebra.LinearAlgebra, - Math.Algebra.Commutative.Monomial, - Math.Algebra.Commutative.MPoly, Math.Algebra.Commutative.GBasis, Math.Algebra.Field.Base, Math.Algebra.Field.Extension, Math.Algebra.Group.PermutationGroup, Math.Algebra.Group.SchreierSims, Math.Algebra.Group.RandomSchreierSims, Math.Algebra.Group.Subquotients, diff --git a/Math/Algebra/Commutative/GBasis.hs b/Math/Algebra/Commutative/GBasis.hs deleted file mode 100644 --- a/Math/Algebra/Commutative/GBasis.hs +++ /dev/null @@ -1,283 +0,0 @@ --- Copyright (c) David Amos, 2008. All rights reserved. - -{-# OPTIONS_GHC -fglasgow-exts #-} - -module Math.Algebra.Commutative.GBasis where - -import Data.List -import qualified Data.Map as M - -import Math.Algebra.Commutative.Monomial -import Math.Algebra.Commutative.MPoly - --- Sources: --- Cox, Little, O'Shea: Ideals, Varieties and Algorithms --- Giovini, Mora, Niesi, Robbiano, Traverso, "One sugar cube please, or Selection strategies in the Buchberger algorithm" - - - -sPoly f g = let h = lcmT (lt f) (lt g) - in h `divT` lt f .* f - h `divT` lt g .* g --- The point about the s-poly is that it cancels out the leading terms of the two polys, exposing their second terms - -isGB fs = all (\h -> h %% fs == 0) (pairWith sPoly fs) - - --- Cox p87 -gb1 fs = gb' fs (pairWith sPoly fs) where - gb' gs (h:hs) = let h' = h %% gs in - if h' == 0 then gb' gs hs else gb' (h':gs) (hs ++ map (sPoly h') gs) - gb' gs [] = gs - --- [f xi xj | xi <- xs, xj <- xs, i < j] -pairWith f (x:xs) = map (f x) xs ++ pairWith f xs -pairWith _ [] = [] - - --- Cox p89-90 -reduce gs = reduce' [] gs where - reduce' gs' (g:gs) | g' == 0 = reduce' gs' gs - | otherwise = reduce' (g':gs') gs - where g' = g %% (gs'++gs) - reduce' gs' [] = reverse $ sort $ map toMonic gs' --- the reverse means that when using an elimination order, the elimination ideal will be at the end - - --- Cox et al p106-7 --- No need to calculate an spoly fi fj if --- 1. the lm fi and lm fj are coprime, or --- 2. there exists some fk, with (i,k) (j,k) already considered, and lm fk divides lcm (lm fi) (lm fj) --- some slight inefficiencies from looking up fi, fj repeatedly -gb2 fs = reduce $ gb' fs (pairs [1..s]) s where - s = length fs - gb' gs ((i,j):ps) t = - let fi = gs!i; fj = gs!j in - if coprimeM (lm fi) (lm fj) || criterion fi fj - -- if lcmM (lm fi) (lm fj) == lm fi * lm fj || criterion fi fj - then gb' gs ps t - else let h = sPoly fi fj %% gs in - if h == 0 then gb' gs ps t else gb' (gs++[h]) (ps ++ [ (i,t+1) | i <- [1..t] ]) (t+1) - where - criterion fi fj = let l = lcmM (lm fi) (lm fj) in any (test l) [1..t] - test l k = k `notElem` [i,j] - && ordpair i k `notElem` ps - && ordpair j k `notElem` ps - && lm (gs!k) `dividesM` l - gb' gs [] _ = gs - - -pairs (x:xs) = map (\y->(x,y)) xs ++ pairs xs -pairs [] = [] - -xs ! i = xs !! (i-1) -- in other words, index the list from 1 not 0 - -ordpair x y | x < y = (x,y) - | otherwise = (y,x) - --- version of gb2 where we eliminate pairs as they're created, rather than as they're processed -gb2b fs = reduce $ gb1' [] fs [] 0 where - gb1' gs (f:fs) ps t = gb1' (gs ++ [f]) fs ps' (t+1) where - ps' = updatePairs gs ps f t - gb1' ls [] ps t = gb2' ls ps t - gb2' gs ((i,j):ps) t = - let h = sPoly (gs!i) (gs!j) %% gs in - if h == 0 - then gb2' gs ps t - else let ps' = updatePairs gs ((i,j):ps) h t in gb2' (gs++[h]) ps' (t+1) - gb2' gs [] _ = gs - updatePairs gs ps f t = - [p | p@(i,j) <- ps, - not (lm f `dividesM` lcmM (lm (gs!i)) (lm (gs!j)))] - ++ [ (i,t+1) | (gi,i) <- zip gs [1..t], - not (coprimeM (lm gi) (lm f)), - not (criterion (lcmM (lm gi) (lm f)) i) ] - where criterion l i = any (`dividesM` l) [lm gk | (gk,k) <- zip gs [1..t], k /= i, ordpair i k `notElem` ps] - - --- Cox et al 108 --- 1. list smallest fs first, as more likely to reduce --- 2. order the pairs with smallest lcm fi fj first ("normal selection strategy") -gb3b fs = - let fs' = sort $ filter (/=0) fs - in reduce $ gb1' [] fs' [] 0 where - gb1' gs (f:fs) ps t = gb1' (gs ++ [f]) fs ps' (t+1) where - ps' = updatePairs gs ps f t - gb1' ls [] ps t = gb2' ls ps t - gb2' gs ((l,(i,j)):ps) t = - let h = sPoly (gs!i) (gs!j) %% gs in - if h == 0 - then gb2' gs ps t - else let ps' = updatePairs gs ((l,(i,j)):ps) h t in gb2' (gs++[h]) ps' (t+1) - gb2' gs [] _ = gs - updatePairs :: (Ord (Monomial ord), Fractional r) => [MPoly ord r] -> [(Monomial ord, (Int,Int))] -> (MPoly ord r) -> Int -> [(Monomial ord, (Int,Int))] - updatePairs gs ps f t = - mergeBy cmpFst - [p | p@(l,(i,j)) <- ps, - not (lm f `dividesM` l)] - $ sortBy cmpFst - [ (l,(i,t+1)) | (gi,i) <- zip gs [1..t], l <- [lcmM (lm gi) (lm f)], - not (coprimeM (lm gi) (lm f)), - not (criterion l i) ] - where criterion l i = any (`dividesM` l) [lm gk | (gk,k) <- zip gs [1..t], k /= i, ordpair i k `notElem` map snd ps] - -cmpFst (a,_) (b,_) = compare a b - -mergeBy cmp (t:ts) (u:us) = - case cmp t u of - LT -> t : mergeBy cmp ts (u:us) - EQ -> t : mergeBy cmp ts (u:us) - GT -> u : mergeBy cmp (t:ts) us -mergeBy _ ts us = ts ++ us -- one of them is null - - --- naive implementation of "sugar strategy" -gb4b fs = - let fs' = sort $ filter (/=0) fs - in reduce $ gb1' [] fs' [] 0 where - gb1' gs (f:fs) ps t = gb1' (gs ++ [f]) fs ps' (t+1) where - ps' = updatePairs gs ps f t - gb1' ls [] ps t = gb2' ls ps t - gb2' gs ((sl,(i,j)):ps) t = - let h = sPoly (gs!i) (gs!j) %% gs in - if h == 0 - then gb2' gs ps t - else let ps' = updatePairs gs ((sl,(i,j)):ps) h t in gb2' (gs++[h]) ps' (t+1) - gb2' gs [] _ = gs - updatePairs :: (Ord (Monomial ord), Fractional r) => - [MPoly ord r] -> [((Int,Monomial ord), (Int,Int))] -> (MPoly ord r) -> Int -> [((Int,Monomial ord), (Int,Int))] - updatePairs gs ps f t = - mergeBy cmpFst - [p | p@((s,l),(i,j)) <- ps, - not (lm f `dividesM` l)] - $ sortBy cmpFst - [ ((s,l),(i,t+1)) | (gi,i) <- zip gs [1..t], l <- [lcmM (lm gi) (lm f)], s <- [sugar gi f l], - not (coprimeM (lm gi) (lm f)), - not (criterion l i) ] - where criterion l i = any (`dividesM` l) [lm gk | (gk,k) <- zip gs [1..t], k /= i, ordpair i k `notElem` map snd ps] - - --- Giovini et al --- The point of sugar is, given fi, fj, to give an upper bound on the degree of sPoly fi fj without having to calculate it --- We can then select by preference pairs with lower sugar, expecting therefore that the s-polys will have lower degree - --- |Given a list of polynomials over a field, return a Groebner basis for the ideal generated by the polynomials -gb :: (Ord (Monomial ord), Fractional k, Ord k) => - [MPoly ord k] -> [MPoly ord k] -gb fs = - -- let fs' = sort $ filter (/=0) fs - let fs' = sort $ map toMonic $ filter (/=0) fs - in reduce $ gb1' [] fs' [] 0 where - gb1' gs (f:fs) ps t = gb1' (gs ++ [f]) fs ps' (t+1) - where ps' = updatePairs gs ps f (t+1) - gb1' ls [] ps t = gb2' ls ps t - gb2' gs (p@(_,(i,j)):ps) t = - if h == 0 - then gb2' gs ps t - else gb2' (gs++[h]) ps' (t+1) - where h = toMonic $ sPoly (gs!i) (gs!j) %% gs - ps' = updatePairs gs (p:ps) h (t+1) - gb2' gs [] _ = gs - updatePairs :: (Ord (Monomial ord), Fractional r) => - [MPoly ord r] -> [((Int,Monomial ord), (Int,Int))] -> (MPoly ord r) -> Int -> [((Int,Monomial ord), (Int,Int))] - updatePairs gs ps gk k = - let newps = [let l = lcmM (lm gi) (lm gk) in ((sugar gi gk l, l), (i,k)) | (gi,i) <- zip gs [1..k-1] ] - ps' = [p | p@((sij,tij),(i,j)) <- ps, ((sik,tik),_) <- [newps ! i], ((sjk,tjk),_) <- [newps ! j], - not ( (tik `properlyDividesM` tij) && (tjk `properlyDividesM` tij) ) ] -- sloppy variant - newps' = discard1 [] newps - newps'' = sortBy cmpSug $ discard2 [] $ sortBy cmpNormal newps' - in mergeBy cmpSug ps' newps'' - where - discard1 ls (r@((_sik,tik),(i,_k)):rs) = - if lm (gs!i) `coprimeM` lm gk - -- then discard [l | l@((_,tjk),_) <- ls, tjk /= tik] [r | r@((_,tjk),_) <- ls, tjk /= tik] - then discard1 (filter (\((_,tjk),_) -> tjk /= tik) ls) (filter (\((_,tjk),_) -> tjk /= tik) rs) - else discard1 (r:ls) rs - discard1 ls [] = ls - discard2 ls (r@((_sik,tik),(i,k)):rs) = discard2 (r:ls) $ filter (\((_sjk,tjk),(j,k')) -> not (k == k' && tik `dividesM` tjk)) rs - discard2 ls [] = ls --- The type annotation on updatePairs appears to be required --- The two calls to toMonic are designed to prevent coefficient explosion, but it is unproven that they are effective - --- sugar of sPoly f g, where h = lcm (lt f) (lt g) --- this is an upper bound on deg (sPoly f g) -sugar f g h = degM h + max (deg f - degM (lm f)) (deg g - degM (lm g)) - -cmpNormal ((s1,t1),(i1,j1)) ((s2,t2),(i2,j2)) = compare (t1,j1) (t2,j2) - -cmpSug ((s1,t1),(i1,j1)) ((s2,t2),(i2,j2)) = compare (s1,t1,j1) (s2,t2,j2) - - --- earlier version of gb3b -gb3 fs = - let gs = sort $ filter (/=0) fs - ps = sortBy cmpFst $ pairWith (\(i,fi) (j,fj) -> (lcmM (lm fi) (lm fj), (i,j)) ) $ zip [1..] gs - in reduce $ gb' gs ps s - where - s = length fs - gb' :: (Ord (Monomial ord), Fractional r) => [MPoly ord r] -> [(Monomial ord, (Int,Int))] -> Int -> [MPoly ord r] - gb' gs ((l,(i,j)):ps) t = - let fi = gs!i; fj = gs!j in - if coprimeM (lm fi) (lm fj) || any (test l) [1..t] - then gb' gs ps t - else let h = sPoly fi fj %% gs in - if h == 0 - then gb' gs ps t - else let ps' = mergeBy cmpFst ps $ sortBy cmpFst $ zip [lcmM (lm h) (lm fi) | fi <- gs] [(i,t+1) | i <- [1..t]] - -- else let ps' = mergeBy cmpFst ps $ zip [lcmM (lm h) (lm fi) | fi <- gs] [(i,t+1) | i <- [1..t]] - in gb' (gs++[h]) ps' (t+1) - where - test l k = k `notElem` [i,j] - && ordpair i k `notElem` map snd ps - && ordpair j k `notElem` map snd ps - && lm (gs!k) `dividesM` l - gb' gs [] _ = gs --- Note that the type annotation on gb' appears to be required. I think this is a bug in the type inference algorithm - - - --- earlier version of gb4b -gb4 fs = - let gs = sort $ filter (/=0) fs - ps = sortBy cmpFst $ pairWith (\(i,fi) (j,fj) -> let l = lcmM (lm fi) (lm fj) in ((sugar fi fj l, l), (i,j)) ) $ zip [1..] gs - in reduce $ gb' gs ps s - where - s = length fs - gb' :: (Ord (Monomial ord), Fractional r) => [MPoly ord r] -> [((Int,Monomial ord), (Int,Int))] -> Int -> [MPoly ord r] - gb' gs (((s,l),(i,j)):ps) t = - let fi = gs!i; fj = gs!j in - if coprimeM (lm fi) (lm fj) || any (test l) [1..t] - then gb' gs ps t - else let h = sPoly fi fj %% gs in - if h == 0 - then gb' gs ps t - else let ps' = mergeBy cmpFst ps $ sortBy cmpFst $ zip [let l = lcmM (lm fi) (lm h) in (sugar fi h l, l) | fi <- gs] [(i,t+1) | i <- [1..t]] - in gb' (gs++[h]) ps' (t+1) - where - test l k = k `notElem` [i,j] - && ordpair i k `notElem` map snd ps - && ordpair j k `notElem` map snd ps - && lm (gs!k) `dividesM` l - gb' gs [] _ = gs --- Note that the type annotation on gb' appears to be required. I think this is a bug in the type inference algorithm - - - -{- -merge (t:ts) (u:us) = - if t <= u - then t : merge ts (u:us) - else u : merge (t:ts) us -merge ts us = ts ++ us -- one of them is null --} - --- OPERATIONS ON IDEALS - --- Cox et al, p181 --- Geometric interpretation: V(I+J) = V(I) `intersect` V(J) -sumI fs gs = gb $ fs ++ gs - --- Cox et al, p183 --- Geometric interpretation: V(I.J) = V(I) `union` V(J) -productI fs gs = gb [f * g | f <- fs, g <- gs] - - diff --git a/Math/Algebra/Commutative/MPoly.hs b/Math/Algebra/Commutative/MPoly.hs deleted file mode 100644 --- a/Math/Algebra/Commutative/MPoly.hs +++ /dev/null @@ -1,210 +0,0 @@ --- Copyright (c) David Amos, 2008. All rights reserved. - -{-# 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 -import Data.List as L -import Control.Arrow (first, second) -import Data.Ratio (denominator) - -import Math.Algebra.Field.Base -import Math.Algebra.Commutative.Monomial - - --- MULTIVARIATE POLYNOMIALS - --- |Type for multivariate polynomials. --- ord is a phantom type defining how terms are ordered, r is the type of the ring we are working over. --- For example, a common choice will be MPoly Grevlex Q, meaning polynomials over Q with the grevlex term ordering -newtype MPoly ord r = MP [(Monomial ord,r)] deriving (Eq) --- deriving instance (Ord (Monomial ord), Ord r) => Ord (MPoly ord r) --- standalone deriving supported from GHC 6.8 - -instance (Ord (Monomial ord), Ord r) => Ord (MPoly ord r) where - compare (MP ts) (MP us) = compare ts us - -instance (Show r, Num r) => Show (MPoly ord r) where - show (MP []) = "0" - show (MP ts) = - let (c:cs) = concatMap showTerm ts - in if c == '+' then cs else c:cs - where showTerm (m,c) = - case show c of - "1" -> "+" ++ show m - "-1" -> "-" ++ show m - cs@(x:_) -> (if x == '-' then cs else '+':cs) ++ (if m == 1 then "" else show m) - - -instance (Ord (Monomial ord), Num r) => Num (MPoly ord r) where - MP ts + MP us = MP (mergeTerms ts us) - negate (MP ts) = MP $ map (second negate) ts - MP ts * MP us = MP $ collect $ sortBy cmpTerm $ [(g*h,c*d) | (g,c) <- ts, (h,d) <- us] - {- - -- The following appears to be slightly slower, perhaps because sortBy is compiled - MP (t@(g,c):ts) * MP (u@(h,d):us) = - let MP vs = MP ts * MP us - in MP $ mergeTerms ((g*h,c*d):vs) $ mergeTerms [(g*h,c*d) | (h,d) <- us] [(g*h,c*d) | (g,c) <- ts] - _ * _ = MP [] - -} - fromInteger 0 = MP [] - fromInteger n = MP [(fromInteger 1, fromInteger n)] - -cmpTerm (a,c) (b,d) = case compare a b of EQ -> EQ; GT -> LT; LT -> GT -- in mpolys we put "larger" terms first - --- inputs in descending order -mergeTerms (t@(g,c):ts) (u@(h,d):us) = - case compare g h of - GT -> t : mergeTerms ts (u:us) - LT -> u : mergeTerms (t:ts) us - EQ -> if e == 0 then mergeTerms ts us else (g,e) : mergeTerms ts us - where e = c + d -mergeTerms ts us = ts ++ us -- one of them is null - -collect (t1@(g,c):t2@(h,d):ts) - | g == h = collect $ (g,c+d):ts - | c == 0 = collect $ t2:ts - | otherwise = t1 : collect (t2:ts) -collect ts = ts - --- Fractional instance so that we can enter fractional coefficients --- Only lets us divide by field elements (with unit monomial), not any other polynomials -instance (Ord (Monomial ord), Fractional r) => Fractional (MPoly ord r) where - recip (MP [(m,c)]) = MP [(recip m, recip c)] - -- recip (MP [(m,c)]) | m == fromInteger 1 = MP [(m, recip c)] - recip _ = error "MPoly.recip: only supported for (non-zero) constants or monomials" - --- |Create a variable with the supplied name. --- By convention, variable names should usually be a single letter followed by none, one or two digits. -var :: String -> MPoly Grevlex Q -var v = MP [(Monomial $ M.singleton v 1, 1)] :: MPoly Grevlex Q - -a, b, c, d, s, t, u, v, w, x, y, z :: MPoly Grevlex Q -a = var "a" -b = var "b" -c = var "c" -d = var "d" -s = var "s" -t = var "t" -u = var "u" -v = var "v" -w = var "w" -x = var "x" -y = var "y" -z = var "z" - -x_ i = var ("x" ++ show i) - -x0, x1, x2, x3 :: MPoly Grevlex Q -x0 = x_ 0 -x1 = x_ 1 -x2 = x_ 2 -x3 = x_ 3 - - --- convertMP :: Ord (Monomial ord') => MPoly ord k -> MPoly ord' k -convertMP (MP ts) = MP $ sortBy cmpTerm $ map (first convertM) ts - --- |Convert a polynomial to lex term ordering -toLex :: MPoly ord k -> MPoly Lex k -toLex = convertMP - --- |Convert a polynomial to glex term ordering -toGlex :: MPoly ord k -> MPoly Glex k -toGlex = convertMP - --- |Convert a polynomial to grevlex term ordering -toGrevlex :: MPoly ord k -> MPoly Grevlex k -toGrevlex = convertMP - -toElim :: MPoly ord k -> MPoly Elim k -toElim = convertMP - - -varLex v = toLex $ var v - -varElim v = toElim $ var v - - --- DIVISION ALGORITHM - -lt (MP (t:ts)) = t - -lm = fst . lt - -deg 0 = -1 -deg (MP ts) = maximum [degM m | (m,c) <- ts] --- the true degree of the polynomial, not the degree of the leading term --- required for sugar strategy when computing Groebner basis - -mulT (m,c) (m',c') = (m*m',c*c') - -divT (m,c) (m',c') = (m/m',c/c') - -dividesT (m,_) (m',_) = dividesM m m' - -properlyDividesT (m,_) (m',_) = dividesM m m' && m /= m' - -lcmT (m,c) (m',c') = (lcmM m m',1) - - -infixl 7 .* -t .* MP ts = MP $ map (mulT t) ts -- preserves term order - - --- 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 = MP [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 (MP (t:ts)) = (MP [t], MP ts) - -infixl 7 %% -f %% gs = r where (_,r) = quotRemMP f gs - --- div and mod by single mpoly -divModMP f g = (q,r) where ([q],r) = quotRemMP f [g] - -divMP f g = q where ([q],_) = quotRemMP f [g] - -modMP f g = r where (_,r) = quotRemMP f [g] - - --- OTHER STUFF - --- injection of field elements into polynomial ring -inject 0 = MP [] -inject c = MP [(fromInteger 1, c)] - -toMonic 0 = 0 -toMonic (MP ts@((_,c):_)) - | c == 1 = MP ts - | otherwise = MP $ map (second (/c)) ts - --- multiply out all denominators -toZ (MP ts) = MP $ map (second (*c)) ts where c = fromInteger $ foldl lcm 1 $ [denominator c | (m,Q c) <- ts] - --- substitute terms for variables in an MPoly --- eg subst [(x,a),(y,a+b),(z,c^2)] (x*y+z) -> a*(a+b)+c^2 -subst vts (MP us) = sum [inject c * substM m | (m,c) <- us] where - substM (Monomial m) = product [substV v ^ i | (v,i) <- M.toList m] - substV v = - let v' = MP [(Monomial $ M.singleton v 1, 1)] in - case L.lookup v' vts of - Just t -> t - Nothing -> v' -- no substitute, so keep as is - -support (MP ts) = [MP [(m,1)] | m <- reverse $ L.sort $ support' ts] - where support' ts = foldl L.union [] [supportM m | (m,c) <- ts] diff --git a/Math/Algebra/Commutative/Monomial.hs b/Math/Algebra/Commutative/Monomial.hs deleted file mode 100644 --- a/Math/Algebra/Commutative/Monomial.hs +++ /dev/null @@ -1,102 +0,0 @@ --- Copyright (c) David Amos, 2008. All rights reserved. - -{-# OPTIONS_GHC -fglasgow-exts #-} - -module Math.Algebra.Commutative.Monomial where - -import qualified Data.Map as M -import Data.List as L -import Control.Arrow - - --- MONOMIAL - -newtype Monomial ord = Monomial (M.Map String Int) deriving (Eq) - -instance Show (Monomial ord) where - show (Monomial a) | M.null a = "1" - | otherwise = concatMap showVar $ M.toList a - where showVar (v,1) = v - showVar (v,i) = v ++ "^" ++ show i - -instance Num (Monomial ord) where - Monomial a * Monomial b = Monomial $ M.filter (/=0) $ M.unionWith (+) a b - -- The filter (/=0) here means we can handle Laurent monomials, and isn't significantly slower - -- Monomial a * Monomial b = Monomial $ M.unionWith (+) a b - fromInteger 1 = Monomial M.empty - -instance Fractional (Monomial ord) where - recip (Monomial m) = Monomial $ M.map negate m - -- can only do the above if (*) is doing filter (/=0) - -- Monomial a / Monomial b = Monomial $ M.filter (/=0) $ M.unionWith (+) a (M.map negate b) - --- |Phantom type representing lex term ordering -data Lex - --- |Phantom type representing glex term ordering -data Glex - --- |Phantom type representing grevlex term ordering -data Grevlex - --- |Phantom type for an elimination term ordering. --- In the ordering, xis come before yjs come before zks, but within the xis, or yjs, or zks, grevlex ordering is used -data Elim - - -diffs a b = M.elems m where Monomial m = a/b --- note that we're guaranteed that all elts of diffs a b are non-zero - -instance Ord (Monomial Lex) where - compare a b = case diffs a b of - [] -> EQ - as -> if head as > 0 then GT else LT - -instance Ord (Monomial Glex) where - compare a b = let ds = diffs a b in - case compare (sum ds) 0 of - GT -> GT - LT -> LT - EQ -> if null ds then EQ else - if head ds > 0 then GT else LT - -instance Ord (Monomial Grevlex) where - compare a b = let ds = diffs a b in - case compare (sum ds) 0 of - GT -> GT - LT -> LT - EQ -> if null ds then EQ else - if last ds < 0 then GT else LT - --- a monomial order for terms in x1,x2,..,y1,y2,..,z1,z2,.. etc --- in which it is grevlex separately on first the xs, then if they're equal, the ys, then if they're equal, the zs, etc -instance Ord (Monomial Elim) where - compare a b = let Monomial m = a/b in - case M.assocs m of - [] -> EQ - (l:s,i):vs -> grevlex $ i : map snd (takeWhile (\(l':_,_) -> l'==l) vs) - where grevlex ds = case compare (sum ds) 0 of - GT -> GT - LT -> LT - EQ -> if last ds < 0 then GT else LT - - -convertM :: Monomial a -> Monomial b -convertM (Monomial x) = Monomial x - - -degM (Monomial m) = sum $ M.elems m - -dividesM (Monomial a) (Monomial b) = M.isSubmapOfBy (<=) a b - -properlyDividesM a b = dividesM a b && a /= b - -lcmM (Monomial a) (Monomial b) = Monomial $ M.unionWith max a b - -gcdM (Monomial a) (Monomial b) = Monomial $ M.intersectionWith min a b - -coprimeM (Monomial a) (Monomial b) = M.null $ M.intersection a b - --- the support of a monomial is the variables it contains -supportM :: Monomial ord -> [Monomial ord] -- type signature needed to say that output has same term order as input -supportM (Monomial m) = [Monomial (M.singleton v 1) | v <- M.keys m] diff --git a/Math/Algebra/LinearAlgebra.hs b/Math/Algebra/LinearAlgebra.hs --- a/Math/Algebra/LinearAlgebra.hs +++ b/Math/Algebra/LinearAlgebra.hs @@ -13,7 +13,7 @@ module Math.Algebra.LinearAlgebra where import qualified Data.List as L -import Math.Algebra.Field.Base -- not actually used in this module +import Math.Core.Field -- not actually used in this module infixr 8 *>, *>> diff --git a/Math/CommutativeAlgebra/GroebnerBasis.hs b/Math/CommutativeAlgebra/GroebnerBasis.hs --- a/Math/CommutativeAlgebra/GroebnerBasis.hs +++ b/Math/CommutativeAlgebra/GroebnerBasis.hs @@ -280,3 +280,97 @@ quotientP fs g = map ( // g ) $ intersectI fs [g] where h // g = let ([u],_) = quotRemMP h [g] in u +-- |@eliminate vs gs@ returns the elimination ideal obtained from the ideal generated by gs by eliminating the variables vs. +eliminate :: (Fractional k, Ord k, MonomialConstructor m, Monomial (m v), Ord (m v)) => + [Vect k (m v)] -> [Vect k (m v)] -> [Vect k (m v)] +eliminate vs gs = let subs = subFst vs in eliminateFst [g `bind` subs | g <- gs] + where subFst :: (Num k, MonomialConstructor m, Eq (m v), Mon (m v)) => + [Vect k (m v)] -> v -> Vect k (Elim2 (m v) (m v)) + subFst vs = (\v -> let v' = var v in if v' `elem` vs then toElimFst v' else toElimSnd v') + +{- +-- !! NOT WORKING +-- |@elimExcept vs gs@ returns the elimination ideal obtained from the ideal generated by gs by eliminating all variables except vs. +elimExcept :: (Fractional k, Ord k, MonomialConstructor m, Monomial (m v), Ord (m v)) => + [Vect k (m v)] -> [Vect k (m v)] -> [Vect k (m v)] +elimExcept vs gs = let subs = subSnd vs in eliminateFst [g `bind` subs | g <- gs] + where subSnd :: (Num k, MonomialConstructor m, Eq (m v), Mon (m v)) => + [Vect k (m v)] -> v -> Vect k (Elim2 (m v) (m v)) + subSnd vs = (\v -> let v' = var v in if v' `elem` vs then toElimSnd v' else toElimFst v') +-} + +-- MONOMIAL BASES FOR QUOTIENT ALGEBRAS + +-- basis for the polynomial ring in variables vs +mbasis vs = mbasis' [1] + where mbasis' ms = ms ++ mbasis' (toSet [v*m | v <- vs, m <- ms]) + +-- |Given variables vs, and a Groebner basis gs, @mbasisQA vs gs@ returns a monomial basis for the quotient algebra k[vs]/\. +-- For example, @mbasisQA [x,y] [x^2+y^2-1]@ returns a monomial basis for k[x,y]/\. +-- In general, the monomial basis is likely to be infinite. +mbasisQA :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) => + [Vect k m] -> [Vect k m] -> [Vect k m] +mbasisQA vs gs = mbasisQA' [1] + where mbasisQA' [] = [] -- the quotient algebra is finite-dimensional + mbasisQA' ms = ms ++ mbasisQA' (toSet [f | v <- vs, m <- ms, let f = v*m, f %% gs == f]) + +-- |Given an ideal I, the leading term ideal lt(I) consists of the leading terms of all elements of I. +-- If I is generated by gs, then @ltIdeal gs@ returns generators for lt(I). +ltIdeal :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) => + [Vect k m] -> [Vect k m] +ltIdeal gs = map (return . lm) $ gb gs + +-- number of monomials of degree i in n variables +numMonomials n i = toInteger (i+n-1) `choose` toInteger (n-1) + +-- |Given variables vs, and a homogeneous ideal gs, @hilbertFunQA vs gs@ returns the Hilbert function for the quotient algebra k[vs]/\. +-- Given an integer i, the Hilbert function returns the number of degree i monomials in a basis for k[vs]/\. +-- For a homogeneous ideal, this number is independent of the monomial ordering used +-- (even though the elements of the monomial basis themselves are dependent on the ordering). +-- +-- If the ideal I is not homogeneous, then R/I is not graded, and the Hilbert function is not well-defined. +-- Specifically, the number of degree i monomials in a basis is likely to depend on which monomial ordering you use. +hilbertFunQA :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) => + [Vect k m] -> [Vect k m] -> Int -> Integer +hilbertFunQA vs gs i = hilbertFunQA' (ltIdeal gs) i + where n = length vs + hilbertFunQA' _ i | i < 0 = 0 + hilbertFunQA' (m:ms) i = hilbertFunQA' ms i - hilbertFunQA' (ms `quotientP` m) (i - deg m) + hilbertFunQA' [] i = numMonomials n i +-- For example, consider k[x,y]/ +-- Under Lex ordering, the monomial basis is 1,y,y^2,y^3,... +-- Under Glex ordering, the monomial basis is 1,x,y,x^2,xy,x^3,x^2y,... +-- So the Hilbert function is not well-defined. +-- Note though that this function does still correctly return the number of degree i monomials for the given monomial ordering + +-- naive implementation which simply counts monomials +hilbertSeriesQA1 vs gs = hilbertSeriesQA1' [1] + where hilbertSeriesQA1' [] = [] -- repeat 0 + hilbertSeriesQA1' ms = length ms : hilbertSeriesQA1' (toSet [f | v <- vs, m <- ms, let f = v*m, f %% gs == f]) + +-- Eisenbud p325, p357 / Schenck p56 +-- This can be made more efficient by choosing which m to recurse on +-- |Given variables vs, and a homogeneous ideal gs, @hilbertSeriesQA vs gs@ returns the Hilbert series for the quotient algebra k[vs]/\. +-- The Hilbert series should be interpreted as a formal power series where the coefficient of t^i is the Hilbert function evaluated at i. +-- That is, the i'th element in the series is the number of degree i monomials in a basis for k[vs]/\. +hilbertSeriesQA :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) => + [Vect k m] -> [Vect k m] -> [Integer] +hilbertSeriesQA vs gs = hilbertSeriesQA' $ ltIdeal gs + where hilbertSeriesQA' (m:ms) = hilbertSeriesQA' ms <-> (replicate (deg m) 0 ++ hilbertSeriesQA' (ms `quotientI` [m])) + hilbertSeriesQA' [] = [numMonomials n i | i <- [0..] ] + n = length vs + (a:as) <-> (b:bs) = (a-b) : (as <-> bs) + as <-> [] = as + [] <-> bs = map negate bs + +-- |For i \>\> 0, the Hilbert function becomes a polynomial in i, called the Hilbert polynomial. +hilbertPolyQA :: (Fractional k, Ord k, Monomial m, Ord m, Algebra k m) => + [Vect k m] -> [Vect k m] -> GlexPoly Q String +hilbertPolyQA vs gs = hilbertPolyQA' (ltIdeal gs) i + where n = toInteger $ length vs + i = glexvar "i" + hilbertPolyQA' [] x = product [ x + fromInteger j | j <- [1..n-1] ] / (fromInteger $ product [1..n-1]) + hilbertPolyQA' (m:ms) x = hilbertPolyQA' ms x - hilbertPolyQA' (ms `quotientP` m) (x - fromIntegral (deg m)) + +-- The dimension of a variety +dim vs gs = 1 + deg (hilbertPolyQA vs gs) diff --git a/Math/CommutativeAlgebra/Polynomial.hs b/Math/CommutativeAlgebra/Polynomial.hs --- a/Math/CommutativeAlgebra/Polynomial.hs +++ b/Math/CommutativeAlgebra/Polynomial.hs @@ -7,6 +7,10 @@ -- -- A monomial ordering is required to specify how monomials are to be ordered. -- The Lex, Glex, and Grevlex monomial orders are defined, with the possibility to add others. +-- +-- In order to make use of this module, some variables must be defined, for example: +-- +-- > [t,u,v,x,y,z] = map glexvar ["t","u","v","x","y","z"] module Math.CommutativeAlgebra.Polynomial where import Math.Core.Field @@ -160,6 +164,7 @@ -- > [x,y,z] = map lexvar ["x","y","z"] lexvar :: v -> LexPoly Q v lexvar v = return $ Lex $ M 1 [(v,1)] +-- lexvar = var instance (Num k, Ord v, Show v) => Algebra k (Lex v) where unit x = x *> return munit @@ -193,6 +198,7 @@ -- > [x,y,z] = map glexvar ["x","y","z"] glexvar :: v -> GlexPoly Q v glexvar v = return $ Glex $ M 1 [(v,1)] +-- glexvar = var instance (Num k, Ord v, Show v) => Algebra k (Glex v) where unit x = x *> return munit @@ -228,6 +234,7 @@ -- > [x,y,z] = map grevlexvar ["x","y","z"] grevlexvar :: v -> GrevlexPoly Q v grevlexvar v = return $ Grevlex $ M 1 [(v,1)] +-- grevlexvar = var instance (Num k, Ord v, Show v) => Algebra k (Grevlex v) where unit x = x *> return munit @@ -268,14 +275,55 @@ -- VARIABLE SUBSTITUTION --- 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 v constraint, --- secondly because Haskell doesn't support type functions. -bind :: (MonomialConstructor m, Num k, Ord a, Show a, Algebra k a) => +-- |Given (Num k, MonomialConstructor m), then Vect k (m v) is the free commutative algebra over v. +-- As such, it is a monad (in the mathematical sense). The following pseudo-code (not legal Haskell) +-- shows how this would work: +-- +-- > instance (Num k, Monomial m) => Monad (\v -> Vect k (m v)) where +-- > return = var +-- > (>>=) = bind +-- +-- bind corresponds to variable substitution, so @v `bind` f@ returns the result of making the substitutions +-- encoded in f into v. +-- +-- Note that the type signature is slightly more general than that required by (>>=). +-- For a monad, we would only require: +-- +-- > bind :: (MonomialConstructor m, Num k, Ord (m v), Show (m v), Algebra k (m v)) => +-- > Vect k (m u) -> (u -> Vect k (m v)) -> Vect k (m v) +-- +-- Instead, the given type signature allows us to substitute in elements of any algebra. +-- This is occasionally useful. + +-- |bind performs variable substitution +bind :: (Num k, MonomialConstructor m, Ord a, Show a, Algebra k a) => Vect k (m v) -> (v -> Vect k a) -> Vect k a -V ts `bind` f = sum [c *> product [f x ^ i | (x,i) <- mindices m] | (m, c) <- ts] +v `bind` f = linear (\m -> product [f x ^ i | (x,i) <- mindices m]) v +-- V ts `bind` f = sum [c *> product [f x ^ i | (x,i) <- mindices m] | (m, c) <- ts] +-- We can't express the Monad instance directly in Haskell, firstly because of the Ord v constraint (? - not used), +-- secondly because Haskell doesn't support type functions. + flipbind f = linear (\m -> product [f x ^ i | (x,i) <- mindices m]) + +-- |Evaluate a polynomial at a point. +-- For example @eval (x^2+y^2) [(x,1),(y,2)]@ evaluates x^2+y^2 at the point (x,y)=(1,2). +eval :: (Num k, MonomialConstructor m, Eq (m v), Show v) => + Vect k (m v) -> [(Vect k (m v), k)] -> k +eval f vs = unwrap $ f `bind` sub + where sub x = case lookup (var x) vs of + Just xval -> xval *> return () + Nothing -> error ("eval: no binding given for " ++ show x) + +-- |Perform variable substitution on a polynomial. +-- For example @subst (x*z-y^2) [(x,u^2),(y,u*v),(z,v^2)]@ performs the substitution x -> u^2, y -> u*v, z -> v^2. +subst :: (Num k, MonomialConstructor m, Eq (m u), Show u, Ord (m v), Show (m v), Algebra k (m v)) => + Vect k (m u) -> [(Vect k (m u), Vect k (m v))] -> Vect k (m v) +subst f vs = f `bind` sub + where sub x = case lookup (var x) vs of + Just xsub -> xsub + Nothing -> error ("eval: no binding given for " ++ show x) +-- The type could be more general than this, but haven't so far found a use case -- DIVISION ALGORITHM FOR POLYNOMIALS diff --git a/Math/Projects/MiniquaternionGeometry.hs b/Math/Projects/MiniquaternionGeometry.hs --- a/Math/Projects/MiniquaternionGeometry.hs +++ b/Math/Projects/MiniquaternionGeometry.hs @@ -18,10 +18,7 @@ import Math.Projects.ChevalleyGroup.Classical -import Test.QuickCheck - - -- Sources: -- Miniquaternion Geometry, Room & Kirkpatrick -- Survey of Non-Desarguesian Planes, Charles Weibel @@ -131,34 +128,6 @@ instance FiniteField J9 where basisFq _ = [1,i,j,k] eltsFq _ = j9 - - --- Near fields - -prop_NearField (a,b,c) = - a+(b+c) == (a+b)+c && -- addition is associative - a+b == b+a && -- addition is commutative - a+0 == a && -- additive identity - a+(-a) == 0 && -- additive inverse - a*(b*c) == (a*b)*c && -- multiplication is associative - a*1 == a && 1*a == a && -- multiplicative identity - (a+b)*c == a*c + b*c && -- right-distributivity - a*0 == 0 - -instance Arbitrary F9 where - arbitrary = do x <- arbitrary :: Gen Int - return (f9 !! (x `mod` 9)) - -instance Arbitrary J9 where - arbitrary = do x <- arbitrary :: Gen Int - return (j9 !! (x `mod` 9)) - -prop_NearFieldF9 (a,b,c) = prop_NearField (a,b,c) where - types = (a,b,c) :: (F9,F9,F9) - -prop_NearFieldJ9 (a,b,c) = prop_NearField (a,b,c) where - types = (a,b,c) :: (J9,J9,J9) - -- PROJECTIVE PLANES diff --git a/Math/Test/TCommutativeAlgebra.hs b/Math/Test/TCommutativeAlgebra.hs deleted file mode 100644 --- a/Math/Test/TCommutativeAlgebra.hs +++ /dev/null @@ -1,95 +0,0 @@ --- Copyright (c) David Amos, 2008. All rights reserved. - -{-# LANGUAGE FlexibleInstances #-} - -module Math.Test.TCommutativeAlgebra where - -import Math.Algebra.Field.Base -import Math.Algebra.Commutative.Monomial -import Math.Algebra.Commutative.MPoly -import Math.Algebra.Commutative.GBasis - -import Test.QuickCheck - --- > quickCheck prop_CommRingMPoly --- > verboseCheck prop_ComRingMPoly -- to see what input data is being used - --- Commutative Ring (with 1) -prop_CommRing (a,b,c) = - a+(b+c) == (a+b)+c && -- addition is associative - a+b == b+a && -- addition is commutative - a+0 == a && -- additive identity - a+(-a) == 0 && -- additive inverse - a*(b*c) == (a*b)*c && -- multiplication is associative - a*b == b*a && -- multiplication is commutative - a*1 == a && -- multiplicative identity - a*(b+c) == a*b + a*c -- distributivity - -monomial is = product $ zipWith (^) (map x_ [1..]) (map (max 0) is) - --- mpoly :: [(Integer,[Int])] -> MPoly Grevlex Q -mpoly ais = sum [fromInteger a * monomial is | (a,is) <- ais] - -{- --- can take a long time to run, probably because of the test for associativity of multiplication -prop_CommRingMPoly (ais,bjs,cks) = prop_CommRing (f,g,h) where - f = mpoly ais - g = mpoly bjs - h = mpoly cks - types = (ais,bjs,cks) :: ( [(Integer,[Int])], [(Integer,[Int])], [(Integer,[Int])] ) --} - -instance Arbitrary (MPoly Grevlex Q) where - -- arbitrary = do ais <- arbitrary :: Gen [(Integer,[Int])] - arbitrary = do ais <- sized $ \n -> resize (n `div` 2) arbitrary :: Gen [(Integer,[Int])] - return (mpoly ais) - -- coarbitrary = undefined -- !! only required if we want to test functions over the type - -prop_CommRingMPoly (f,g,h) = prop_CommRing (f,g,h) where - types = (f,g,h) :: (MPoly Grevlex Q, MPoly Grevlex Q, MPoly Grevlex Q) - - --- Sources for tests: --- [IVA] - Cox, Little, O'Shea: Ideals, Varieties and Algorithms --- [UAG] - Cox, Little, O'Shea: Using Algebraic Geometry - - -test = and [ - gb (map toGlex [x*z-y^2,x^3-z^2]) == map toGlex [y^6-z^5,x*y^4-z^4,x^2*y^2-z^3,x^3-z^2,x*z-y^2], -- IVA p93 - gb (map toLex [x^2+y^2+z^2-1,x^2+z^2-y,x-z]) == map toLex [x-z,y-2*z^2,z^4+1/2*z^2-1/4], -- IVA p94 - gb (map toLex [x^2+y^2+z^2-1,x*y*z-1]) == map toLex [x+y^3*z+y*z^3-y*z,y^4*z^2+y^2*z^4-y^2*z^2+1], -- IVA p116 - gb [x*y+z-x*z,x^2-z,2*x^3-x^2*y*z-1] == [z^4-3*z^3-4*y*z+2*z^2-y+2*z-2,y*z^2+2*y*z-2*z^2+1,y^2-2*y*z+z^2-z,x+y-z] -- Grevlex, UAG p50-1 - ] - - - -{- -http://www.cs.amherst.edu/~dac/iva.html -states that IVA, 2nd ed, 5th printing (the one I have) has a production error causing many +s and -s to appear incorrectly - -This explains the following misprints I've found: -p117: -gb (map toLex [x*y-4,y^2-(x^3-1)]) --> [x-1/16y^4-1/16y^2,y^5+y^3-64] -IVA p117 claims it should be -y^3 in the second poly -But my answer is clearly correct, by looking at the reduction sequence for x*y-4 -x*y-4 -> 1/16(y^5+y^3)-4 -> 0 - x-1/16(y^4+y^2) y^5+y^3-64 -By contrast, reducing over their set clearly stops at 1/8y^3 - -gb (map toLex [x-t-u,y-t^2-2*t*u,z-t^3-3*t^2*u]) -The answer I get has some sign differences compared to IVA p127 --} - -{- -The code has no trouble chomping through some of the examples that took a long time in the Sugar paper, eg -gb [x+y+z+t+u, x*y+y*z+z*t+t*u+u*x, x*y*z+y*z*t+z*t*u+t*u*x+u*x*y, x*y*z*t+y*z*t*u+z*t*u*x+t*u*x*y+u*x*y*z, x*y*z*t*u-1] -gb $ map toLex [x+y+z+t+u, x*y+y*z+z*t+t*u+u*x, x*y*z+y*z*t+z*t*u+t*u*x+u*x*y, x*y*z*t+y*z*t*u+z*t*u*x+t*u*x*y+u*x*y*z, x*y*z*t*u-1] -gb [w^31-w^6-w-x, w^8-y, w^10-z] -gb $ map toLex [w^31-w^6-w-x, w^8-y, w^10-z] - -However, for some reason, the code gets indigestion on the following -gb $ map toLex [y*(1+x^2)^4 - 2*(5+19*x^2-45*x^4+x^6-4*x^8), z*(1+x^2)^4-2*(x+51*x^3+3*x^5+17*x^7)] - -(For comparison, the v1 implementation of gbasis can manage, even though its performance on the sugar examples is only comparable) --} diff --git a/Math/Test/TCommutativeAlgebra/TGroebnerBasis.hs b/Math/Test/TCommutativeAlgebra/TGroebnerBasis.hs --- a/Math/Test/TCommutativeAlgebra/TGroebnerBasis.hs +++ b/Math/Test/TCommutativeAlgebra/TGroebnerBasis.hs @@ -20,18 +20,22 @@ testlistLexGb, testlistGlexGb, testlistIntersectI, - testlistQuotientI + testlistQuotientI, + testlistEliminate, + -- testlistElimExcept, + testlistHilbertPolyQA ] -data Var = X | Y | Z deriving (Eq,Ord) +data Var = X | Y | Z | W deriving (Eq,Ord) instance Show Var where show X = "x" show Y = "y" show Z = "z" + show W = "w" -[x,y,z] = map grevlexvar [X,Y,Z] +[x,y,z,w] = map grevlexvar [X,Y,Z,W] testcaseGb desc input output = TestCase (assertEqual desc output (gb input)) @@ -73,3 +77,37 @@ testcaseQuotientI "[x*z, y*z] [z]" [x*z, y*z] [z] [x,y], -- IVA p192 testcaseQuotientI "[y^2, z^2] [y*z]" [y^2, z^2] [y*z] [y,z] -- Schenck p56 (in passing) ] + + +testcaseEliminate vs gs gs' = TestCase $ assertEqual "Eliminate" gs' (eliminate vs gs) + +testlistEliminate = + let [t,u,v,x,y,z,x',y',z'] = map glexvar ["t","u","v","x","y","z","x'","y'","z'"] in + TestList [ + testcaseEliminate [x,y,z] [x^2+y^2-z^2,x'-(x+z),y'-y,z'-(z-x)] [x'*z'-y'^2], -- Reid p15 + testcaseEliminate [t] [(t^2+1)*x-2*t, (t^2+1)*y-(t^2-1)] [x^2+y^2-1], + testcaseEliminate [u,v] [x'-u^2,y'-u*v,z'-v^2] [x'*z'-y'^2] -- Reid p16 + ] + +{- +testcaseElimExcept vs gs gs' = TestCase $ assertEqual "Eliminate" gs' (eliminate vs gs) + +testlistElimExcept = + let [t,u,v,x,y,z,x',y',z'] = map glexvar ["t","u","v","x","y","z","x'","y'","z'"] in + TestList [ + testcaseElimExcept [x',y',z'] [x^2+y^2-z^2,x'-(x+z),y'-y,z'-(z-x)] [x'*z'-y'^2], -- Reid p15 + testcaseElimExcept [x,y] [(t^2+1)*x-2*t, (t^2+1)*y-(t^2-1)] [x^2+y^2-1], + testcaseElimExcept [x',y',z'] [x'-u^2,y'-u*v,z'-v^2] [x'*z'-y'^2] -- Reid p16 + ] +-} + +testcaseHilbertPolyQA desc hp vs gs = + TestCase $ assertEqual desc hp (hilbertPolyQA vs gs) + +testlistHilbertPolyQA = + let i = glexvar "i" in + TestList [ + testcaseHilbertPolyQA "hilbertPoly " (3*i+1) [x,y,z,w] [y*z-x*w,z^2-y*w,y^2-x*z] -- Schenck p56-7 + ] + + diff --git a/Math/Test/TCore/TField.hs b/Math/Test/TCore/TField.hs --- a/Math/Test/TCore/TField.hs +++ b/Math/Test/TCore/TField.hs @@ -59,31 +59,32 @@ arbitrary = do {n <- arbitrary; return (f25 !! mod (fromInteger n) 25)} -test = do putStrLn "Testing F2..." - quickCheck (prop_Field :: (F2,F2,F2) -> Bool) - putStrLn "Testing F3..." - quickCheck (prop_Field :: (F3,F3,F3) -> Bool) - putStrLn "Testing F5..." - quickCheck (prop_Field :: (F5,F5,F5) -> Bool) - putStrLn "Testing F7..." - quickCheck (prop_Field :: (F7,F7,F7) -> Bool) - putStrLn "Testing F11..." - quickCheck (prop_Field :: (F11,F11,F11) -> Bool) - putStrLn "Testing F13..." - quickCheck (prop_Field :: (F13,F13,F13) -> Bool) - putStrLn "Testing F17..." - quickCheck (prop_Field :: (F17,F17,F17) -> Bool) - putStrLn "Testing F19..." - quickCheck (prop_Field :: (F19,F19,F19) -> Bool) - putStrLn "Testing F23..." - quickCheck (prop_Field :: (F23,F23,F23) -> Bool) - putStrLn "Testing F4..." - quickCheck (prop_Field :: (F4,F4,F4) -> Bool) - putStrLn "Testing F8..." - quickCheck (prop_Field :: (F8,F8,F8) -> Bool) - putStrLn "Testing F9..." - quickCheck (prop_Field :: (F9,F9,F9) -> Bool) - putStrLn "Testing F16..." - quickCheck (prop_Field :: (F16,F16,F16) -> Bool) - putStrLn "Testing F25..." - quickCheck (prop_Field :: (F25,F25,F25) -> Bool) +quickCheckField = + do putStrLn "Testing F2..." + quickCheck (prop_Field :: (F2,F2,F2) -> Bool) + putStrLn "Testing F3..." + quickCheck (prop_Field :: (F3,F3,F3) -> Bool) + putStrLn "Testing F5..." + quickCheck (prop_Field :: (F5,F5,F5) -> Bool) + putStrLn "Testing F7..." + quickCheck (prop_Field :: (F7,F7,F7) -> Bool) + putStrLn "Testing F11..." + quickCheck (prop_Field :: (F11,F11,F11) -> Bool) + putStrLn "Testing F13..." + quickCheck (prop_Field :: (F13,F13,F13) -> Bool) + putStrLn "Testing F17..." + quickCheck (prop_Field :: (F17,F17,F17) -> Bool) + putStrLn "Testing F19..." + quickCheck (prop_Field :: (F19,F19,F19) -> Bool) + putStrLn "Testing F23..." + quickCheck (prop_Field :: (F23,F23,F23) -> Bool) + putStrLn "Testing F4..." + quickCheck (prop_Field :: (F4,F4,F4) -> Bool) + putStrLn "Testing F8..." + quickCheck (prop_Field :: (F8,F8,F8) -> Bool) + putStrLn "Testing F9..." + quickCheck (prop_Field :: (F9,F9,F9) -> Bool) + putStrLn "Testing F16..." + quickCheck (prop_Field :: (F16,F16,F16) -> Bool) + putStrLn "Testing F25..." + quickCheck (prop_Field :: (F25,F25,F25) -> Bool) diff --git a/Math/Test/TGraph.hs b/Math/Test/TGraph.hs --- a/Math/Test/TGraph.hs +++ b/Math/Test/TGraph.hs @@ -50,7 +50,8 @@ map is2ArcTransitive [c 7, q 3, G.to1n coxeterGraph] ++ map is3ArcTransitive [c 7, G.to1n petersen] ++ map (not . is3ArcTransitive) [q 3] ++ - [isArcTransitive (j v k i) | v <- [3..5], k <- [1..v `div` 2], i <- [0..k] ] ++ -- [AGT] p60 + -- [isArcTransitive (j v k i) | v <- [3..5], k <- [1..v `div` 2], i <- [0..k] ] ++ -- [AGT] p60 + -- !! j 4 2 0 is not connected, so this test now gives error. Not sure how it passed before [is2ArcTransitive (j (2*k+1) k 0) | k <- [1..2] ] ++ [isDistanceTransitive (j v k (k-1)) | v <- [3..5], k <- [1..v `div` 2] ] ++ -- [AGT] p75 [isDistanceTransitive (j (2*k+1) k 0) | k <- [1..2] ] ++ diff --git a/Math/Test/TProjects/TMiniquaternionGeometry.hs b/Math/Test/TProjects/TMiniquaternionGeometry.hs new file mode 100644 --- /dev/null +++ b/Math/Test/TProjects/TMiniquaternionGeometry.hs @@ -0,0 +1,52 @@ +-- Copyright (c) David Amos, 2009-2011. All rights reserved. + +module Math.Test.TProjects.TMiniquaternionGeometry where + +{- +import qualified Data.List as L + +import Math.Common.ListSet as LS + +import Math.Algebra.Field.Base +import Math.Combinatorics.FiniteGeometry (pnf, ispnf, orderPGL) +import Math.Combinatorics.Graph (combinationsOf) +import Math.Combinatorics.GraphAuts +import Math.Algebra.Group.PermutationGroup hiding (order) +import qualified Math.Algebra.Group.SchreierSims as SS +import Math.Algebra.Group.RandomSchreierSims +import Math.Combinatorics.Design as D +import Math.Algebra.LinearAlgebra -- ( (<.>), (<+>) ) + +import Math.Projects.ChevalleyGroup.Classical +-} + +import Test.QuickCheck +import Math.Projects.MiniquaternionGeometry + + +-- Near fields + +prop_NearField (a,b,c) = + a+(b+c) == (a+b)+c && -- addition is associative + a+b == b+a && -- addition is commutative + a+0 == a && -- additive identity + a+(-a) == 0 && -- additive inverse + a*(b*c) == (a*b)*c && -- multiplication is associative + a*1 == a && 1*a == a && -- multiplicative identity + (a+b)*c == a*c + b*c && -- right-distributivity + a*0 == 0 + +instance Arbitrary F9 where + arbitrary = do x <- arbitrary :: Gen Int + return (f9 !! (x `mod` 9)) + +instance Arbitrary J9 where + arbitrary = do x <- arbitrary :: Gen Int + return (j9 !! (x `mod` 9)) + +prop_NearFieldF9 (a,b,c) = prop_NearField (a,b,c) where + types = (a,b,c) :: (F9,F9,F9) + +prop_NearFieldJ9 (a,b,c) = prop_NearField (a,b,c) where + types = (a,b,c) :: (J9,J9,J9) + diff --git a/Math/Test/TestAll.hs b/Math/Test/TestAll.hs --- a/Math/Test/TestAll.hs +++ b/Math/Test/TestAll.hs @@ -5,12 +5,23 @@ import Math.Test.TPermutationGroup import Math.Test.TSubquotients import Math.Test.TFiniteGeometry -import Math.Test.TCommutativeAlgebra import Math.Test.TNonCommutativeAlgebra import Math.Test.TField import Math.Test.TRootSystem +import Math.Test.TCore.TField + +import Math.Test.TCombinatorics.TDigraph +import Math.Test.TCombinatorics.TIncidenceAlgebra +import Math.Test.TCombinatorics.TMatroid +import Math.Test.TCombinatorics.TPoset +import Math.Test.TCommutativeAlgebra.TGroebnerBasis +import Math.Test.TProjects.TMiniquaternionGeometry + + + import Test.QuickCheck +import Test.HUnit testall = and [Math.Test.TGraph.test @@ -18,13 +29,22 @@ ,Math.Test.TPermutationGroup.test ,Math.Test.TSubquotients.test ,Math.Test.TFiniteGeometry.test - ,Math.Test.TCommutativeAlgebra.test ,Math.Test.TField.test ,Math.Test.TRootSystem.test ] quickCheckAll = do - quickCheck prop_CommRingMPoly quickCheck prop_NonCommRingNPoly quickCheck prop_GroupPerm + quickCheckField + quickCheck prop_NearFieldF9 + quickCheck prop_NearFieldJ9 + +hunitAll = runTestTT $ TestList [ + testlistDigraph, + testlistIncidenceAlgebra, + testlistMatroid, + testlistPoset, + testlistGroebnerBasis + ]