diff --git a/Math/Combinat.hs b/Math/Combinat.hs --- a/Math/Combinat.hs +++ b/Math/Combinat.hs @@ -1,17 +1,24 @@ -- | A collection of functions to generate combinatorial -- objects like partitions, compositions, permutations, --- Young tableaux, various trees, etc. +-- Young tableaux, various trees, etc etc. -- -- -- The long-term goals are -- --- (1) to be efficient; +-- (1) generate most of the standard structures; +-- +-- (2) while being efficient; -- --- (2) to be able to enumerate the structures --- with constant memory usage. +-- (3) to be able to enumerate the structures +-- with constant memory usage; -- --- The short-term goal is to generate +-- (4) and to be able to randomly sample from them. +-- +-- (5) finally, be a repository of algorithms +-- +-- +-- The short-term goal is simply to generate -- many interesting structures. -- -- @@ -40,6 +47,7 @@ , module Math.Combinat.Permutations , module Math.Combinat.Tableaux , module Math.Combinat.Trees + , module Math.Combinat.LatticePaths , module Math.Combinat.Graphviz ) where @@ -51,4 +59,5 @@ import Math.Combinat.Permutations import Math.Combinat.Tableaux import Math.Combinat.Trees +import Math.Combinat.LatticePaths import Math.Combinat.Graphviz diff --git a/Math/Combinat/Compositions.hs b/Math/Combinat/Compositions.hs --- a/Math/Combinat/Compositions.hs +++ b/Math/Combinat/Compositions.hs @@ -6,10 +6,19 @@ module Math.Combinat.Compositions where -import Math.Combinat.Numbers (factorial,binomial) +-------------------------------------------------------------------------------- -------------------------------------------------------- +import System.Random +import Math.Combinat.Sets ( randomChoice ) +import Math.Combinat.Numbers ( factorial , binomial ) +import Math.Combinat.Helper + +-------------------------------------------------------------------------------- +-- * generating all compositions + +-- | A /composition/ of an integer @n@ into @k@ parts is an ordered @k@-tuple of nonnegative (sometimes positive) integers +-- whose sum is @n@. type Composition = [Int] -- | Compositions fitting into a given shape and having a given degree. @@ -41,7 +50,7 @@ allCompositions' :: [Int] -> [[Composition]] allCompositions' shape = map (compositions' shape) [0..d] where d = sum shape --- | Compositions of a given length. +-- | Nonnegative compositions of a given length. compositions :: Integral a => a -- ^ length @@ -72,4 +81,29 @@ countCompositions1 :: Integral a => a -> a -> Integer countCompositions1 len d = countCompositions len (d-len) -------------------------------------------------------- +-------------------------------------------------------------------------------- +-- * random compositions + +-- | @randomComposition k n@ returns a uniformly random composition +-- of the number @n@ as an (ordered) sum of @k@ /nonnegative/ numbers +randomComposition :: RandomGen g => Int -> Int -> g -> ([Int],g) +randomComposition k n g0 = + if k<1 || n<0 + then error "randomComposition: k should be positive, and n should be nonnegative" + else (comp, g1) + where + (cs,g1) = randomChoice (k-1) (n+k-1) g0 + comp = pairsWith (\x y -> y-x-1) (0 : cs ++ [n+k]) + +-- | @randomComposition1 k n@ returns a uniformly random composition +-- of the number @n@ as an (ordered) sum of @k@ /positive/ numbers +randomComposition1 :: RandomGen g => Int -> Int -> g -> ([Int],g) +randomComposition1 k n g0 = + if k<1 || n a +unGen g = case g of + Gen x -> x + Inv x -> x + -- | A /word/, describing (non-uniquely) an element of a group. -- The identity element is represented (among others) by the empty word. type Word a = [Generator a] -------------------------------------------------------------------------------- + +-- | Generators are shown as small letters: @a@, @b@, @c@, ... +-- and their inverses are shown as capital letters, so @A=a^-1@, @B=b^-1@, etc. +showGen :: Generator Int -> Char +showGen (Gen i) = chr (96+i) +showGen (Inv i) = chr (64+i) + +showWord :: Word Int -> String +showWord = map showGen + +-------------------------------------------------------------------------------- instance Functor Generator where fmap f g = case g of @@ -64,6 +88,45 @@ go 0 = [[]] go n = [ x:xs | xs <- go (n-1) , x <- elems ] elems = [ Gen a | a<-[1..g] ] + +-------------------------------------------------------------------------------- +-- * Random words + +-- | A random group generator (or its inverse) between @1@ and @g@ +randomGenerator + :: RandomGen g + => Int -- ^ @g@ = number of generators + -> g -> (Generator Int, g) +randomGenerator d g0 = (gen,g2) where + (b,g1) = random g0 + (k,g2) = randomR (1,d) g1 + gen = if b then Gen k else Inv k + +-- | A random group generator (but never its inverse) between @1@ and @g@ +randomGeneratorNoInv + :: RandomGen g + => Int -- ^ @g@ = number of generators + -> g -> (Generator Int, g) +randomGeneratorNoInv d g0 = (Gen k,g1) where + (k,g1) = randomR (1,d) g0 + +-- | A random word of length @n@ using @g@ generators (or their inverses) +randomWord + :: RandomGen g + => Int -- ^ @g@ = number of generators + -> Int -- ^ @n@ = length of the word + -> g -> (Word Int, g) +randomWord d n g0 = (word,g1) where + (g1,word) = mapAccumL (\g _ -> swap (randomGenerator d g)) g0 [1..n] + +-- | A random word of length @n@ using @g@ generators (but not their inverses) +randomWordNoInv + :: RandomGen g + => Int -- ^ @g@ = number of generators + -> Int -- ^ @n@ = length of the word + -> g -> (Word Int, g) +randomWordNoInv d n g0 = (word,g1) where + (g1,word) = mapAccumL (\g _ -> swap (randomGeneratorNoInv d g)) g0 [1..n] -------------------------------------------------------------------------------- -- * The free group on @g@ generators diff --git a/Math/Combinat/Helper.hs b/Math/Combinat/Helper.hs --- a/Math/Combinat/Helper.hs +++ b/Math/Combinat/Helper.hs @@ -1,26 +1,57 @@ +-- | Miscellaneous helper functions + module Math.Combinat.Helper where +-------------------------------------------------------------------------------- + import Control.Monad import Data.List import Data.Ord -import qualified Data.Set as Set +import Data.Set (Set) ; import qualified Data.Set as Set +import Data.Map (Map) ; import qualified Data.Map as Map + import Debug.Trace -------------------------------------------------------------------------------- +-- * debugging debug :: Show a => a -> b -> b debug x y = trace ("-- " ++ show x ++ "\n") y -{-# SPECIALIZE swap :: (a,a) -> (a,a) #-} +-------------------------------------------------------------------------------- +-- * pairs + +{-# SPECIALIZE swap :: (a ,a ) -> (a ,a ) #-} {-# SPECIALIZE swap :: (Int,Int) -> (Int,Int) #-} 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 + go _ = [] + +pairsWith :: (a -> a -> b) -> [a] -> [b] +pairsWith f = go where + go (x:xs@(y:_)) = f x y : go xs + go _ = [] + +-------------------------------------------------------------------------------- +-- * equality and ordering + equating :: Eq b => (a -> b) -> a -> a -> Bool equating f x y = (f x == f y) @@ -43,6 +74,7 @@ | otherwise = x : worker (Set.insert x s) xs -------------------------------------------------------------------------------- +-- * first \/ last -- | The boolean argument will @True@ only for the last element mapWithLast :: (Bool -> a -> b) -> [a] -> [b] @@ -60,7 +92,7 @@ go b (x : xs) = f b False x : go False xs -------------------------------------------------------------------------------- --- helpers for ASCII drawing +-- * helpers for ASCII drawing -- | extend lines with spaces so that they have the same line mkLinesUniformWidth :: [String] -> [String] @@ -85,18 +117,25 @@ vConcatLines = concat -------------------------------------------------------------------------------- +-- * counting -- helps testing the random rutines count :: Eq a => a -> [a] -> Int count x xs = length $ filter (==x) xs +histogram :: (Eq a, Ord a) => [a] -> [(a,Int)] +histogram xs = Map.toList table where + table = Map.fromListWith (+) [ (x,1) | x<-xs ] + -------------------------------------------------------------------------------- +-- * maybe fromJust :: Maybe a -> a fromJust (Just x) = x fromJust Nothing = error "fromJust: Nothing" -------------------------------------------------------------------------------- +-- * bool intToBool :: Int -> Bool intToBool 0 = False @@ -108,6 +147,7 @@ boolToInt True = 1 -------------------------------------------------------------------------------- +-- * iteration -- iterated function application nest :: Int -> (a -> a) -> a -> a diff --git a/Math/Combinat/LatticePaths.hs b/Math/Combinat/LatticePaths.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/LatticePaths.hs @@ -0,0 +1,199 @@ + +-- | Dyck paths, lattice paths, etc + +{-# LANGUAGE BangPatterns #-} +module Math.Combinat.LatticePaths where + +-------------------------------------------------------------------------------- + +import Math.Combinat.Numbers +import Math.Combinat.Trees.Binary + +-------------------------------------------------------------------------------- +-- * Types + +-- | A step in a lattice path +data Step + = UpStep -- ^ the step @(1,1)@ + | DownStep -- ^ the step @(1,-1)@ + deriving (Eq,Ord,Show) + +-- | A lattice path is a path using only the allowed steps, never going below the zero level line @y=0@. +-- +-- Note that if you rotate such a path by 45 degrees counterclockwise, +-- you get a path which uses only the steps @(1,0)@ and @(0,1)@, and stays +-- above the main diagonal (hence the name, we just use a different convention). +-- +type LatticePath = [Step] + +-------------------------------------------------------------------------------- + +-- | A lattice path is called \"valid\", if it never goes below the @y=0@ line. +isValidPath :: LatticePath -> Bool +isValidPath = 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 + go !h !y [] = h + go !h !y (t:ts) = case t of + UpStep -> go (max h (y+1)) (y+1) ts + DownStep -> go h (y-1) ts + +-- | Endpoint of a lattice path, which starts from @(0,0)@. +pathEndpoint :: LatticePath -> (Int,Int) +pathEndpoint = go 0 0 where + go !x !y [] = (x,y) + go !x !y (t:ts) = case t of + UpStep -> go (x+1) (y+1) ts + DownStep -> go (x+1) (y-1) ts + +-- | Returns the coordinates of the path (excluding the starting point @(0,0)@, but including +-- the endpoint +pathCoordinates :: LatticePath -> [(Int,Int)] +pathCoordinates = go 0 0 where + go _ _ [] = [] + go !x !y (t:ts) = let x' = x + 1 + y' = case t of { UpStep -> y+1 ; DownStep -> y-1 } + in (x',y') : go x' y' ts + +-- | 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 +pathNumberOfZeroTouches = pathNumberOfTouches' 0 + +-- | Number of points on the path which touch the level line at height @h@ +-- (excluding the starting point @(0,0)@, but including the endpoint). +pathNumberOfTouches' + :: Int -- ^ @h@ = the touch level + -> LatticePath -> Int +pathNumberOfTouches' h = go 0 0 0 where + go !cnt _ _ [] = cnt + go !cnt !x !y (t:ts) = let y' = case t of { UpStep -> y+1 ; DownStep -> y-1 } + cnt' = if y'==h then cnt+1 else cnt + in go cnt' (x+1) y' ts + + +-------------------------------------------------------------------------------- +-- * Dyck paths + +-- | @dyckPaths m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@. +-- +-- Remark: Dyck paths are obviously in bijection with nested parentheses, and thus +-- also with binary trees. +-- +dyckPaths :: Int -> [LatticePath] +dyckPaths = map (map f) . nestedParentheses where + f p = case p of { LeftParen -> UpStep ; RightParen -> 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 + +-------------------------------------------------------------------------------- +-- * Bounded Dyck paths + +-- | @boundedDyckPaths h m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ whose height is at most @h@. +-- Synonym for 'boundedDyckPathsNaive'. +-- +boundedDyckPaths + :: Int -- ^ @h@ = maximum height + -> Int -- ^ @m@ = half-length + -> [LatticePath] +boundedDyckPaths = boundedDyckPathsNaive + +-- | @boundedDyckPathsNaive h m@ lists all Dyck paths from @(0,0)@ to @(2m,0)@ whose height is at most @h@. +-- +-- > sort (boundedDyckPaths h m) == sort [ p | p <- dyckPaths m , pathHeight p <= h ] +-- > sort (boundedDyckPaths m m) == sort (dyckPaths m) +-- +-- Naive recursive algorithm, resulting order is pretty ad-hoc. +-- +boundedDyckPathsNaive + :: Int -- ^ @h@ = maximum height + -> Int -- ^ @m@ = half-length + -> [LatticePath] +boundedDyckPathsNaive = worker where + worker !h !m + | h<0 = [] + | m<0 = [] + | m==0 = [[]] + | h<=0 = [] + | otherwise = as ++ bs + where + bracket p = UpStep : p ++ [DownStep] + as = [ bracket p | p <- boundedDyckPaths (h-1) (m-1) ] + bs = [ bracket p ++ q | k <- [1..m-1] , p <- boundedDyckPaths (h-1) (k-1) , q <- boundedDyckPaths h (m-k) ] + +-------------------------------------------------------------------------------- +-- * More general lattice paths + +-- | All lattice paths from @(0,0)@ to @(x,y)@. Clearly empty unless @x-y@ is even. +-- Synonym for 'latticePathsNaive' +-- +latticePaths :: (Int,Int) -> [LatticePath] +latticePaths = latticePathsNaive + +-- | All lattice paths from @(0,0)@ to @(x,y)@. Clearly empty unless @x-y@ is even. +-- +-- Note that +-- +-- > sort (dyckPaths n) == sort (latticePaths (0,2*n)) +-- +-- Naive recursive algorithm, resulting order is pretty ad-hoc. +-- +latticePathsNaive :: (Int,Int) -> [LatticePath] +latticePathsNaive (x,y) = worker x y where + worker !x !y + | odd (x-y) = [] + | x<0 = [] + | y<0 = [] + | y==0 = dyckPaths (div x 2) + | x==1 && y==1 = [[UpStep]] + | otherwise = as ++ bs + where + bracket p = UpStep : p ++ [DownStep] + as = [ UpStep : p | p <- worker (x-1) (y-1) ] + bs = [ bracket p ++ q | k <- [1..(div x 2)] , p <- dyckPaths (k-1) , q <- worker (x-2*k) y ] + +-------------------------------------------------------------------------------- +-- * Zero-level touches + +-- | @touchingDyckPathsNaive 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 + :: Int -- ^ @k@ = number of touches + -> Int -- ^ @m@ = half-length + -> [LatticePath] +touchingDyckPaths = touchingDyckPathsNaive + + +-- | @touchingDyckPathsNaive 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). +-- +-- > sort (touchingDyckPathsNaive k m) == sort [ p | p <- dyckPaths m , pathNumberOfZeroTouches p == k ] +-- +-- Naive recursive algorithm, resulting order is pretty ad-hoc. +-- +touchingDyckPathsNaive + :: Int -- ^ @k@ = number of touches + -> Int -- ^ @m@ = half-length + -> [LatticePath] +touchingDyckPathsNaive = worker where + worker !k !m + | m == 0 = if k==0 then [[]] else [] + | k <= 0 = [] + | m < 0 = [] + | k == 1 = [ bracket p | p <- dyckPaths (m-1) ] + | otherwise = [ bracket p ++ q | l <- [1..m-1] , p <- dyckPaths (l-1) , q <- worker (k-1) (m-l) ] + where + bracket p = UpStep : p ++ [DownStep] + +-------------------------------------------------------------------------------- + diff --git a/Math/Combinat/Numbers.hs b/Math/Combinat/Numbers.hs --- a/Math/Combinat/Numbers.hs +++ b/Math/Combinat/Numbers.hs @@ -10,6 +10,8 @@ import Data.Array +import Math.Combinat.Helper ( sum' ) + -------------------------------------------------------------------------------- -- | @(-1)^k@ @@ -102,6 +104,9 @@ -- | (Signed) Stirling numbers of the first kind. OEIS:A008275. -- This function uses 'signedStirling1stArray', so it shouldn't be used -- to compute /many/ Stirling numbers. +-- +-- Argument order: @signedStirling1st n k@ +-- signedStirling1st :: Integral a => a -> a -> Integer signedStirling1st n k | k < 1 = 0 @@ -114,6 +119,9 @@ -- | Stirling numbers of the second kind. OEIS:A008277. -- This function uses an explicit formula. +-- +-- Argument order: @stirling2nd n k@ +-- stirling2nd :: Integral a => a -> a -> Integer stirling2nd n k | k < 1 = 0 @@ -138,5 +146,32 @@ / toRational (k+1) -------------------------------------------------------------------------------- +-- * Bell numbers + +-- | Bell numbers (Sloane's A000110) from B(0) up to B(n). B(0)=B(1)=1, B(2)=2, etc. +-- +-- The Bell numbers count the number of /set partitions/ of a set of size @n@ +-- +-- See +-- +bellNumbersArray :: Integral a => a -> Array Int Integer +bellNumbersArray nn = arr where + arr = array (0::Int,n) kvs + n = fromIntegral nn :: Int + kvs = (0,1) : [ (k, f k) | k<-[1..n] ] + f n = sum' [ binomial (n-1) k * arr ! k | k<-[0..n-1] ] + +-- | The n-th Bell number B(n), using the Stirling numbers of the second kind. +-- This may be slower than using 'bellNumbersArray'. +bellNumber :: Integral a => a -> Integer +bellNumber nn + | n < 0 = error "bellNumber: expecting a nonnegative index" + | n == 0 = 1 + | otherwise = sum' [ stirling2nd n k | k<-[1..n] ] + where + n = fromIntegral nn :: Int + +-------------------------------------------------------------------------------- + diff --git a/Math/Combinat/Partitions.hs b/Math/Combinat/Partitions.hs --- a/Math/Combinat/Partitions.hs +++ b/Math/Combinat/Partitions.hs @@ -1,9 +1,10 @@ --- | Partitions. Partitions are nonincreasing sequences of positive integers. +-- | 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. - +-- module Math.Combinat.Partitions ( -- * Type and basic stuff Partition @@ -33,6 +34,13 @@ , allPartitions , countAllPartitions' , countAllPartitions + -- * Ferrer diagrams + , printFerrerDiagram + , ferrerDiagram + , ferrerDiagramEnglishNotation + , ferrerDiagramFrenchNotation + , ferrerDiagramEnglishNotation' + , ferrerDiagramFrenchNotation' -- * Paritions of multisets, vector partitions , partitionMultiset , IntVector @@ -50,7 +58,7 @@ -------------------------------------------------------------------------------- --- | The additional invariant enforced here is that partitions +-- | A partition of an integer. The additional invariant enforced here is that partitions -- are monotone decreasing sequences of positive integers. newtype Partition = Partition [Int] deriving (Eq,Ord,Show,Read) @@ -62,7 +70,7 @@ toPartitionUnsafe :: [Int] -> Partition toPartitionUnsafe = Partition --- | Checks whether the input is a partition. See the note at 'isPartition'! +-- | Checks whether the input is an integer partition. See the note at 'isPartition'! toPartition :: [Int] -> Partition toPartition xs = if isPartition xs then toPartitionUnsafe xs @@ -119,7 +127,7 @@ _elements :: [Int] -> [(Int,Int)] _elements shape = [ (i,j) | (i,l) <- zip [1..] shape, j<-[1..l] ] --- | Computes the number of \"automorphisms\" of a given partition. +-- | Computes the number of \"automorphisms\" of a given integer partition. countAutomorphisms :: Partition -> Integer countAutomorphisms = _countAutomorphisms . fromPartition @@ -128,7 +136,7 @@ --------------------------------------------------------------------------------- --- | Partitions of d, fitting into a given rectangle, as lists. +-- | Integer partitions of @d@, fitting into a given rectangle, as lists. _partitions' :: (Int,Int) -- ^ (height,width) -> Int -- ^ d @@ -153,24 +161,24 @@ countPartitions' (h,w) d = sum [ countPartitions' (i,w-1) (d-i) | i <- [1..min d h] ] --- | Partitions of d, as lists +-- | Partitions of @d@, as lists _partitions :: Int -> [[Int]] _partitions d = _partitions' (d,d) d --- | Partitions of d. +-- | Partitions of @d@. partitions :: Int -> [Partition] partitions d = partitions' (d,d) d countPartitions :: Int -> Integer countPartitions d = countPartitions' (d,d) d --- | All partitions fitting into a given rectangle. +-- | All integer partitions fitting into a given rectangle. allPartitions' :: (Int,Int) -- ^ (height,width) -> [[Partition]] allPartitions' (h,w) = [ partitions' (h,w) i | i <- [0..d] ] where d = h*w --- | All partitions up to a given degree. +-- | All integer partitions up to a given degree (that is, all integer partitions whose sum is less or equal to @d@) allPartitions :: Int -> [[Partition]] allPartitions d = [ partitions i | i <- [0..d] ] @@ -182,6 +190,32 @@ countAllPartitions :: Int -> Integer countAllPartitions d = sum [ countPartitions i | i <- [0..d] ] + +-------------------------------------------------------------------------------- +-- * Ferrer diagrams + +printFerrerDiagram :: Partition -> IO () +printFerrerDiagram = putStrLn . ferrerDiagram + +-- | Synonym for 'ferrerDiagramEnglishNotation' +ferrerDiagram :: Partition -> String +ferrerDiagram = ferrerDiagramEnglishNotation + +ferrerDiagramEnglishNotation :: Partition -> String +ferrerDiagramEnglishNotation = ferrerDiagramEnglishNotation' '@' + +ferrerDiagramFrenchNotation :: Partition -> String +ferrerDiagramFrenchNotation = ferrerDiagramFrenchNotation' '@' + +ferrerDiagramEnglishNotation' :: Char -> Partition -> String +ferrerDiagramEnglishNotation' ch part = unlines (map f ys) where + ys = fromPartition part + f n = replicate n ch + +ferrerDiagramFrenchNotation' :: Char -> Partition -> String +ferrerDiagramFrenchNotation' ch part = unlines (map f ys) where + ys = reverse $ fromPartition $ dualPartition part + f n = replicate n ch -------------------------------------------------------------------------------- diff --git a/Math/Combinat/Permutations.hs b/Math/Combinat/Permutations.hs --- a/Math/Combinat/Permutations.hs +++ b/Math/Combinat/Permutations.hs @@ -51,7 +51,7 @@ #ifdef QUICKCHECK -- * QuickCheck , checkAll -#endif QUICKCHECK +#endif ) where @@ -66,6 +66,7 @@ import Data.Array import Data.Array.ST +import Data.Array.Unsafe import Math.Combinat.Helper import Math.Combinat.Numbers (factorial,binomial) diff --git a/Math/Combinat/Sets.hs b/Math/Combinat/Sets.hs --- a/Math/Combinat/Sets.hs +++ b/Math/Combinat/Sets.hs @@ -3,31 +3,56 @@ module Math.Combinat.Sets ( - choose + -- * choices + choose , choose_ , combine , compose + -- * tensor products , tuplesFromList , listTensor - -- + -- * sublists , kSublists , sublists , countKSublists , countSublists + -- * random choice + , randomChoice ) where -import Math.Combinat.Numbers (binomial) +-------------------------------------------------------------------------------- +import Data.List ( sort , mapAccumL ) +import System.Random + +-- import Data.Map (Map) +-- import qualified Data.Map as Map + +import Math.Combinat.Numbers ( binomial ) +import Math.Combinat.Helper ( swap ) + -------------------------------------------------------------------------------- +-- * choices -- | All possible ways to choose @k@ elements from a list, without --- repetitions. \"Antisymmetric power\" for lists. Synonym for "kSublists". +-- repetitions. \"Antisymmetric power\" for lists. Synonym for 'kSublists'. choose :: Int -> [a] -> [[a]] choose 0 _ = [[]] choose k [] = [] choose k (x:xs) = map (x:) (choose (k-1) xs) ++ choose k xs +-- | @choose_ k n@ returns all possible ways of choosing @k@ disjoint elements from @[1..n]@ +-- +-- > choose_ k n == choose k [1..n] +-- +choose_ :: Int -> Int -> [[Int]] +choose_ k n = if n<0 || k<0 + then error "choose_: n and k should nonnegative" + else if k>n || k<0 + then [] + else choose k [1..n] + -- | All possible ways to choose @k@ elements from a list, /with repetitions/. --- \"Symmetric power\" for lists. See also "Math.Combinat.Combinations". +-- \"Symmetric power\" for lists. See also "Math.Combinat.Compositions". -- TODO: better name? combine :: Int -> [a] -> [[a]] combine 0 _ = [[]] @@ -38,6 +63,9 @@ compose :: Int -> [a] -> [[a]] compose = combine +-------------------------------------------------------------------------------- +-- * tensor products + -- | \"Tensor power\" for lists. Special case of 'listTensor': -- -- > tuplesFromList k xs == listTensor (replicate k xs) @@ -58,6 +86,7 @@ --listTensor (xs:xss) = [ y:ys | y <- xs, ys <- listTensor xss ] -------------------------------------------------------------------------------- +-- * sublists -- | Sublists of a list having given number of elements. kSublists :: Int -> [a] -> [[a]] @@ -79,3 +108,43 @@ countSublists n = 2 ^ n -------------------------------------------------------------------------------- +-- * random choice + +-- | @randomChoice k n@ returns a uniformly random choice of @k@ elements from the set @[1..n]@ +-- +-- Example: +-- +-- > do +-- > cs <- replicateM 10000 (getStdRandom (randomChoice 3 7)) +-- > mapM_ print $ histogram cs +-- +randomChoice :: RandomGen g => Int -> Int -> g -> ([Int],g) +randomChoice k n g0 = + if k>n || k<0 + then error "randomChoice: k out of range" + else (makeChoiceFromIndices as, g1) + where + -- choose numbers from [1..n], [1..n-1], [1..n-2] etc + (g1,as) = mapAccumL (\g m -> swap (randomR (1,m) g)) g0 [n,n-1..n-k+1] + +-- | From a list of $k$ numbers, where the first is in the interval @[1..n]@, +-- the second in @[1..n-1]@, the third in @[1..n-2]@, we create a size @k@ subset of @n@. +makeChoiceFromIndices :: [Int] -> [Int] +makeChoiceFromIndices = sort . go [] where + + go :: [Int] -> [Int] -> [Int] + go acc (b:bs) = b' : go (insert b' acc) bs where b' = skip b acc + go _ [] = [] + + -- skip over the already occupied positions. Second argument should be a sorted list + skip :: Int -> [Int] -> Int + skip x (y:ys) = if x>=y then skip (x+1) ys else x + skip x [] = x + + -- insert into a sorted list + insert :: Int -> [Int] -> [Int] + insert x (y:ys) = if x<=y then x:y:ys else y : insert x ys + insert x [] = [x] + +-------------------------------------------------------------------------------- + diff --git a/Math/Combinat/Trees/Binary.hs b/Math/Combinat/Trees/Binary.hs --- a/Math/Combinat/Trees/Binary.hs +++ b/Math/Combinat/Trees/Binary.hs @@ -50,6 +50,7 @@ import Data.Array import Data.Array.ST +import Data.Array.Unsafe import Data.List import Data.Tree (Tree(..),Forest(..)) @@ -207,7 +208,11 @@ -------------------------------------------------------------------------------- -- * Nested parentheses --- | Synonym for 'fasc4A_algorithm_P'. +-- | Generates all sequences of nested parentheses of length @2n@ in +-- lexigraphic order. +-- +-- Synonym for 'fasc4A_algorithm_P'. +-- nestedParentheses :: Int -> [[Paren]] nestedParentheses = fasc4A_algorithm_P diff --git a/combinat.cabal b/combinat.cabal --- a/combinat.cabal +++ b/combinat.cabal @@ -1,5 +1,5 @@ Name: combinat -Version: 0.2.5.0 +Version: 0.2.6.0 Synopsis: Generation of various combinatorial objects. Description: A collection of functions to generate (and if there is a formula, count) combinatorial objects like partitions, @@ -13,7 +13,7 @@ Homepage: http://code.haskell.org/~bkomuves/ Stability: Experimental Category: Math -Tested-With: GHC == 7.4.2 +Tested-With: GHC == 7.8.3 Cabal-Version: >= 1.6 Build-Type: Simple @@ -21,27 +21,22 @@ Description: Compile with the QuickCheck tests. default: False -Flag splitBase - Description: Choose the new smaller, split-up base package. - Flag base4 Description: Base v4 Library - if flag(splitBase) - if flag(base4) - Build-Depends: base >= 4 && < 5, array, containers, random, transformers - cpp-options: -DBASE_VERSION=4 - else - Build-Depends: base >= 3 && < 4, array, containers, random, transformers - cpp-options: -DBASE_VERSION=3 - if flag(withQuickCheck) - Build-Depends: QuickCheck - else - Build-Depends: base < 3 - cpp-options: -DBASE_VERSION=2 + if flag(base4) + Build-Depends: base >= 4 && < 5, array >= 0.4, containers, random, transformers + cpp-options: -DBASE_VERSION=4 + else + Build-Depends: base >= 3 && < 4, array >= 0.4, containers, random, transformers + cpp-options: -DBASE_VERSION=3 + if flag(withQuickCheck) + Build-Depends: QuickCheck + + Exposed-Modules: Math.Combinat Math.Combinat.Numbers Math.Combinat.Numbers.Series @@ -56,10 +51,10 @@ Math.Combinat.Trees Math.Combinat.Trees.Binary Math.Combinat.Trees.Nary + Math.Combinat.LatticePaths Math.Combinat.FreeGroups Math.Combinat.Graphviz - - Other-Modules: Math.Combinat.Helper + Math.Combinat.Helper Extensions: CPP, MultiParamTypeClasses, ScopedTypeVariables, GeneralizedNewtypeDeriving