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

-- |A module providing a type for non-commutative polynomials.
module Math.Algebra.NonCommutative.NCPoly where

import Data.List as L
import Math.Algebra.Field.Base

-- (NON-COMMUTATIVE) MONOMIALS

newtype Monomial v = M [v] deriving (Eq)

instance Ord v => Ord (Monomial v) where
compare (M xs) (M ys) = compare (length xs,xs) (length ys,ys)
-- Glex ordering

instance (Eq v, Show v) => Show (Monomial v) where
show (M xs) | null xs = "1"
| otherwise = concatMap showPower (L.group xs)
where showPower [v] = showVar v
showPower vs@(v:_) = showVar v ++ "^" ++ show (length vs)
showVar v = filter (/= '"') (show v)
-- Taken from NonComMonomial - why don't we just use it directly

instance (Eq v, Show v) => Num (Monomial v) where
M xs * M ys = M (xs ++ ys)
fromInteger 1 = M []

-- try to find l, r such that a = lbr
divM (M a) (M b) = divM' [] a where
divM' ls (r:rs) =
if b `L.isPrefixOf` (r:rs)
then Just (M \$ reverse ls, M \$ drop (length b) (r:rs))
else divM' (r:ls) rs
divM' _ [] = Nothing

-- (NON-COMMUTATIVE) POLYNOMIALS

newtype NPoly r v = NP [(Monomial v,r)] deriving (Eq)

instance (Ord r, Ord v) => Ord (NPoly r v) where
compare (NP ts) (NP us) = compare ts us

instance (Show r, Eq v, Show v) => Show (NPoly r v) where
show (NP []) = "0"
show (NP ts) =
let (c:cs) = concatMap showTerm ts
in if c == '+' then cs else c:cs
where showTerm (m,a) =
case show a of
"1" -> "+" ++ show m
"-1" -> "-" ++ show m
-- cs@(x:_) -> (if x == '-' then cs else '+':cs) ++ (if m == 1 then "" else show m)
cs -> showCoeff cs ++ (if m == 1 then "" else show m)
showCoeff (c:cs) = if any (`elem` ['+','-']) cs
then "+(" ++ c:cs ++ ")"
else if c == '-' then c:cs else '+':c:cs

instance (Eq r, Num r, Ord v, Show v) => Num (NPoly r v) where
NP ts + NP us = NP (mergeTerms ts us)
negate (NP ts) = NP \$ map (\(m,c) -> (m,-c)) ts
NP ts * NP us = NP \$ collect \$ L.sortBy cmpTerm \$ [(g*h,c*d) | (g,c) <- ts, (h,d) <- us]
fromInteger 0 = NP []
fromInteger n = NP [(fromInteger 1, fromInteger n)]

cmpTerm (a,c) (b,d) = case compare a b of EQ -> EQ; GT -> LT; LT -> GT -- in mpolys we put "larger" terms first

-- inputs in descending order
mergeTerms (t@(g,c):ts) (u@(h,d):us) =
case cmpTerm t u of
LT -> t : mergeTerms ts (u:us)
GT -> u : mergeTerms (t:ts) us
EQ -> if e == 0 then mergeTerms ts us else (g,e) : mergeTerms ts us
where e = c + d
mergeTerms ts us = ts ++ us -- one of them is null

collect (t1@(g,c):t2@(h,d):ts)
| g == h = collect \$ (g,c+d):ts
| c == 0  = collect \$ t2:ts
| otherwise = t1 : collect (t2:ts)
collect ts = ts

-- Fractional instance so that we can enter fractional coefficients
-- Only lets us divide by field elements (with unit monomial), not any other polynomials
instance (Eq k, Fractional k, Ord v, Show v) => Fractional (NPoly k v) where
recip (NP [(1,c)]) = NP [(1, recip c)]
recip _ = error "NPoly.recip: only supported for (non-zero) constants"

-- SOME VARIABLES (INDETERMINATES)
-- The idea is that you define your own type of indeterminates as required, along the same lines as this

data Var = X | Y | Z deriving (Eq,Ord)

instance Show Var where
show X = "x"
show Y = "y"
show Z = "z"

-- |Create a non-commutative variable for use in forming non-commutative polynomials.
-- For example, we could define x = var "x", y = var "y". Then x*y /= y*x.
var :: (Num k) => v -> NPoly k v
var v = NP [(M [v], 1)]

