LaurentMPoly.hs

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

module Math.Projects.KnotTheory.LaurentMPoly where

import qualified Data.Map as M
import Data.List as L

import Math.Algebra.Field.Base


-- It would be possible to refactor this and other code into common code for semigroup rings
-- But there are enough small fiddly differences that it's easier not to


-- LAURENT MONOMIALS

newtype LaurentMonomial = LM (M.Map String Q) deriving (Eq)
-- We allow exponents to be rationals, because we want to support rings such as Z[q^1/2,q^-1/2] in connection with Hecke algebras

-- the multidegree - not sure how meaningful this is if we have negative indices too
degLM (LM m) = sum $ M.elems m

-- Glex ordering
instance Ord LaurentMonomial where
    compare a b = let ds = M.elems m where LM m = a/b in
                  case compare (sum ds) 0 of
                  GT -> GT
                  LT -> LT
                  EQ -> if null ds then EQ else
                        if head ds > 0 then GT else LT

instance Show LaurentMonomial where
    show (LM a) | M.null a = "1"
                | otherwise = concatMap showVar $ M.toList a
                where showVar (v,1) = v
                      showVar (v,i) = v ++ "^" ++ show i

instance Num LaurentMonomial where
    LM a * LM b = LM $ M.filter (/=0) $ M.unionWith (+) a b
    fromInteger 1 = LM M.empty

instance Fractional (LaurentMonomial) where
    recip (LM m) = LM $ M.map negate m


-- untested
-- numeratorLM (LM a) = LM $ M.filter (>0) a

denominatorLM (LM a) = recip $ LM $ M.filter (<0) a

-- not valid for arguments with negative exponents (because 0 won't trump -i)
lcmLM (LM a) (LM b) = LM $ M.unionWith max a b

-- not tested -- for arguments with non-zero denominators
-- gcdLM (LM a) (LM b) = LM $ M.intersectionWith min a b

divLM a b = let LM c = a/b in if all (>=0) (M.elems c) then Just (LM c) else Nothing


-- LAURENT POLYNOMIALS

newtype LaurentMPoly r = LP [(LaurentMonomial,r)] deriving (Eq,Ord)

instance Show r => Show (LaurentMPoly r) where
    show (LP []) = "0"
    show (LP ts) =
        let (c:cs) = concatMap showTerm (reverse ts) -- we show Laurent polys with smallest terms first
        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)
                  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
              -- we don't attempt sign reversal within brackets in case we have expressions like t^-1 inside the brackets

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

cmpTerm (a,c) (b,d) = case compare a b of EQ -> EQ; GT -> LT; LT -> GT
-- we have to put largest terms first so that quotRem works

-- 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 single terms, not any other polynomials
instance (Eq r, Fractional r) => Fractional (LaurentMPoly r) where
    recip (LP [(m,c)]) = LP [(recip m, recip c)]
    recip _ = error "LaurentMPoly.recip: only supported for (non-zero) constants or monomials"


lm (LP ((m,c):ts)) = m
lc (LP ((m,c):ts)) = c
lt (LP ((m,c):ts)) = LP [(m,c)]

quotRemLP f g
    | g == 0 = error "quotRemLP: division by zero"
    | denominatorLP f /= 1 || denominatorLP g /= 1 = error "quotRemLP: negative exponents"
    | otherwise = quotRemLP' f (0,0)
    where
    quotRemLP' 0 (q,r) = (q,r)
    quotRemLP' h (q,r) =
        case lm h `divLM` lm g of
        Just m -> let t = LP [(m, lc h / lc g)]
                  in quotRemLP' (h-t*g) (q+t,r)
        Nothing -> let lth = lt h -- can't reduce lt h, so add it to the remainder and try to reduce the remaining terms
                   in quotRemLP' (h-lth) (q, r+lth)


-- g must be a binomial without negative exponents - eg i^2+1
reduceLP f g@(LP [_,_]) =
    let fn = f * fd
        fd = denominatorLP f
        (_,rn) = quotRemLP fn g
        (_,rd) = quotRemLP fd g
    in rn / rd


var v = LP [(LM $ M.singleton v 1, 1)]

t = var "t" :: LaurentMPoly Q
x = var "x" :: LaurentMPoly Q
y = var "y" :: LaurentMPoly Q
z = var "z" :: LaurentMPoly Q


denominatorLP (LP ts) = LP [(m',1)] where
    m' = foldl lcmLM 1 [denominatorLM m | (m,c) <- ts]

{-
-- not tested for terms with non-zero denominator
gcdTermsLP (LP ts) = LP [(m',1)] where
    m' = foldl gcdLM 1 [m | (m,c) <- ts]
-}

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

sqrtvar v = LP [(LM $ M.singleton v (1/2), 1)]

-- substitute terms for variables in an MPoly
-- eg subst [(x,a),(y,a+b),(z,c^2)] (x*y+z) -> a*(a+b)+c^2
subst vts (LP us) = sum [inject c * substM m | (m,c) <- us] where
    substM (LM m) = product [substV v ^^^ i | (v,i) <- M.toList m]
    substV v =
        let v' = var v in
        case L.lookup v' vts of
        Just t -> t
        Nothing -> v' -- no substitute, so keep as is

f ^^^ i | denominatorQ i == 1 = f ^^ numeratorQ i -- exponent is an integer
        | otherwise = case f of
                      LP [(LM m,1)] -> LP [(LM $ M.map (*i) m ,1)]-- base is a monomial
                      otherwise     -> error ("(^^^): Cannot calculate " ++ show f ++ " ^^^ " ++ show i)
{-
-- halve all indices - useful when we really want to be working over k[t^1/2,t^-1/2]
halfExponents (LP ts) =
    if any odd (concatMap (\(LM m,c) -> M.elems m) ts)
    then error ("halfExponents: " ++ show (LP ts))
    else LP $ map (\(LM m, c) -> (LM $ M.filter (/=0) $ M.map (`div` 2) m, c)) ts

halfExponents' f@(LP ts) =
    let f'@(LP us) = LP $ map (\(LM m, c) -> (LM $ M.map (`div` 2) m, c)) ts
    in if any (==0) (concatMap (\(LM m,c) -> M.elems m) us)
       then Left f
       else Right f'

quarterExponents' f@(LP ts) =
    let f'@(LP us) = LP $ map (\(LM m, c) -> (LM $ M.map (`div` 4) m, c)) ts
    in if any (==0) (concatMap (\(LM m,c) -> M.elems m) us)
       then Left f
       else Right f'
-}