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

{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}

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

import Prelude hiding ( (*>) )

import Math.Algebra.Field.Base hiding (powers)
import Math.Algebras.VectorSpace
import Math.Algebras.TensorProduct
import Math.Algebras.Structures
import qualified Data.List as L


data NonComMonomial v = NCM Int [v] deriving (NonComMonomial v -> NonComMonomial v -> Bool
(NonComMonomial v -> NonComMonomial v -> Bool)
-> (NonComMonomial v -> NonComMonomial v -> Bool)
-> Eq (NonComMonomial v)
forall v. Eq v => NonComMonomial v -> NonComMonomial v -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: NonComMonomial v -> NonComMonomial v -> Bool
$c/= :: forall v. Eq v => NonComMonomial v -> NonComMonomial v -> Bool
== :: NonComMonomial v -> NonComMonomial v -> Bool
$c== :: forall v. Eq v => NonComMonomial v -> NonComMonomial v -> Bool
Eq)

instance Ord v => Ord (NonComMonomial v) where
    compare :: NonComMonomial v -> NonComMonomial v -> Ordering
compare (NCM lx :: Int
lx xs :: [v]
xs) (NCM ly :: Int
ly ys :: [v]
ys) = (Int, [v]) -> (Int, [v]) -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (-Int
lx, [v]
xs) (-Int
ly, [v]
ys)
-- ie Glex ordering

instance (Eq v, Show v) => Show (NonComMonomial v) where
    show :: NonComMonomial v -> String
show (NCM _ []) = "1"
    show (NCM _ vs :: [v]
vs) = ([v] -> String) -> [[v]] -> String
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap [v] -> String
forall a. Show a => [a] -> String
showPower ([v] -> [[v]]
forall a. Eq a => [a] -> [[a]]
L.group [v]
vs)
        where showPower :: [a] -> String
showPower [v :: a
v] = a -> String
forall a. Show a => a -> String
showVar a
v
              showPower vs :: [a]
vs@(v :: a
v:_) = a -> String
forall a. Show a => a -> String
showVar a
v String -> ShowS
forall a. [a] -> [a] -> [a]
++ "^" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show ([a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
vs)
              showVar :: a -> String
showVar v :: a
v = (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
v)

instance Mon (NonComMonomial v) where
    munit :: NonComMonomial v
munit = Int -> [v] -> NonComMonomial v
forall v. Int -> [v] -> NonComMonomial v
NCM 0 []
    mmult :: NonComMonomial v -> NonComMonomial v -> NonComMonomial v
mmult (NCM i :: Int
i xs :: [v]
xs) (NCM j :: Int
j ys :: [v]
ys) = Int -> [v] -> NonComMonomial v
forall v. Int -> [v] -> NonComMonomial v
NCM (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
j) ([v]
xs[v] -> [v] -> [v]
forall a. [a] -> [a] -> [a]
++[v]
ys)

instance (Eq k, Num k, Ord v) => Algebra k (NonComMonomial v) where
    unit :: k -> Vect k (NonComMonomial v)
unit 0 = Vect k (NonComMonomial v)
forall k b. Vect k b
zerov -- V []
    unit x :: k
x = [(NonComMonomial v, k)] -> Vect k (NonComMonomial v)
forall k b. [(b, k)] -> Vect k b
V [(NonComMonomial v
forall m. Mon m => m
munit,k
x)]
    mult :: Vect k (Tensor (NonComMonomial v) (NonComMonomial v))
-> Vect k (NonComMonomial v)
mult = Vect k (NonComMonomial v) -> Vect k (NonComMonomial v)
forall k b. (Eq k, Num k, Ord b) => Vect k b -> Vect k b
nf (Vect k (NonComMonomial v) -> Vect k (NonComMonomial v))
-> (Vect k (Tensor (NonComMonomial v) (NonComMonomial v))
    -> Vect k (NonComMonomial v))
-> Vect k (Tensor (NonComMonomial v) (NonComMonomial v))
-> Vect k (NonComMonomial v)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Tensor (NonComMonomial v) (NonComMonomial v) -> NonComMonomial v)
-> Vect k (Tensor (NonComMonomial v) (NonComMonomial v))
-> Vect k (NonComMonomial v)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\(a :: NonComMonomial v
a,b :: NonComMonomial v
b) -> NonComMonomial v
a NonComMonomial v -> NonComMonomial v -> NonComMonomial v
forall m. Mon m => m -> m -> m
`mmult` NonComMonomial v
b)

