```-- 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]
```