-- 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)
data Lex
data Glex
data Grevlex
data Elim -- a term order for elimination
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]