{-
-- This is the monoid algebra for non-commutative monomials (which is the free monoid)
instance (Num k, Ord v) => Algebra k (NonComMonomial v) where
    unit 0 = zero -- V []
    unit x = V [(munit,x)] where munit = NCM 0 []
    mult (V ts) = nf $ fmap (\(a,b) -> a `mmult` b) (V ts)
        where mmult (NCM lu us) (NCM lv vs) = NCM (lu+lv) (us++vs)
    -- mult (V ts) = nf $ V [(a `mmult` b, x) | (T a b, x) <- ts]
-}

{-
-- 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 (NonComMonomial v) where
    counit (V ts) = sum [x | (m,x) <- ts] -- trace
    comult = fmap (\m -> (m,m))
-}


class Monomial m where
    var :: v -> Vect Q (m v)
    powers :: Eq v => m v -> [(v,Int)]
-- why do we need "powers"??

V ts :: [(m t, k)]
ts bind :: Vect k (m t) -> (t -> Vect k b) -> Vect k b
`bind` f :: t -> 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 [t -> Vect k b
f t
x Vect k b -> Int -> Vect k b
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
i | (x :: t
x,i :: Int
i) <- m t -> [(t, Int)]
forall (m :: * -> *) v. (Monomial m, Eq v) => m v -> [(v, Int)]
powers m t
m] | (m :: m t
m, c :: k
c) <- [(m t, k)]
ts] 
-- flipbind f = linear (\m -> product [f x ^ i | (x,i) <- powers m])

instance Monomial NonComMonomial where
    var :: v -> Vect Q (NonComMonomial v)
var v :: v
v = [(NonComMonomial v, Q)] -> Vect Q (NonComMonomial v)
forall k b. [(b, k)] -> Vect k b
V [(Int -> [v] -> NonComMonomial v
forall v. Int -> [v] -> NonComMonomial v
NCM 1 [v
v],1)]
    powers :: NonComMonomial v -> [(v, Int)]
powers (NCM _ vs :: [v]
vs) = ([v] -> (v, Int)) -> [[v]] -> [(v, Int)]
forall a b. (a -> b) -> [a] -> [b]
map [v] -> (v, Int)
forall a. [a] -> (a, Int)
power ([v] -> [[v]]
forall a. Eq a => [a] -> [[a]]
L.group [v]
vs)
        where power :: [a] -> (a, Int)
