diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,27 @@
+Copyright (c) Brent Yorgey 2009
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. Neither the name of the author nor the names of other contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGE.
diff --git a/MathObj/FactoredRational.hs b/MathObj/FactoredRational.hs
new file mode 100644
--- /dev/null
+++ b/MathObj/FactoredRational.hs
@@ -0,0 +1,196 @@
+-- | A representation of rational numbers as lists of prime powers,
+--   allowing efficient representation, multiplication and division of
+--   large numbers, especially of the sort occurring in combinatorial
+--   computations.
+-- 
+--   The module also includes a method for generating factorials in
+--   factored form directly, and for computing Euler's totient and
+--   generating all divisors of factored integers.
+module MathObj.FactoredRational 
+    ( -- * Type
+      T
+
+      -- * Utilities
+    , factorial
+    , eulerPhi
+    , divisors
+ 
+    ) where
+
+import qualified Algebra.Additive as Additive
+import qualified Algebra.Ring as Ring
+import qualified Algebra.Field as Field
+import qualified Algebra.IntegralDomain as Integral
+
+import qualified Algebra.ZeroTestable as ZeroTestable
+import qualified Algebra.Real as Real
+import qualified Algebra.ToRational as ToRational
+import qualified Algebra.RealIntegral as RealIntegral
+import qualified Algebra.ToInteger as ToInteger
+
+import Data.Numbers.Primes
+
+import PreludeBase 
+import NumericPrelude
+
+-- Represent rational numbers by their prime factorizations.
+-- Perhaps this should use a sparse representation instead, using a Map from 
+-- primes to powers?  Well, that should be easy enough to change later.
+
+-- | The type of factored rationals.
+--
+--   Instances are provided for Eq, Ord, Additive, Ring, ZeroTestable,
+--   Real, ToRational, Integral, RealIntegral, ToInteger, and Field.
+--
+--   Note that currently, addition is performed on factored rationals
+--   by converting them to normal rationals, performing the addition,
+--   and factoring.  This could probably be made more efficient by
+--   finding a common denominator, pulling out common factors from the
+--   numerators, and performing the addition and factoring only on the
+--   relatively prime parts.
+data T = FQZero            -- ^ zero
+       | FQ Bool [Integer] -- ^ prime exponents with sign bit, True = negative.
+
+-- XXX this ought to be improved.
+instance Show T where
+  show FQZero = "0"
+  show (FQ True pows) = "(-1)" ++ showPows pows
+  show (FQ False pows) = showPows pows
+
+showPows :: [Integer] -> String
+showPows pows = concat $ zipWith showPow primes pows
+  where showPow p 0 = ""
+        showPow p 1 = "(" ++ show p ++ ")"
+        showPow p n = "(" ++ show p ++ "^" ++ show n ++ ")"
+
+instance Additive.C T where
+  zero = FQZero
+  FQZero + a = a
+  a + FQZero = a
+  x + y = fromRational' (toRational x + toRational y) 
+
+  negate FQZero   = FQZero
+  negate (FQ s e) = FQ (not s) e
+
+instance Ring.C T where
+  FQZero * _ = FQZero
+  _ * FQZero = FQZero
+  (FQ s1 e1) * (FQ s2 e2) = FQ (s1 /= s2) (zipWithExt 0 0 (+) e1 e2)
+
+  fromInteger 0 = FQZero
+  fromInteger n | n < 0     = FQ True (factor (negate n))
+                | otherwise = FQ False (factor n)
+
+  _ ^ 0 = one
+  FQZero   ^ _ = FQZero
+  (FQ s e) ^ n = FQ s (map (*n) e)
+
+-- | Zip two lists together with a combining function, using default
+--   values to extend the lists if one is shorter than the other.
+zipWithExt :: a -> b -> (a -> b -> c) -> [a] -> [b] -> [c]
+zipWithExt da db f = zipWithExt'
+  where zipWithExt' []     bs     = zipWith f (repeat da) bs
+        zipWithExt' as     []     = zipWith f as (repeat db)
+        zipWithExt' (a:as) (b:bs) = f a b : zipWithExt' as bs
+
+-- | A simple factoring method. 
+--
+--   We should probably just depend on another module with some
+--   dedicated, efficient factoring code written by someone really
+--   smart, but this simple method works OK for now.
+--
+--   Precondition: argument is positive.
+factor :: Integer -> [Integer]
+factor n = factor' n primes
+  where
+    factor' 1 _ = []
+    factor' n (p:ps) = let (k,n') = logRem n p
+                       in  k : factor' n' ps
+   
+-- | @logRem n p@ computes (k,n'), where k is the highest power of p
+--   that divides n, and n' = n `div` p^k.
+logRem :: Integer -> Integer -> (Integer, Integer)
+logRem = logRem' 0
+  where logRem' k n p | n `mod` p == 0 = logRem' (k+1) (n `div` p) p
+                      | otherwise = (k,n)
+
+instance ZeroTestable.C T where 
+  isZero FQZero = True
+  isZero _      = False
+
+instance Eq T where
+  FQZero == FQZero   = True
+  FQZero == (FQ _ _) = False
+  (FQ _ _) == FQZero = False
+  (FQ s1 e1) == (FQ s2 e2) = s1 == s2 && dropZeros e1 == dropZeros e2
+    where dropZeros = reverse . dropWhile (==0) . reverse
+
+instance Ord T where
+  compare FQZero FQZero = EQ
+  compare FQZero (FQ False _) = LT
+  compare FQZero (FQ True  _) = GT
+  compare (FQ False _) FQZero = GT
+  compare (FQ True  _) FQZero = LT
+  compare (FQ False _) (FQ True  _)   = GT
+  compare (FQ True  _) (FQ False _)   = LT
+  compare fq1 fq2 = compare (toRational fq1) (toRational fq2)
+
+instance Real.C T where
+  abs FQZero   = FQZero
+  abs (FQ _ e) = FQ False e
+
+instance ToRational.C T where
+  toRational FQZero   = 0
+  toRational (FQ s e) = (if s then negate else id)
+                      . product
+                      . zipWith (^-) (map (%1) primes)
+                      $ e
+
+instance Integral.C T where
+  divMod a b =
+    if isZero b
+      then (undefined,a)
+      else (a/b,0)
+
+instance RealIntegral.C T 
+  -- default definition is fine
+
+instance ToInteger.C T where
+  toInteger FQZero = 0
+  toInteger (FQ s e) | any (<0) e = error "non-integer in FactoredRational.toInteger"
+                     | otherwise  = (if s then negate else id)
+                                  . product
+                                  . zipWith (^) primes
+                                  $ e
+
+instance Field.C T where
+  recip FQZero = error "division by zero"
+  recip (FQ s e) = FQ s (map negate e)
+
+-- | Efficiently compute n! directly as a factored rational.
+factorial :: Integer -> T
+factorial 0 = one
+factorial 1 = one
+factorial n = FQ False (takeWhile (>0) . map (factorialFactors n) $ primes)
+
+-- | @factorialFactors n p@ computes the power of prime p in the
+--   factorization of n!.
+factorialFactors :: Integer -> Integer -> Integer
+factorialFactors n p = sum 
+                     . takeWhile (>0)
+                     . map (n `div`)
+                     $ iterate (*p) p
+
+-- | Compute Euler's totient function (@eulerPhi n@ is the number of
+--   integers less than and relatively prime to n).  Only makes sense
+--   for (nonnegative) integers.
+eulerPhi :: T -> Integer
+eulerPhi FQZero = 1
+eulerPhi (FQ _ pows) = product $ zipWith phiP primes pows
+  where phiP _ 0 = 1
+        phiP p a = p^(a-1) * (p-1)
+
+-- | List of the divisors of n.  Only makes sense for integers.
+divisors :: T -> [T]
+divisors FQZero = [1]
+divisors (FQ b pows) = map (FQ b) $ mapM (enumFromTo 0) pows
diff --git a/MathObj/Monomial.hs b/MathObj/Monomial.hs
new file mode 100644
--- /dev/null
+++ b/MathObj/Monomial.hs
@@ -0,0 +1,137 @@
+{-# LANGUAGE PatternGuards #-}
+
+-- | Monomials in a countably infinite set of variables x1, x2, x3, ...
+module MathObj.Monomial
+    ( -- * Type
+      T(..)
+
+      -- * Creating monomials
+    , constant
+    , x
+
+      -- * Utility functions
+    , degree
+    , pDegree
+    , scaleMon
+
+    ) 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 Algebra.Field as Field
+
+import qualified Data.Map as M
+import Data.Ord (comparing)
+import Control.Arrow ((***))
+import Data.List (sort, intercalate)
+
+import NumericPrelude
+import PreludeBase
+
+-- | A monomial is a map from variable indices to integer powers,
+--   paired with a (polymorphic) coefficient.  Note that negative
+--   integer powers are handled just fine, so monomials form a field.
+--
+--   Instances are provided for Eq, Ord, ZeroTestable, Additive, Ring,
+--   Differential, and Field.  Note that adding two monomials only
+--   makes sense if they have matching variables and exponents.  The
+--   Differential instance represents partial differentiation with
+--   respect to x1.
+--
+--   The Ord instance for monomials orders them first by permutation
+--   degree, then by largest variable index (largest first), then by
+--   exponent (largest first).  This may seem a bit odd, but in fact
+--   reflects the use of these monomials to implement cycle index
+--   series, where this ordering corresponds nicely to generation
+--   of integer partitions. To make the library more general we could
+--   parameterize monomials by the desired ordering.
+data T a = Cons { coeff  :: a 
+                , powers :: M.Map Integer Integer
+                }
+
+instance (ZeroTestable.C a, Ring.C a, Eq a, Show a) => Show (T a) where
+  show (Cons a pows) | isZero a  = "0"
+                     | a == 1    = showVars pows
+                     | a == (-1) = "-" ++ showVars pows
+                     | otherwise = show a ++ " " ++ showVars pows
+
+showVars :: M.Map Integer Integer -> String
+showVars m = intercalate " " $ concatMap showVar (M.toList m)
+  where showVar (_,0) = []
+        showVar (v,1) = ["x" ++ show v]
+        showVar (v,p) = ["x" ++ show v ++ "^" ++ show p]
+
+-- | The degree of a monomial is the sum of its exponents.
+degree :: T a -> Integer
+degree (Cons _ m) = M.fold (+) 0 m
+
+-- | The \"partition degree\" of a monomial is the sum of the products
+--   of each variable index with its exponent.  For example, x1^3 x2^2
+--   x4^3 has partition degree 1*3 + 2*2 + 4*3 = 19.  The terminology
+--   comes from the fact that, for example, we can view x1^3 x2^2 x4^3
+--   as corresponding to an integer partition of 19 (namely, 1 + 1 + 1
+--   + 2 + 2 + 4 + 4 + 4).
+pDegree :: T a -> Integer
+pDegree (Cons _ m) = sum . map (uncurry (*)) . M.assocs $ m
+
+-- | Create a constant monomial.
+constant :: a -> T a
+constant a = Cons a M.empty
+
+-- | Create the monomial xn for a given n.
+x :: (Ring.C a) => Integer -> T a
+x n = Cons Ring.one (M.singleton n 1)
+
+-- | Scale all the variable subscripts by a constant.  Useful for
+--   operations like plethyistic substitution or Mobius inversion.
+scaleMon :: Integer -> T a -> T a
+scaleMon n (Cons a m) = Cons a (M.mapKeys (n*) m)
+
+instance Eq (T a) where
+  (Cons _ m1) == (Cons _ m2) = m1 == m2
+
+instance Ord (T a) where
+  compare m1 m2
+    | d1 < d2   = LT
+    | d1 > d2   = GT
+    | otherwise = comparing q m1 m2
+    where d1 = pDegree m1
+          d2 = pDegree m2
+          q  = map Rev . reverse . sort . M.assocs . powers
+
+newtype Rev a = Rev { getRev :: a }
+  deriving Eq
+instance Ord a => Ord (Rev a) where
+  compare (Rev a) (Rev b) = compare b a
+
+instance (ZeroTestable.C a) => ZeroTestable.C (T a) where
+  isZero (Cons a _) = isZero a
+  
+instance (Additive.C a, ZeroTestable.C a) => Additive.C (T a) where
+  zero = Cons zero M.empty
+  negate (Cons a m) = Cons (negate a) m
+
+  -- precondition: m1 == m2
+  (Cons a1 m1) + (Cons a2 _m2) | isZero s  = Cons s M.empty
+                               | otherwise = Cons s m1
+                               where s = a1 + a2
+
+instance (Ring.C a, ZeroTestable.C a) => Ring.C (T a) where
+  fromInteger n = Cons (fromInteger n) M.empty
+  (Cons a1 m1) * (Cons a2 m2) = Cons (a1*a2) 
+                                     (M.filterWithKey (\_ p -> not (isZero p)) $
+                                        M.unionWith (+) m1 m2
+                                     )
+
+-- Partial differentiation with respect to x1.
+instance (ZeroTestable.C a, Ring.C a) => Differential.C (T a) where
+  differentiate (Cons a m) 
+    | Just 1 <- M.lookup 1 m = Cons a M.empty
+    | Just p <- M.lookup 1 m = Cons (a*fromInteger p) (M.adjust (subtract 1) 1 m)
+    | otherwise              = Cons 0 M.empty
+
+instance (ZeroTestable.C a, Field.C a, Eq a) => Field.C (T a) where
+  recip (Cons 0 _)    = error "Monomial.recip: division by zero"
+  recip (Cons a pows) = Cons (recip a) (M.map negate pows)
diff --git a/MathObj/MultiVarPolynomial.hs b/MathObj/MultiVarPolynomial.hs
new file mode 100644
--- /dev/null
+++ b/MathObj/MultiVarPolynomial.hs
@@ -0,0 +1,156 @@
+-- | 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
+
+    ) 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 (+) xs ys
+
+-- | Merge two sorted lists, with a combining function for elements
+--   that are equal.
+merge :: Ord a => (a -> a -> a) -> [a] -> [a] -> [a]
+merge _ [] ys = ys
+merge _ xs [] = xs
+merge f xxs@(x:xs) yys@(y:ys) | x < y     = x : merge f xs yys
+                              | x > y     = y : merge f xxs ys
+                              | otherwise = (f x y) : merge f xs ys
+
+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
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/np-extras.cabal b/np-extras.cabal
new file mode 100644
--- /dev/null
+++ b/np-extras.cabal
@@ -0,0 +1,22 @@
+name:           np-extras
+version:        0.1
+license:        BSD3
+license-file:   LICENSE
+build-type:     Simple
+cabal-version:  >= 1.2.3
+tested-with:    GHC == 6.10.3
+author:         Brent Yorgey
+maintainer:     Brent Yorgey <byorgey@cis.upenn.edu>
+category:       Math
+synopsis:       NumericPrelude extras
+description:    Various extras to extend the NumericPrelude, including
+                multivariate polynomials and factored rationals.
+
+Library
+  build-depends: base >= 3.0 && < 4.2, numeric-prelude >= 0.1.1 && < 0.2,
+                 primes >= 0.1.1 && < 0.2, containers >= 0.2 && < 0.3
+  exposed-modules:
+    MathObj.FactoredRational
+    MathObj.Monomial
+    MathObj.MultiVarPolynomial
+  extensions: NoImplicitPrelude