combinat 0.2.6.1 → 0.2.6.2

raw patch · 8 files changed

+487/−31 lines, 8 filesPVP: minor bump suggested

API additions: PVP suggests at least a minor version bump

API changes (from Hackage documentation)

Files

• Math/Combinat/Helper.hs +11/−10
• Math/Combinat/LatticePaths.hs +106/−4
• Math/Combinat/Numbers.hs +7/−5
• Math/Combinat/Partitions.hs +62/−9
• Math/Combinat/Partitions/NonCrossing.hs +195/−0
• Math/Combinat/Partitions/Set.hs +99/−0
• Math/Combinat/Trees/Binary.hs +3/−1
• combinat.cabal +4/−2

@@ -29,16 +29,6 @@ swap :: (a,b) -> (b,a) swap (x,y) = (y,x) ------------------------------------------------------------------------------------ * lists--{-# SPECIALIZE sum' :: [Int]     -> Int     #-}-{-# SPECIALIZE sum' :: [Integer] -> Integer #-}-sum' :: Num a => [a] -> a-sum' = foldl' (+) 0----------------------------------------------------------------------------------- pairs :: [a] -> [(a,a)] pairs = go where   go (x:xs@(y:_)) = (x,y) : go xs@@ -50,6 +40,14 @@   go _            = []  --------------------------------------------------------------------------------+-- * lists++{-# SPECIALIZE sum' :: [Int]     -> Int     #-}+{-# SPECIALIZE sum' :: [Integer] -> Integer #-}+sum' :: Num a => [a] -> a+sum' = foldl' (+) 0++-------------------------------------------------------------------------------- -- * equality and ordering   equating :: Eq b => (a -> b) -> a -> a -> Bool@@ -62,6 +60,9 @@  reverseCompare :: Ord a => a -> a -> Ordering reverseCompare x y = reverseOrdering $ compare x y++reverseSort :: Ord a => [a] -> [a]+reverseSort = sortBy reverseCompare  groupSortBy :: (Eq b, Ord b) => (a -> b) -> [a] -> [[a]] groupSortBy f = groupBy (equating f) . sortBy (comparing f) 

@@ -6,6 +6,8 @@ 
 --------------------------------------------------------------------------------
 
+import System.Random
+
 import Math.Combinat.Numbers
 import Math.Combinat.Trees.Binary
 
@@ -36,6 +38,14 @@                  in  if y'<0 then False 
                              else go y' ts
 
+-- | A Dyck path is a lattice path whose last point lies on the @y=0@ line
+isDyckPath :: LatticePath -> Bool
+isDyckPath = go 0 where
+  go !y []     = y==0
+  go !y (t:ts) = let y' = case t of { UpStep -> y+1 ; DownStep -> y-1 }
+                 in  if y'<0 then False 
+                             else go y' ts
+
 -- | Maximal height of a lattice path
 pathHeight :: LatticePath -> Int
 pathHeight = go 0 0 where
@@ -53,7 +63,7 @@     DownStep -> go (x+1) (y-1) ts
 
 -- | Returns the coordinates of the path (excluding the starting point @(0,0)@, but including
--- the endpoint
+-- the endpoint)
 pathCoordinates :: LatticePath -> [(Int,Int)]
 pathCoordinates = go 0 0 where
   go _  _  []     = []
@@ -61,6 +71,13 @@                         y' = case t of { UpStep -> y+1 ; DownStep -> y-1 }
                     in  (x',y') : go x' y' ts
 
+-- | Number of peaks of a path (excluding the endpoint)
+pathNumberOfPeaks :: LatticePath -> Int
+pathNumberOfPeaks = go 0 where
+  go !k (x:xs@(y:_)) = go (if x==UpStep && y==DownStep then k+1 else k) xs
+  go !k [x] = k
+  go !k [ ] = k
+
 -- | Number of points on the path which touch the @y=0@ zero level line
 -- (excluding the starting point @(0,0)@, but including the endpoint; that is, for Dyck paths it this is always positive!).
 pathNumberOfZeroTouches :: LatticePath -> Int
@@ -86,14 +103,41 @@ -- Remark: Dyck paths are obviously in bijection with nested parentheses, and thus
 -- also with binary trees.
 --
+-- Order is reverse lexicographical:
+--
+-- > sort (dyckPaths m) == reverse (dyckPaths m)
+-- 
 dyckPaths :: Int -> [LatticePath]
-dyckPaths = map (map f) . nestedParentheses where
-  f p = case p of { LeftParen -> UpStep ; RightParen -> DownStep }
+dyckPaths = map nestedParensToDyckPath . nestedParentheses 
 
+-- | @dyckPaths m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@. 
+--
+-- > sort (dyckPathsNaive m) == sort (dyckPaths m) 
+--  
+-- Naive recursive algorithm, order is ad-hoc
+--
+dyckPathsNaive :: Int -> [LatticePath]
+dyckPathsNaive = worker where
+  worker  0 = [[]]
+  worker  m = as ++ bs where
+    as = [ bracket p      | p <- worker (m-1) ] 
+    bs = [ bracket p ++ q | k <- [1..m-1] , p <- worker (k-1) , q <- worker (m-k) ]
+  bracket p = UpStep : p ++ [DownStep]
+
 -- | The number of Dyck paths from @(0,0)@ to @(2m,0)@ is simply the m\'th Catalan number.
 countDyckPaths :: Int -> Integer
 countDyckPaths m = catalan m
 
+-- | The trivial bijection
+nestedParensToDyckPath :: [Paren] -> LatticePath
+nestedParensToDyckPath = map f where
+  f p = case p of { LeftParen -> UpStep ; RightParen -> DownStep }
+
+-- | The trivial bijection in the other direction
+dyckPathToNestedParens :: LatticePath -> [Paren]
+dyckPathToNestedParens = map g where
+  g s = case s of { UpStep -> LeftParen ; DownStep -> RightParen }
+
 --------------------------------------------------------------------------------
 -- * Bounded Dyck paths
 
@@ -163,7 +207,7 @@ --------------------------------------------------------------------------------
 -- * Zero-level touches
 
--- | @touchingDyckPathsNaive k m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ which touch the 
+-- | @touchingDyckPaths k m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ which touch the 
 -- zero level line @y=0@ exactly @k@ times (excluding the starting point, but including the endpoint;
 -- thus, @k@ should be positive). Synonym for 'touchingDyckPathsNaive'.
 touchingDyckPaths
@@ -194,6 +238,64 @@     | otherwise = [ bracket p ++ q | l <- [1..m-1] , p <- dyckPaths (l-1) , q <- worker (k-1) (m-l) ]
     where
       bracket p = UpStep : p ++ [DownStep] 
+
+--------------------------------------------------------------------------------
+-- * Dyck paths with given number of peaks
+
+-- | @peakingDyckPaths k m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ with exactly @k@ peaks.
+--
+-- Synonym for 'peakingDyckPathsNaive'
+--
+peakingDyckPaths
+  :: Int      -- ^ @k@ = number of peaks
+  -> Int      -- ^ @m@ = half-length
+  -> [LatticePath]
+peakingDyckPaths = peakingDyckPathsNaive 
+
+-- | @peakingDyckPathsNaive k m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ with exactly @k@ peaks.
+--
+-- > sort (peakingDyckPathsNaive k m) = sort [ p | p <- dyckPaths m , pathNumberOfPeaks p == k ]
+--  
+-- Naive recursive algorithm, resulting order is pretty ad-hoc.
+--
+peakingDyckPathsNaive 
+  :: Int      -- ^ @k@ = number of peaks
+  -> Int      -- ^ @m@ = half-length
+  -> [LatticePath]
+peakingDyckPathsNaive = worker where
+  worker !k !m
+    | m == 0    = if k==0 then [[]] else []       
+    | k <= 0    = []
+    | m <  0    = []
+    | k == 1    = [ singlePeak m ] 
+    | otherwise = as ++ bs ++ cs
+    where
+      as = [ bracket p      |                                 p <- worker k (m-1)                           ]
+      bs = [ smallHill ++ q |                                                       q <- worker (k-1) (m-1) ]
+      cs = [ bracket p ++ q | l <- [2..m-1] , a <- [1..k-1] , p <- worker a (l-1) , q <- worker (k-a) (m-l) ]
+      smallHill     = [ UpStep , DownStep ]
+      singlePeak !m = replicate m UpStep ++ replicate m DownStep 
+      bracket p = UpStep : p ++ [DownStep] 
+
+-- | Dyck paths of length @2m@ with @k@ peaks are counted by the Narayana numbers @N(m,k) = \binom{m}{k} \binom{m}{k-1} / m@
+countPeakingDyckPaths
+  :: Int      -- ^ @k@ = number of peaks
+  -> Int      -- ^ @m@ = half-length
+  -> Integer
+countPeakingDyckPaths k m 
+  | m == 0    = if k==0 then 1 else 0
+  | k <= 0    = 0
+  | m <  0    = 0
+  | k == 1    = 1
+  | otherwise = div (binomial m k * binomial m (k-1)) (fromIntegral m)
+
+--------------------------------------------------------------------------------
+-- * Random lattice paths
+
+-- | A uniformly random Dyck path of length @2m@
+randomDyckPath :: RandomGen g => Int -> g -> (LatticePath,g)
+randomDyckPath m g0 = (nestedParensToDyckPath parens, g1) where
+  (parens,g1) = randomNestedParentheses m g0
 
 --------------------------------------------------------------------------------
 

@@ -109,9 +109,10 @@ -- signedStirling1st :: Integral a => a -> a -> Integer signedStirling1st n k -  | k < 1     = 0-  | k > n     = 0-  | otherwise = signedStirling1stArray n ! (fromIntegral k)+  | k==0 && n==0 = 1+  | k < 1        = 0+  | k > n        = 0+  | otherwise    = signedStirling1stArray n ! (fromIntegral k)  -- | (Unsigned) Stirling numbers of the first kind. See 'signedStirling1st'. unsignedStirling1st :: Integral a => a -> a -> Integer@@ -124,8 +125,9 @@ -- stirling2nd :: Integral a => a -> a -> Integer stirling2nd n k -  | k < 1     = 0-  | k > n     = 0+  | k==0 && n==0 = 1+  | k < 1        = 0+  | k > n        = 0   | otherwise = sum xs `div` factorial k where       xs = [ paritySign (k-i) * binomial k i * (fromIntegral i)^n | i<-[0..k] ] 

@@ -2,9 +2,14 @@ -- | Partitions of integers and multisets.  -- Integer partitions are nonincreasing sequences of positive integers. ----- See also ---   Donald E. Knuth: The Art of Computer Programming, vol 4, pre-fascicle 3B.+-- See: --+--  * Donald E. Knuth: The Art of Computer Programming, vol 4, pre-fascicle 3B.+--+--  * <http://en.wikipedia.org/wiki/Partition_(number_theory)>+--++{-# LANGUAGE BangPatterns #-} module Math.Combinat.Partitions   ( -- * Type and basic stuff     Partition@@ -23,6 +28,7 @@   , _elements   , countAutomorphisms   , _countAutomorphisms+  , HasNumberOfParts(..)     -- * Generation   , partitions'     , _partitions' @@ -30,10 +36,12 @@   , partitions   , _partitions   , countPartitions+  , partitionsWithKParts   , allPartitions'     , allPartitions    , countAllPartitions'   , countAllPartitions+  , countPartitionsWithKParts     -- * Ferrer diagrams   , printFerrerDiagram    , ferrerDiagram @@ -103,7 +111,7 @@ -- | The weight of the partition  --   (that is, the sum of the corresponding sequence). weight :: Partition -> Int-weight (Partition part) = sum part+weight (Partition part) = sum' part  -- | The dual (or conjugate) partition. dualPartition :: Partition -> Partition@@ -133,18 +141,26 @@  _countAutomorphisms :: [Int] -> Integer _countAutomorphisms = multinomial . map length . group- + --------------------------------------------------------------------------------- +class HasNumberOfParts p where+  numberOfParts :: p -> Int++instance HasNumberOfParts Partition where+  numberOfParts (Partition p) = length p+  +---------------------------------------------------------------------------------+ -- | Integer partitions of @d@, fitting into a given rectangle, as lists. _partitions'    :: (Int,Int)     -- ^ (height,width)   -> Int           -- ^ d   -> [[Int]]         _partitions' _ 0 = [[]] -_partitions' (0,_) d = if d==0 then [[]] else []-_partitions' (_,0) d = if d==0 then [[]] else []-_partitions' (h,w) d = +_partitions' ( 0 , _) d = if d==0 then [[]] else []+_partitions' ( _ , 0) d = if d==0 then [[]] else []+_partitions' (!h ,!w) d =    [ i:xs | i <- [1..min d h] , xs <- _partitions' (i,w-1) (d-i) ]  -- | Partitions of d, fitting into a given rectangle. The order is again lexicographic.@@ -189,9 +205,46 @@   --sum [ countPartitions' (h,w) i | i <- [0..d] ] where d = h*w  countAllPartitions :: Int -> Integer-countAllPartitions d = sum [ countPartitions i | i <- [0..d] ]+countAllPartitions d = sum' [ countPartitions i | i <- [0..d] ]  --------------------------------------------------------------------------------+-- partitions with given number of parts++-- | Lists partitions of @n@ into @k@ parts.+--+-- > sort (partitionsWithKParts k n) == sort [ p | p <- partitions n , numberOfParts p == k ]+--+-- Naive recursive algorithm.+--+partitionsWithKParts +  :: Int    -- ^ @k@ = number of parts+  -> Int    -- ^ @n@ = the integer we partition+  -> [Partition]+partitionsWithKParts k n = map Partition $ go n k n where+{-+  h = max height+  k = number of parts+  n = integer+-}+  go !h !k !n +    | k <  0     = []+    | k == 0     = if h>=0 && n==0 then [[] ] else []+    | k == 1     = if h>=n && n>=1 then [[n]] else []+    | otherwise  = [ a:p | a <- [1..(min h (n-k+1))] , p <- go a (k-1) (n-a) ]++countPartitionsWithKParts +  :: Int    -- ^ @k@ = number of parts+  -> Int    -- ^ @n@ = the integer we partition+  -> Integer+countPartitionsWithKParts k n = go n k n where+  go !h !k !n +    | k <  0     = 0+    | k == 0     = if h>=0 && n==0 then 1 else 0+    | k == 1     = if h>=n && n>=1 then 1 else 0+    | otherwise  = sum' [ go a (k-1) (n-a) | a<-[1..(min h (n-k+1))] ]+++-------------------------------------------------------------------------------- -- * Ferrer diagrams  printFerrerDiagram :: Partition -> IO ()@@ -269,7 +322,7 @@           reverse first ++            (c,u,v-1) :             [ (c,u,u) | (c,u,_) <- reverse second ] -      _ -> error "should not happen"+      _ -> error "fasc3B_algorithm_M: should not happen"            to_vector cuvs =      accumArray (flip const) 0 (1,m)

@@ -0,0 +1,195 @@+
+-- | Non-crossing partitions.
+--
+-- See eg. <http://en.wikipedia.org/wiki/Noncrossing_partition>
+--
+
+{-# LANGUAGE BangPatterns #-}
+module Math.Combinat.Partitions.NonCrossing where
+
+--------------------------------------------------------------------------------
+
+import Control.Applicative
+
+import Data.List
+import Data.Ord
+
+import System.Random
+
+import Math.Combinat.Numbers
+import Math.Combinat.LatticePaths
+import Math.Combinat.Helper
+import Math.Combinat.Partitions.Set
+import Math.Combinat.Partitions ( HasNumberOfParts(..) )
+
+--------------------------------------------------------------------------------
+-- * The type of non-crossing partitions
+
+-- | A non-crossing partition of the set @[1..n]@ in standard form: 
+-- entries decreasing in each block  and blocks listed in increasing order of their first entries.
+newtype NonCrossing = NonCrossing [[Int]] deriving (Eq,Ord,Show,Read)
+
+-- | Checks whether a set partition is noncrossing.
+--
+-- Implementation method: we convert to a Dyck path and then back again, and finally compare. 
+-- Probably not very efficient, but should be better than a naive check for crosses...)
+--
+_isNonCrossing :: [[Int]] -> Bool
+_isNonCrossing zzs0 = _isNonCrossingUnsafe (_standardizeNonCrossing zzs0)
+
+-- | Warning: This function assumes the standard ordering!
+_isNonCrossingUnsafe :: [[Int]] -> Bool
+_isNonCrossingUnsafe zzs = 
+  case _nonCrossingPartitionToDyckPathMaybe zzs of
+    Nothing   -> False
+    Just dyck -> case dyckPathToNonCrossingPartitionMaybe dyck of
+      Nothing                -> False
+      Just (NonCrossing yys) -> yys == zzs
+
+-- | Convert to standard form: entries decreasing in each block 
+-- and blocks listed in increasing order of their first entries.
+_standardizeNonCrossing :: [[Int]] -> [[Int]]
+_standardizeNonCrossing = sortBy (comparing myhead) . map reverseSort where
+  myhead xs = case xs of
+    (x:xs) -> x
+    []     -> error "_standardizeNonCrossing: empty subset"
+
+fromNonCrossing :: NonCrossing -> [[Int]]
+fromNonCrossing (NonCrossing xs) = xs
+
+toNonCrossingUnsafe :: [[Int]] -> NonCrossing
+toNonCrossingUnsafe = NonCrossing
+
+-- | Throws an error if the input is not a non-crossing partition
+toNonCrossing :: [[Int]] -> NonCrossing
+toNonCrossing xxs = case toNonCrossingMaybe xxs of
+  Just nc -> nc
+  Nothing -> error "toNonCrossing: not a non-crossing partition"
+
+toNonCrossingMaybe :: [[Int]] -> Maybe NonCrossing
+toNonCrossingMaybe xxs0 = 
+  if _isNonCrossingUnsafe xxs
+    then Just $ NonCrossing xxs
+    else Nothing
+  where 
+    xxs = _standardizeNonCrossing xxs0
+
+-- | If a set partition is actually non-crossing, then we can convert it
+setPartitionToNonCrossing :: SetPartition -> Maybe NonCrossing
+setPartitionToNonCrossing (SetPartition zzs0) =
+  if _isNonCrossingUnsafe zzs
+    then Just $ NonCrossing zzs
+    else Nothing
+  where
+    zzs = _standardizeNonCrossing zzs0
+
+instance HasNumberOfParts NonCrossing where
+  numberOfParts (NonCrossing p) = length p
+
+--------------------------------------------------------------------------------
+-- * Bijection to Dyck paths
+
+-- | Bijection between Dyck paths and noncrossing partitions
+--
+-- Based on: David Callan: /Sets, Lists and Noncrossing Partitions/
+--
+-- Fails if the input is not a Dyck path.
+dyckPathToNonCrossingPartition :: LatticePath -> NonCrossing
+dyckPathToNonCrossingPartition = NonCrossing . go 0 [] [] [] where
+  go :: Int -> [Int] -> [Int] -> [[Int]] -> LatticePath -> [[Int]] 
+  go !cnt stack small big path =
+    case path of
+      (x:xs) -> case x of 
+        UpStep   -> let cnt' = cnt + 1 in case xs of
+          (y:ys)   -> case y of
+            UpStep   -> go cnt' (cnt':stack) small                  big  xs  
+            DownStep -> go cnt' (cnt':stack) []    (reverse small : big) xs
+          []       -> error "dyckPathToNonCrossingPartition: last step is an UpStep (thus input was not a Dyck path)"
+        DownStep -> case stack of
+          (k:ks)   -> go cnt ks (k:small) big xs
+          []       -> error "dyckPathToNonCrossingPartition: empty stack, shouldn't happen (thus input was not a Dyck path)"
+      [] -> tail $ reverse (reverse small : big)
+
+-- | Safe version of 'dyckPathToNonCrossingPartition'
+dyckPathToNonCrossingPartitionMaybe :: LatticePath -> Maybe NonCrossing
+dyckPathToNonCrossingPartitionMaybe = fmap NonCrossing . go 0 [] [] [] where
+  go :: Int -> [Int] -> [Int] -> [[Int]] -> LatticePath -> Maybe [[Int]] 
+  go !cnt stack small big path =
+    case path of
+      (x:xs) -> case x of 
+        UpStep   -> let cnt' = cnt + 1 in case xs of
+          (y:ys)   -> case y of
+            UpStep   -> go cnt' (cnt':stack) small                  big  xs  
+            DownStep -> go cnt' (cnt':stack) []    (reverse small : big) xs
+          []       -> Nothing
+        DownStep -> case stack of
+          (k:ks)   -> go cnt ks (k:small) big xs
+          []       -> Nothing
+      [] -> Just $ tail $ reverse (reverse small : big)
+
+-- | The inverse bijection (should never fail proper 'NonCrossing'-s)
+nonCrossingPartitionToDyckPath :: NonCrossing -> LatticePath
+nonCrossingPartitionToDyckPath (NonCrossing zzs) = go 0 zzs where
+  go !k (ys@(y:_):yys) = replicate (y-k) UpStep ++ replicate (length ys) DownStep ++ go y yys
+  go !k []             = []
+  go _  _              = error "nonCrossingPartitionToDyckPath: shouldnt't happen"
+
+-- | Safe version 'nonCrossingPartitionToDyckPath'
+_nonCrossingPartitionToDyckPathMaybe :: [[Int]] -> Maybe LatticePath
+_nonCrossingPartitionToDyckPathMaybe = go 0 where
+  go !k (ys@(y:_):yys) = fmap (\zs -> replicate (y-k) UpStep ++ replicate (length ys) DownStep ++ zs) (go y yys)
+  go !k []             = Just []
+  go _  _              = Nothing
+
+--------------------------------------------------------------------------------
+
+{- 
+-- this should be mapped to NonCrossing [[3],[5,4,2],[7,6,1],[9,8]]
+testpath = [u,u,u,d,u,u,d,d,d,u,u,d,d,d,u,u,d,d] where
+  u = UpStep
+  d = DownStep
+
+testnc = NonCrossing [[3],[5,4,2],[7,6,1],[9,8]]
+-}
+
+--------------------------------------------------------------------------------
+-- * Generating non-crossing partitions
+
+-- | Lists all non-crossing partitions of @[1..n]@
+--
+-- Equivalent to (but orders of magnitude faster than) filtering out the non-crossing ones:
+--
+-- > (sort $ catMaybes $ map setPartitionToNonCrossing $ setPartitions n) == sort (nonCrossingPartitions n)
+--
+nonCrossingPartitions :: Int -> [NonCrossing]
+nonCrossingPartitions = map dyckPathToNonCrossingPartition . dyckPaths
+
+-- | Lists all non-crossing partitions of @[1..n]@ into @k@ parts.
+--
+-- > sort (nonCrossingPartitionsWithKParts k n) == sort [ p | p <- nonCrossingPartitions n , numberOfParts p == k ]
+--
+nonCrossingPartitionsWithKParts 
+  :: Int   -- ^ @k@ = number of parts 
+  -> Int   -- ^ @n@ = size of the set
+  -> [NonCrossing]
+nonCrossingPartitionsWithKParts k n = map dyckPathToNonCrossingPartition $ peakingDyckPaths k n
+
+-- | Non-crossing partitions are counted by the Catalan numbers
+countNonCrossingPartitions :: Int -> Integer
+countNonCrossingPartitions = countDyckPaths
+
+-- | Non-crossing partitions with @k@ parts are counted by the Naranaya numbers
+countNonCrossingPartitionsWithKParts 
+  :: Int   -- ^ @k@ = number of parts 
+  -> Int   -- ^ @n@ = size of the set
+  -> Integer
+countNonCrossingPartitionsWithKParts = countPeakingDyckPaths
+
+--------------------------------------------------------------------------------
+
+-- | Uniformly random non-crossing partition
+randomNonCrossingPartition :: RandomGen g => Int -> g -> (NonCrossing,g)
+randomNonCrossingPartition n g0 = (dyckPathToNonCrossingPartition dyck, g1) where
+  (dyck,g1) = randomDyckPath n g0
+
+--------------------------------------------------------------------------------