power vs :: [a]
vs@(v :: a
v:_) = (a
v,[a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
vs)


type NCPoly v = Vect Q (NonComMonomial v)

{-
x,y,z :: NCPoly String
x = var "x"
y = var "y"
z = var "z"
-}


-- DIVISION

class DivisionBasis m where
    divM :: m -> m -> Maybe (m,m)
    -- divM a b tries to find l, r such that a = lbr
{-
    findOverlap :: m -> m -> Maybe (m,m,m)
    -- given two monomials f g, find if possible a,b,c with f=ab g=bc
-}

instance Eq v => DivisionBasis (NonComMonomial v) where
    divM :: NonComMonomial v
-> NonComMonomial v -> Maybe (NonComMonomial v, NonComMonomial v)
divM (NCM _ a :: [v]
a) (NCM _ b :: [v]
b) = [v] -> [v] -> Maybe (NonComMonomial v, NonComMonomial v)
divM' [] [v]
a where
        divM' :: [v] -> [v] -> Maybe (NonComMonomial v, NonComMonomial v)
divM' ls :: [v]
ls (r :: v
r:rs :: [v]
rs) =
            if [v]
b [v] -> [v] -> Bool
forall a. Eq a => [a] -> [a] -> Bool
`L.isPrefixOf` (v
rv -> [v] -> [v]
forall a. a -> [a] -> [a]
:[v]
rs)
            then (NonComMonomial v, NonComMonomial v)
-> Maybe (NonComMonomial v, NonComMonomial v)
forall a. a -> Maybe a
Just ([v] -> NonComMonomial v
forall v. [v] -> NonComMonomial v
ncm ([v] -> NonComMonomial v) -> [v] -> NonComMonomial v
forall a b. (a -> b) -> a -> b
$ [v] -> [v]
forall a. [a] -> [a]
reverse [v]
ls, [v] -> NonComMonomial v
forall v. [v] -> NonComMonomial v
ncm ([v] -> NonComMonomial v) -> [v] -> NonComMonomial v
forall a b. (a -> b) -> a -> b
$ Int -> [v] -> [v]
forall a. Int -> [a] -> [a]
drop ([v] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [v]
b) (v
rv -> [v] -> [v]
forall a. a -> [a] -> [a]
:[v]
rs))
            else [v] -> [v] -> Maybe (NonComMonomial v, NonComMonomial v)
divM' (v
rv -> [v] -> [v]
forall a. a -> [a] -> [a]
:[v]
ls) [v]
rs
        divM' _ [] = Maybe (NonComMonomial v, NonComMonomial v)
forall a. Maybe a
Nothing
{-
    findOverlap (NCM _ xs) (NCM _ ys) = findOverlap' [] xs ys where
        findOverlap' as [] cs = Nothing -- (reverse as, [], cs)
        findOverlap' as (b:bs) cs =
            if (b:bs) `L.isPrefixOf` cs
            then Just (ncm $ reverse as, ncm $ b:bs, ncm $ drop (length (b:bs)) cs)
            else findOverlap' (b:as) bs cs
-}
ncm :: [v] -> NonComMonomial v
ncm xs :: [v]
xs = Int -> [v] -> NonComMonomial v
forall v. Int -> [v] -> NonComMonomial v
NCM ([v] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [v]
xs) [v]
xs

lm :: Vect k b -> b
lm (V ((m :: b
m,c :: k
c):ts :: [(b, k)]
ts)) = b
m
lc :: Vect k b -> k
lc (V ((m :: b
m,c :: k
c):ts :: [(b, k)]
ts)) = k
c
lt :: Vect k b -> Vect k b
lt (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]

-- given f, gs, find ls, rs, f' such that f = sum (zipWith3 (*) ls gs rs) + f', with f' not divisible by any g
quotRemNP :: Vect k m -> [Vect k m] -> ([(Vect k m, Vect k m)], Vect k m)
quotRemNP f :: Vect k m
f gs :: [Vect k m]
gs | (Vect k m -> Bool) -> [Vect k m] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Vect k m -> Vect k m -> Bool
forall a. Eq a => a -> a -> Bool
/=0) [Vect k m]
gs = Vect k m
-> ([(Vect k m, Vect k m)], Vect k m)
-> ([(Vect k m, Vect k m)], Vect k m)
quotRemNP' Vect k m
f (Int -> (Vect k m, Vect k m) -> [(Vect k m, Vect k m)]
forall a. Int -> a -> [a]
replicate Int
n (0,0), 0)
               | Bool
otherwise = String -> ([(Vect k m, Vect k m)], Vect k m)
forall a. HasCallStack => String -> a
error "quotRemNP: division by zero"
    where
    n :: Int
n = [Vect k m] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Vect k m]
gs
    quotRemNP' :: Vect k m
-> ([(Vect k m, Vect k m)], Vect k m)
-> ([(Vect k m, Vect k m)], Vect k m)
quotRemNP' 0 (lrs :: [(Vect k m, Vect k m)]
lrs,f' :: Vect k m
f') = ([(Vect k m, Vect k m)]
lrs,Vect k m
f')
    quotRemNP' h :: Vect k m
