{- | Count and create combinatorial objects. Also see 'combinat' package. -} module Combinatorics ( permute, permuteFast, permuteShare, permuteRep, choose, variateRep, variate, tuples, partitions, rectifications, setPartitions, pairPartitions, chooseUnrank, chooseUnrankMaybe, chooseRank, factorial, binomial, binomialSeq, binomialGen, binomialSeqGen, multinomial, factorials, binomials, catalanNumber, catalanNumbers, derangementNumber, derangementNumbers, setPartitionNumbers, surjectiveMappingNumber, surjectiveMappingNumbers, fibonacciNumber, fibonacciNumbers, ) where import qualified PowerSeries import qualified Combinatorics.Private as Comb import Data.Function.HT (nest, ) import Data.Maybe.HT (toMaybe, ) import Data.Tuple.HT (mapFst, ) import qualified Data.List.Match as Match import Data.List.HT (mapAdjacent, removeEach, ) import Data.List (genericIndex, ) import Control.Monad (liftM2, ) {- $setup >>> import qualified Combinatorics as Comb >>> import qualified Test.QuickCheck as QC >>> import Test.Utility (equalFuncList, equalFuncList2) >>> >>> import Control.Applicative (liftA2, (<$>)) >>> import qualified Data.List.Match as Match >>> import qualified Data.List.Key as Key >>> import qualified Data.List as List >>> import qualified Data.Set as Set >>> import Data.Tuple.HT (uncurry3) >>> import Data.List.HT (allEqual, isAscending) >>> import Data.List (sort, nub) >>> import Data.Eq.HT (equating) >>> >>> genPermuteRep :: Int -> QC.Gen [(Char, Int)] >>> genPermuteRep n = do >>> xns <- QC.listOf $ liftA2 (,) QC.arbitrary $ QC.choose (0,n) >>> return $ Match.take (takeWhile (<=n) $ scanl1 (+) $ map snd xns) xns >>> >>> genVariate :: QC.Gen [Char] >>> genVariate = take 7 <$> QC.arbitrary >>> >>> genBinomial :: QC.Gen (Integer, Integer) >>> genBinomial = do >>> n <- QC.choose (0,100) >>> k <- QC.choose (0,n) >>> return (n,k) >>> >>> genChooseIndex :: QC.Gen (Integer, Integer, Integer) >>> genChooseIndex = do >>> n <- QC.choose (0,25) >>> k <- QC.choose (0,n) >>> i <- QC.choose (0, Comb.binomial n k - 1) >>> return (n,k,i) -} {-* Generate compositions from a list of elements. -} -- several functions for permutation -- cf. Equation.hs {- | Generate list of all permutations of the input list. The list is sorted lexicographically. >>> Comb.permute "abc" ["abc","acb","bac","bca","cab","cba"] >>> Comb.permute "aabc" ["aabc","aacb","abac","abca","acab","acba","aabc","aacb","abac","abca","acab","acba","baac","baca","baac","baca","bcaa","bcaa","caab","caba","caab","caba","cbaa","cbaa"] prop> QC.forAll (take 6 <$> QC.arbitrary :: QC.Gen [Int]) $ \xs -> allEqual $ map (\p -> sort (p xs)) $ Comb.permute : Comb.permuteFast : Comb.permuteShare : [] -} permute :: [a] -> [[a]] permute = Comb.permuteRec {- | Generate list of all permutations of the input list. It is not lexicographically sorted. It is slightly faster and consumes less memory than the lexicographical ordering 'permute'. -} permuteFast :: [a] -> [[a]] permuteFast x = permuteFastStep [] x [] {- | Each element of (allcycles x) has a different element at the front. Iterate cycling on the tail elements of each element list of (allcycles x). -} permuteFastStep :: [a] -> [a] -> [[a]] -> [[a]] permuteFastStep suffix [] tl = suffix:tl permuteFastStep suffix x tl = foldr (\c -> permuteFastStep (head c : suffix) (tail c)) tl (allCycles x) {- | All permutations share as much suffixes as possible. The reversed permutations are sorted lexicographically. -} permuteShare :: [a] -> [[a]] permuteShare x = map fst $ -- map (\(y,[]) -> y) $ -- safer but inefficient nest (length x) (concatMap permuteShareStep) [([], x)] permuteShareStep :: ([a], [a]) -> [([a], [a])] permuteShareStep (perm,todo) = map (mapFst (:perm)) (removeEach todo) {- | >>> Comb.permuteRep [('a',2), ('b',1), ('c',1)] ["aabc","aacb","abac","abca","acab","acba","baac","baca","bcaa","caab","caba","cbaa"] prop> QC.forAll (genPermuteRep 7) $ \xs -> let perms = Comb.permuteRep $ Key.nub fst xs in perms == nub perms prop> QC.forAll (genPermuteRep 10) $ \xs -> let perms = Comb.permuteRep $ Key.nub fst xs in List.sort perms == Set.toList (Set.fromList perms) prop> QC.forAll (genPermuteRep 10) $ isAscending . Comb.permuteRep . Key.nub fst . sort prop> QC.forAll (QC.choose (0,10)) $ \n k -> Comb.choose n k == Comb.permuteRep [(False, n-k), (True, k)] -} permuteRep :: [(a,Int)] -> [[a]] permuteRep = Comb.permuteRep {- | >>> map (map (\b -> if b then 'x' else '.')) $ Comb.choose 5 3 ["..xxx",".x.xx",".xx.x",".xxx.","x..xx","x.x.x","x.xx.","xx..x","xx.x.","xxx.."] >>> map (map (\b -> if b then 'x' else '.')) $ Comb.choose 3 5 [] prop> QC.forAll (QC.choose (0,10)) $ \n k -> all (\x -> n == length x && k == length (filter id x)) (Comb.choose n k) -} choose :: Int -> Int -> [[Bool]] choose = Comb.chooseRec {- | Generate all choices of n elements out of the list x with repetitions. \"variation\" seems to be used historically, but I like it more than \"k-permutation\". >>> Comb.variateRep 2 "abc" ["aa","ab","ac","ba","bb","bc","ca","cb","cc"] -} variateRep :: Int -> [a] -> [[a]] variateRep = Comb.variateRep {- | Generate all choices of n elements out of the list x without repetitions. >>> Comb.variate 2 "abc" ["ab","ac","ba","bc","ca","cb"] >>> Comb.variate 2 "abcd" ["ab","ac","ad","ba","bc","bd","ca","cb","cd","da","db","dc"] >>> Comb.variate 3 "abcd" ["abc","abd","acb","acd","adb","adc","bac","bad","bca","bcd","bda","bdc","cab","cad","cba","cbd","cda","cdb","dab","dac","dba","dbc","dca","dcb"] prop> QC.forAll genVariate $ \xs -> Comb.variate (length xs) xs == Comb.permute xs prop> \xs -> equating (take 1000) (Comb.variate (length xs) xs) (Comb.permute (xs::String)) -} variate :: Int -> [a] -> [[a]] variate = Comb.variateRec {- | Generate all choices of n elements out of the list x respecting the order in x and without repetitions. >>> Comb.tuples 2 "abc" ["ab","ac","bc"] >>> Comb.tuples 2 "abcd" ["ab","ac","ad","bc","bd","cd"] >>> Comb.tuples 3 "abcd" ["abc","abd","acd","bcd"] -} tuples :: Int -> [a] -> [[a]] tuples = Comb.tuplesRec {- | >>> Comb.partitions "abc" [("abc",""),("bc","a"),("ac","b"),("c","ab"),("ab","c"),("b","ac"),("a","bc"),("","abc")] prop> QC.forAll genVariate $ \xs -> length (Comb.partitions xs) == 2 ^ length xs -} partitions :: [a] -> [([a],[a])] partitions = foldr (\x -> concatMap (\(lxs,rxs) -> [(x:lxs,rxs), (lxs,x:rxs)])) [([],[])] {- | Number of possibilities arising in rectification of a predicate in deductive database theory. Stefan Brass, \"Logische Programmierung und deduktive Datenbanken\", 2007, page 7-60 This is isomorphic to the partition of @n@-element sets into @k@ non-empty subsets. >>> Comb.rectifications 4 "abc" ["aabc","abac","abbc","abca","abcb","abcc"] >>> map (length . uncurry Comb.rectifications) $ do x<-[0..10]; y<-[0..x]; return (x,[1..y::Int]) [1,0,1,0,1,1,0,1,3,1,0,1,7,6,1,0,1,15,25,10,1,0,1,31,90,65,15,1,0,1,63,301,350,140,21,1,0,1,127,966,1701,1050,266,28,1,0,1,255,3025,7770,6951,2646,462,36,1,0,1,511,9330,34105,42525,22827,5880,750,45,1] prop> QC.forAll (QC.choose (0,7)) $ \k xs -> isAscending . Comb.rectifications k . nub . sort $ (xs::String) -} rectifications :: Int -> [a] -> [[a]] rectifications = let recourse _ 0 xt = if null xt then [[]] else [] recourse ys n xt = let n1 = pred n in liftM2 (:) ys (recourse ys n1 xt) ++ case xt of [] -> [] (x:xs) -> map (x:) (recourse (ys++[x]) n1 xs) in recourse [] {- | Their number is @k^n@. -} {- setPartitionsEmpty :: Int -> [a] -> [[[a]]] setPartitionsEmpty k = let recourse [] = [replicate k []] recourse (x:xs) = map (\(ys0,y,ys1) -> ys0 ++ [x:y] ++ ys1) $ concatMap splitEverywhere (recourse xs) {- do xs1 <- recourse xs (ys0,y,ys1) <- splitEverywhere xs1 return (ys0 ++ [x:y] ++ ys1) -} in recourse -} -- TestMe: isAscending . Comb.setPartitions k . nub . sort setPartitions :: Int -> [a] -> [[[a]]] setPartitions 0 xs = if null xs then [[]] else [ ] setPartitions _ [] = [] setPartitions 1 xs = [[xs]] -- unnecessary for correctness, but useful for efficiency setPartitions k (x:xs) = do (rest, choosen) <- partitions xs part <- setPartitions (pred k) rest return ((x:choosen) : part) {- | All ways of separating a list of terms into pairs. All partitions are given in a canonical form, sorted lexicographically. The canonical form is: The list of pairs is ordered with respect to the first pair members, and the elements in each pair are ordered. The order is implied by the order in the input list. -} pairPartitions :: [a] -> [[(a,a)]] pairPartitions xs = if odd (length xs) then [] else let go (y:ys) = do (z,zs) <- removeEach ys fmap ((y,z):) $ go zs go [] = [[]] in go xs {-* Rank and unrank combinatorial objects. -} {- | @chooseUnrank n k i == choose n k !! i@ prop> QC.forAll (QC.choose (0,10)) $ \n k -> map (Comb.chooseUnrank n k) [0 .. Comb.binomial n k - 1] == Comb.choose n k prop> QC.forAll genChooseIndex $ \(n,k,i) -> Comb.chooseRank (Comb.chooseUnrank n k i) == (n, k, i) prop> \bs -> uncurry3 Comb.chooseUnrank (Comb.chooseRank bs :: (Integer, Integer, Integer)) == bs -} chooseUnrank :: Integral a => a -> a -> a -> [Bool] chooseUnrank = Comb.chooseUnrankRec chooseUnrankMaybe :: Int -> Int -> Int -> Maybe [Bool] chooseUnrankMaybe n k i = toMaybe (0 <= i && i < binomial n k) (chooseUnrank n k i) -- error ("chooseUnrank: out of range " ++ show (n, k, i)) {- | -} chooseRank :: Integral a => [Bool] -> (a, a, a) chooseRank = foldl (\(n,k0,i0) (bins,b) -> let (k1,i1) = if b then (succ k0, i0 + genericIndex (bins++[0]) k1) else (k0,i0) in (succ n, k1, i1)) (0,0,0) . zip binomials . reverse {-* Generate complete lists of combinatorial numbers. -} {- | prop> QC.forAll (take 8 <$> QC.arbitrary) $ \xs -> length (Comb.permute xs) == Comb.factorial (length (xs::String)) prop> QC.forAll (take 6 <$> QC.arbitrary) $ \xs -> sum (map sum (Comb.permute xs)) == sum xs * Comb.factorial (length xs) -} factorial :: Integral a => a -> a factorial n = product [1..n] {- | Pascal's triangle containing the binomial coefficients. prop> QC.forAll (QC.choose (0,12)) $ \n k -> length (Comb.choose n k) == Comb.binomial n k prop> QC.forAll genBinomial $ \(n,k) -> let (q, r) = divMod (Comb.factorial n) (Comb.factorial k * Comb.factorial (n-k)) in r == 0 && Comb.binomial n k == q prop> QC.forAll (take 16 <$> QC.arbitrary) $ \xs k -> length (Comb.tuples k xs) == Comb.binomial (length (xs::String)) k -} binomial :: Integral a => a -> a -> a binomial = Comb.binomial binomialSeq :: Integral a => a -> [a] binomialSeq = Comb.binomialSeq binomialGen :: (Integral a, Fractional b) => b -> a -> b binomialGen n k = genericIndex (binomialSeqGen n) k binomialSeqGen :: (Fractional b) => b -> [b] binomialSeqGen n = scanl (\acc (num,den) -> acc*num / den) 1 (zip (iterate (subtract 1) n) (iterate (1+) 1)) {- | prop> QC.forAll (genPermuteRep 10) $ \xs -> length (Comb.permuteRep xs) == Comb.multinomial (map snd xs) prop> QC.forAll (QC.listOf $ QC.choose (0,300::Integer)) $ \xs -> Comb.multinomial xs == Comb.multinomial (sort xs) -} multinomial :: Integral a => [a] -> a multinomial = product . mapAdjacent binomial . scanr1 (+) {-* Generate complete lists of factorial numbers. -} {- | prop> equalFuncList Comb.factorial Comb.factorials 1000 -} factorials :: Num a => [a] factorials = Comb.factorials {-| Pascal's triangle containing the binomial coefficients. Only efficient if a prefix of all rows is required. It is not efficient for picking particular rows or even particular elements. prop> equalFuncList2 Comb.binomial Comb.binomials 100 -} binomials :: Num a => [[a]] binomials = Comb.binomials {- | @catalanNumber n@ computes the number of binary trees with @n@ nodes. -} catalanNumber :: Integer -> Integer catalanNumber n = case divMod (binomial (2*n) n) (n+1) of (c,0) -> c _ -> error "catalanNumber: Integer implementation broken" {- | Compute the sequence of Catalan numbers by recurrence identity. It is @catalanNumbers !! n == catalanNumber n@ prop> equalFuncList Comb.catalanNumber Comb.catalanNumbers 1000 -} catalanNumbers :: Num a => [a] catalanNumbers = let xs = 1 : PowerSeries.mul xs xs in xs derangementNumber :: Integer -> Integer derangementNumber n = sum (scanl (*) ((-1) ^ mod n 2) [-n,1-n..(-1)]) {- | Number of fix-point-free permutations with @n@ elements. prop> equalFuncList Comb.derangementNumber Comb.derangementNumbers 1000 -} derangementNumbers :: Num a => [a] derangementNumbers = Comb.derangementNumbersPS0 -- generation of all possibilities and computation of their number should be in different modules {- | Number of partitions of an @n@ element set into @k@ non-empty subsets. Known as Stirling numbers . prop> QC.forAll (QC.choose (0,10000)) $ \k -> QC.forAll (take 7 <$> QC.arbitrary) $ \xs -> length (Comb.setPartitions k xs) == (Comb.setPartitionNumbers !! length (xs::String) ++ repeat 0) !! k prop> QC.forAll (QC.choose (0,7)) $ \k xs -> length (Comb.rectifications k xs) == (Comb.setPartitionNumbers !! k ++ repeat 0) !! length (xs::String) -} setPartitionNumbers :: Num a => [[a]] setPartitionNumbers = Comb.setPartitionNumbers {- | @surjectiveMappingNumber n k@ computes the number of surjective mappings from a @n@ element set to a @k@ element set. -} surjectiveMappingNumber :: Integer -> Integer -> Integer surjectiveMappingNumber n k = foldl subtract 0 $ zipWith (*) (map (^n) [0..]) (binomialSeq k) {- | prop> equalFuncList2 Comb.surjectiveMappingNumber Comb.surjectiveMappingNumbers 20 -} surjectiveMappingNumbers :: Num a => [[a]] surjectiveMappingNumbers = Comb.surjectiveMappingNumbersPS {- | Multiply two Fibonacci matrices, that is matrices of the form > /F[n-1] F[n] \ > \F[n] F[n+1]/ -} fiboMul :: (Integer,Integer,Integer) -> (Integer,Integer,Integer) -> (Integer,Integer,Integer) fiboMul (f0,f1,f2) (g0,g1,g2) = let h0 = f0*g0 + f1*g1 h1 = f0*g1 + f1*g2 -- h1 = f1*g0 + f2*g1 h2 = f1*g1 + f2*g2 in (h0,h1,h2) {- Fast computation using matrix power of > /0 1\ > \1 1/ Hard-coded fast power with integer exponent. Better use a generic algorithm. -} fibonacciNumber :: Integer -> Integer fibonacciNumber x = let aux 0 = (1,0,1) aux (-1) = (-1,1,0) aux n = let (m,r) = divMod n 2 f = aux m f2 = fiboMul f f in if r==0 then f2 else fiboMul (0,1,1) f2 (_,y,_) = aux x in y {- | Number of possibilities to compose a 2 x n rectangle of n bricks. > ||| |-- --| > ||| |-- --| prop> equalFuncList Comb.fibonacciNumber Comb.fibonacciNumbers 10000 -} fibonacciNumbers :: [Integer] fibonacciNumbers = let xs = 0 : ys ys = 1 : zipWith (+) xs ys in xs {- * Auxiliary functions -} {- candidates for Useful -} {- | Create a list of all possible rotations of the input list. -} allCycles :: [a] -> [[a]] allCycles x = Match.take x (map (Match.take x) (iterate tail (cycle x)))