-- 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 (Monomial v -> Monomial v -> Bool
(Monomial v -> Monomial v -> Bool)
-> (Monomial v -> Monomial v -> Bool) -> Eq (Monomial v)
forall v. Eq v => Monomial v -> Monomial v -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Monomial v -> Monomial v -> Bool
$c/= :: forall v. Eq v => Monomial v -> Monomial v -> Bool
== :: Monomial v -> Monomial v -> Bool
$c== :: forall v. Eq v => Monomial v -> Monomial v -> Bool
Eq)

instance Ord v => Ord (Monomial v) where
    compare :: Monomial v -> Monomial v -> Ordering
compare (M xs :: [v]
xs) (M ys :: [v]
ys) = (Int, [v]) -> (Int, [v]) -> Ordering
forall a. Ord a => a -> a -> Ordering
compare ([v] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [v]
xs,[v]
xs) ([v] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [v]
ys,[v]
ys)
-- Glex ordering


instance (Eq v, Show v) => Show (Monomial v) where
    show :: Monomial v -> String
show (M xs :: [v]
xs) | [v] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [v]
xs = "1"
                | Bool
otherwise = ([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]
xs)
        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)
-- Taken from NonComMonomial - why don't we just use it directly


instance (Eq v, Show v) => Num (Monomial v) where
    M xs :: [v]
xs * :: Monomial v -> Monomial v -> Monomial v
* M ys :: [v]
ys = [v] -> Monomial v
forall v. [v] -> Monomial v
M ([v]
xs [v] -> [v] -> [v]
forall a. [a] -> [a] -> [a]
++ [v]
ys)
    fromInteger :: Integer -> Monomial v
fromInteger 1 = [v] -> Monomial v
forall v. [v] -> Monomial v
M []

-- try to find l, r such that a = lbr

divM :: Monomial v -> Monomial v -> Maybe (Monomial v, Monomial v)
divM (M a :: [v]
a) (M b :: [v]
b) = [v] -> [v] -> Maybe (Monomial v, Monomial v)
divM' [] [v]
a where
    divM' :: [v] -> [v] -> Maybe (Monomial v, Monomial 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 (Monomial v, Monomial v) -> Maybe (Monomial v, Monomial v)
forall a. a -> Maybe a
Just ([v] -> Monomial v
forall v. [v] -> Monomial v
M ([v] -> Monomial v) -> [v] -> Monomial v
forall a b. (a -> b) -> a -> b
$ [v] -> [v]
forall a. [a] -> [a]
reverse [v]
ls, [v] -> Monomial v
forall v. [v] -> Monomial v
M ([v] -> Monomial v) -> [v] -> Monomial 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 (Monomial v, Monomial v)
divM' (v
rv -> [v] -> [v]
forall a. a -> [a] -> [a]
:[v]
ls) [v]
rs
    divM' _ [] = Maybe (Monomial v, Monomial v)
forall a. Maybe a
Nothing


-- (NON-COMMUTATIVE) POLYNOMIALS


newtype NPoly r v = NP [(Monomial v,r)] deriving (NPoly r v -> NPoly r v -> Bool
(NPoly r v -> NPoly r v -> Bool)
-> (NPoly r v -> NPoly r v -> Bool) -> Eq (NPoly r v)
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
forall r v. (Eq v, Eq r) => NPoly r v -> NPoly r v -> Bool
/= :: NPoly r v -> NPoly r v -> Bool
$c/= :: forall r v. (Eq v, Eq r) => NPoly r v -> NPoly r v -> Bool
== :: NPoly r v -> NPoly r v -> Bool
$c== :: forall r v. (Eq v, Eq r) => NPoly r v -> NPoly r v -> Bool
Eq)

instance (Ord r, Ord v) => Ord (NPoly r v) where
    compare :: NPoly r v -> NPoly r v -> Ordering
compare (NP ts :: [(Monomial v, r)]
ts) (NP us :: [(Monomial v, r)]
us) = [(Monomial v, r)] -> [(Monomial v, r)] -> Ordering
forall a. Ord a => a -> a -> Ordering
compare [(Monomial v, r)]
ts [(Monomial v, r)]
us

