Polynomial.hs

-- Copyright (c) 2011, David Amos. All rights reserved.

{-# LANGUAGE GeneralizedNewtypeDeriving, MultiParamTypeClasses, FlexibleInstances #-}

-- |A module defining the algebra of commutative polynomials over a field k.
-- Polynomials are represented as the free k-vector space with the monomials as basis.
-- 
-- 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.
module Math.CommutativeAlgebra.Polynomial where

import Math.Core.Field
import Math.Algebras.VectorSpace
import Math.Algebras.TensorProduct
import Math.Algebras.Structures

-- |In order to work with monomials, we need to be able to multiply them and divide them.
-- Multiplication is defined by the Mon (monoid) class. Division is defined in this class.
-- The functions here are primarily intended for internal use only.
class (Eq m, Show m, Mon m) => Monomial m where
    -- mdivmaybe :: m -> m -> Maybe m
    mdivides :: m -> m -> Bool
    mdiv :: m -> m -> m
    mgcd :: m -> m -> m
    mlcm :: m -> m -> m
    mcoprime :: m -> m -> Bool
    mdeg :: m -> Int

-- mlcm m1 m2 = let m = mgcd m1 m2 in mmult m1 (mdiv m2 m)

-- mdividesne m1 m2 = m1 /= m2 && mdivides m1 m2

mproperlydivides m1 m2 = m1 /= m2 && mdivides m1 m2


-- |We want to be able to construct monomials over any set of variables that we choose.
-- Although we will often use String as the type of our variables,
-- it is useful to define polymorphic types for monomials.
class MonomialConstructor m where
    mvar :: v -> m v
    mindices :: m v -> [(v,Int)]

-- |@var v@ creates a variable in the vector space of polynomials.
-- For example, if we want to work in Q[x,y,z], we might define:
--
-- > [x,y,z] = map var ["x","y","z"] :: [GlexPoly Q String]
--
-- Notice that, in general, it is necessary to provide a type annotation so that
-- the compiler knows which field and which term order to use.
var :: (Num k, MonomialConstructor m) => v -> Vect k (m v)
var = return . mvar

-- class MonomialOrder m where
--     isGraded :: m -> Bool


-- MONOMIALS

-- |The underlying implementation of monomials in variables of type v. Most often, we will be interested in MonImpl String,
-- with the variable \"x\" represented by M 1 [(\"x\",1)]. However, other types can be used instead.
--
-- No Ord instance is defined for MonImpl v, so it cannot be used as the basis for a free vector space of polynomials.
-- Instead, several different newtype wrappers are provided, corresponding to different monomial orderings.
data MonImpl v = M Int [(v,Int)] deriving (Eq)
-- The initial Int is the degree of the monomial. Storing it speeds up equality tests and comparisons

instance Show v => Show (MonImpl v) where
    show (M _ []) = "1"
    show (M _ 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

instance (Ord v) => Mon (MonImpl v) where
    munit = M 0 []
    mmult (M si xis) (M sj yjs) = M (si+sj) $ addmerge xis yjs

instance (Ord v, Show v) => Monomial (MonImpl v) where
{-
    mdivmaybe m1 m2 =
        let m@(M s xis) = mdiv m1 m2
        in if s < 0 || any (< 0) (map snd xis)
           then Nothing
           else Just m
-}
    mdivides (M si xis) (M sj yjs) = si <= sj && mdivides' xis yjs where
        mdivides' ((x,i):xis) ((y,j):yjs) =
            case compare x y of
            LT -> False
            GT -> mdivides' ((x,i):xis) yjs
            EQ -> if i<=j then mdivides' xis yjs else False
        mdivides' [] _ = True
        mdivides' _ [] = False
    mdiv (M si xis) (M sj yjs) = M (si-sj) $ addmerge xis $ map (\(y,j) -> (y,-j)) yjs
    -- we don't check that the result has no negative indices
    mgcd (M _ xis) (M _ yjs) = mgcd' 0 [] xis yjs
        where mgcd' s zks ((x,i):xis) ((y,j):yjs) =
                  case compare x y of
                  LT -> mgcd' s zks xis ((y,j):yjs)
                  GT -> mgcd' s zks ((x,i):xis) yjs
                  EQ -> let k = min i j in mgcd' (s+k) ((x,k):zks) xis yjs
              mgcd' s zks _ _ = M s (reverse zks)
    mlcm (M si xis) (M sj yjs) = mlcm' 0 [] xis yjs
        where mlcm' s zks ((x,i):xis) ((y,j):yjs) =
                  case compare x y of
                  LT -> mlcm' (s+i) ((x,i):zks) xis ((y,j):yjs)
                  GT -> mlcm' (s+j) ((y,j):zks) ((x,i):xis) yjs
                  EQ -> let k = max i j in mlcm' (s+k) ((x,k):zks) xis yjs
              mlcm' s zks xis yjs = let zks' = xis ++ yjs; s' = sum (map snd zks') -- either xis or yjs is null
                                    in M (s+s') (reverse zks ++ zks')
    mcoprime (M _ xis) (M _ yjs) = mcoprime' xis yjs
        where mcoprime' ((x,i):xis) ((y,j):yjs) =
                  case compare x y of
                  LT -> mcoprime' xis ((y,j):yjs)
                  GT -> mcoprime' ((x,i):xis) yjs
                  EQ -> False
              mcoprime' _ _ = True
    -- mcoprime m1 m2 = mgcd m1 m2 == munit
    mdeg (M s _) = s

instance MonomialConstructor MonImpl where
    mvar v = M 1 [(v,1)]
    mindices (M si xis) = xis


-- LEX ORDER

-- |A type representing monomials with Lex ordering.
--
-- Lex stands for lexicographic ordering.
-- For example, in Lex ordering, monomials up to degree two would be ordered as follows: x^2+xy+xz+x+y^2+yz+y+z^2+z+1.
newtype Lex v = Lex (MonImpl v) deriving (Eq, Mon, Monomial, MonomialConstructor) -- GeneralizedNewtypeDeriving

instance Show v => Show (Lex v) where
    show (Lex m) = show m

instance Ord v => Ord (Lex v) where
    compare (Lex (M si xis)) (Lex (M sj yjs)) = compare' xis yjs
        where compare' ((x,i):xis) ((y,j):yjs) =
                  case compare x y of
                  LT -> LT
                  GT -> GT
                  EQ -> case compare i j of
                        LT -> GT
                        GT -> LT
                        EQ -> compare' xis yjs
              compare' [] [] = EQ
              compare' _ [] = LT
              compare' [] _ = GT
        -- unfortunately we can't use the following, because we want [] sorted after everything, not before
        -- compare [(x,-i) | (x,i) <- xis] [(y,-j) | (y,j) <- yjs]

-- instance MonomialOrder Lex where isGraded _ = False

-- |A type representing polynomials with Lex term ordering.
type LexPoly k v = Vect k (Lex v)

-- |@lexvar v@ creates a variable in the algebra of commutative polynomials over Q with Lex term ordering.
-- It is provided as a shortcut, to avoid having to provide a type annotation, as with @var@.
-- For example, the following code creates variables called x, y and z:
--
-- > [x,y,z] = map lexvar ["x","y","z"]
lexvar :: v -> LexPoly Q v
lexvar v = return $ Lex $ M 1 [(v,1)]

instance (Num k, Ord v, Show v) => Algebra k (Lex v) where
    unit x = x *> return munit
    mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy


-- GLEX ORDER

-- |A type representing monomials with Glex ordering.
--
-- Glex stands for graded lexicographic. Thus monomials are ordered first by degree, then by lexicographic order.
-- For example, in Glex ordering, monomials up to degree two would be ordered as follows: x^2+xy+xz+y^2+yz+z^2+x+y+z+1.
newtype Glex v = Glex (MonImpl v) deriving (Eq, Mon, Monomial, MonomialConstructor) -- GeneralizedNewtypeDeriving

instance Show v => Show (Glex v) where
    show (Glex m) = show m

instance Ord v => Ord (Glex v) where
    compare (Glex (M si xis)) (Glex (M sj yjs)) =
        compare (-si, [(x,-i) | (x,i) <- xis]) (-sj, [(y,-j) | (y,j) <- yjs])

-- instance MonomialOrder Glex where isGraded _ = True

-- |A type representing polynomials with Glex term ordering.
type GlexPoly k v = Vect k (Glex v)

-- |@glexvar v@ creates a variable in the algebra of commutative polynomials over Q with Glex term ordering.
-- It is provided as a shortcut, to avoid having to provide a type annotation, as with @var@.
-- For example, the following code creates variables called x, y and z:
--
-- > [x,y,z] = map glexvar ["x","y","z"]
glexvar :: v -> GlexPoly Q v
glexvar v = return $ Glex $ M 1 [(v,1)]

instance (Num k, Ord v, Show v) => Algebra k (Glex v) where
    unit x = x *> return munit
    mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy


-- GREVLEX ORDER

-- |A type representing monomials with Grevlex ordering.
--
-- Grevlex stands for graded reverse lexicographic. Thus monomials are ordered first by degree, then by reverse lexicographic order.
-- For example, in Grevlex ordering, monomials up to degree two would be ordered as follows: x^2+xy+y^2+xz+yz+z^2+x+y+z+1.
--
-- In general, Grevlex leads to the smallest Groebner bases.
newtype Grevlex v = Grevlex (MonImpl v) deriving (Eq, Mon, Monomial, MonomialConstructor) -- GeneralizedNewtypeDeriving

instance Show v => Show (Grevlex v) where
    show (Grevlex m) = show m

instance Ord v => Ord (Grevlex v) where
    compare (Grevlex (M si xis)) (Grevlex (M sj yjs)) =
        compare (-si, reverse xis) (-sj, reverse yjs)

-- instance MonomialOrder Grevlex where isGraded _ = True

-- |A type representing polynomials with Grevlex term ordering.
type GrevlexPoly k v = Vect k (Grevlex v)

-- |@grevlexvar v@ creates a variable in the algebra of commutative polynomials over Q with Grevlex term ordering.
-- It is provided as a shortcut, to avoid having to provide a type annotation, as with @var@.
-- For example, the following code creates variables called x, y and z:
--
-- > [x,y,z] = map grevlexvar ["x","y","z"]
grevlexvar :: v -> GrevlexPoly Q v
grevlexvar v = return $ Grevlex $ M 1 [(v,1)]

instance (Num k, Ord v, Show v) => Algebra k (Grevlex v) where
    unit x = x *> return munit
    mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy


-- ELIMINATION ORDER

data Elim2 a b = Elim2 !a !b deriving (Eq)

instance (Ord a, Ord b) => Ord (Elim2 a b) where
    compare (Elim2 a1 b1) (Elim2 a2 b2) = compare (a1,b1) (a2,b2)

instance (Show a, Show b) => Show (Elim2 a b) where
    show (Elim2 ma mb) = case (show ma, show mb) of
                       ("1","1") -> "1"
                       (ma',"1") -> ma'
                       ("1",mb') -> mb'
                       (ma',mb') -> ma' ++ mb'

instance (Mon a, Mon b) => Mon (Elim2 a b) where
    munit = Elim2 munit munit
    mmult (Elim2 a1 b1) (Elim2 a2 b2) = Elim2 (mmult a1 a2) (mmult b1 b2)

instance (Monomial a, Monomial b) => Monomial (Elim2 a b) where
    -- mdivmaybe :: Ord v => m v -> m v -> Maybe (m v)
    mdivides (Elim2 a1 b1) (Elim2 a2 b2) = mdivides a1 a2 && mdivides b1 b2
    mdiv (Elim2 a1 b1) (Elim2 a2 b2) = Elim2 (mdiv a1 a2) (mdiv b1 b2)
    mgcd (Elim2 a1 b1) (Elim2 a2 b2) = Elim2 (mgcd a1 a2) (mgcd b1 b2)
    mlcm (Elim2 a1 b1) (Elim2 a2 b2) = Elim2 (mlcm a1 a2) (mlcm b1 b2)
    mcoprime (Elim2 a1 b1) (Elim2 a2 b2) = mcoprime a1 a2 && mcoprime b1 b2
    mdeg (Elim2 a b) = mdeg a + mdeg b

instance (Num k, Ord a, Mon a, Ord b, Mon b) => Algebra k (Elim2 a b) where
    unit x = x *> return munit
    mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy


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

flipbind f = linear (\m -> product [f x ^ i | (x,i) <- mindices m])


-- DIVISION ALGORITHM FOR POLYNOMIALS

lt (V (t:ts)) = t -- leading term
lm = fst . lt     -- leading monomial
lc = snd . lt     -- leading coefficient

-- deg :: (Num k, Monomial m, MonomialOrder m) => Vect k m -> Int
deg (V []) = -1
deg f = maximum $ [mdeg m | (m,c) <- terms f]
{-
deg f | isGraded (lm f) = mdeg (lm f)
      | otherwise = maximum $ [mdeg m | (m,c) <- terms f]
-}
-- the true degree of the polynomial, not the degree of the leading term
-- required for sugar strategy when computing Groebner basis

toMonic 0 = 0
toMonic f = (1 / lc f) *> f

-- tdivmaybe (m1,x1) (m2,x2) = fmap (\m -> (m,x1/x2)) $ mdivmaybe m1 m2

tdivides (m1,x1) (m2,x2) = mdivides m1 m2

tdiv (m1,x1) (m2,x2) = (mdiv m1 m2, x1/x2)

tgcd (m1,_) (m2,_) = (mgcd m1 m2, 1)
-- tlcm (m1,_) (m2,_) = (mlcm m1 m2, 1)

tmult (m,c) (m',c') = (mmult m m',c*c')

infixl 7 *->
t *-> V ts = V $ map (tmult 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 `tdivides` lt h
        then let t = V [lt h `tdiv` 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)

{-
-- version using mdivmaybe - surprisingly, this seems to be slower
quotRemMP2 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) =
        case tdivmaybe (lt h) (lt g) of
        Nothing -> divisionStep h (gs,u:us',us,r)
        Just t -> let h' = h - V [t] * g
                      u' = u + V [t]
                  in quotRemMP' h' (reverse us' ++ u':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 %%

-- |@f %% gs@ is the reduction of a polynomial f with respect to a list of polynomials gs.
-- In the case where the gs are a Groebner basis for an ideal I,
-- then @f %% gs@ is the equivalence class representative of f in R/I,
-- and is zero if and only if f is in I.
(%%) :: (Fractional k, Monomial m, Ord m, Algebra k m) =>
     Vect k m -> [Vect k m] -> Vect k m
f %% gs = r where (_,r) = quotRemMP f gs


-- |As a convenience, a partial instance of Fractional is defined for polynomials.
-- The instance is well-defined only for scalars, and gives an error if used on other values.
-- The purpose of this is to allow entry of fractional scalars, in expressions such as @x/2@.
-- On the other hand, an expression such as @2/x@ will return an error.
instance (Fractional k, Monomial m, Ord m, Algebra k m) => Fractional (Vect k m) where
    recip (V [(m,c)]) | m == munit = V [(m,1/c)]
                      | otherwise = error "Polynomial recip: only defined for scalars"
    fromRational x = V [(munit, fromRational x)]