combinat 0.1 → 0.2

raw patch · 7 files changed

+373/−93 lines, 7 filesdep +randomdep ~arraysetup-changedPVP ok

version bump matches the API change (PVP)

Dependencies added: random

Dependency ranges changed: array

API changes (from Hackage documentation)

Files

• Math/Combinat.hs +6/−1
• Math/Combinat/Helper.hs +10/−1
• Math/Combinat/Permutations.hs +195/−21
• Math/Combinat/Sets.hs +35/−0
• Math/Combinat/Trees.hs +118/−62
• Setup.lhs +2/−3
• combinat.cabal +7/−5

@@ -20,10 +20,14 @@ --  * underscore prefix: use plain lists instead of other types with  --    enforced invariants; --+--  * \"random\" prefix: generates random objects +--    (typically with uniform distribution); +-- --  * \"count\" prefix: counting functions.  module Math.Combinat -  ( module Math.Combinat.Tuples+  ( module Math.Combinat.Sets+  , module Math.Combinat.Tuples   , module Math.Combinat.Combinations   , module Math.Combinat.Partitions   , module Math.Combinat.Permutations@@ -31,6 +35,7 @@   , module Math.Combinat.Trees   ) where +import Math.Combinat.Sets import Math.Combinat.Tuples import Math.Combinat.Combinations import Math.Combinat.Partitions