x = var X :: NPoly Q Var
y = var Y :: NPoly Q Var
z = var Z :: NPoly Q Var

-- DIVISION ALGORITHM

lm (NP ((m,c):ts)) = m
lc (NP ((m,c):ts)) = c
lt (NP (t:ts)) = NP [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 f gs | all (/=0) gs = quotRemNP' f (replicate n (0,0), 0)
| otherwise = error "quotRemNP: division by zero"
where
n = length gs
quotRemNP' 0 (lrs,f') = (lrs,f')
quotRemNP' h (lrs,f') = divisionStep h (gs,[],lrs,f')
divisionStep h (g:gs, lrs', (l,r):lrs, f') =
case lm h `divM` lm g of
Just (l',r') -> let l'' = NP [(l',lc h / lc g)]
r'' = NP [(r',1)]
h' = h - l'' * g * r''
in quotRemNP' h' (reverse lrs' ++ (l+l'',r+r''):lrs, f')
Nothing -> divisionStep h (gs,(l,r):lrs',lrs,f')
divisionStep h ([],lrs',[],f') =
let lth = lt h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms
in quotRemNP' (h-lth) (reverse lrs', f'+lth)

-- It is only marginally (5-10%) more space/time efficient not to track the (lazily unevaluated) factors
remNP f gs | all (/=0) gs = remNP' f 0
-- let result = remNP' f 0 in if result == remNP2 f gs then result else error ("remNP2 " ++ show f ++ " " ++ show gs)
| otherwise = error "remNP: division by zero"
where
n = length gs
remNP' 0 f' = f'
remNP' h f' = divisionStep h gs f'
divisionStep h (g:gs) f' =
case lm h `divM` lm g of
Just (l',r') -> let l'' = NP [(l',lc h / lc g)]
r'' = NP [(r',1)]
h' = h - l'' * g * r''
in remNP' h' f'
Nothing -> divisionStep h gs f'
divisionStep h [] f' =
let lth = lt h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms
in remNP' (h-lth) (f'+lth)

infixl 7 %%
-- f %% gs = r where (_,r) = quotRemNP f gs
f %% gs = remNP f gs

-- !! Not sure if the following is valid
-- The idea is to avoid dividing by lc g, because sometimes our coefficient ring is not a field
-- Passes all the knot theory tests
-- However, it may be that if we ever get a non-invertible element at the front, we are in trouble anyway
remNP2 f gs | all (/=0) gs = remNP' f 0
| otherwise = error "remNP: division by zero"
where
n = length gs
remNP' 0 f' = f'
remNP' h f' = divisionStep h gs f'
divisionStep h (g:gs) f' =
case lm h `divM` lm g of
Just (l',r') -> let l'' = NP [(l',1)] -- NP [(l',lc h / lc g)]
r'' = NP [(r',1)]
lcg = inject (lc g)
lch = inject (lc h)
-- h' = h - l'' * g * r''
h' = lcg * h - lch * l'' * g * r''
in remNP' h' (lcg * f') -- must multiply f' by lcg too (otherwise get incorrect results, eg tlBasis 4)
Nothing -> divisionStep h gs f'
divisionStep h [] f' =
let lth = lt h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms
in remNP' (h-lth) (f'+lth)

-- OTHER STUFF

toMonic 0 = 0
toMonic (NP ts@((_,c):_))
| c == 1 = NP ts
| otherwise = NP \$ map (\(m,d)->(m,d/c)) ts

-- injection of field elements into polynomial ring
inject 0 = NP []
inject c = NP [(fromInteger 1, c)]

-- substitute terms for variables in an NPoly
-- eg subst [(x,a),(y,a+b),(z,c^2)] (x*y+z) -> a*(a+b)+c^2
subst vts (NP us) = sum [inject c * substM m | (m,c) <- us] where
substM (M xs) = product [substV x | x <- xs]
substV v =
let v' = NP [(M [v], 1)] in
case L.lookup v' vts of
Just t -> t
Nothing -> error ("subst: no substitute supplied for " ++ show v')

-- INVERTIBLE
-- To support algebras which have invertible elements

class Invertible a where
inv :: a -> a

x ^- k = inv x ^ k```