-- Copyright (c) 2010, 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 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 (Eq) instance Ord v => Ord (NonComMonomial v) where compare (NCM lx xs) (NCM ly ys) = compare (-lx, xs) (-ly, ys) -- ie Glex ordering instance (Eq v, Show v) => Show (NonComMonomial v) where show (NCM _ []) = "1" show (NCM _ vs) = concatMap showPower (L.group vs) where showPower [v] = showVar v showPower vs@(v:_) = showVar v ++ "^" ++ show (length vs) showVar v = filter (/= '"') (show v) instance Mon (NonComMonomial v) where munit = NCM 0 [] mmult (NCM i xs) (NCM j ys) = NCM (i+j) (xs++ys) instance (Num k, Ord v) => Algebra k (NonComMonomial v) where unit 0 = zero -- V [] unit x = V [(munit,x)] mult = nf . fmap (\(a,b) -> a `mmult` 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 `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 NonComMonomial where var v = V [(NCM 1 [v],1)] powers (NCM _ vs) = map power (L.group vs) where power vs@(v:_) = (v,length 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 (NCM _ a) (NCM _ b) = divM' [] a where divM' ls (r:rs) = if b `L.isPrefixOf` (r:rs) then Just (ncm $ reverse ls, ncm $ drop (length b) (r:rs)) else divM' (r:ls) rs divM' _ [] = 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 xs = NCM (length xs) xs lm (V ((m,c):ts)) = m lc (V ((m,c):ts)) = c lt (V (t:ts)) = V [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'' = V [(l',lc h / lc g)] r'' = V [(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 | 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'' = V [(l',lc h / lc g)] r'' = V [(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