@@ -4,7 +4,16 @@ import Debug.Trace  debug :: Show a => a -> b -> b-debug x y = trace (show x) y+debug x y = trace ("-- " ++ show x ++ "\n") y++{-# SPECIALIZE swap :: (a,a) -> (a,a) #-}+{-# SPECIALIZE swap :: (Int,Int) -> (Int,Int) #-}+swap :: (a,b) -> (b,a)+swap (x,y) = (y,x)++-- helps testing the random rutines +count :: Eq a => a -> [a] -> Int+count x xs = length $ filter (==x) xs  fromJust :: Maybe a -> a fromJust (Just x) = x

@@ -3,57 +3,231 @@ --   Donald E. Knuth: The Art of Computer Programming, vol 4, pre-fascicle 2B. -- {-# OPTIONS_GHC -fno-warn-name-shadowing #-}-module Math.Combinat.Permutations where+{-# LANGUAGE ScopedTypeVariables #-}+module Math.Combinat.Permutations +  ( -- * Types+    Permutation+  , fromPermutation+  , toPermutationUnsafe+  , isPermutation+  , toPermutation+  , permutationSize+    -- * Permutation groups+  , permute+  , permuteList+  , multiply+  , inverse+    -- * Simple permutations+  , permutations+  , _permutations+  , permutationsNaive+  , _permutationsNaive+  , countPermutations+    -- * Random permutations+  , randomPermutation+  , _randomPermutation+  , randomCyclicPermutation+  , _randomCyclicPermutation+  , randomPermutationDurstenfeld+  , randomCyclicPermutationSattolo+    -- * Multisets+  , permuteMultiset+  , countPermuteMultiset+  , fasc2B_algorithm_L+  ) +  where -import Data.List+import Control.Monad+import Control.Monad.ST++import Data.List hiding (permutations) import Data.Array+import Data.Array.ST  import Math.Combinat.Helper +import System.Random+ --------------------------------------------------------{- -- * Types --- | Standard notation for permutations+-- | Standard notation for permutations. Internally it is an array of the integers @[1..n]@.  newtype Permutation = Permutation (Array Int Int) deriving (Eq,Ord,Show,Read) +{- -- | Disjoint cycle notation for permutations newtype DisjCycles  = DisjCycles [[Int]] deriving (Eq,Ord,Show,Read) -} +fromPermutation :: Permutation -> [Int]+fromPermutation (Permutation ar) = elems ar++-- | Assumes that the input is a permutation of the numbers @[1..n]@.+toPermutationUnsafe :: [Int] -> Permutation+toPermutationUnsafe xs = Permutation perm where+  n = length xs+  perm = listArray (1,n) xs++-- | Checks whether the input is a permutation of the numbers @[1..n]@.+isPermutation :: [Int] -> Bool+isPermutation xs = (ar!0 == 0) && and [ ar!j == 1 | j<-[1..n] ] where+  n = length xs+  -- the zero index is an unidiomatic hack+  ar = accumArray (+) 0 (0,n) $ map f xs +  f :: Int -> (Int,Int)+  f j = if j<1 || j>n then (0,1) else (j,1)++-- | Checks the input.+toPermutation :: [Int] -> Permutation+toPermutation xs = if isPermutation xs +  then toPermutationUnsafe xs+  else error "toPermutation: not a permutation"++-- | Returns @n@, where the input is a permutation of the numbers @[1..n]@+permutationSize :: Permutation -> Int+permutationSize (Permutation ar) = snd $ bounds ar+ -------------------------------------------------------+-- * Permutation groups+    +-- | Action of a permutation on a set. If our permutation is +-- encoded with the sequence @[p1,p2,...,pn]@, then in the+-- two-line notation we have+--+-- > ( 1  2  3  ... n  )+-- > ( p1 p2 p3 ... pn )+--+-- We adopt the convention that permutations act /on the left/ +-- (as opposed to Knuth, where they act on the right).+-- Thus, +-- +-- > permute pi1 (permute pi2 set) == permute (pi1 `multiply` pi2) set+--   +-- The second argument should be an array with bounds @(1,n)@.+-- The function checks the array bounds.+permute :: Permutation -> Array Int a -> Array Int a    +permute pi@(Permutation perm) ar = +  if (a==1) && (b==n) +    then listArray (1,n) [ ar!(perm!i) | i <- [1..n] ] +    else error "permute: array bounds do not match"+  where+    (_,n) = bounds perm  +    (a,b) = bounds ar   ++-- | The list should be of length @n@.+permuteList :: Permutation -> [a] -> [a]    +permuteList perm xs = elems $ permute perm $ listArray (1,n) xs where+  n = permutationSize perm++-- | Multiplies two permutations together. See 'permute' for our+-- conventions.  +multiply :: Permutation -> Permutation -> Permutation+multiply pi1@(Permutation perm1) (Permutation perm2) = +  if (n==m) +    then Permutation result+    else error "multiply: permutations of different sets"  +  where+	  (_,n) = bounds perm1+	  (_,m) = bounds perm2    +	  result = permute pi1 perm2    +  +infixr 7 `multiply`  +    +-- | The inverse permutation+inverse :: Permutation -> Permutation    +inverse (Permutation perm1) = Permutation result+  where+    result = array (1,n) $ map swap $ assocs perm1+    (_,n) = bounds perm1++------------------------------------------------------- -- * Permutations of distinct elements --- | Permutations of [1..n] in lexicographic order, naive algorithm.-_permutations :: Int -> [[Int]]  -_permutations 0 = [[]]-_permutations 1 = [[1]]-_permutations n = helper [1..n] where+-- | A synonym for 'permutationsNaive'+permutations :: Int -> [Permutation]+permutations = permutationsNaive++_permutations :: Int -> [[Int]]+_permutations = _permutationsNaive++-- | Permutations of @[1..n]@ in lexicographic order, naive algorithm.+permutationsNaive :: Int -> [Permutation]+permutationsNaive n = map toPermutationUnsafe $ _permutations n ++_permutationsNaive :: Int -> [[Int]]  +_permutationsNaive 0 = [[]]+_permutationsNaive 1 = [[1]]+_permutationsNaive n = helper [1..n] where   helper [] = [[]]   helper xs = [ i : ys | i <- xs , ys <- helper (xs `minus` i) ]   minus [] _ = []   minus (x:xs) i = if x < i then x : minus xs i else xs--{--permutations :: Int -> [Permutation]-permutations n = map toPermutationUnsafe $ _permutations n --}-+           -- | # = n! countPermutations :: Int -> Integer countPermutations = factorial  -------------------------------------------------------+-- * Random permutations++-- | A synonym for 'randomPermutationDurstenfeld'.+randomPermutation :: RandomGen g => Int -> g -> (Permutation,g)+randomPermutation = randomPermutationDurstenfeld++_randomPermutation :: RandomGen g => Int -> g -> ([Int],g)+_randomPermutation n rndgen = (fromPermutation perm, rndgen') where+  (perm, rndgen') = randomPermutationDurstenfeld n rndgen ++-- | A synonym for 'randomCyclicPermutationSattolo'.+randomCyclicPermutation :: RandomGen g => Int -> g -> (Permutation,g)+randomCyclicPermutation = randomCyclicPermutationSattolo++_randomCyclicPermutation :: RandomGen g => Int -> g -> ([Int],g)+_randomCyclicPermutation n rndgen = (fromPermutation perm, rndgen') where+  (perm, rndgen') = randomCyclicPermutationSattolo n rndgen ++-- | Generates a uniformly random permutation of @[1..n]@.+-- Durstenfeld's algorithm (see <http://en.wikipedia.org/wiki/Knuth_shuffle>).+randomPermutationDurstenfeld :: RandomGen g => Int -> g -> (Permutation,g)+randomPermutationDurstenfeld = randomPermutationDurstenfeldSattolo False++-- | Generates a uniformly random /cyclic/ permutation of @[1..n]@.+-- Sattolo's algorithm (see <http://en.wikipedia.org/wiki/Knuth_shuffle>).+randomCyclicPermutationSattolo :: RandomGen g => Int -> g -> (Permutation,g)+randomCyclicPermutationSattolo = randomPermutationDurstenfeldSattolo True++randomPermutationDurstenfeldSattolo :: RandomGen g => Bool -> Int -> g -> (Permutation,g)+randomPermutationDurstenfeldSattolo isSattolo n rnd = res where+  res = runST $ do+    ar <- newArray_ (1,n) +    forM_ [1..n] $ \i -> writeArray ar i i+    rnd' <- worker n (if isSattolo then n-1 else n) rnd ar +    perm <- unsafeFreeze ar+    return (Permutation perm, rnd')+  worker :: RandomGen g => Int -> Int -> g -> STUArray s Int Int -> ST s g +  worker n m rnd ar = +    if n==1 +      then return rnd +      else do+        let (k,rnd') = randomR (1,m) rnd+        when (k /= n) $ do+          y <- readArray ar k +          z <- readArray ar n+          writeArray ar n y+          writeArray ar k z+        worker (n-1) (m-1) rnd' ar ++------------------------------------------------------- -- * Permutations of a multiset --- | Generates all permutations of a multiset. ---   The order is lexicographic.  -permute :: (Eq a, Ord a) => [a] -> [[a]] -permute = fasc2B_algorithm_L+-- | Generates all permutations of a multiset.  +--   The order is lexicographic. A synonym for 'fasc2B_algorithm_L'+permuteMultiset :: (Eq a, Ord a) => [a] -> [[a]] +permuteMultiset = fasc2B_algorithm_L  -- | # = \\frac { (\sum_i n_i) ! } { \\prod_i (n_i !) }    -countPermute :: (Eq a, Ord a) => [a] -> Integer-countPermute xs = factorial n `div` product [ factorial (length z) | z <- group ys ] +countPermuteMultiset :: (Eq a, Ord a) => [a] -> Integer+countPermuteMultiset xs = factorial n `div` product [ factorial (length z) | z <- group ys ]    where     ys = sort xs     n = length xs

@@ -0,0 +1,35 @@++-- | Subsets. ++module Math.Combinat.Sets +  ( kSublists+  , sublists+  , countKSublists+  , countSublists+  ) +  where++import Math.Combinat.Helper++-------------------------------------------------------++-- | Sublists of a list having given number of elements.+kSublists :: Int -> [a] -> [[a]]+kSublists 0 _  = [[]]+kSublists k [] = []+kSublists k (x:xs) = map (x:) (kSublists (k-1) xs) ++ kSublists k xs  ++-- | @# = \binom { n } { k }@.+countKSublists :: Int -> Int -> Integer+countKSublists k n = binomial (fromIntegral n) (fromIntegral k)++-- | All sublists of a list.+sublists :: [a] -> [[a]]+sublists [] = [[]]+sublists (x:xs) = map (x:) (sublists xs) ++ sublists xs ++-- | @# = 2^n@.+countSublists :: Int -> Integer+countSublists n = 2 ^ n++-------------------------------------------------------

@@ -6,6 +6,8 @@   ( -- * Types     BinTree(..)   , leaf+  , BinTree'(..)+  , forgetNodeDecorations   , module Data.Tree    , Paren(..)   , parenthesesToString@@ -21,22 +23,38 @@   , binaryTreeToNestedParentheses     -- * Nested parentheses   , nestedParentheses +  , randomNestedParentheses+  , nthNestedParentheses+  , countNestedParentheses   , fasc4A_algorithm_P+  , fasc4A_algorithm_W+  , fasc4A_algorithm_U     -- * Binary trees   , binaryTrees   , countBinaryTrees   , binaryTreesNaive+  , randomBinaryTree+  , fasc4A_algorithm_R   )    where +import Control.Monad+import Control.Monad.ST++import Data.Array+import Data.Array.ST+ import Data.List import Data.Tree (Tree(..),Forest(..)) +import System.Random+ import Math.Combinat.Helper  ------------------------------------------------------- -- * Types +-- | A binary tree with leaves decorated with type @a@. data BinTree a   = Branch (BinTree a) (BinTree a)   | Leaf a@@ -45,6 +63,22 @@ leaf :: BinTree () leaf = Leaf () +-- | A binary tree with leaves and internal nodes decorated +-- with types @a@ and @b@, respectively.+data BinTree' a b+  = Branch' (BinTree' a b) b (BinTree' a b)+  | Leaf' a+  deriving (Eq,Ord,Show,Read)++forgetNodeDecorations :: BinTree' a b -> BinTree a+forgetNodeDecorations (Branch' left _ right) = +  Branch (forgetNodeDecorations left) (forgetNodeDecorations right)+forgetNodeDecorations (Leaf' decor) = Leaf decor +  +instance Functor BinTree where+  fmap f (Branch left right) = Branch (fmap f left) (fmap f right)+  fmap f (Leaf x) = Leaf (f x)+     -------------------------------------------------------  data Paren = LeftParen | RightParen deriving (Eq,Ord,Show,Read)@@ -136,6 +170,17 @@ nestedParentheses :: Int -> [[Paren]] nestedParentheses = fasc4A_algorithm_P +-- | Synonym for 'fasc4A_algorithm_W'.+randomNestedParentheses :: RandomGen g => Int -> g -> ([Paren],g)+randomNestedParentheses = fasc4A_algorithm_W++-- | Synonym for 'fasc4A_algorithm_U'.+nthNestedParentheses :: Int -> Integer -> [Paren]+nthNestedParentheses = fasc4A_algorithm_U++countNestedParentheses :: Int -> Integer+countNestedParentheses = countBinaryTrees+ -- | Generates all sequences of nested parentheses of length 2n. -- Order is lexigraphic (when right parentheses are considered  -- smaller then left ones).@@ -169,11 +214,47 @@ 	      _ -> findj ( lls, [] ) ( reverse rs ++ xs , ys)  	    RightParen -> Just ( reverse ys ++ xs ++ reverse (LeftParen:rs) ++ ls , [] )     +-- | Generates a uniformly random sequence of nested parentheses of length 2n.    +-- Based on \"Algorithm W\" in Knuth.+fasc4A_algorithm_W :: RandomGen g => Int -> g -> ([Paren],g)+fasc4A_algorithm_W n' rnd = worker (rnd,n,n,[]) where+  n = fromIntegral n' :: Integer  +  -- the numbers we use are of order n^2, so for n >> 2^16 +  -- on a 32 bit machine, we need big integers.+  worker :: RandomGen g => (g,Integer,Integer,[Paren]) -> ([Paren],g)+  worker (rnd,_,0,parens) = (parens,rnd)+  worker (rnd,p,q,parens) = +    if x<(q+1)*(q-p) +      then worker (rnd' , p   , q-1 , LeftParen :parens)+      else worker (rnd' , p-1 , q   , RightParen:parens)+    where +      (x,rnd') = randomR ( 0 , (q+p)*(q-p+1)-1 ) rnd +-- | Nth sequence of nested parentheses of length 2n. +-- The order is the same as in 'fasc4A_algorithm_P'.+-- Based on \"Algorithm U\" in Knuth.+fasc4A_algorithm_U +  :: Int               -- ^ n+  -> Integer           -- ^ N; should satisfy 1 <= N <= C(n) +  -> [Paren]+fasc4A_algorithm_U n' bign0 = reverse $ worker (bign0,c0,n,n,[]) where+  n = fromIntegral n' :: Integer+  c0 = foldl f 1 [2..n]  +  f c p = ((4*p-2)*c) `div` (p+1) +  worker :: (Integer,Integer,Integer,Integer,[Paren]) -> [Paren]+  worker (_   ,_,_,0,parens) = parens+  worker (bign,c,p,q,parens) = +    if bign <= c' +      then worker (bign    , c'   , p   , q-1 , RightParen:parens)+      else worker (bign-c' , c-c' , p-1 , q   , LeftParen :parens)+    where+      c' = ((q+1)*(q-p)*c) `div` ((q+p)*(q-p+1))+   ------------------------------------------------------- -- * Binary trees  -- | Generates all binary trees with n nodes. +--   At the moment just a synonym for 'binaryTreesNaive'. binaryTrees :: Int -> [BinTree ()] binaryTrees = binaryTreesNaive @@ -182,7 +263,7 @@ -- This is also the counting function for forests and nested parentheses. countBinaryTrees :: Int -> Integer countBinaryTrees n = binomial (2*n) n `div` (1 + fromIntegral n)-+     -- | Generates all binary trees with n nodes. The naive algorithm. binaryTreesNaive :: Int -> [BinTree ()] binaryTreesNaive 0 = [ leaf ]@@ -193,67 +274,42 @@   , r <- binaryTreesNaive (n-1-i)    ] ------ binary tree zipper--data Ctx a-  = Top -  | L (Ctx a) (BinTree a)-  | R (BinTree a) (Ctx a) --type Loc a = (BinTree a, Ctx a)--left :: Loc a -> Loc a-left (Branch l r , c) = (l , L c r)-left (Leaf _ , _) = error "left: Leaf"--right :: Loc a -> Loc a-right (Branch l r , c) = (r , R l c)-right (Leaf _ , _) = error "right: Leaf"- -top :: BinTree a -> Loc a-top t = (t, Top)- -up :: Loc a -> Loc a-up (t, L c r) = (Branch t r, c)-up (t, R l c) = (Branch l t, c)-up (t, Top  ) = error "up: top"--upmost :: Loc a -> Loc a-upmost l@(t, Top) = l-upmost l = upmost (up l)- -modify :: (BinTree a -> BinTree a) -> Loc a -> Loc a-modify f (t, c) = (f t, c)---------{---- | Generates all binary trees with n nodes.--- Based on \"Algorithm B\" in Knuth, uses tree zippers.-fasc4A_algorithm_B	:: Int -> [BinTree ()]-fasc4A_algorithm_B 0 = [ leaf ]-fasc4A_algorithm_B n = unfold1 next start where-  start = nest n (\t -> Branch t leaf) leaf+-- | Generates an uniformly random binary tree, using 'fasc4A_algorithm_R'.+randomBinaryTree :: RandomGen g => Int -> g -> (BinTree (), g)+randomBinaryTree n rnd = (tree,rnd') where+  (decorated,rnd') = fasc4A_algorithm_R n rnd      +  tree = fmap (const ()) $ forgetNodeDecorations decorated -  killLeft  (Branch _ r) = Branch leaf r-  killRight (Branch l _) = Branch l leaf-  -  next t = case findj (top t) of-    Nothing -> Nothing-    Just locj@(s,c) -> case findk (top s) of-      lock@(u,Top) -> Just $ promote (modify killLeft locj ) lock -      lock@(u,_  ) -> Just $ promote locj (modify killRight lock)+-- | Grows a uniformly random binary tree. +-- \"Algorithm R\" (Remy's procudere) in Knuth.+-- Nodes are decorated with odd numbers, leaves with even numbers (from the+-- set @[0..2n]@). Uses mutable arrays internally.+fasc4A_algorithm_R :: RandomGen g => Int -> g -> (BinTree' Int Int, g)+fasc4A_algorithm_R n0 rnd = res where+  res = runST $ do+    ar <- newArray (0,2*n0) 0+    rnd' <- worker rnd 1 ar+    links <- unsafeFreeze ar+    return (toTree links, rnd')+  toTree links = f (links!0) where+    f i = if odd i +      then Branch' (f $ links!i) i (f $ links!(i+1)) +      else Leaf' i  +  worker :: RandomGen g => g -> Int -> STUArray s Int Int -> ST s g+  worker rnd n ar = do +    if n > n0+      then return rnd+      else do+        writeArray ar (n2-b)   n2+        lk <- readArray ar k+        writeArray ar (n2-1+b) lk+        writeArray ar k        (n2-1)+        worker rnd' (n+1) ar      +    where  +      n2 = n+n+      (x,rnd') = randomR (0,4*n-3) rnd+      (k,b) = x `divMod` 2       -  findj :: Loc () -> Maybe (Loc ())-  findj (Branch (Leaf _) t , c) = case t of-    Branch l r -> findj $ left (Branch t leaf , c) -    Leaf _ -> Nothing-  findj loc@(Leaf _ , c) = Just loc--  findk :: Loc () -> Loc ()-  findk loc@( Branch l (Leaf _) , _) = loc-  findk loc@( Branch l r , _) = findk (right loc)+-------------------------------------------------------         -  promote :: Loc () -> Loc () -> BinTree ()-  promote locj lock = undefined--}+

@@ -1,4 +1,3 @@-#! /usr/bin/env runhaskell
- 
-> import Distribution.Simple
+#! /usr/bin/env runhaskell+> import Distribution.Simple > main = defaultMain

@@ -1,5 +1,5 @@ Name:                combinat-Version:             0.1+Version:             0.2 Synopsis:            Generation of various combinatorial objects. Description:         A collection of functions to generate combinatorial                      objects like partitions, combinations, permutations,@@ -9,10 +9,9 @@ Author:              Balazs Komuves Copyright:           (c) 2008 Balazs Komuves Maintainer:          bkomuves (plus) hackage (at) gmail (dot) com-Stability:           Unstable---Portability:         Portable+Stability:           Experimental Category:            Math-Tested-With:         GHC == 6.8.3+Tested-With:         GHC == 6.10.1 Cabal-Version:       >= 1.2 Build-Type:          Simple @@ -21,11 +20,12 @@  Library   if flag(splitBase)-    Build-Depends:       base >= 3, array, containers+    Build-Depends:       base >= 3, array, containers, random   else     Build-Depends:       base <  3    Exposed-Modules:     Math.Combinat, +                       Math.Combinat.Sets,                        Math.Combinat.Tuples,                         Math.Combinat.Combinations,                        Math.Combinat.Partitions,@@ -34,6 +34,8 @@                        Math.Combinat.Trees      Other-Modules:       Math.Combinat.Helper++  Extensions:          MultiParamTypeClasses, ScopedTypeVariables    Hs-Source-Dirs:      .