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

{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}

-- |A module defining the algebra of commutative polynomials over a field k
module Math.Algebras.Commutative where

import Math.Algebra.Field.Base hiding (powers)
import Math.Algebras.VectorSpace
import Math.Algebras.TensorProduct
import Math.Algebras.Structures

-- GLEX MONOMIALS

data GlexMonomial v = Glex Int [(v,Int)] deriving (Eq)
-- The initial Int is the degree of the monomial. Storing it speeds up equality tests and comparisons

-- type GlexMonomialS = GlexMonomial String

instance Ord v => Ord (GlexMonomial v) where
compare (Glex si xis) (Glex sj yjs) = compare (-si, xis) (-sj, yjs)

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

{-
-- GlexMonomial is a functor and a monad
-- However, this isn't all that much use, and to make proper use of it we'd need a "nf" function
-- So leaving this commented out

-- map one basis to another
instance Functor GlexMonomial where
fmap f (Glex si xis) = Glex si [(f x, i) | (x,i) <- xis]
-- Note that as we can't assume the Ord instance, we would need to call "nf" afterwards

-- GlexMonomial is the free commutative monoid, and hence a monad
return x = Glex 1 [(x,1)]
(Glex _ xis) >>= f = let parts = [(i, sj, yjs) | (x,i) <- xis, let Glex sj yjs = f x]
in Glex (sum [i*sj | (i,sj,_) <- parts])
(concatMap (\(i,_,yjs)->map (\(y,j)->(y,i*j)) yjs) parts)
-- this isn't really much use - it's variable substitution, but we're only allowed to substitute monomials for each var
-- Note that as we can't assume the Ord instance, we would need to call "nf" afterwards
-}

-- This is the monoid algebra for commutative monomials (which are the free commutative monoid)
instance (Num k, Ord v) => Algebra k (GlexMonomial v) where
unit 0 = zero -- V []
unit x = V [(munit,x)] where munit = Glex 0 []
mult (V ts) = nf \$ fmap (\(a,b) -> a `mmult` b) (V ts)
where mmult (Glex si xis) (Glex sj yjs) = Glex (si+sj) \$ addmerge xis yjs

{-
-- This is just the Set Coalgebra, so better to use a generic instance
-- Also, not used anywhere. Hence commented out
instance Num k => Coalgebra k (GlexMonomial v) where
counit (V ts) = sum [x | (m,x) <- ts] -- trace
comult = fmap (\m -> T m m) -- diagonal
-- comult (V ts) = V [(T m m, x) | (m, x) <- ts]
-}
type GlexPoly k v = Vect k (GlexMonomial v)

glexVar v = V [(Glex 1 [(v,1)], 1)]

class Monomial m where
var :: v -> Vect Q (m v)
powers :: m v -> [(v,Int)]

-- |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 b constraint,
-- secondly because Haskell doesn't support type functions.
bind :: (Monomial m, Num k, Ord b, Show b, Algebra k b) =>
Vect k (m v) -> (v -> Vect k b) -> Vect k b
V ts `bind` f = sum [c `smultL` product [f x ^ i | (x,i) <- powers m] | (m, c) <- ts]
-- V ts `bind` f = sum [product [f x ^ i | (x,i) <- powers m] * unit c | (m, c) <- ts]

instance Monomial GlexMonomial where
var = glexVar
powers (Glex _ xis) = xis

-- DIVISION

lt (V (t:ts)) = t

class DivisionBasis b where
dividesB :: b -> b -> Bool
divB :: b -> b -> b

dividesT (b1,x1) (b2,x2) = dividesB b1 b2
divT (b1,x1) (b2,x2) = (divB b1 b2, x1/x2)

-- 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 `dividesT` lt h
then let t = V [lt h `divT` 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)

infixl 7 %%

(%%) :: (Fractional k, Ord b, Show b, Algebra k b, DivisionBasis b)
=> Vect k b -> [Vect k b] -> Vect k b
f %% gs = r where (_,r) = quotRemMP f gs

instance Ord v => DivisionBasis (GlexMonomial v) where
dividesB (Glex si xis) (Glex sj yjs) = si <= sj && dividesB' xis yjs where
dividesB' ((x,i):xis) ((y,j):yjs) =
case compare x y of
LT -> False
GT -> dividesB' ((x,i):xis) yjs
EQ -> if i<=j then dividesB' xis yjs else False
dividesB' [] _ = True
dividesB' _ [] = False
divB (Glex si xis) (Glex sj yjs) = Glex (si-sj) \$ divB' xis yjs where
divB' ((x,i):xis) ((y,j):yjs) =
case compare x y of
LT -> (x,i) : divB' xis ((y,j):yjs)
EQ -> if i == j then divB' xis yjs else (x,i-j) : divB' xis yjs -- we don't bother to check i > j
GT -> error "divB'" -- (y,-j) : divB' ((x,i):xis) yjs
divB' xis [] = xis
divB' [] yjs = error "divB'"

{-
-- Need to thread this through Maybe properly, so perhaps use do notation
divB2 (Glex si xis) (Glex sj yjs)
| si < sj = Nothing
| otherwise = case divB' xis yjs of
Nothing -> Nothing
Just zks -> Glex (si-sj) zks
where divB' ((x,i):xis) ((y,j):yjs) =
case compare x y of
LT -> (x,i) : divB' xis ((y,j):yjs)
EQ -> case compare i j of
LT -> Nothing
EQ -> divB' xis yjs
GT -> (x,i-j) : divB' xis yjs
GT -> Nothing
-}
-- !! could change divB to return Maybe, and avoid need for dividesB

```