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

    {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}

    -- |A module defining the algebra of commutative polynomials over a field k.
    --
    -- Most users should probably use Math.CommutativeAlgebra.Polynomial instead, which is basically the same thing
    -- but more fully-featured. This module will probably be deprecated at some point, but remains for now because
    -- it has a simpler implementation which may be more helpful for people wanting to understand the code.
    module Math.Algebras.Commutative where

    import Prelude hiding ( (*>) )

    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 (GlexMonomial v -> GlexMonomial v -> Bool
    (GlexMonomial v -> GlexMonomial v -> Bool)
    -> (GlexMonomial v -> GlexMonomial v -> Bool)
    -> Eq (GlexMonomial v)
    forall v. Eq v => GlexMonomial v -> GlexMonomial v -> Bool
    forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
    /= :: GlexMonomial v -> GlexMonomial v -> Bool
    $c/= :: forall v. Eq v => GlexMonomial v -> GlexMonomial v -> Bool
    == :: GlexMonomial v -> GlexMonomial v -> Bool
    $c== :: forall v. Eq v => GlexMonomial v -> GlexMonomial v -> Bool
    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 :: GlexMonomial v -> GlexMonomial v -> Ordering
    compare (Glex si :: Int
    si xis :: [(v, Int)]
    xis) (Glex sj :: Int
    sj yjs :: [(v, Int)]
    yjs) = (Int, [(v, Int)]) -> (Int, [(v, Int)]) -> Ordering
    forall a. Ord a => a -> a -> Ordering
    compare (-Int
    si, [(v
    x,-Int
    i) | (x :: v
    x,i :: Int
    i) <- [(v, Int)]
    xis]) (-Int
    sj, [(v
    y,-Int
    j) | (y :: v
    y,j :: Int
    j) <- [(v, Int)]
    yjs])

    instance Show v => Show (GlexMonomial v) where
        show :: GlexMonomial v -> String
    show (Glex _ []) = "1"
        show (Glex _ xis :: [(v, Int)]
    xis) = ((v, Int) -> String) -> [(v, Int)] -> String
    forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
    concatMap (\(x :: v
    x,i :: Int
    i) -> if Int
    iInt -> Int -> Bool
    forall a. Eq a => a -> a -> Bool
    ==1 then v -> String
    forall a. Show a => a -> String
    showVar v
    x else v -> String
    forall a. Show a => a -> String
    showVar v
    x String -> ShowS
    forall a. [a] -> [a] -> [a]
    ++ "^" String -> ShowS
    forall a. [a] -> [a] -> [a]
    ++ Int -> String
    forall a. Show a => a -> String
    show Int
    i) [(v, Int)]
    xis
            where showVar :: a -> String
    showVar x :: a
    x = (Char -> Bool) -> ShowS
    forall a. (a -> Bool) -> [a] -> [a]
    filter ( Char -> Char -> Bool
    forall a. Eq a => a -> a -> Bool
    /= '"' ) (a -> String
    forall a. Show a => a -> String
    show a
    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

    instance Applicative GlexMonomial where
        pure = return
        (<*>) = ap

    -- 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 :: k -> Vect k (GlexMonomial v)
    unit x :: k
    x = k
    x k -> Vect k (GlexMonomial v) -> Vect k (GlexMonomial v)
    forall k b. (Eq k, Num k) => k -> Vect k b -> Vect k b
    *> GlexMonomial v -> Vect k (GlexMonomial v)
    forall (m :: * -> *) a. Monad m => a -> m a
    return GlexMonomial v
    forall v. GlexMonomial v
    munit
            where munit :: GlexMonomial v
    munit = Int -> [(v, Int)] -> GlexMonomial v
    forall v. Int -> [(v, Int)] -> GlexMonomial v
    Glex 0 []
        mult :: Vect k (Tensor (GlexMonomial v) (GlexMonomial v))
    -> Vect k (GlexMonomial v)
    mult xy :: Vect k (Tensor (GlexMonomial v) (GlexMonomial v))
    xy = Vect k (GlexMonomial v) -> Vect k (GlexMonomial v)
    forall k b. (Eq k, Num k, Ord b) => Vect k b -> Vect k b
    nf (Vect k (GlexMonomial v) -> Vect k (GlexMonomial v))
    -> Vect k (GlexMonomial v) -> Vect k (GlexMonomial v)
    forall a b. (a -> b) -> a -> b
    $ (Tensor (GlexMonomial v) (GlexMonomial v) -> GlexMonomial v)
    -> Vect k (Tensor (GlexMonomial v) (GlexMonomial v))
    -> Vect k (GlexMonomial v)
    forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
    fmap (\(a :: GlexMonomial v
    a,b :: GlexMonomial v
    b) -> GlexMonomial v
    a GlexMonomial v -> GlexMonomial v -> GlexMonomial v
    forall v.
    Ord v =>
    GlexMonomial v -> GlexMonomial v -> GlexMonomial v
    `mmult` GlexMonomial v
    b) Vect k (Tensor (GlexMonomial v) (GlexMonomial v))
    xy
            where mmult :: GlexMonomial v -> GlexMonomial v -> GlexMonomial v
    mmult (Glex si :: Int
    si xis :: [(v, Int)]
    xis) (Glex sj :: Int
    sj yjs :: [(v, Int)]
    yjs) = Int -> [(v, Int)] -> GlexMonomial v
    forall v. Int -> [(v, Int)] -> GlexMonomial v
    Glex (Int
    siInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    sj) ([(v, Int)] -> GlexMonomial v) -> [(v, Int)] -> GlexMonomial v
    forall a b. (a -> b) -> a -> b
    $ [(v, Int)] -> [(v, Int)] -> [(v, Int)]
    forall a b.
    (Ord a, Num b, Eq b) =>
    [(a, b)] -> [(a, b)] -> [(a, b)]
    addmerge [(v, Int)]
    xis [(v, Int)]
    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 :: Vect k (GlexMonomial v) -> k
    counit = Vect k () -> k
    forall k. Num k => Vect k () -> k
    unwrap (Vect k () -> k)
    -> (Vect k (GlexMonomial v) -> Vect k ())
    -> Vect k (GlexMonomial v)
    -> k
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vect k () -> Vect k ()
    forall k b. (Eq k, Num k, Ord b) => Vect k b -> Vect k b
    nf (Vect k () -> Vect k ())
    -> (Vect k (GlexMonomial v) -> Vect k ())
    -> Vect k (GlexMonomial v)
    -> Vect k ()
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (GlexMonomial v -> ()) -> Vect k (GlexMonomial v) -> Vect k ()
    forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
    fmap (\m :: GlexMonomial v
    m -> () )  -- trace
        -- counit (V ts) = sum [x | (m,x) <- ts]  -- trace
        comult :: Vect k (GlexMonomial v)
    -> Vect k (Tensor (GlexMonomial v) (GlexMonomial v))
    comult = (GlexMonomial v -> Tensor (GlexMonomial v) (GlexMonomial v))
    -> Vect k (GlexMonomial v)
    -> Vect k (Tensor (GlexMonomial v) (GlexMonomial v))
    forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
    fmap (\m :: GlexMonomial v
    m -> (GlexMonomial v
    m,GlexMonomial v
    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 -> GlexPoly k v
    glexVar v :: v
    v = [(GlexMonomial v, k)] -> GlexPoly k v
    forall k b. [(b, k)] -> Vect k b
    V [(Int -> [(v, Int)] -> GlexMonomial v
    forall v. Int -> [(v, Int)] -> GlexMonomial v
    Glex 1 [(v
    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 :: [(m v, k)]
    ts bind :: Vect k (m v) -> (v -> Vect k b) -> Vect k b
    `bind` f :: v -> Vect k b
    f = [Vect k b] -> Vect k b
    forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
    sum [k
    c k -> Vect k b -> Vect k b
    forall k b. (Eq k, Num k) => k -> Vect k b -> Vect k b
    *> [Vect k b] -> Vect k b
    forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
    product [v -> Vect k b
    f v
    x Vect k b -> Int -> Vect k b
    forall a b. (Num a, Integral b) => a -> b -> a
    ^ Int
    i | (x :: v
    x,i :: Int
    i) <- m v -> [(v, Int)]
    forall (m :: * -> *) v. Monomial m => m v -> [(v, Int)]
    powers m v
    m] | (m :: m v
    m, c :: k
    c) <- [(m v, k)]
    ts]
    -- flipbind f = linear (\m -> product [f x ^ i | (x,i) <- powers m])

    instance Monomial GlexMonomial where
        var :: v -> Vect Q (GlexMonomial v)
    var = v -> Vect Q (GlexMonomial v)
    forall k v. Num k => v -> GlexPoly k v
    glexVar
        powers :: GlexMonomial v -> [(v, Int)]
    powers (Glex _ xis :: [(v, Int)]
    xis) = [(v, Int)]
    xis

    -- DIVISION

    lt :: Vect k b -> (b, k)
    lt (V (t :: (b, k)
    t:ts :: [(b, k)]
    ts)) = (b, k)
    t

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

    dividesT :: (b, b) -> (b, b) -> Bool
    dividesT (b1 :: b
    b1,x1 :: b
    x1) (b2 :: b
    b2,x2 :: b
    x2) = b -> b -> Bool
    forall b. DivisionBasis b => b -> b -> Bool
    dividesB b
    b1 b
    b2
    divT :: (a, b) -> (a, b) -> (a, b)
    divT (b1 :: a
    b1,x1 :: b
    x1) (b2 :: a
    b2,x2 :: b
    x2) = (a -> a -> a
    forall b. DivisionBasis b => b -> b -> b
    divB a
    b1 a
    b2, b
    x1b -> b -> b
    forall a. Fractional a => a -> a -> a
    /b
    x2)

    -- given f, gs, find as, r such that f = sum (zipWith (*) as gs) + r, with r not divisible by any g
    quotRemMP :: Vect b b -> [Vect b b] -> ([Vect b b], Vect b b)
    quotRemMP f :: Vect b b
    f gs :: [Vect b b]
    gs = Vect b b -> ([Vect b b], Vect b b) -> ([Vect b b], Vect b b)
    quotRemMP' Vect b b
    f (Int -> Vect b b -> [Vect b b]
    forall a. Int -> a -> [a]
    replicate Int
    n 0, 0) where
        n :: Int
    n = [Vect b b] -> Int
    forall (t :: * -> *) a. Foldable t => t a -> Int
    length [Vect b b]
    gs
        quotRemMP' :: Vect b b -> ([Vect b b], Vect b b) -> ([Vect b b], Vect b b)
    quotRemMP' 0 (us :: [Vect b b]
    us,r :: Vect b b
    r) = ([Vect b b]
    us,Vect b b
    r)
        quotRemMP' h :: Vect b b
    h (us :: [Vect b b]
    us,r :: Vect b b
    r) = Vect b b
    -> ([Vect b b], [Vect b b], [Vect b b], Vect b b)
    -> ([Vect b b], Vect b b)
    divisionStep Vect b b
    h ([Vect b b]
    gs,[],[Vect b b]
    us,Vect b b
    r)
        divisionStep :: Vect b b
    -> ([Vect b b], [Vect b b], [Vect b b], Vect b b)
    -> ([Vect b b], Vect b b)
    divisionStep h :: Vect b b
    h (g :: Vect b b
    g:gs :: [Vect b b]
    gs,us' :: [Vect b b]
    us',u :: Vect b b
    u:us :: [Vect b b]
    us,r :: Vect b b
    r) =
            if Vect b b -> (b, b)
    forall k b. Vect k b -> (b, k)
    lt Vect b b
    g (b, b) -> (b, b) -> Bool
    forall b b b. DivisionBasis b => (b, b) -> (b, b) -> Bool
    `dividesT` Vect b b -> (b, b)
    forall k b. Vect k b -> (b, k)
    lt Vect b b
    h
            then let t :: Vect b b
    t = [(b, b)] -> Vect b b
    forall k b. [(b, k)] -> Vect k b
    V [Vect b b -> (b, b)
    forall k b. Vect k b -> (b, k)
    lt Vect b b
    h (b, b) -> (b, b) -> (b, b)
    forall a b.
    (DivisionBasis a, Fractional b) =>
    (a, b) -> (a, b) -> (a, b)
    `divT` Vect b b -> (b, b)
    forall k b. Vect k b -> (b, k)
    lt Vect b b
    g]
                     h' :: Vect b b
    h' = Vect b b
    h Vect b b -> Vect b b -> Vect b b
    forall a. Num a => a -> a -> a
    - Vect b b
    tVect b b -> Vect b b -> Vect b b
    forall a. Num a => a -> a -> a
    *Vect b b
    g
                     u' :: Vect b b
    u' = Vect b b
    uVect b b -> Vect b b -> Vect b b
    forall a. Num a => a -> a -> a
    +Vect b b
    t
                 in Vect b b -> ([Vect b b], Vect b b) -> ([Vect b b], Vect b b)
    quotRemMP' Vect b b
    h' ([Vect b b] -> [Vect b b]
    forall a. [a] -> [a]
    reverse [Vect b b]
    us' [Vect b b] -> [Vect b b] -> [Vect b b]
    forall a. [a] -> [a] -> [a]
    ++ Vect b b
    u'Vect b b -> [Vect b b] -> [Vect b b]
    forall a. a -> [a] -> [a]
    :[Vect b b]
    us, Vect b b
    r)
            else Vect b b
    -> ([Vect b b], [Vect b b], [Vect b b], Vect b b)
    -> ([Vect b b], Vect b b)
    divisionStep Vect b b
    h ([Vect b b]
    gs,Vect b b
    uVect b b -> [Vect b b] -> [Vect b b]
    forall a. a -> [a] -> [a]
    :[Vect b b]
    us',[Vect b b]
    us,Vect b b
    r)
        divisionStep h :: Vect b b
    h ([],us' :: [Vect b b]
    us',[],r :: Vect b b
    r) =
            let (lth :: Vect b b
    lth,h' :: Vect b b
    h') = Vect b b -> (Vect b b, Vect b b)
    forall k b. Vect k b -> (Vect k b, Vect k b)
    splitlt Vect b b
    h
            in Vect b b -> ([Vect b b], Vect b b) -> ([Vect b b], Vect b b)
    quotRemMP' Vect b b
    h' ([Vect b b] -> [Vect b b]
    forall a. [a] -> [a]
    reverse [Vect b b]
    us', Vect b b
    rVect b b -> Vect b b -> Vect b b
    forall a. Num a => a -> a -> a
    +Vect b b
    lth)
        splitlt :: Vect k b -> (Vect k b, Vect k b)
    splitlt (V (t :: (b, k)
    t:ts :: [(b, k)]
    ts)) = ([(b, k)] -> Vect k b
    forall k b. [(b, k)] -> Vect k b
    V [(b, k)
    t], [(b, k)] -> Vect k b
    forall k b. [(b, k)] -> Vect k b
    V [(b, k)]
    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 :: Vect k b
    f %% :: Vect k b -> [Vect k b] -> Vect k b
    %% gs :: [Vect k b]
    gs = Vect k b
    r where (_,r :: Vect k b
    r) = Vect k b -> [Vect k b] -> ([Vect k b], Vect k b)
    forall b b.
    (DivisionBasis b, Fractional b, Eq b, Ord b, Show b,
     Algebra b b) =>
    Vect b b -> [Vect b b] -> ([Vect b b], Vect b b)
    quotRemMP Vect k b
    f [Vect k b]
    gs

    instance Ord v => DivisionBasis (GlexMonomial v) where
        dividesB :: GlexMonomial v -> GlexMonomial v -> Bool
    dividesB (Glex si :: Int
    si xis :: [(v, Int)]
    xis) (Glex sj :: Int
    sj yjs :: [(v, Int)]
    yjs) = Int
    si Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    <= Int
    sj Bool -> Bool -> Bool
    && [(v, Int)] -> [(v, Int)] -> Bool
    forall a a. (Ord a, Ord a) => [(a, a)] -> [(a, a)] -> Bool
    dividesB' [(v, Int)]
    xis [(v, Int)]
    yjs where
            dividesB' :: [(a, a)] -> [(a, a)] -> Bool
    dividesB' ((x :: a
    x,i :: a
    i):xis :: [(a, a)]
    xis) ((y :: a
    y,j :: a
    j):yjs :: [(a, a)]
    yjs) =
                case a -> a -> Ordering
    forall a. Ord a => a -> a -> Ordering
    compare a
    x a
    y of
                LT -> Bool
    False
                GT -> [(a, a)] -> [(a, a)] -> Bool
    dividesB' ((a
    x,a
    i)(a, a) -> [(a, a)] -> [(a, a)]
    forall a. a -> [a] -> [a]
    :[(a, a)]
    xis) [(a, a)]
    yjs
                EQ -> if a
    ia -> a -> Bool
    forall a. Ord a => a -> a -> Bool
    <=a
    j then [(a, a)] -> [(a, a)] -> Bool
    dividesB' [(a, a)]
    xis [(a, a)]
    yjs else Bool
    False
            dividesB' [] _ = Bool
    True
            dividesB' _ [] = Bool
    False
        divB :: GlexMonomial v -> GlexMonomial v -> GlexMonomial v
    divB (Glex si :: Int
    si xis :: [(v, Int)]
    xis) (Glex sj :: Int
    sj yjs :: [(v, Int)]
    yjs) = Int -> [(v, Int)] -> GlexMonomial v
    forall v. Int -> [(v, Int)] -> GlexMonomial v
    Glex (Int
    siInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    sj) ([(v, Int)] -> GlexMonomial v) -> [(v, Int)] -> GlexMonomial v
    forall a b. (a -> b) -> a -> b
    $ [(v, Int)] -> [(v, Int)] -> [(v, Int)]
    forall a b.
    (Ord a, Num b, Eq b) =>
    [(a, b)] -> [(a, b)] -> [(a, b)]
    divB' [(v, Int)]
    xis [(v, Int)]
    yjs where
            divB' :: [(a, b)] -> [(a, b)] -> [(a, b)]
    divB' ((x :: a
    x,i :: b
    i):xis :: [(a, b)]
    xis) ((y :: a
    y,j :: b
    j):yjs :: [(a, b)]
    yjs) =
                case a -> a -> Ordering
    forall a. Ord a => a -> a -> Ordering
    compare a
    x a
    y of
                LT -> (a
    x,b
    i) (a, b) -> [(a, b)] -> [(a, b)]
    forall a. a -> [a] -> [a]
    : [(a, b)] -> [(a, b)] -> [(a, b)]
    divB' [(a, b)]
    xis ((a
    y,b
    j)(a, b) -> [(a, b)] -> [(a, b)]
    forall a. a -> [a] -> [a]
    :[(a, b)]
    yjs)
                EQ -> if b
    i b -> b -> Bool
    forall a. Eq a => a -> a -> Bool
    == b
    j then [(a, b)] -> [(a, b)] -> [(a, b)]
    divB' [(a, b)]
    xis [(a, b)]
    yjs else (a
    x,b
    ib -> b -> b
    forall a. Num a => a -> a -> a
    -b
    j) (a, b) -> [(a, b)] -> [(a, b)]
    forall a. a -> [a] -> [a]
    : [(a, b)] -> [(a, b)] -> [(a, b)]
    divB' [(a, b)]
    xis [(a, b)]
    yjs -- we don't bother to check i > j
                GT -> String -> [(a, b)]
    forall a. HasCallStack => String -> a
    error "divB'" -- (y,-j) : divB' ((x,i):xis) yjs
            divB' xis :: [(a, b)]
    xis [] = [(a, b)]
    xis
            divB' [] yjs :: [(a, b)]
    yjs = String -> [(a, b)]
    forall a. HasCallStack => String -> a
    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