diff --git a/Math/Combinat.hs b/Math/Combinat.hs --- a/Math/Combinat.hs +++ b/Math/Combinat.hs @@ -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 diff --git a/Math/Combinat/Helper.hs b/Math/Combinat/Helper.hs --- a/Math/Combinat/Helper.hs +++ b/Math/Combinat/Helper.hs @@ -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 diff --git a/Math/Combinat/Permutations.hs b/Math/Combinat/Permutations.hs --- a/Math/Combinat/Permutations.hs +++ b/Math/Combinat/Permutations.hs @@ -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 ). +randomPermutationDurstenfeld :: RandomGen g => Int -> g -> (Permutation,g) +randomPermutationDurstenfeld = randomPermutationDurstenfeldSattolo False + +-- | Generates a uniformly random /cyclic/ permutation of @[1..n]@. +-- Sattolo's algorithm (see ). +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 diff --git a/Math/Combinat/Sets.hs b/Math/Combinat/Sets.hs new file mode 100644 --- /dev/null +++ b/Math/Combinat/Sets.hs @@ -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 + +------------------------------------------------------- diff --git a/Math/Combinat/Trees.hs b/Math/Combinat/Trees.hs --- a/Math/Combinat/Trees.hs +++ b/Math/Combinat/Trees.hs @@ -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 --} + diff --git a/Setup.lhs b/Setup.lhs --- a/Setup.lhs +++ b/Setup.lhs @@ -1,4 +1,3 @@ -#! /usr/bin/env runhaskell - -> import Distribution.Simple +#! /usr/bin/env runhaskell +> import Distribution.Simple > main = defaultMain diff --git a/combinat.cabal b/combinat.cabal --- a/combinat.cabal +++ b/combinat.cabal @@ -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: .