np-extras-0.2: MathObj/MultiVarPolynomial.hs

-- | Polynomials in a countably infinite set of variables x1, x2, x3, ...
module MathObj.MultiVarPolynomial 
    ( -- * Type
      T(..)

      -- * Constructing polynomials
    , fromMonomials
    , lift0
    , lift1
    , lift2
    , x
    , constant

      -- * Operations

    , compose

    , merge

    ) where

import qualified Algebra.Additive as Additive
import qualified Algebra.Ring as Ring
import qualified Algebra.ZeroTestable as ZeroTestable
import qualified Algebra.Differential as Differential

import qualified MathObj.Monomial as Mon

import qualified Data.Map as M

import NumericPrelude
import PreludeBase

-- | A polynomial is just a list of monomials, construed as their sum.
--   We maintain the invariant that polynomials are always sorted by
--   the ordering on monomials defined in "MathObj.Monomial": first by
--   partition degree, then by largest variable index (decreasing),
--   then by exponent of the highest-index variable (decreasing).
--   This works out nicely for operations on cycle index series.
--
--   Instances are provided for Additive, Ring, Differential
--   (partial differentiation with respect to x1), and Show.
newtype T a = Cons [Mon.T a]

instance (ZeroTestable.C a, Ring.C a, Ord a, Show a) => Show (T a) where
  show (Cons []) = "0"
  show (Cons (m:ms)) = show m ++ concatMap showMon ms
    where showMon m | Mon.coeff m < 0 = " - " ++ show (negate m)
                    | otherwise       = " + " ++ show m

{-# INLINE fromMonomials #-}
fromMonomials :: [Mon.T a] -> T a
fromMonomials = lift0

{-# INLINE lift0 #-}
lift0 :: [Mon.T a] -> T a
lift0 = Cons

{-# INLINE lift1 #-}
lift1 :: ([Mon.T a] -> [Mon.T a]) -> (T a -> T a)
lift1 f (Cons xs) = Cons (f xs)

{-# INLINE lift2 #-}
lift2 :: ([Mon.T a] -> [Mon.T a] -> [Mon.T a]) -> (T a -> T a -> T a)
lift2 f (Cons xs) (Cons ys) = Cons (f xs ys)

-- | Create the polynomial xn for a given n.
x :: (Ring.C a) => Integer -> T a
x n = fromMonomials [Mon.x n]

-- | Create a constant polynomial.
constant :: a -> T a
constant a = fromMonomials [Mon.constant a]

-- | Add two polynomials.  We assume that they are already sorted, so
--   that addition works on infinite polynomials.
add :: (Ord a, Additive.C a) => [a] -> [a] -> [a]
add xs ys = merge True (+) xs ys

-- | Merge two sorted lists, with a flag specifying whether to keep
--   singletons, and a combining function for elements that are equal.
merge :: Ord a => Bool -> (a -> a -> a) -> [a] -> [a] -> [a]
merge True  _ [] ys = ys
merge False _ [] _  = []
merge True  _ xs [] = xs
merge False _ _  [] = []
merge b f xxs@(x:xs) yys@(y:ys) | x < y     = if' b (x:) id $ merge b f xs yys
                                | x > y     = if' b (y:) id $ merge b f xxs ys
                                | otherwise = f x y : merge b f xs ys
  where if' True x _ = x
        if' False _ y = y

instance (Additive.C a, ZeroTestable.C a) => Additive.C (T a) where
  zero   = fromMonomials []
  negate = lift1 $ map negate
  (+)    = lift2 add

-- | Multiply two (sorted) polynomials.
mul :: (Ring.C a, Ord a) => [a] -> [a] -> [a]
mul [] _ = []
mul _ [] = []
mul (x:xs) (y:ys) = x*y : add (map (x*) ys) (mul xs (y:ys))

instance (Ring.C a, ZeroTestable.C a) => Ring.C (T a) where
  fromInteger n = fromMonomials [fromInteger n]
  (*) = lift2 mul

-- Partial differentiation with respect to x1.
instance (ZeroTestable.C a, Ring.C a) => Differential.C (T a) where
  differentiate = lift1 $ filter (not . isZero) . map Differential.differentiate

-- | Plethyistic substitution: F o G = F(G(x1,x2,x3...),
--   G(x2,x4,x6...), G(x3,x6,x9...), ...)  See Bergeron, Labelle, and
--   Leroux, \"Combinatorial Species and Tree-Like Structures\",
--   p. 43.
compose :: (Ring.C a, ZeroTestable.C a) => T a -> T a -> T a
compose (Cons []) _ = Cons []
compose (Cons (x:_)) (Cons []) = Cons [x]
compose (Cons xs) yys@(Cons (y:ys)) 
  | Mon.degree y == 0 = error "MultiVarPolynomial.compose: inner series must not have a constant term."
  | otherwise = comp xs yys

-- | We need to be careful to make sure this is suitably
--   lazy. For example, this works for finite polynomials:
--
-- > comp ms p = sum . map (substMon p) $ ms
--
--   but not for infinite ones!
--
--   This is accomplished by calling a recursive helper function
--   taking as an extra argument a running sum containing only terms
--   with partition degree greater than or equal to the most recently
--   processed monomial.  Plethyistically substituting a polynomial
--   (with no constant term) into a monomial of partition degree d
--   produces a polynomial with all terms of partition degree >= d, so
--   when we encounter a monomial with partition degree d, we know we
--   are done with all terms in the running sum of lesser partition
--   degree.
--
--   Precondition: the second argument has no constant term.
comp :: (Ring.C a, ZeroTestable.C a) => [Mon.T a] -> T a -> T a
comp ms p = comp' 0 ms
  where -- comp' :: T a -> [Mon.T a] -> T a
        comp' part []     = part
        comp' part (m:ms) = lift2 (++) done $ comp' (rest + substMon p m) ms
          where (done,rest) = splitPoly ((< Mon.pDegree m) . Mon.pDegree) part

-- | Plethyistic substitution of a polynomial into a monomial.
substMon :: (ZeroTestable.C a, Ring.C a) => T a -> Mon.T a -> T a
substMon poly m
  = (constant (Mon.coeff m) *)
  . M.foldWithKey (\sub pow -> (*) (scalePoly sub poly ^pow)) 1 
  $ Mon.powers m
  
-- | @scalePoly n Z@ changes Z(x_1, x_2, x_3, ...) into Z(x_n, x_2n, x_3n, ...)
scalePoly :: Integer -> T a -> T a
scalePoly n = lift1 $ map (Mon.scaleMon n)

-- | Split a polynomial into two pieces based on a predicate.
splitPoly :: (Mon.T a -> Bool) -> T a -> (T a, T a)
splitPoly p (Cons xs) = (Cons ys, Cons zs)
  where (ys, zs) = span p xs