{-# 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)]