instance (Show r, Eq v, Show v) => Show (NPoly r v) where
    show :: NPoly r v -> String
show (NP []) = "0"
    show (NP ts :: [(Monomial v, r)]
ts) =
        let (c :: Char
c:cs :: String
cs) = ((Monomial v, r) -> String) -> [(Monomial v, r)] -> String
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (Monomial v, r) -> String
forall a a. (Show a, Show a, Eq a, Num a) => (a, a) -> String
showTerm [(Monomial v, r)]
ts
        in if Char
c Char -> Char -> Bool
forall a. Eq a => a -> a -> Bool
== '+' then String
cs else Char
cChar -> ShowS
forall a. a -> [a] -> [a]
:String
cs
        where showTerm :: (a, a) -> String
showTerm (m :: a
m,a :: a
a) =
                  case a -> String
forall a. Show a => a -> String
show a
a of
                  "1" -> "+" String -> ShowS
forall a. [a] -> [a] -> [a]
++ a -> String
forall a. Show a => a -> String
show a
m
                  "-1" -> "-" String -> ShowS
forall a. [a] -> [a] -> [a]
++ a -> String
forall a. Show a => a -> String
show a
m
                  -- cs@(x:_) -> (if x == '-' then cs else '+':cs) ++ (if m == 1 then "" else show m)

                  cs :: String
cs -> ShowS
showCoeff String
cs String -> ShowS
forall a. [a] -> [a] -> [a]
++ (if a
m a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== 1 then "" else a -> String
forall a. Show a => a -> String
show a
m)
              showCoeff :: ShowS
