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