-- 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, [(x,-i) | (x,i) <- xis]) (-sj, [(y,-j) | (y,j) <- 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
instance Monad GlexMonomial where
    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 (Eq k, Num k, Ord v) => Algebra k (GlexMonomial v) where
    unit x = x *> return munit
        where munit = Glex 0 []
    mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy
        where mmult (Glex si xis) (Glex sj yjs) = Glex (si+sj) $ addmerge xis yjs


-- GlexPoly can be given the set coalgebra structure, which is compatible with the monoid algebra structure
instance (Eq k, Num k) => Coalgebra k (GlexMonomial v) where
    counit = unwrap . nf . fmap (\m -> () )  -- trace
    -- counit (V ts) = sum [x | (m,x) <- ts]  -- trace
    comult = fmap (\m -> (m,m) )             -- diagonal

type GlexPoly k v = Vect k (GlexMonomial v)


-- |glexVar creates a variable in the algebra of commutative polynomials with Glex term ordering.
-- For example, the following code creates variables called x, y and z:
--
-- > [x,y,z] = map glexVar ["x","y","z"] :: GlexPoly Q String
glexVar :: (Num k) => v -> GlexPoly k 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, Eq k, 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 *> product [f x ^ i | (x,i) <- powers m] | (m, c) <- ts] 
-- flipbind f = linear (\m -> product [f x ^ i | (x,i) <- powers m])


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

-- |(%%) reduces a polynomial with respect to a list of polynomials.
(%%) :: (Eq k, 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