showCoeff (c :: Char
c:cs :: String
cs) = if (Char -> Bool) -> String -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Char -> String -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` ['+','-']) String
cs
                                 then "+(" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Char
cChar -> ShowS
forall a. a -> [a] -> [a]
:String
cs String -> ShowS
forall a. [a] -> [a] -> [a]
++ ")"
                                 else if Char
c Char -> Char -> Bool
forall a. Eq a => a -> a -> Bool
== '-' then Char
cChar -> ShowS
forall a. a -> [a] -> [a]
:String
cs else '+'Char -> ShowS
forall a. a -> [a] -> [a]
:Char
cChar -> ShowS
forall a. a -> [a] -> [a]
:String
cs

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

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


-- inputs in descending order

mergeTerms :: [(a, b)] -> [(a, b)] -> [(a, b)]
mergeTerms (t :: (a, b)
t@(g :: a
g,c :: b
c):ts :: [(a, b)]
ts) (u :: (a, b)
u@(h :: a
h,d :: b
d):us :: [(a, b)]
us) =
    case (a, b) -> (a, b) -> Ordering
forall a b b. Ord a => (a, b) -> (a, b) -> Ordering
cmpTerm (a, b)
t (a, b)
u of
    LT -> (a, b)
t (a, b) -> [(a, b)] -> [(a, b)]
forall a. a -> [a] -> [a]
: [(a, b)] -> [(a, b)] -> [(a, b)]
mergeTerms [(a, b)]
ts ((a, b)
u(a, b) -> [(a, b)] -> [(a, b)]
forall a. a -> [a] -> [a]
:[(a, b)]
us)
    GT -> (a, b)
u (a, b) -> [(a, b)] -> [(a, b)]
forall a. a -> [a] -> [a]
: [(a, b)] -> [(a, b)] -> [(a, b)]
mergeTerms ((a, b)
t(a, b) -> [(a, b)] -> [(a, b)]
forall a. a -> [a] -> [a]
:[(a, b)]
ts) [(a, b)]
us
    EQ -> if b
e b -> b -> Bool
forall a. Eq a => a -> a -> Bool
== 0 then [(a, b)] -> [(a, b)] -> [(a, b)]
mergeTerms [(a, b)]
ts [(a, b)]
us else (a
g,b
e) (a, b) -> [(a, b)] -> [(a, b)]
forall a. a -> [a] -> [a]
: [(a, b)] -> [(a, b)] -> [(a, b)]
mergeTerms [(a, b)]
ts [(a, b)]
us
    where e :: b
e = b
c b -> b -> b
forall a. Num a => a -> a -> a
+ b
d
mergeTerms ts :: [(a, b)]
ts us :: [(a, b)]
us = [(a, b)]
ts [(a, b)] -> [(a, b)] -> [(a, b)]
forall a. [a] -> [a] -> [a]
++ [(a, b)]
us -- one of them is null


collect :: [(a, a)] -> [(a, a)]
collect (t1 :: (a, a)
t1@(g :: a
g,c :: a
c):t2 :: (a, a)
t2@(h :: a
h,d :: a
d):ts :: [(a, a)]
ts)
    | a
g a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
h = [(a, a)] -> [(a, a)]
collect ([(a, a)] -> [(a, a)]) -> [(a, a)] -> [(a, a)]
forall a b. (a -> b) -> a -> b
$ (a
g,a
ca -> a -> a
forall a. Num a => a -> a -> a
+a
d)(a, a) -> [(a, a)] -> [(a, a)]
forall a. a -> [a] -> [a]
:[(a, a)]
ts
    | a
c a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== 0  = [(a, a)] -> [(a, a)]
collect ([(a, a)] -> [(a, a)]) -> [(a, a)] -> [(a, a)]
forall a b. (a -> b) -> a -> b
$ (a, a)
t2(a, a) -> [(a, a)] -> [(a, a)]
forall a. a -> [a] -> [a]
:[(a, a)]
ts
    | Bool
otherwise = (a, a)
t1 (a, a) -> [(a, a)] -> [(a, a)]
forall a. a -> [a] -> [a]
: [(a, a)] -> [(a, a)]
collect ((a, a)
t2(a, a) -> [(a, a)] -> [(a, a)]
forall a. a -> [a] -> [a]
:[(a, a)]
ts)
collect ts :: [(a, a)]
ts = [(a, a)]
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 :: NPoly k v -> NPoly k v
recip (NP [(1,c :: k
c)]) = [(Monomial v, k)] -> NPoly k v
forall r v. [(Monomial v, r)] -> NPoly r v
NP [(1, k -> k
forall a. Fractional a => a -> a
recip k
c)]
    recip _ = String -> NPoly k v
forall a. HasCallStack => String -> a
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 (Var -> Var -> Bool
(Var -> Var -> Bool) -> (Var -> Var -> Bool) -> Eq Var
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Var -> Var -> Bool
$c/= :: Var -> Var -> Bool
== :: Var -> Var -> Bool
$c== :: Var -> Var -> Bool
Eq,Eq Var
Eq Var =>
(Var -> Var -> Ordering)
-> (Var -> Var -> Bool)
-> (Var -> Var -> Bool)
-> (Var -> Var -> Bool)
-> (Var -> Var -> Bool)
-> (Var -> Var -> Var)
-> (Var -> Var -> Var)
-> Ord Var
Var -> Var -> Bool
Var -> Var -> Ordering
Var -> Var -> Var
forall a.
Eq a =>
(a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: Var -> Var -> Var
$cmin :: Var -> Var -> Var
max :: Var -> Var -> Var
$cmax :: Var -> Var -> Var
>= :: Var -> Var -> Bool
$c>= :: Var -> Var -> Bool
> :: Var -> Var -> Bool
$c> :: Var -> Var -> Bool
<= :: Var -> Var -> Bool
$c<= :: Var -> Var -> Bool
< :: Var -> Var -> Bool
$c< :: Var -> Var -> Bool
compare :: Var -> Var -> Ordering
$ccompare :: Var -> Var -> Ordering
$cp1Ord :: Eq Var
Ord)

instance Show Var where
    show :: Var -> String
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 -> NPoly k v
var v :: v
v = [(Monomial v, k)] -> NPoly k v
forall r v. [(Monomial v, r)] -> NPoly r v
NP [([v] -> Monomial v
forall v. [v] -> Monomial v
M [v
v], 1)]

x :: NPoly Q Var
x = Var -> NPoly Q Var
forall k v. Num k => v -> NPoly k v
var Var
X :: NPoly Q Var
y :: NPoly Q Var
y = Var -> NPoly Q Var
forall k v. Num k => v -> NPoly k v
var Var
Y :: NPoly Q Var
z :: NPoly Q Var
z = Var -> NPoly Q Var
forall k v. Num k => v -> NPoly k v
var Var
Z :: NPoly Q Var


-- DIVISION ALGORITHM


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

        in NPoly r v
-> ([(NPoly r v, NPoly r v)], NPoly r v)
-> ([(NPoly r v, NPoly r v)], NPoly r v)
quotRemNP' (NPoly r v
hNPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
-NPoly r v
lth) ([(NPoly r v, NPoly r v)] -> [(NPoly r v, NPoly r v)]
forall a. [a] -> [a]
reverse [(NPoly r v, NPoly r v)]
lrs', NPoly r v
f'NPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
+NPoly r v
lth)

-- It is only marginally (5-10%) more space/time efficient not to track the (lazily unevaluated) factors

remNP :: NPoly r v -> [NPoly r v] -> NPoly r v
remNP f :: NPoly r v
f gs :: [NPoly r v]
gs | (NPoly r v -> Bool) -> [NPoly r v] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (NPoly r v -> NPoly r v -> Bool
forall a. Eq a => a -> a -> Bool
/=0) [NPoly r v]
gs = NPoly r v -> NPoly r v -> NPoly r v
remNP' NPoly r v
f 0
-- let result = remNP' f 0 in if result == remNP2 f gs then result else error ("remNP2 " ++ show f ++ " " ++ show gs)

           | Bool
otherwise = String -> NPoly r v
forall a. HasCallStack => String -> a
error "remNP: division by zero"
    where
    n :: Int
n = [NPoly r v] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [NPoly r v]
gs
    remNP' :: NPoly r v -> NPoly r v -> NPoly r v
remNP' 0 f' :: NPoly r v
f' = NPoly r v
f'
    remNP' h :: NPoly r v
h f' :: NPoly r v
f' = NPoly r v -> [NPoly r v] -> NPoly r v -> NPoly r v
divisionStep NPoly r v
h [NPoly r v]
gs NPoly r v
f'
    divisionStep :: NPoly r v -> [NPoly r v] -> NPoly r v -> NPoly r v
divisionStep h :: NPoly r v
h (g:gs) f' :: NPoly r v
f' =
        case NPoly r v -> Monomial v
forall r v. NPoly r v -> Monomial v
lm NPoly r v
h Monomial v -> Monomial v -> Maybe (Monomial v, Monomial v)
forall v.
Eq v =>
Monomial v -> Monomial v -> Maybe (Monomial v, Monomial v)
`divM` NPoly r v -> Monomial v
forall r v. NPoly r v -> Monomial v
lm NPoly r v
g of
        Just (l' :: Monomial v
l',r' :: Monomial v
r') -> let l'' :: NPoly r v
l'' = [(Monomial v, r)] -> NPoly r v
forall r v. [(Monomial v, r)] -> NPoly r v
NP [(Monomial v
l',NPoly r v -> r
forall r v. NPoly r v -> r
lc NPoly r v
h r -> r -> r
forall a. Fractional a => a -> a -> a
/ NPoly r v -> r
forall r v. NPoly r v -> r
lc NPoly r v
g)]
                            r'' :: NPoly r v
r'' = [(Monomial v, r)] -> NPoly r v
forall r v. [(Monomial v, r)] -> NPoly r v
NP [(Monomial v
r',1)]
                            h' :: NPoly r v
h' = NPoly r v
h NPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
- NPoly r v
l'' NPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
* NPoly r v
g NPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
* NPoly r v
r''
                        in NPoly r v -> NPoly r v -> NPoly r v
remNP' NPoly r v
h' NPoly r v
f'
        Nothing -> NPoly r v -> [NPoly r v] -> NPoly r v -> NPoly r v
divisionStep NPoly r v
h [NPoly r v]
gs NPoly r v
f'
    divisionStep h :: NPoly r v
h [] f' :: NPoly r v
f' =
        let lth :: NPoly r v
lth = NPoly r v -> NPoly r v
forall r v. NPoly r v -> NPoly r v
lt NPoly r v
h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms

        in NPoly r v -> NPoly r v -> NPoly r v
remNP' (NPoly r v
hNPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
-NPoly r v
lth) (NPoly r v
f'NPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
+NPoly r v
lth)

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

f :: NPoly r v
f %% :: NPoly r v -> [NPoly r v] -> NPoly r v
%% gs :: [NPoly r v]
gs = NPoly r v -> [NPoly r v] -> NPoly r v
forall v r.
(Fractional r, Ord v, Show v, Eq r) =>
NPoly r v -> [NPoly r v] -> NPoly r v
remNP NPoly r v
f [NPoly r v]
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 :: NPoly r v -> [NPoly r v] -> NPoly r v
remNP2 f :: NPoly r v
f gs :: [NPoly r v]
gs | (NPoly r v -> Bool) -> [NPoly r v] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (NPoly r v -> NPoly r v -> Bool
forall a. Eq a => a -> a -> Bool
/=0) [NPoly r v]
gs = NPoly r v -> NPoly r v -> NPoly r v
remNP' NPoly r v
f 0
           | Bool
otherwise = String -> NPoly r v
forall a. HasCallStack => String -> a
error "remNP: division by zero"
    where
    n :: Int
n = [NPoly r v] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [NPoly r v]
gs
    remNP' :: NPoly r v -> NPoly r v -> NPoly r v
remNP' 0 f' :: NPoly r v
f' = NPoly r v
f'
    remNP' h :: NPoly r v
h f' :: NPoly r v
f' = NPoly r v -> [NPoly r v] -> NPoly r v -> NPoly r v
divisionStep NPoly r v
h [NPoly r v]
gs NPoly r v
f'
    divisionStep :: NPoly r v -> [NPoly r v] -> NPoly r v -> NPoly r v
divisionStep h :: NPoly r v
h (g:gs) f' :: NPoly r v
f' =
        case NPoly r v -> Monomial v
forall r v. NPoly r v -> Monomial v
lm NPoly r v
h Monomial v -> Monomial v -> Maybe (Monomial v, Monomial v)
forall v.
Eq v =>
Monomial v -> Monomial v -> Maybe (Monomial v, Monomial v)
`divM` NPoly r v -> Monomial v
forall r v. NPoly r v -> Monomial v
lm NPoly r v
g of
        Just (l' :: Monomial v
l',r' :: Monomial v
r') -> let l'' :: NPoly r v
l'' = [(Monomial v, r)] -> NPoly r v
forall r v. [(Monomial v, r)] -> NPoly r v
NP [(Monomial v
l',1)] -- NP [(l',lc h / lc g)]

                            r'' :: NPoly r v
r'' = [(Monomial v, r)] -> NPoly r v
forall r v. [(Monomial v, r)] -> NPoly r v
NP [(Monomial v
r',1)]
                            lcg :: NPoly r v
lcg = r -> NPoly r v
forall r v. (Num r, Eq r, Eq v, Show v) => r -> NPoly r v
inject (NPoly r v -> r
forall r v. NPoly r v -> r
lc NPoly r v
g)
                            lch :: NPoly r v
lch = r -> NPoly r v
forall r v. (Num r, Eq r, Eq v, Show v) => r -> NPoly r v
inject (NPoly r v -> r
forall r v. NPoly r v -> r
lc NPoly r v
h)
                            -- h' = h - l'' * g * r''

                            h' :: NPoly r v
h' = NPoly r v
lcg NPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
* NPoly r v
h NPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
- NPoly r v
lch NPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
* NPoly r v
l'' NPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
* NPoly r v
g NPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
* NPoly r v
r''
                        in NPoly r v -> NPoly r v -> NPoly r v
remNP' NPoly r v
h' (NPoly r v
lcg NPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
* NPoly r v
f') -- must multiply f' by lcg too (otherwise get incorrect results, eg tlBasis 4)

        Nothing -> NPoly r v -> [NPoly r v] -> NPoly r v -> NPoly r v
divisionStep NPoly r v
h [NPoly r v]
gs NPoly r v
f'
    divisionStep h :: NPoly r v
h [] f' :: NPoly r v
f' =
        let lth :: NPoly r v
lth = NPoly r v -> NPoly r v
forall r v. NPoly r v -> NPoly r v
lt NPoly r v
h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms

        in NPoly r v -> NPoly r v -> NPoly r v
remNP' (NPoly r v
hNPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
-NPoly r v
lth) (NPoly r v
f'NPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
+NPoly r v
lth)


-- OTHER STUFF


toMonic :: NPoly r v -> NPoly r v
toMonic 0 = 0
toMonic (NP ts :: [(Monomial v, r)]
ts@((_,c :: r
c):_))
    | r
c r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== 1 = [(Monomial v, r)] -> NPoly r v
forall r v. [(Monomial v, r)] -> NPoly r v
NP [(Monomial v, r)]
ts
    | Bool
otherwise = [(Monomial v, r)] -> NPoly r v
forall r v. [(Monomial v, r)] -> NPoly r v
NP ([(Monomial v, r)] -> NPoly r v) -> [(Monomial v, r)] -> NPoly r v
forall a b. (a -> b) -> a -> b
$ ((Monomial v, r) -> (Monomial v, r))
-> [(Monomial v, r)] -> [(Monomial v, r)]
forall a b. (a -> b) -> [a] -> [b]
map (\(m :: Monomial v
m,d :: r
d)->(Monomial v
m,r
dr -> r -> r
forall a. Fractional a => a -> a -> a
/r
c)) [(Monomial v, r)]
ts

-- injection of field elements into polynomial ring

inject :: r -> NPoly r v
inject 0 = [(Monomial v, r)] -> NPoly r v
forall r v. [(Monomial v, r)] -> NPoly r v
NP []
inject c :: r
c = [(Monomial v, r)] -> NPoly r v
forall r v. [(Monomial v, r)] -> NPoly r v
NP [(Integer -> Monomial v
forall a. Num a => Integer -> a
fromInteger 1, r
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 :: [(NPoly r v, NPoly r v)] -> NPoly r v -> NPoly r v
subst vts :: [(NPoly r v, NPoly r v)]
vts (NP us :: [(Monomial v, r)]
us) = [NPoly r v] -> NPoly r v
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [r -> NPoly r v
forall r v. (Num r, Eq r, Eq v, Show v) => r -> NPoly r v
inject r
c NPoly r v -> NPoly r v -> NPoly r v
forall a. Num a => a -> a -> a
* Monomial v -> NPoly r v
substM Monomial v
m | (m :: Monomial v
m,c :: r
c) <- [(Monomial v, r)]
us] where
    substM :: Monomial v -> NPoly r v
substM (M xs :: [v]
xs) = [NPoly r v] -> NPoly r v
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [v -> NPoly r v
substV v
x | v
x <- [v]
xs]
    substV :: v -> NPoly r v
substV v :: v
v =
        let v' :: NPoly r v
v' = [(Monomial v, r)] -> NPoly r v
forall r v. [(Monomial v, r)] -> NPoly r v
NP [([v] -> Monomial v
forall v. [v] -> Monomial v
M [v
v], 1)] in
        case NPoly r v -> [(NPoly r v, NPoly r v)] -> Maybe (NPoly r v)
forall a b. Eq a => a -> [(a, b)] -> Maybe b
L.lookup NPoly r v
v' [(NPoly r v, NPoly r v)]
vts of
        Just t :: NPoly r v
t -> NPoly r v
t
        Nothing -> String -> NPoly r v
forall a. HasCallStack => String -> a
error ("subst: no substitute supplied for " String -> ShowS
forall a. [a] -> [a] -> [a]
++ NPoly r v -> String
forall a. Show a => a -> String
show NPoly r v
v')


-- INVERTIBLE

-- To support algebras which have invertible elements


class Invertible a where
    inv :: a -> a

x :: a
x ^- :: a -> b -> a
^- k :: b
k = a -> a
forall a. Invertible a => a -> a
inv a
x a -> b -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ b
k