Monomial.hs

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