module Math.Algebra.Commutative.MPoly where
import qualified Data.Map as M
import Data.List as L
import Control.Arrow (first, second)
import Data.Ratio (denominator)
import Math.Algebra.Field.Base
import Math.Algebra.Commutative.Monomial
newtype MPoly ord r = MP [(Monomial ord,r)] deriving (Eq)
instance (Ord (Monomial ord), Ord r) => Ord (MPoly ord r) where
compare (MP ts) (MP us) = compare ts us
instance (Show r, Num r) => Show (MPoly ord r) where
show (MP []) = "0"
show (MP ts) =
let (c:cs) = concatMap showTerm ts
in if c == '+' then cs else c:cs
where showTerm (m,c) =
case show c of
"1" -> "+" ++ show m
"-1" -> "-" ++ show m
cs@(x:_) -> (if x == '-' then cs else '+':cs) ++ (if m == 1 then "" else show m)
instance (Ord (Monomial ord), Num r) => Num (MPoly ord r) where
MP ts + MP us = MP (mergeTerms ts us)
negate (MP ts) = MP $ map (second negate) ts
MP ts * MP us = MP $ collect $ sortBy cmpTerm $ [(g*h,c*d) | (g,c) <- ts, (h,d) <- us]
fromInteger 0 = MP []
fromInteger n = MP [(fromInteger 1, fromInteger n)]
cmpTerm (a,c) (b,d) = case compare a b of EQ -> EQ; GT -> LT; LT -> GT
mergeTerms (t@(g,c):ts) (u@(h,d):us) =
case compare g h of
GT -> t : mergeTerms ts (u:us)
LT -> 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
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
instance (Ord (Monomial ord), Fractional r) => Fractional (MPoly ord r) where
recip (MP [(m,c)]) = MP [(recip m, recip c)]
recip _ = error "MPoly.recip: only supported for (non-zero) constants or monomials"
var :: String -> MPoly Grevlex Q
var v = MP [(Monomial $ M.singleton v 1, 1)] :: MPoly Grevlex Q
a, b, c, d, s, t, u, v, w, x, y, z :: MPoly Grevlex Q
a = var "a"
b = var "b"
c = var "c"
d = var "d"
s = var "s"
t = var "t"
u = var "u"
v = var "v"
w = var "w"
x = var "x"
y = var "y"
z = var "z"
x_ i = var ("x" ++ show i)
x0, x1, x2, x3 :: MPoly Grevlex Q
x0 = x_ 0
x1 = x_ 1
x2 = x_ 2
x3 = x_ 3
convertMP (MP ts) = MP $ sortBy cmpTerm $ map (first convertM) ts
toLex :: MPoly ord k -> MPoly Lex k
toLex = convertMP
toGlex :: MPoly ord k -> MPoly Glex k
toGlex = convertMP
toGrevlex :: MPoly ord k -> MPoly Grevlex k
toGrevlex = convertMP
toElim :: MPoly ord k -> MPoly Elim k
toElim = convertMP
varLex v = toLex $ var v
varElim v = toElim $ var v
lt (MP (t:ts)) = t
lm = fst . lt
deg 0 = 1
deg (MP ts) = maximum [degM m | (m,c) <- ts]
mulT (m,c) (m',c') = (m*m',c*c')
divT (m,c) (m',c') = (m/m',c/c')
dividesT (m,_) (m',_) = dividesM m m'
properlyDividesT (m,_) (m',_) = dividesM m m' && m /= m'
lcmT (m,c) (m',c') = (lcmM m m',1)
infixl 7 .*
t .* MP ts = MP $ map (mulT t) ts
quotRemMP f gs = quotRemMP' f (replicate n 0, 0) where
n = length gs
quotRemMP' 0 (us,r) = (us,r)
quotRemMP' h (us,r) = divisionStep h (gs,[],us,r)
divisionStep h (g:gs,us',u:us,r) =
if lt g `dividesT` lt h
then let t = MP [lt h `divT` lt g]
h' = h t*g
u' = u+t
in quotRemMP' h' (reverse us' ++ u':us, r)
else divisionStep h (gs,u:us',us,r)
divisionStep h ([],us',[],r) =
let (lth,h') = splitlt h
in quotRemMP' h' (reverse us', r+lth)
splitlt (MP (t:ts)) = (MP [t], MP ts)
infixl 7 %%
f %% gs = r where (_,r) = quotRemMP f gs
divModMP f g = (q,r) where ([q],r) = quotRemMP f [g]
divMP f g = q where ([q],_) = quotRemMP f [g]
modMP f g = r where (_,r) = quotRemMP f [g]
inject 0 = MP []
inject c = MP [(fromInteger 1, c)]
toMonic 0 = 0
toMonic (MP ts@((_,c):_))
| c == 1 = MP ts
| otherwise = MP $ map (second (/c)) ts
toZ (MP ts) = MP $ map (second (*c)) ts where c = fromInteger $ foldl lcm 1 $ [denominator c | (m,Q c) <- ts]
subst vts (MP us) = sum [inject c * substM m | (m,c) <- us] where
substM (Monomial m) = product [substV v ^ i | (v,i) <- M.toList m]
substV v =
let v' = MP [(Monomial $ M.singleton v 1, 1)] in
case L.lookup v' vts of
Just t -> t
Nothing -> v'
support (MP ts) = [MP [(m,1)] | m <- reverse $ L.sort $ support' ts]
where support' ts = foldl L.union [] [supportM m | (m,c) <- ts]