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