h (lrs :: [(Vect k m, Vect k m)]
lrs,f' :: Vect k m
f') = Vect k m
-> ([Vect k m], [(Vect k m, Vect k m)], [(Vect k m, Vect k m)],
    Vect k m)
-> ([(Vect k m, Vect k m)], Vect k m)
divisionStep Vect k m
h ([Vect k m]
gs,[],[(Vect k m, Vect k m)]
lrs,Vect k m
f')
    divisionStep :: Vect k m
-> ([Vect k m], [(Vect k m, Vect k m)], [(Vect k m, Vect k m)],
    Vect k m)
-> ([(Vect k m, Vect k m)], Vect k m)
divisionStep h :: Vect k m
h (g :: Vect k m
g:gs :: [Vect k m]
gs, lrs' :: [(Vect k m, Vect k m)]
lrs', (l :: Vect k m
l,r :: Vect k m
r):lrs :: [(Vect k m, Vect k m)]
lrs, f' :: Vect k m
f') =
        case Vect k m -> m
forall k b. Vect k b -> b
lm Vect k m
h m -> m -> Maybe (m, m)
forall m. DivisionBasis m => m -> m -> Maybe (m, m)
`divM` Vect k m -> m
forall k b. Vect k b -> b
lm Vect k m
g of
        Just (l' :: m
l',r' :: m
r') -> let l'' :: Vect k m
l'' = [(m, k)] -> Vect k m
forall k b. [(b, k)] -> Vect k b
V [(m
l',Vect k m -> k
forall k b. Vect k b -> k
lc Vect k m
h k -> k -> k
forall a. Fractional a => a -> a -> a
/ Vect k m -> k
forall k b. Vect k b -> k
lc Vect k m
g)]
                            r'' :: Vect k m
r'' = [(m, k)] -> Vect k m
forall k b. [(b, k)] -> Vect k b
V [(m
r',1)]
                            h' :: Vect k m
h' = Vect k m
h Vect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
- Vect k m
l'' Vect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
* Vect k m
g Vect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
* Vect k m
r''
                        in Vect k m
-> ([(Vect k m, Vect k m)], Vect k m)
-> ([(Vect k m, Vect k m)], Vect k m)
quotRemNP' Vect k m
h' ([(Vect k m, Vect k m)] -> [(Vect k m, Vect k m)]
forall a. [a] -> [a]
reverse [(Vect k m, Vect k m)]
lrs' [(Vect k m, Vect k m)]
-> [(Vect k m, Vect k m)] -> [(Vect k m, Vect k m)]
forall a. [a] -> [a] -> [a]
++ (Vect k m
lVect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
+Vect k m
l'',Vect k m
rVect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
+Vect k m
r'')(Vect k m, Vect k m)
-> [(Vect k m, Vect k m)] -> [(Vect k m, Vect k m)]
forall a. a -> [a] -> [a]
:[(Vect k m, Vect k m)]
lrs, Vect k m
f')
        Nothing -> Vect k m
-> ([Vect k m], [(Vect k m, Vect k m)], [(Vect k m, Vect k m)],
    Vect k m)
-> ([(Vect k m, Vect k m)], Vect k m)
divisionStep Vect k m
h ([Vect k m]
gs,(Vect k m
l,Vect k m
r)(Vect k m, Vect k m)
-> [(Vect k m, Vect k m)] -> [(Vect k m, Vect k m)]
forall a. a -> [a] -> [a]
:[(Vect k m, Vect k m)]
lrs',[(Vect k m, Vect k m)]
lrs,Vect k m
f')
    divisionStep h :: Vect k m
h ([],lrs' :: [(Vect k m, Vect k m)]
lrs',[],f' :: Vect k m
f') =
        let lth :: Vect k m
lth = Vect k m -> Vect k m
forall k b. Vect k b -> Vect k b
lt Vect k m
h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms
        in Vect k m
-> ([(Vect k m, Vect k m)], Vect k m)
-> ([(Vect k m, Vect k m)], Vect k m)
quotRemNP' (Vect k m
hVect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
-Vect k m
lth) ([(Vect k m, Vect k m)] -> [(Vect k m, Vect k m)]
forall a. [a] -> [a]
reverse [(Vect k m, Vect k m)]
lrs', Vect k m
f'Vect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
+Vect k m
lth)

-- It is only marginally (5-10%) more space/time efficient not to track the (lazily unevaluated) factors
remNP :: Vect k m -> [Vect k m] -> Vect k m
remNP f :: Vect k m
f gs :: [Vect k m]
gs | (Vect k m -> Bool) -> [Vect k m] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Vect k m -> Vect k m -> Bool
forall a. Eq a => a -> a -> Bool
/=0) [Vect k m]
gs = Vect k m -> Vect k m -> Vect k m
remNP' Vect k m
f 0
           | Bool
otherwise = String -> Vect k m
forall a. HasCallStack => String -> a
error "remNP: division by zero"
    where
    n :: Int
n = [Vect k m] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Vect k m]
gs
    remNP' :: Vect k m -> Vect k m -> Vect k m
remNP' 0 f' :: Vect k m
f' = Vect k m
f'
    remNP' h :: Vect k m
h f' :: Vect k m
f' = Vect k m -> [Vect k m] -> Vect k m -> Vect k m
divisionStep Vect k m
h [Vect k m]
gs Vect k m
f'
    divisionStep :: Vect k m -> [Vect k m] -> Vect k m -> Vect k m
divisionStep h :: Vect k m
h (g:gs) f' :: Vect k m
f' =
        case Vect k m -> m
forall k b. Vect k b -> b
lm Vect k m
h m -> m -> Maybe (m, m)
forall m. DivisionBasis m => m -> m -> Maybe (m, m)
`divM` Vect k m -> m
forall k b. Vect k b -> b
lm Vect k m
g of
        Just (l' :: m
l',r' :: m
r') -> let l'' :: Vect k m
l'' = [(m, k)] -> Vect k m
forall k b. [(b, k)] -> Vect k b
V [(m
l',Vect k m -> k
forall k b. Vect k b -> k
lc Vect k m
h k -> k -> k
forall a. Fractional a => a -> a -> a
/ Vect k m -> k
forall k b. Vect k b -> k
lc Vect k m
g)]
                            r'' :: Vect k m
r'' = [(m, k)] -> Vect k m
forall k b. [(b, k)] -> Vect k b
V [(m
r',1)]
                            h' :: Vect k m
h' = Vect k m
h Vect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
- Vect k m
l'' Vect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
* Vect k m
g Vect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
* Vect k m
r''
                        in Vect k m -> Vect k m -> Vect k m
remNP' Vect k m
h' Vect k m
f'
        Nothing -> Vect k m -> [Vect k m] -> Vect k m -> Vect k m
divisionStep Vect k m
h [Vect k m]
gs Vect k m
f'
    divisionStep h :: Vect k m
h [] f' :: Vect k m
f' =
        let lth :: Vect k m
lth = Vect k m -> Vect k m
forall k b. Vect k b -> Vect k b
lt Vect k m
h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms
        in Vect k m -> Vect k m -> Vect k m
remNP' (Vect k m
hVect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
-Vect k m
lth) (Vect k m
f'Vect k m -> Vect k m -> Vect k m
forall a. Num a => a -> a -> a
+Vect k m
lth)

infixl 7 %%
-- f %% gs = r where (_,r) = quotRemNP f gs
f :: Vect k m
f %% :: Vect k m -> [Vect k m] -> Vect k m
%% gs :: [Vect k m]
gs = Vect k m -> [Vect k m] -> Vect k m
forall k m.
(DivisionBasis m, Fractional k, Eq k, Ord m, Show m,
 Algebra k m) =>
Vect k m -> [Vect k m] -> Vect k m
remNP Vect k m
f [Vect k m]
gs