np-extras (empty) → 0.1

raw patch · 6 files changed

+540/−0 lines, 6 filesdep +basedep +containersdep +numeric-preludesetup-changed

Dependencies added: base, containers, numeric-prelude, primes

Files

• LICENSE +27/−0
• MathObj/FactoredRational.hs +196/−0
• MathObj/Monomial.hs +137/−0
• MathObj/MultiVarPolynomial.hs +156/−0
• Setup.hs +2/−0
• np-extras.cabal +22/−0

@@ -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.

@@ -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

@@ -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)

@@ -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

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -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