@@ -0,0 +1,99 @@+
+-- | Set partitions.
+--
+-- See eg. <http://en.wikipedia.org/wiki/Partition_of_a_set>
+-- 
+
+{-# LANGUAGE BangPatterns #-}
+module Math.Combinat.Partitions.Set where
+
+--------------------------------------------------------------------------------
+
+import Data.List
+import Data.Ord
+
+import System.Random
+
+import Math.Combinat.Sets
+import Math.Combinat.Numbers
+import Math.Combinat.Helper
+import Math.Combinat.Partitions ( HasNumberOfParts(..) )
+
+--------------------------------------------------------------------------------
+-- * The type of set partitions
+
+-- | A partition of the set @[1..n]@ (in standard order)
+newtype SetPartition = SetPartition [[Int]] deriving (Eq,Ord,Show,Read)
+
+_standardizeSetPartition :: [[Int]] -> [[Int]]
+_standardizeSetPartition = sortBy (comparing myhead) . map sort where
+  myhead xs = case xs of
+    (x:xs) -> x
+    []     -> error "_standardizeSetPartition: empty subset"
+
+fromSetPartition :: SetPartition -> [[Int]]
+fromSetPartition (SetPartition zzs) = zzs
+
+toSetPartitionUnsafe :: [[Int]] -> SetPartition
+toSetPartitionUnsafe = SetPartition
+
+toSetPartition :: [[Int]] -> SetPartition
+toSetPartition zzs = if _isSetPartition zzs
+  then SetPartition (_standardizeSetPartition zzs)
+  else error "toSetPartition: not a set partition"
+
+_isSetPartition :: [[Int]] -> Bool
+_isSetPartition zzs = sort (concat zzs) == [1..n] where 
+  n = sum' (map length zzs)
+
+instance HasNumberOfParts SetPartition where
+  numberOfParts (SetPartition p) = length p
+
+--------------------------------------------------------------------------------
+-- * Generating set partitions
+
+-- | Synonym for 'setPartitionsNaive'
+setPartitions :: Int -> [SetPartition]
+setPartitions = setPartitionsNaive
+
+-- | Synonym for 'setPartitionsWithKPartsNaive'
+--
+-- > sort (setPartitionsWithKParts k n) == sort [ p | p <- setPartitions n , numberOfParts p == k ]
+-- 
+setPartitionsWithKParts   
+  :: Int    -- ^ @k@ = number of parts
+  -> Int    -- ^ @n@ = size of the set
+  -> [SetPartition]
+setPartitionsWithKParts = setPartitionsWithKPartsNaive
+
+-- | List all set partitions of @[1..n]@, naive algorithm
+setPartitionsNaive :: Int -> [SetPartition]
+setPartitionsNaive n = map (SetPartition . _standardizeSetPartition) $ go [1..n] where
+  go :: [Int] -> [[[Int]]]
+  go []     = [[]]
+  go (z:zs) = [ s : rest | k <- [1..n] , s0 <- choose (k-1) zs , let s = z:s0 , rest <- go (zs \\ s) ]
+
+-- | Set partitions of the set @[1..n]@ into @k@ parts
+setPartitionsWithKPartsNaive 
+  :: Int    -- ^ @k@ = number of parts
+  -> Int    -- ^ @n@ = size of the set
+  -> [SetPartition]
+setPartitionsWithKPartsNaive k n = map (SetPartition . _standardizeSetPartition) $ go k [1..n] where
+  go :: Int -> [Int] -> [[[Int]]]
+  go !k []     = if k==0 then [[]] else []
+  go  1 zs     = [[zs]]
+  go !k (z:zs) = [ s : rest | l <- [1..n] , s0 <- choose (l-1) zs , let s = z:s0 , rest <- go (k-1) (zs \\ s) ]
+
+
+-- | Set partitions are counted by the Bell numbers
+countSetPartitions :: Int -> Integer
+countSetPartitions = bellNumber 
+
+-- | Set partitions of size @k@ are counted by the Stirling numbers of second kind
+countSetPartitionsWithKParts 
+  :: Int    -- ^ @k@ = number of parts
+  -> Int    -- ^ @n@ = size of the set
+  -> Integer
+countSetPartitionsWithKParts k n = stirling2nd n k
+
+--------------------------------------------------------------------------------

@@ -228,7 +228,7 @@ countNestedParentheses = countBinaryTrees  -- | Generates all sequences of nested parentheses of length 2n.--- Order is lexigraphic (when right parentheses are considered +-- Order is lexicographical (when right parentheses are considered  -- smaller then left ones). -- Based on \"Algorithm P\" in Knuth, but less efficient because of -- the \"idiomatic\" code.@@ -249,6 +249,7 @@ 	      LeftParen  ->  	        {- debug ("---",reverse ls,l,r,rs) $ -} 	        findj ( lls , [] ) ( reverse (RightParen:rs) , [] ) +  next _ = error "fasc4A_algorithm_P: fatal error shouldn't happen"    findj :: ([Paren],[Paren]) -> ([Paren],[Paren]) -> Maybe ([Paren],[Paren])   findj ( [] , _ ) _ = Nothing@@ -259,6 +260,7 @@ 	      (a:_:as) -> findj ( ls, RightParen:rs ) ( as , LeftParen:a:ys ) 	      _ -> findj ( lls, [] ) ( reverse rs ++ xs , ys)  	    RightParen -> Just ( reverse ys ++ xs ++ reverse (LeftParen:rs) ++ ls , [] )+  findj _ _ = error "fasc4A_algorithm_P: fatal error shouldn't happen"      -- | Generates a uniformly random sequence of nested parentheses of length 2n.     -- Based on \"Algorithm W\" in Knuth.

@@ -1,5 +1,5 @@ Name:                combinat-Version:             0.2.6.1+Version:             0.2.6.2 Synopsis:            Generation of various combinatorial objects. Description:         A collection of functions to generate (and if there is                       a formula, count) combinatorial objects like partitions, @@ -45,6 +45,8 @@                        Math.Combinat.Tuples                         Math.Combinat.Compositions                        Math.Combinat.Partitions+                       Math.Combinat.Partitions.Set+                       Math.Combinat.Partitions.NonCrossing                        Math.Combinat.Permutations                        Math.Combinat.Tableaux                        Math.Combinat.Tableaux.Kostka@@ -57,7 +59,7 @@                        Math.Combinat.Helper    Extensions:          CPP, MultiParamTypeClasses, ScopedTypeVariables, -                       GeneralizedNewtypeDeriving +                       GeneralizedNewtypeDeriving, BangPatterns     Hs-Source-Dirs:      .